# On the Probability of Magnus Carlsen reaching 2900

Sohan Bendre
Chess Ed
Shiva Maharaj
Chess Ed
Nick Polson111Nick Polson is Professor at Chicago Booth<EMAIL_ADDRESS>
Booth School of Business
University of Chicago
Vadim Sokolov222Vadim Sokolov is an Assistant Professor at Operations Research
at George Mason University<EMAIL_ADDRESS>
Department of Systems Engineering
and Operations Research
George Mason University

###### Abstract

How likely is it that Magnus Carlsen will achieve an Elo rating of $2900$?
This has been a goal of Magnus and is of great current interest to the chess
community. Our paper uses probabilistic methods to address this question. The
probabilistic properties of Elo’s rating system have long been studied and we
provide an application of such methods. By applying a Brownian motion model of
Stern as a simple tool we provide answers. Our research also has fundamental
bearing on the choice of the $K$-factor used in Elo’s system for GrandMaster
(GM) chess play. Finally, we conclude with a discussion of policy issues
involved with the choice of $K$-factor.

Key Words: Elo, Chess, GrandMaster, $K$-factor, Maxwell-Boltzmann, Magnus
Carlsen

##  1 Introduction

Recently, chess grandmaster Magnus Carlsen announced his newest goal, to reach
a 2900 FIDE rating. Carlsen (2022, Personal Communication) saw the problem as
follows: he realized that he needed a _”hot streak”_ to achieve his goal. We
analyze Magnus’ hot streak at the beginning of 2019 where his rating jumped
from 2835 to 2881 in a matter of months to analyze how well he would have to
play to achieve his goal of 2900 with a reasonable success. This was also one
of his main questions, and he was also curious as to how probabilistic methods
could help in his understanding of how his Elo Elo (1978) rating changes. We
directly address two problems. First, how likely is it that Magnus Carlsen
will achieve an Elo rating of $2900$ and within what time period? Second, what
is the best possible strategy for Magnus to achieve his goal?

Many argue that current implementation of Elo’s system is the reason that
Carlsen has such a difficult time increasing his rating. Many believe that the
problem lies in the value of the $K$-factor. For grandmaster (GM) play,
$K=10$, the lowest possible $K$-value for a FIDE rated chess player. The
$K$-factor governs the variability in rating changes. Carlsen (2022, personal
communication) suggested that a $K=15$ would give him a reasonable chance of
achieving his goal. Moreover, he suggested that he thought he was fairly rated
by the system at the current moment and a $K=20$ would be too large for GM
play.

As of August 2022, Carlsen’s official rating is 2860, far above 2nd ranked
Ding Liren with a rating of 2806. At first sight, it might seem that Carlsen’s
goal is very achievable to the average chess player, whose rating might change
by 50 points or more in a single tournament. But for Carlsen it has proven to
be quite difficult, as his rating hasn’t fluctuated by more than 20 points in
the past two years. In reality, it seems virtually impossible that Carlsen can
ever reach 2900, due to the constraints of the current chess rating system.
Specifically, the choice of $K$-factor used in Elo’s rating system for chess
grandmaster (GM) play is only $K=10$, which governs how many points he can
gain with a win. Intuitively, as he is already the best player in the world
(by a wide margin) if he draws (which is common in GM play) he will lose
points and, of course, any loss will lead to a large setback. For example,
after an incredible performance at the Tata Steel tournament ($9.5/13$) he
only gained $3$ points only to lose this in a loss against a much lower rated
player.

###  1.1 Elo Rating System in Chess

The basic Elo (Elo, 1978) probability model for assessing the probability of
$A$’s play against $B$ is

$P(A)=W(R_{A}-R_{B}),$

where $P(A)$ is the probability of $A$ wining and $R_{A}$ is the strength of
$A$.

The function $W$ maps the difference in scores into $(0,1)$ interval and is
symmetric, monotonically increasing. Further, it approaches one as difference
in ratings grows. The Elo rating system is used by the World Chess Federation
(FIDE) and numerous other chess organizations to determine how many points
need to be assigned to each player after a game. The Elo rating system is a
statistical model for estimating the expected score of a player after a game.
Elo’s system was adopted by FIDE in $1970$ and has been in use ever since. A
few other scoring systems have been recently proposed, including the
TrueSkill, created by Microsoft Research Dangauthier et al. (2007), the Glicko
rating system created by statistician Mark Glickman Glickman and Jones (1999),
and the Chessmetrics rating system created by Jeff Sonas. Nevertheless, Elo’s
system continues to be the prevalent chess rating system given its large
success.

In Elo’s system the difference of $G$ points is equivalent the odds ratio of
$K:1$. The $K$-factor used by FIDE is equal to 10, and the logistic growth
rate $G$ is $400$. A difference of $400$ means higher scored player is likely
to win $10$ out $11$ games. Hence, the FIDE Elo formula for calculating the
odds of a game between $A$ and $B$ with corresponding ratings of $R_{A}$ and
$R_{B}$ is as follows

$\mathrm{Odds}(A~{}\mathrm{beats}~{}B)=\dfrac{P(A)}{P(B)}=\dfrac{P(A)}{1-P(A)}=K^{\frac{R_{A}-R_{B}}{400}},$

where $P(A)$ is the probability that player Re-expressing this gives

$P(A)=\dfrac{1}{1+10^{-\frac{R_{A}-R_{B}}{400}}}.$

After a game, the Elo rating of each player is updated to adjust for the
observed outcome of win $(1)$, loss $(0)$ or draw $(0.5)$. The Elo rating of
the player A is updated as follows

$R_{A}^{+}=R_{A}+K(S_{A}-P(A))\;\;{\rm where}\;\;S_{A}\in\\{0,1,1/2\\}$

is the result of the game. At the end of a tournament, adjustments for all the
games will be applied to update the post tournament rating. Our approach
builds on the Elo ratings analysis of Aldous (2017).

Elo’s system is considered to be efficient and become the standard in many
games of play. For example, it was shown that when it was applied to outcomes
of tennis matches, it performed poorly due to the inefficiency how points get
assigned to ATP (men) and WTA (women) players Williams et al. (2021). However,
it is still a topic of debate whether this system can be improved. Our view is
that a discussion of the $K$-factor and the implication for player’s ratings
is the more pressing question. Section 4 provides a further discussion.

To a large extent, as Elo’s model makes clear, the larger values of $K$ will
allow for greater fluctuation in rating. Due to this, it is quite difficult
for Carlsen’s rating to increase by a large amount, but it also does ensure
that his rating won’t immensely decrease either. If the $K$ factor was
increased to allow greater rating changes, theoretically if Carlsen won all
his games he may be able to reach the $2900$ threshold. However, drawing or
losing a game would have a negative impact on his rating, because he is so
much more highly rated than all other players. For example, with the current
$K$ factor, if Carlsen drew against the 2nd-ranked player in the world, Ding
Liren (rating 2806), his rating would decrease by about 0.83 points. If
instead he lost to Ding Liren, his rating would drop by a dramatic $5.83$
points. However, increasing $K$ to $20$ (too large in Carlsen’s view), a draw
would cause his rating to drop about $1.65$ points and a loss would mean a
drop of about $11.65$ points. Considering that draws at the grandmaster level
are occurring $70-75$% of the time nowadays, the situation seems tough for
Carlsen any way you look at it due to the dramatic increase in competition.
Reaching such a goal would be one of the greatest sporting achievements.

The rest of the paper is outlined as follows. Section 1.2 provides a brief
review of Elo ratings and how they pertain to chess. Section 2 provides a
simple Brownian motion model for assessing the two fundamental questions of
how likely is it that Magnus Carlsen will reach 2900, and at what level of
play would he have to achieve to have a reasonable chance of achieving $2900$?
To do this we use the Brownian motion model for sports scores as originally
developed by Stern (1991) and extended by a number of authors, Feng et al.
(2016) for EPL, Polson and Stern (2015) for implied volatility of a sports
game. Then we provide an empirical analysis of Magnus Carlsen’s chess games
and changes in Elo ratings. In particular, we analyze Magnus’ hot streak at
the beginning of 2019. We use this to analyze how well he would have to play
to achieve his goal of 2900 with a reasonable success. Finally, Section 3
concludes with a discussion of the choice of $K$-factors in chess.

###  1.2 Elo Ratings and Changing Strengths

The key feature of the rating is that it needs to change in time and actually
reflect current strength of the player. The Elo rating system has several
important advantages

  * •

Dynamic with a rating that gets updated after each game or tournament

  * •

Implicitly weight recent games higher, and work as exponential smoother

  * •

The updating model is straightforward to understand and implement.

However, there is a still a question of how accurate the rating is and does it
truly reflect a current ability of a player. Aldous (2017) provides several
asymptotic theoretical results. However, it is also noted that asymptotic
results are not relevant. He shows that, the typical error of predicted win
probabilities will not be substantially less than 10%, regardless of number of
matches played.

Glickman and Jones (1999) also showed empirically that the quality of
prediction differs for player of different rank. He used maximum likelihood
estimator (MLE) to identify the adjusting factor $\alpha$

$P(A)=\dfrac{1}{1+10^{-\alpha\frac{R_{A}-R_{B}}{400}}}.$

He found out that for all groups the best estimate of $\alpha<1$, meaning that
the growth rate has to be more than 400. This suggests, that rating need to be
updated “faster”. This also can be achieved by increasing the K-value.

Figure 1: Probability curve for different values of $\alpha$

Glickman and Jones (1999) estimated that $\alpha=0.59$ for player with the
rating between 1400 and 1600, and it is $0.95$ for players with rating between
2200 and 2700. Thus, the Elo performs most poorly on the middle-range players,
that happen to be the largest group

Figure 2: Histogram of Elo ratings of player registered with FIDE as of August
2022

The other important feature of the Elo system is the exponential decay of
increase in probabilities. In other words, the probability of win (and thus
adjustment to the Elo rating) increases exponentially slow as the difference
in ratings goes up. Meaning, then the gap between Carlsen and the number two
is large, it is very hard for Carlsen to improve the rating. To show this fact
empirically, we plot the histogram of Carlsen’s rating changes (per game) over
the last two years (total of 110 classical games)

Figure 3: Histogram of Magnus Carlsen’s change in rating in 2020-2022 period

As we can see from Figure 3, that the histogram is skewed to the left.

##  2 The Probability that Magnus Carlsen will make $2900$

In order to address the question of how likely is it for Magnus to make $2900$
we need a probability model. Stern (1988) provides a Brownian motion model for
addressing the probabilities of winning a sports contest. Feng et al. (2016)
extend this model to a Skellam process for the outcome of EPL football matches
which have the same distribution profile as chess $(win,draw,loss)$. Whilst
this leads to more accurate assessments, it has the caveat of
interpretability. Stern’s Brownian motion model is easily interpretable and
can be view as a limit of the underlying process which due to small changes in
the Elo ratings will be an excellent approximation.

Figure 1 plots his rating changes for the whole period. Figure 2 for the a
recent sample of $110$ classical games (the last two years) which we will use
to estimate model parameters. Table 1 takes the last $110$ games of Magnus
Carlsen and shows his empirical $(win,draw,loss)$ probabilities and his
empirical gains/loses in his Elo ratings. We denote these by
$(\hat{p}_{1},\hat{p}_{2},\hat{p}_{3})$ and
$(\hat{e}_{1},\hat{e}_{2},\hat{e}_{3})$, respectively. $\hat{\cdot}$’s are
used to indicate that they are estimates. Notice that as Carlsen is the best
player in the world, he loses Elo rating points when he draws even though this
is the most likely outcome. This skewness in the distribution is what makes
his goal of $2900$ extremely hard. Appendix A gives specific dates, results
and rating changes.

Figure 3 plots his rating changes for the hot streak period at the beginning
of 2019. Table 2 provides the same set of statistics, but for his hot streak
in $2019$. Notice the difference in the estimates. His extremely high win
percentage (hot streak) allowed him to gain over $50$ points in a short period
of time.

##### Brownian Motion Model.

Let $X_{t}$ denote Magnus Carlsen’s current rating at time $t$. Let $X_{0}$
denote the Elo rating at the beginning of the time period. We assume that
$X_{0}$ is a fair assessment at the beginning of the period and the Elo system
is in equilibrium. Let $\mu$ denote the instantaneous skill level over the
period $(0,t)$. This might result in a change in performance due to extra
effort. Stern (1991) shows that changes in outcome approximately follow a
Brownian motion model of the form

$X_{t}=X_{0}+\mu t+\sigma B_{t}$

where $\sigma$ denotes the volatility of the Elo rating and $B_{t}$ is a
standard Brownian motion. In particular, $B_{0}=0$ and $B_{t}\sim N(0,t)$.

Let $X_{MagnusCarlsen}$ denote the random variable that measures the changes
in Magnus Carlsen’s ratings. Then we can estimate mean, $\hat{\mu}_{X}$ and
standard deviation, $\hat{\sigma}_{X}$, via the formulas

$\hat{\mu}_{X}=\sum_{i=1}^{3}\hat{e}_{i}\hat{p}_{i}\;\;{\rm
and}\;\;\hat{\sigma}_{X}^{2}=\sum_{i=1}^{3}\hat{e}_{i}^{2}\hat{p}_{i}-\hat{\mu}_{X}^{2}$

summed over the three possibilities of $.(win,draw,loss)$. Put simply, we
estimate mean and standard deviation from empirical averages.

Following Stern, then we can calculate the desired probability—-how likely is
it for Magnus to reach his goal of $3000$, starting at $X_{0}=2860$ with his
current skill and volatility of his Elo rating from the last two years, namely

$\displaystyle\mathbb{P}(X_{T}>3000|X_{0}=2860)$
$\displaystyle=\mathbb{P}\left(Z>\frac{36-\hat{\mu}_{X}T}{\hat{\sigma}_{X}\sqrt{T}}\right)$
$\displaystyle\mathbb{P}(X_{T}>3000|X_{0}=2860)$
$\displaystyle=1-\Phi\left(\frac{40-\hat{\mu}_{X}T}{\hat{\sigma}_{X}\sqrt{T}}\right)$

where $\Phi$ denotes the normal cumulative distribution function.

Table 3 shows how likely on a one year and three year basis. Under current
circumstances, as intuitively thought, it is high unlikely he will be able to
make his goal with an estimated probability $\hat{p}=1/40$. This increases
over the three year period, but again is small, $\hat{p}=1/10$.

There are two further questions we can address with our model

###  2.1 Hot Streak of 2019. What level of play is required to achieve $2900$?

Let $p$ denote the probability of achieving the goal. Let $\mu$ denote the
required level of skill to achieve this goal Following Polson and Stern (2018)
who define the implied volatility of a sports game, we can use Elo’s formula
to infer an implied skill level (a.k.a. performance) to achieve a desired
probability of success of achieving their goal.

From the Stern calculation, with $X_{T}-X_{0}=40$ where $X_{T}$ is the desired
goal

$1-p=\Phi\left(\frac{X_{T}-X_{0}-\mu_{X}T}{\hat{\sigma}_{X}\sqrt{T}}\right)$
(1)

This can be inverted to address the level of play, denoted by $\mu_{imp}$ for
implied level of play, required to guarantee the given $p$.

$\mu_{imp}=\frac{1}{T}(X_{T}-X_{0})-\frac{1}{\sqrt{T}}\Phi^{-1}\left(1-p\right)\hat{\sigma}_{X}$

This is an instantaneous level of play whose advantage compounds with $T$.

Table 4 shows the probabilities if he plays like $2019$ and also that if he
wants a 50/50 chance in the next year of achieving his goal, he’ll have to
play at a level of $\mu=$.

###  2.2 Simulated Paths

We simulate two paths for Magnus Carlsen. First simulation is based on his
most recent performance during the 2020-2022 period. Our data set contains 110
classical games from this period, starting from Tata Steel Masters 2020 and
ending with World Chess Olympiad 2022.

Figure 4(a) shows the histogram of the ratings of Magnus Carlsen’s opponents
during this period. Figure 4(b) compares the predicted outcome by the Elo
system (blue dots) and the actual outcome (black dots)

|
---|---
(a) Opponent Ratings | (b) Game Outcomes
Figure 4: Opponent Ratings and Game Outcomes for 2020-2022 period.

Table 1 summarizes the data set the observed changes in rating after each
game. We have a total of 110 games in this period with average change in
rating being -0.11 with standard deviation of 2.67.

N | $\mu$ | $\sigma$ | Min | Max
---|---|---|---|---
110 | -0.11 | 2.67 | -7.44 | 4.5
Table 1: Summary statistic for change in rating per game during the 2020-2022
period.

We simulate 2000 times the 200 future games by Magnus Carlsen. First, we
assume that the performance (likely outcome) and the mix of opponents is the
same as in 2020-2022 period.

|
---|---
(a) How many games till 2900 | (b) Simulated Trajectories
Figure 5: Left panel: Histogram of when the rating of 2900 is first reached
over the simulated 1000 games. Right Panel: Several Simulated Trajectories
over the next 200 games. The solid line shows the average across 2000
simulations for this specific game

In this simulation. Magnus reached the rating in 90 simulated trajectories out
of 2000. If Magnus Carlsen continues showing the same performance as he did
during the 2020-2022 period, he has 4.5% chance of reaching 2900.

####  2.2.1 Magnus’ hot streak in 2019

Now, we perform same simulations using the data from 2019’s Magnus’ hot streak
period. During this period the rating went from 2835 to 2872 and peaked in
August at 2882. The Figure 6 shows the rating for each of the months in 2019.

Figure 6: Elo rating for every month of 2019.

Table 2 summarizes the data set the observed changes in rating after each
game. We have a total of 78 games in this period with average change in rating
being 0.48 with standard deviation of 2.42.

N | $\mu$ | $\sigma$ | Min | Max
---|---|---|---|---
78 | 0.48 | 2.42 | -2.8 | 4.7
Table 2: Summary statistic for change in rating per game during the 2019
period.

Figure 7 shows the histogram of the ratings of Magnus Carlsen’s opponents
during this period.

Figure 7: Opponent Ratings

No we simulate 2000 times the 200 future games by Magnus Carlsen using the
2019 data. Again, we assume that the performance (likely outcome) and the mix
of opponents is the same as in 2019 period.

|
---|---
(a) How many games till 2900 | (b) Simulated Trajectories
Figure 8: Left panel: Histogram of when the rating of 2900 is first reached
over the simulated 1000 games. Right Panel: Several Simulated Trajectories
over the next 200 games. The solid line shows the average across 2000
simulations for this specific game

In this simulation. Magnus reached the rating in 1600 simulated trajectories
out of 2000. If Magnus Carlsen continues showing the same performance as he
did during the 2019 period, he has 80% chance of reaching 2900.

###  2.3 What happens if we increase the $K$-factor?

Increasing the $K$-factor, say to $K=15$, essentially changes the volatility
in our model. From equation (1), we see that this increases the probability of
a tail probability, namely the probability of getting to $2900$. The other way
of increasing your probability is to increase your skill level, $\mu$, as
described above. Here, we assume that the skill level does not change and the
average number of wins is the same, but consider a scenario, when K-factor is
15, rather than 10, and using 2020-2022 data:

|
---|---
(a) How many games till 2900 | (b) Simulated Trajectories
Figure 9: Left panel: Histogram of when the rating of 2900 is first reached
over the simulated 1000 games. Right Panel: Several Simulated Trajectories
over the next 200 games. The solid line shows the average across 2000
simulations for this specific game

If Magnus Carlsen continues showing the same performance as he did during the
2019 period, and K-factor is 15, he has 18% chance of reaching 2900.

When K-factor is 15, and we are using 2019 data:

|
---|---
(a) How many games till 2900 | (b) Simulated Trajectories
Figure 10: Left panel: Histogram of when the rating of 2900 is first reached
over the simulated 1000 games. Right Panel: Several Simulated Trajectories
over the next 200 games. The solid line shows the average across 2000
simulations for this specific game

If Magnus Carlsen continues showing the same performance as he did during the
2019 period, and K-factor is 15, he has 95% chance of reaching 2900.

Table 3 summarizes the findings our simulation studies and shows the
probability of reaching the rating of 2900 under two different assumptions
about the skill level (2019 performance vs 2020-2022 performance) and the
K-factor (10 vs 15).

| K=10 | K=15
---|---|---
2020-2022 Data | 4.5% | 18%
2019 Data | 80% | 95%
Table 3: Chances of reaching 2900 in the next 200 games under different future
performance assumptions (2019 vs 2020-2022) and different K factors (10 vs 15)

##  3 Discussion

The problem of lack of volatility in Elo system has been previously discussed
in the literature. For example, Aldous (2017) talks about scaling the Elo
rating. This is an important policy issues, which we believe has to be careful
analyzed using statistical and probabilistic models. Further, as shown in the
Figure below, the calibration of $K$ needs to be done differently for
different skill levels.

Figure 11: Comparing sample distribution of results for different rating
groups. Credit: Viswanath (2016)

In this paper we have demonstrated the effect of changing K-factor on the
probability of Carlsen reaching rating of 2900 in the next 200 games.
Specifically, we have shown that it goes from 4.5% (under the current system)
to 18% when K-factor of 15 is used.

There are a number of ares for future research. $K$-factor can be calibrated
from the Brier score. As the $K$-factor is equivalent to a probabilistic
prediction from the logistic regression duality. Hence, can check for all
grandmasters and then impute optimal $K$. Can also use Elo’s Maxwell-Boltzmann
distribution and see how far in the tail Magnus’ performance is.

Much discussion about whether the $K$-factor should be changed. John Nunn
(2009) asked for ”proof” of why it would be a good idea citing that the Elo
system has worked well for $50$ years. We agree that the Elo system has worked
perfectly well in chess and other sports Feng et al. (2016); Polson and Stern
(2015); Stern (1991), however, a subjective choice of $K=10$ has implications
for how high the current best player can achieve his goals. Our research
outlines a path to answer Nunn’s question and to influence policy.

## References

  * Aldous (2017) David Aldous.  Elo ratings and the sports model: A neglected topic in applied probability?  _Statistical science_ , 32(4):616–629, 2017.
  * Dangauthier et al. (2007) Pierre Dangauthier, Ralf Herbrich, Tom Minka, and Thore Graepel.  TrueSkill Through Time: Revisiting the History of Chess.  In _Advances in Neural Information Processing Systems_ , volume 20. Curran Associates, Inc., 2007.
  * Elo (1978) Arpad E Elo.  _The rating of chessplayers, past and present_.  Arco Pub., 1978.
  * Feng et al. (2016) Guanhao Feng, Nicholas G. Polson, and Jianeng Xu.  The Market for English Premier League (EPL) Odds.  _Journal of Quantitative Analysis in Sports_ , 12(4), January 2016.
  * Glickman and Jones (1999) Mark E Glickman and Albyn C Jones.  Rating the chess rating system.  _Chance_ , 12:21–28, 1999.
  * John Nunn (2009) John Nunn.  On the K-factor: Show me the proof!  https://en.chessbase.com/post/nunn-on-the-k-factor-show-me-the-proof-, April 2009\.
  * Polson and Stern (2015) Nicholas Polson and Hal Stern.  The implied volatility of a sports game.  _Journal of Quantitative Analysis in Sports_ , 11, January 2015.
  * Stern (1991) Hal Stern.  On the Probability of Winning a Football Game.  _The American Statistician_ , 45(3):179–183, 1991\.  ISSN 0003-1305.
  * Viswanath (2016) Ganesh Viswanath.  Elo Rating System: How underrated are the kids?, April 2016.
  * Williams et al. (2021) Leighton Vaughan Williams, Chunping Liu, Lerato Dixon, and Hannah Gerrard.  How well do elo-based ratings predict professional tennis matches?  _Journal of Quantitative Analysis in Sports_ , 17(2):91–105, 2021.


# Axion-assisted Resonance Oscillation Rescues the Dodelson-Widrow Mechanism

Shu-Yuan Guo<EMAIL_ADDRESS>Xuewen Liu<EMAIL_ADDRESS>Bin Zhu
<EMAIL_ADDRESS>Department of Physics, Yantai University, Yantai
264005, China

###### Abstract

The $\rm{keV}$ scale sterile neutrino was a qualified candidate for dark
matter particles in the Dodelson-Widrow mechanism. But the mixing angle,
needed to provide enough amount of dark matter, is in contradiction with the
astrophysical observations. To alleviate such tension, we introduce an
effective interaction, i.e.
$g_{a}(\phi/\Lambda)\partial_{\mu}a\overline{\nu_{\alpha}}\gamma^{\mu}\gamma_{5}\nu_{\alpha}$,
among Standard Model neutrino $\nu_{\alpha}$, axion $a$, and singlet $\phi$.
The axial-vector interaction form is determined by the axion shift symmetry,
and the singlet $\phi$ with dynamically varied vacuum expectation value is
introduced to reinforce the axial-vector coupling strength and evade the
stringent neutrino oscillation constraints. The effective potential generated
by the new interaction could cancel the SM counterpart, resulting in an
enhanced converting probability between SM neutrino and sterile neutrino.
Hence, the production rate of sterile neutrinos can be substantially enlarged
with smaller mixing compared to the DW mechanism.

##  I Introduction

The mystery surrounding dark matter (DM) production and neutrino mass
generation is intriguing. The seesaw mechanism, which introduces an extra
neutrino $\nu_{s}$, is the minimal and natural framework for solving these two
issues. After the electroweak symmetry breaking, the Dirac mass terms
$yv\bar{\nu}_{\alpha}\nu_{s}+\rm{h.c.}$ will inevitably result in a small
mixing between the active and sterile neutrinos,
$\nu_{4}=\cos\theta\nu_{s}+\sin\theta\nu_{\alpha}$, which makes the sterile
neutrino an ideal dark matter candidate in the Dodelson-Widrow (DW) mechanism
Dodelson and Widrow (1994). This mass mixing $\theta$ is too tiny in the early
universe to thermalize the sterile neutrinos with the Standard Model (SM)
bath. Thus there were no initial distributions for sterile neutrino. However,
by the collisions of active neutrinos with other SM particles, even a slight
amount of mixing can produce a sizable population of sterile neutrinos. The
production process will end once the SM neutrinos decoupled from the thermal
bath and the relic abundance of sterile neutrinos is attained with $\theta\sim
10^{-6}$ and $m_{\nu_{s}}\sim 100\mathrm{keV}$ Abazajian (2017). Unfortunately
the active-sterile neutrino mixing is already in tension with searches for DM
decaying into monochromatic X-rays in this minimal setup Watson  _et al._
(2012); Horiuchi  _et al._ (2014); Perez  _et al._ (2017); Dessert  _et al._
(2020); Ng  _et al._ (2019). Together with the bounds from the observation on
the dwarf galaxies Tremaine and Gunn (1979); Boyarsky  _et al._ (2009); Merle
_et al._ (2016); Abazajian (2017), the parameter space has been almost ruled
out.

Given these challenges, less minimal hypotheses that would increase the rate
of sterile neutrinos and reduce the restrictions imposed by astrophysical
observations have been proposed. The existence of a significant lepton-number
asymmetry in the early universe to enlarge the sterile neutrino production
rate resonantly is a common new physics approach Shi and Fuller (1999) to
alleviate the above tension. The problem with this approach is that it can not
explain why the lepton asymmetry is significantly larger than the baryon
asymmetry, and get a strong constraint from Big Bang Nucleosynthesis.
Including additional self-interaction between active neutrinos is another
elegant approach De Gouvêa  _et al._ (2020); Kelly  _et al._ (2020, 2021) to
erasing the tension. It helps efficiently produce sterile neutrino DM and
keeps neutrinos in thermal equilibrium with themselves for a longer time. The
neutrinophilic scalar or vector is responsible for mediating the self-
interaction, which results in a connection between production rate and dark
radiation $\Delta N_{\mathrm{eff}}$. A similar approach to generating the
appropriate relic abundance of $\nu_{s}$ is the freeze-in production via Higgs
singlet decay or vector boson decay processes Kusenko (2006); Petraki and
Kusenko (2008); Roland  _et al._ (2015); Drewes  _et al._ (2017); Shuve and
Yavin (2014). Ref. Berlin and Hooper (2017) proposed a solution by introducing
interactions with axion-like particles, which breaks the connection of the
mixing angles that responsible for sterile neutrino production in the early
universe and for the decay of sterile neutrino at late times. The size of the
mixing decreases as the energy density of axion-like particles diluted over
the cosmological expansion, thus avoiding the X-ray bounds. Even though all
above proposals can explain the sterile neutrino relic abundance, they go far
beyond the minimal configuration. The resonant oscillation seems to be the
best option if we insist on preserving the grace of the DW mechanism.

It is demonstrated that the presence of a primordial condensate of the dark
photon can alter the dispersion relation of the neutrinos in the early
universe Alonso-Álvarez and Cline (2021). A greater active neutrino self-
energy allows for a level crossover with its heavier sterile counterpart,
enabling resonant oscillations between the two states. We notice that the
underlying principle of the approach is that a dark photon originates from the
misalignment mechanism where it has to couple to the curvature R with extra
substantial, highly tuned couplings. But these couplings trigger perturbative
unitarity violations in longitudinal photon-graviton scattering at low
energiesArias  _et al._ (2012). This straightforward and elegant approach also
faces problems.

We point out that the axion can be used to generate resonant oscillation, as
it naturally originates from the misalignment mechanismLinde (1988); Preskill
_et al._ (1983). The minimal interaction between axion and active neutrino
$g_{a}\partial_{\mu}a\overline{\nu_{\alpha}}\gamma^{\mu}\gamma_{5}\nu_{\alpha}$
suffers from stringent constraints, where the magnitude of $g_{a}$ required
for resonance to exist contradicts the neutrino oscillation in the current
universe. In this letter, we propose a brand-new method to generate resonant
oscillations in the early universe and remove the neutrino oscillation
constraints in the late universe through the axion assisted DW mechanism
(denoted as aDW hereafter), with interaction
$g_{a}(\phi/\Lambda)\partial_{\mu}a\overline{\nu_{\alpha}}\gamma^{\mu}\gamma_{5}\nu_{\alpha}$.
The main idea is that the vacuum expectation value of the singlet $\phi$ is
not a constant but a function of temperature, where it becomes significant in
the early universe and vanishes in the late universe. In contrast to the
minimal axion assisted DW mechanism, the ensuing resonant oscillations enable
an observable relic abundance associated with significantly reduced portal
couplings $g_{a}$, evading the astrophysical constraints.

##  II Production of sterile neutrino

Sterile neutrino is among the simplest solutions which generate the nonzero
neutrino masses confirmed by oscillation experiments. The mixing between it
and SM neutrino will then naturally give birth to the initial abundance of
sterile neutrino. An accurate calculation of sterile neutrino production
requires solving the evolution equation in the matrix density formalism, while
in the collision-dominated period (which it’s indeed the case of the
production peak of sterile neutrino, i.e. $T\sim 0.1\leavevmode\nobreak\
\rm{GeV}$), the quasi-classical Boltzmann equation is a good approximation
Abazajian  _et al._ (2001),

$\displaystyle\left(\frac{\partial}{\partial t}-Hp\frac{\partial}{\partial
p}\right)f_{\nu_{s}}(p,t)$ $\displaystyle\simeq$
$\displaystyle\frac{\Gamma_{\alpha}(p)}{2}\langle
P_{m}(\nu_{\alpha}\leftrightarrow\nu_{s};p,t)\rangle f_{\nu_{\alpha}}(p,t),$

here $f_{\nu_{s}}$ and $f_{\nu_{\alpha}}$ stand for the distribution functions
of sterile neutrino and SM neutrino, $\Gamma_{\alpha}$ is the reaction rate
where SM neutrino takes part in. $\langle P_{m}\rangle$ is the averaged
converting probability between sterile neutrino and SM neutrino, expressed as
Volkas and Wong (2000); Lee  _et al._ (2000)

$\langle
P_{m}(\nu_{\alpha}\leftrightarrow\nu_{s};p,t)\rangle\simeq\frac{1}{2}\frac{\Delta^{2}\sin^{2}2\theta}{\Delta^{2}\sin^{2}2\theta+\Gamma_{\alpha}^{2}/4+(\Delta\cos
2\theta-V)^{2}},$ (2)

with $\Delta=m_{s}^{2}/2p$, the effective potential is shown as $V$, which
characterizes the effects of SM neutrinos passing through the thermal
background. The effective potential consists of contributions from SM
electroweak interactions and any other new physics. For DW mechanism, there is
only the SM electroweak contribution, which is approximated as
$V_{\rm{SM}}\simeq-14\pi/(45\alpha)(3-s_{w}^{2})s_{w}^{2}G_{F}^{2}pT^{4}\equiv-
C_{W}^{2}pT^{4}$ Notzold and Raffelt (1988); Alonso-Álvarez and Cline (2021)
for temperature between the mass of corresponding charged lepton and
electroweak phase transition. The $\alpha$ stands for the fine structure
constant, $s_{w}$ is the sine of Weinberg angle, and $G_{F}$ is the Fermi
constant.

To obtain the relic abundance of sterile neutrinos, one can derive the
equation for $r\equiv n_{s}/n_{\alpha}$, with $n_{i}$ denoting the number
density of sterile (SM) neutrino when $i=s$ $(\alpha)$. The resulting equation
reads

$\frac{dr}{d\ln a}=\frac{\gamma}{H}+r\frac{d\ln g_{s}^{\star}}{d\ln a}.$ (3)

Here entropy density conservation has been considered, and $g_{s}^{\ast}$ is
the degrees of freedom with respect to entropy density. $\gamma$ on the right
hand side is defined as

$\gamma\equiv\frac{1}{n_{a}}\int\frac{d^{3}p}{(2\pi)^{3}}\Gamma\langle
P_{m}\rangle f_{\nu_{a}}.$ (4)

For constant $g_{s}^{\ast}$, $\gamma/H$ stands for the relative number of
sterile neutrino, compared to SM neutrino, in each log-interval of $T$. Thus
in the simple DW mechanism there are two free parameters, i.e., the mixing
angle $\theta$ and the sterile neutrino mass $m_{s}$. The relic abundance of
sterile neutrino will be determined once these two parameters are fixed. In
Figure 1, we show the variation of $\gamma/H$ with temperature during a period
of $[10^{-3},10]\leavevmode\nobreak\ \rm{GeV}$ in the DW mechanism (the red
line). We choose the mixing angle as $\sin^{2}2\theta=10^{-12}$, and the
sterile neutrino mass is selected as $m_{s}=10\leavevmode\nobreak\ \rm{keV}$.
The production peak reaches $\gamma/H\sim 10^{-6}$ at temperature $\sim
0.2\leavevmode\nobreak\ \rm{GeV}$.

Such a framework is simple enough to produce sterile neutrino, just from SM
neutrino oscillation. However this mechanism is facing crises from various
experiments. As the sterile neutrino are much heavier than SM neutrino, decay
process of $\nu_{s}\to\nu_{\alpha}\gamma$ is inevitable, the decay rate is Pal
and Wolfenstein (1982)

$\Gamma(\nu_{s}\to\nu_{a}\gamma)\simeq\sin^{2}\theta\frac{9G_{F}^{2}\alpha
m_{s}^{5}}{2^{11}\pi^{4}}.$ (5)

Hence the observation of X-ray Watson  _et al._ (2012); Horiuchi  _et al._
(2014); Perez  _et al._ (2017); Dessert  _et al._ (2020); Ng  _et al._ (2019)
can put strong bounds on the mixing angle. Moreover, for fermionic DM, there
is a lower mass limit due to the requirement that the phase-space density of
DM does not exceed that of the degenerate Fermi gas Tremaine and Gunn (1979).
A bound of $m_{s}>2\leavevmode\nobreak\ \rm{keV}$ is set from the analysis of
DM phase space distribution in dwarf galaxies Boyarsky  _et al._ (2009); Merle
_et al._ (2016); Abazajian (2017). As can be seen in Figure 2, the X-ray
observation has excluded almost the whole parameter space in DW mechanism.
Considering the constraints of dwarf galaxies, there is no room for the DW
mechanism to survive.

##  III The axion assisted DW mechanism

We consider a new interaction of SM neutrino in an effective operator form

$-\mathcal{L}_{a}=g_{a}\left(\frac{\phi}{\Lambda}\right)\partial_{\mu}a\leavevmode\nobreak\
\overline{\nu_{\alpha}}\gamma^{\mu}\gamma_{5}\nu_{\alpha},$ (6)

including the axion field $a$ with mass $m_{a}$, and a singlet $\phi$. The
axial-vector interaction form is fixed by the axion shift symmetry. The
inclusion of $\phi$ is to avoid severe constraints of neutrino oscillation
observations. Scattering between axion and SM neutrino will modify the vacuum
dispersion relation of neutrino, hence generating a Mikheyev-Smirnov-
Wolfenstein matter effect Wolfenstein (1978); Mikheyev and Smirnov (1985);
Mikheev and Smirnov (1986). In later discussions we can see that, without the
inclusion of $\phi$, the coupling strength $g_{a}$ needed to give the correct
DM relic density is in strong contradiction with the oscillation observations.
The coupling strength matrix $g_{a}$ has a negative mass dimension $[M]^{-1}$.
$\Lambda$ is the ultra-violate cutoff scale. Due to the large occupation
number nature, the axion fields can be approximated in a classical wave form
of $a(t,x)=a_{0}\cos(m_{a}t-\vec{p}\cdot\vec{x})$, with
$a_{0}=\sqrt{2\rho_{a}}/m_{a}$ and $\rho_{a}$ being the energy density of
axions. In our consideration, we treat sterile neutrino as a dominant
component of dark matter, while for the axion we set $\rho_{a}$ as $1\%$ of
the dark matter energy density.

Containing the singlet $\phi$ in the effective interaction is an important
ingredient in our setup. $\phi$ is a morphon field Croon  _et al._ (2022)
which develops a dynamically varied vacuum expectation value by interacting
with Ricci scalar $R$ at the early stage of the universe,

$V(\phi)=\frac{1}{2}m_{\phi}^{2}\phi^{2}-\frac{1}{2}\xi
R\phi^{2}+\frac{1}{4}\lambda_{\phi}\phi^{4}.$ (7)

The coupling constant $\xi$ quantifies the interaction strength between the
morphon field and curvature scalar. One common choice is $\xi=1/6$, which
comes from the requirement of conformal invariance Sonego and Faraoni (1993);
Dowker (1971). A comprehensive analysis of the temperature dependence of
Ricci-scalar is carried out in Caldwell and Gubser (2013) for the time range
from inflationary epoch to the present day. In our calculation, we consider
the evolution of Ricci scalar from $T\sim 1$ MeV to $T\sim 10$ GeV. Thus, at
high temperatures the scalar field develops a nonzero vacuum expectation
value,

$\displaystyle
v_{\phi}(T)=\left\\{\begin{array}[]{cl}\sqrt{\frac{1}{\lambda_{\phi}}\left(\xi
R(T)-m_{\phi}^{2}\right)}\,,&\xi R(T)>m_{\phi}^{2}\\\\[11.0pt]
0&\text{otherwise}\,.\end{array}\right.$ (10)

The scalar vacuum expectation value decreases with $T$ until it vanishes at
low temperatures. This is the main feature that it helps generate resonant
oscillations in the early universe and evade constraints in the late universe,
such as the neutrino oscillation bounds.

Now the additional contribution to the effective potential of SM neutrino can
be obtained from the new introduced effective interaction in Equation 6, in
the first-order approximation as

$V_{a}=g_{a}\frac{v_{\phi}(T)}{\Lambda}\left(-\partial_{0}a+\vec{\partial}a\cdot\vec{p}_{\nu}/|\vec{p}_{\nu}|\right).$
(11)

The derivative on space part could be neglected as the velocity of axion in
the Milky Way is of $\mathcal{O}(10^{-3})c$. Then the effective potential is
further approximated as

$V_{a}\approx g_{a}\frac{v_{\phi}(T)}{\Lambda}\sqrt{2\rho_{a}}\sin(m_{a}t),$
(12)

where the plane-wave form of the axion field has been taken.

We have mentioned above that the SM weak interactions contribute a negative
thermal potential $V_{\rm{SM}}=-C_{W}^{2}pT^{4}$. While the new contribution
could be either positive or negative, depending on the axion mass and time (or
temperature). Therefore, in principle there exists cancellation between the
contributions from SM and the new interaction. It would substantially enlarge
the averaged converting probability $\langle P_{m}\rangle$, and resonantly
produce the sterile neutrino dark matter. As a result, the mixing angle could
be much smaller than that in the original DW mechanism. This can help to
release the strong tension against the X-ray observations. The cancellation
condition is explicitly written as

$\Delta\cos 2\theta-V_{\rm{SM}}-V_{a}=0.$ (13)

We firstly give a rough estimation to show the incapability of the minimal
interaction
$g_{a}\partial_{\mu}a\overline{\nu_{\alpha}}\gamma^{\mu}\gamma_{5}\nu_{\alpha}$.
The former two terms in Equation 13 have different dependencies on
temperature, i.e. $\Delta\cos 2\theta\propto T^{-1}$ and $V_{\rm{SM}}\propto
T^{5}$ with setting $p=T$ for simplicity. Thus the lowest value of the sum is
estimated as $\sim 4.1(m_{s}/10\leavevmode\nobreak\
\rm{keV})^{5/3}\leavevmode\nobreak\ \rm{eV}$. The new potential is in a form
of $g_{a}\sqrt{2\rho_{a}}\sin(m_{a}t)$, without including the singlet $\phi$.
We find that the coupling strength $g_{a}$ needed for $V_{a}$, to have a
considerable cancellation with the former two terms, is much above the upper
limits from the oscillation experiments. For instance, $g_{a}$ should be no
less than $\sim 10^{4}\leavevmode\nobreak\ \rm{eV}^{-1}$ for the configuration
we considered about dark matter components. However it’s much larger than the
upper bounds from oscillation experiments Huang and Nath (2018). Taking
IceCube as an example, the upper bounds of $g_{a}$ lie in $[8\times
10^{-9},5\times 10^{-6}]\leavevmode\nobreak\ \rm{eV}^{-1}$ for the axion
masses in $[10^{-22},10^{-2}]\leavevmode\nobreak\ \rm{eV}$ Huang and Nath
(2018). However, the oscillation limitations on $g_{a}$ would not appear in
our case, since the introduced $\phi$ has a temperature-dependent vacuum
expectation value. The selected $m_{\phi}$, i.e.
$m_{\phi}=10^{-15}\leavevmode\nobreak\ \rm{eV}$, would make the vacuum
expectation value of $\phi$ to be vanished at around $5\leavevmode\nobreak\
\rm{MeV}$. Therefore the introduced effective interaction has no effects on
the current oscillation observations. Another simple interaction form between
SM neutrino and scalar is $\phi\bar{\nu}\nu$. From Ref. Ge and Parke (2019) we
know it only results in a neutrino mass correction, and is irrelevant to the
resonance oscillation.

Figure 1: Production of sterile neutrino in the DW and aDW mechanisms. The red
and green lines stand for the variation of $\gamma/H$ with temperature $T$ in
DW and aDW mechanisms, respectively. The parameters are set as
$\sin^{2}2\theta=10^{-12}$, $m_{s}=10\leavevmode\nobreak\ \rm{keV}$,
$m_{\phi}=10^{-15}\leavevmode\nobreak\ \rm{eV}$, $\xi=1/6$,
$m_{a}=10^{-10}\leavevmode\nobreak\ \rm{eV}$, and a combination of
$g_{a}/(\sqrt{\lambda_{\phi}}\Lambda)$ is set to $10^{14}\leavevmode\nobreak\
\rm{eV}^{-2}$. The blue and purple lines show relative size of the SM
potential and the new potential to $\Delta\cos 2\theta$.

To reveal the resonance effects, we depict the variation of the two
potentials, as functions of temperature in Figure 1. In practice, we normalize
these two terms by $\Delta\cos 2\theta$, and momentum is set to $p=T$. The
purple line indicates the relative magnitude of SM potential, and the blue
line stands for the relative magnitude of the new potential $V_{a}$. The black
dashed line represents the normalization factor $\Delta\cos 2\theta$.

Without considering large deviations between the sterile neutrino momentum and
temperature, it can be seen from Figure 1 that $\Delta\cos 2\theta$ is
dominant in the region of $T\lesssim 0.2\leavevmode\nobreak\ \rm{GeV}$, while
the SM potential $|V_{\rm{SM}}|$ is dominant in the region of $T\gtrsim
0.8\leavevmode\nobreak\ \rm{GeV}$. The new potential $V_{a}$ oscillates with
temperature and terminates at $T\sim 5\leavevmode\nobreak\ \rm{MeV}$ because
the vacuum expectation value vanishes at lower temperatures. As the
temperature increases, the oscillation frequency becomes smaller since the
relation between time and temperature is $t\propto T^{-2}$. The maximum value
that can be reached by the oscillated $V_{a}$ gradually approaches and exceeds
$\Delta\cos 2\theta$. Therefore, as the temperature increases, the
cancellation between the three terms is expected to be more pronounced, until
$|V_{\rm{SM}}|$ increases so that SM contribution becomes dominant.

In Figure 1, we also show the dependence of $\gamma/H$ on the temperature for
both DW (red line) and aDW (green line) mechanisms. The aDW scenario will go
back to DW for the high and low temperatures. For the intermediate
temperature, the cancellation leads to an enhancement on $\gamma/H$ in the aDW
mechanism, whose magnitude is up to $\sim 6$ orders higher in the comparison
of the two scenarios. Thus the production of sterile neutrino dark matter is
expected to be greatly enhanced for aDW.

Base on the resonance effects, the sterile neutrino could be produced with a
much smaller mixing angle. In Figure 2, we show the correct relic density in
($m_{s}$, $\sin^{2}2\theta$) plane. The green line stands for the available
parameter space in aDW scenario, and the red line represents the traditional
DW mechanism. The most stringent constraints come from the X-ray observation
and phase space distribution of dwarf galaxies, which are shown by the brown
and purple excluded regions. The most parameter space that predicted by DW
mechanism have been excluded. In our scenario containing the new interactions,
the values of the mixing angle can be much lower comparatively, meanwhile
providing the correct amount of relic abundance. For the parameters we have
selected, sterile neutrinos with mass in $[2,12]\leavevmode\nobreak\ \rm{keV}$
can satisfy all the astrophysical constraints. In Figure 2, we also show the
detecting capabilities of two forthcoming experiments, i.e. Athena Neronov and
Malyshev (2016) and eXTP Zhong  _et al._ (2020); Malyshev  _et al._ (2020). We
see that they could detect some parameter region that is more stringent than
the X-ray constraints, and can be used to test the feasibility of our
scenario.

Figure 2: Correct relic density of sterile neutrino DM as functions of mass
and mixing angle. The red line represents the proper parameter space in the
traditional Dodelson-Widrow mechanism (DW), while the green line stands for
that in the aDW mechanism (aDW). The choices of model parameters are
$m_{\phi}=10^{-15}\leavevmode\nobreak\ \rm{eV}$, $\xi=1/6$,
$m_{a}=10^{-10}\leavevmode\nobreak\ \rm{eV}$, and
$g_{a}/(\sqrt{\lambda_{\phi}}\Lambda)=10^{14}\leavevmode\nobreak\
\rm{eV}^{-2}$. Constraints from X-ray observation are depicted in brown
region, and lower bound on mass of sterile neutrino from dwarf galaxies is
shown by purple region. The detecting capabilities of two forthcoming
experiments, i.e. Athena and eXTP, are showed in orange and blue dashed lines
respectively.

##  IV Conclusion

Sterile neutrino dark matter, generated via the mixing with SM neutrino in DW
mechanism, suffers severe conflict between its relic density and the
astrophysical observations. X-ray constraints excluded the most parameter
space of this scenario. To release this tension, we have introduced an
effective axial-vector interaction between axion, SM neutrino, and a singlet
in this work. The new interactions bring a new contribution to SM neutrino
effective potential, then the cancellation occurs between this new
contribution and SM one. The sterile neutrino thus can be resonantly produced
with a significant enhancement factor. Around the typical DW-peak temperature,
the resonance enhancement could produce sterile neutrino with an amount of
several orders higher than DW mechanism. Such that the much smaller mixing
angles can satisfy the stringent astrophysical constraints and provide the
right abundance of the relic DM at the same time. An inevitable consequence of
new interaction with neutrino is the modification of the neutrino oscillation.
Since the vacuum expectation value of the new scalar dynamically varies in
temperature, the new effective potential would vanish at the current
temperature. So it’s safe to ignore it in our case.

## Acknowledgements

We thank Guo-yuan Huang for helpful discussions. This work was supported by
the National Natural Science Foundation of China under Grants Nos. 12005180,
12275232, by the Natural Science Foundation of Shandong Province under Grants
No. ZR2020QA083, and by the Project of Shandong Province Higher Educational
Science and Technology Program under Grants No. 2019KJJ007.

## References

  * Dodelson and Widrow (1994) S. Dodelson and L. M. Widrow, Phys. Rev. Lett. 72, 17 (1994), arXiv:hep-ph/9303287 .
  * Abazajian (2017) K. N. Abazajian, Phys. Rept. 711-712, 1 (2017), arXiv:1705.01837 [hep-ph] .
  * Watson  _et al._ (2012) C. R. Watson, Z.-Y. Li, and N. K. Polley, JCAP 03, 018 (2012), arXiv:1111.4217 [astro-ph.CO] .
  * Horiuchi  _et al._ (2014) S. Horiuchi, P. J. Humphrey, J. Onorbe, K. N. Abazajian, M. Kaplinghat, and S. Garrison-Kimmel, Phys. Rev. D 89, 025017 (2014), arXiv:1311.0282 [astro-ph.CO] .
  * Perez  _et al._ (2017) K. Perez, K. C. Y. Ng, J. F. Beacom, C. Hersh, S. Horiuchi, and R. Krivonos, Phys. Rev. D 95, 123002 (2017), arXiv:1609.00667 [astro-ph.HE] .
  * Dessert  _et al._ (2020) C. Dessert, N. L. Rodd, and B. R. Safdi, Science 367, 1465 (2020), arXiv:1812.06976 [astro-ph.CO] .
  * Ng  _et al._ (2019) K. C. Y. Ng, B. M. Roach, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos, and D. R. Wik, Phys. Rev. D 99, 083005 (2019), arXiv:1901.01262 [astro-ph.HE] .
  * Tremaine and Gunn (1979) S. Tremaine and J. E. Gunn, Phys. Rev. Lett. 42, 407 (1979).
  * Boyarsky  _et al._ (2009) A. Boyarsky, O. Ruchayskiy, and D. Iakubovskyi, JCAP 03, 005 (2009), arXiv:0808.3902 [hep-ph] .
  * Merle  _et al._ (2016) A. Merle, A. Schneider, and M. Totzauer, JCAP 04, 003 (2016), arXiv:1512.05369 [hep-ph] .
  * Shi and Fuller (1999) X.-D. Shi and G. M. Fuller, Phys. Rev. Lett. 82, 2832 (1999), arXiv:astro-ph/9810076 .
  * De Gouvêa  _et al._ (2020) A. De Gouvêa, M. Sen, W. Tangarife, and Y. Zhang, Phys. Rev. Lett. 124, 081802 (2020), arXiv:1910.04901 [hep-ph] .
  * Kelly  _et al._ (2020) K. J. Kelly, M. Sen, W. Tangarife, and Y. Zhang, Phys. Rev. D 101, 115031 (2020), arXiv:2005.03681 [hep-ph] .
  * Kelly  _et al._ (2021) K. J. Kelly, M. Sen, and Y. Zhang, Phys. Rev. Lett. 127, 041101 (2021), arXiv:2011.02487 [hep-ph] .
  * Kusenko (2006) A. Kusenko, Phys. Rev. Lett. 97, 241301 (2006), arXiv:hep-ph/0609081 .
  * Petraki and Kusenko (2008) K. Petraki and A. Kusenko, Phys. Rev. D 77, 065014 (2008), arXiv:0711.4646 [hep-ph] .
  * Roland  _et al._ (2015) S. B. Roland, B. Shakya, and J. D. Wells, Phys. Rev. D 92, 113009 (2015), arXiv:1412.4791 [hep-ph] .
  * Drewes  _et al._ (2017) M. Drewes _et al._ , JCAP 01, 025 (2017), arXiv:1602.04816 [hep-ph] .
  * Shuve and Yavin (2014) B. Shuve and I. Yavin, Phys. Rev. D 89, 113004 (2014), arXiv:1403.2727 [hep-ph] .
  * Berlin and Hooper (2017) A. Berlin and D. Hooper, Phys. Rev. D 95, 075017 (2017), arXiv:1610.03849 [hep-ph] .
  * Alonso-Álvarez and Cline (2021) G. Alonso-Álvarez and J. M. Cline, JCAP 10, 041 (2021), arXiv:2107.07524 [hep-ph] .
  * Arias  _et al._ (2012) P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, JCAP 06, 013 (2012), arXiv:1201.5902 [hep-ph] .
  * Linde (1988) A. D. Linde, Phys. Lett. B 201, 437 (1988).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abazajian  _et al._ (2001) K. Abazajian, G. M. Fuller, and M. Patel, Phys. Rev. D 64, 023501 (2001), arXiv:astro-ph/0101524 .
  * Volkas and Wong (2000) R. R. Volkas and Y. Y. Y. Wong, Phys. Rev. D 62, 093024 (2000), arXiv:hep-ph/0007185 .
  * Lee  _et al._ (2000) K. S. M. Lee, R. R. Volkas, and Y. Y. Y. Wong, Phys. Rev. D 62, 093025 (2000), arXiv:hep-ph/0007186 .
  * Notzold and Raffelt (1988) D. Notzold and G. Raffelt, Nucl. Phys. B 307, 924 (1988).
  * Pal and Wolfenstein (1982) P. B. Pal and L. Wolfenstein, Phys. Rev. D 25, 766 (1982).
  * Wolfenstein (1978) L. Wolfenstein, Phys. Rev. D 17, 2369 (1978).
  * Mikheyev and Smirnov (1985) S. P. Mikheyev and A. Y. Smirnov, Sov. J. Nucl. Phys. 42, 913 (1985).
  * Mikheev and Smirnov (1986) S. P. Mikheev and A. Y. Smirnov, Nuovo Cim. C 9, 17 (1986).
  * Croon  _et al._ (2022) D. Croon, G. Elor, R. Houtz, H. Murayama, and G. White, Phys. Rev. D 105, L061303 (2022), arXiv:2012.15284 [hep-ph] .
  * Sonego and Faraoni (1993) S. Sonego and V. Faraoni, Class. Quant. Grav. 10, 1185 (1993).
  * Dowker (1971) J. S. Dowker, Annals Phys. 62, 361 (1971).
  * Caldwell and Gubser (2013) R. R. Caldwell and S. S. Gubser, Phys. Rev. D 87, 063523 (2013), arXiv:1302.1201 [astro-ph.CO] .
  * Huang and Nath (2018) G.-Y. Huang and N. Nath, Eur. Phys. J. C 78, 922 (2018), arXiv:1809.01111 [hep-ph] .
  * Ge and Parke (2019) S.-F. Ge and S. J. Parke, Phys. Rev. Lett. 122, 211801 (2019), arXiv:1812.08376 [hep-ph] .
  * Neronov and Malyshev (2016) A. Neronov and D. Malyshev, Phys. Rev. D 93, 063518 (2016), arXiv:1509.02758 [astro-ph.HE] .
  * Zhong  _et al._ (2020) D. Zhong, M. Valli, and K. N. Abazajian, Phys. Rev. D 102, 083008 (2020), arXiv:2003.00148 [astro-ph.HE] .
  * Malyshev  _et al._ (2020) D. Malyshev, C. Thorpe-Morgan, A. Santangelo, J. Jochum, and S.-N. Zhang, Phys. Rev. D 101, 123009 (2020), arXiv:2001.07014 [astro-ph.HE] .