# Signatures of Ultralight Dark Matter in Neutrino Oscillation Experiments

Abhish Dev<EMAIL_ADDRESS>Maryland Center for Fundamental Physics, Department
of Physics, University of Maryland, College Park,MD 20742-4111 USA  Pedro A.N.
Machado<EMAIL_ADDRESS>Theoretical Physics Department, Fermilab, P.O. Box
500, Batavia, IL 60510, USA  Pablo Martínez-Miravé<EMAIL_ADDRESS>Departament de Física Teòrica, Universitat de València, 46100, Burjassot,
Spain  Institut de Física Corpuscular CSIC-Universitat de València, 46980
Paterna, Spain

###### Abstract

We study how neutrino oscillations could probe the existence of ultralight
bosonic dark matter. Three distinct signatures on neutrino oscillations are
identified, depending on the mass of the dark matter and the specific
experimental setup. These are time modulation signals, oscillation probability
distortions due to fast modulations, and fast varying matter effects. We
provide all the necessary information to perform a bottom-up, model-
independent experimental analysis to probe such scenarios. Using the future
DUNE experiment as an example, we estimate its sensitivity to ultralight
scalar dark matter. Our results could be easily used by any other oscillation
experiment.

††preprint: FERMILAB-PUB-20-260-T

##  I Introduction

The nature of dark matter remains one of the greatest unknowns in particle
physics. Although there is plenty of data that corroborates the existence of
dark matter, its nature, mass, and interactions with the standard model are
widely unknown Tanabashi  _et al._ (2018). One possibility is when the mass of
the dark matter (DM) is much below the electronvolt scale. For these masses,
the dark matter De Broglie wavelength becomes macroscopic and could even
compare to the size of solar systems and galaxies. This requires the
ultralight dark matter to be heavier than $10^{-22}$ eV. Due to Fermi-Dirac
statistics, phase space considerations today can be used to better constrain
the mass of fermionic dark matter to be above 100 eV to the keV scale (also
known as the Tremaine-Gunn bound) Tremaine and Gunn (1979); Madsen (1990,
1991); Boyarsky  _et al._ (2009).

Ultralight bosonic fields are common and can arise as pseudo-Nambu Goldstones
from the spontaneous symmetry breaking of an approximate global symmetry such
as the QCD axion Peccei and Quinn (1977a, b); Weinberg (1978); Wilczek (1978).
In addition to being ubiquitous, they could also provide dynamical solutions
to finetuning problems. For example, the QCD axion makes the QCD theta angle
$\theta_{QCD}$ dynamical and relaxes it to zero Vafa and Witten (1984). These
fields can also constitute ultralight dark matter and result in unique fuzzy
dark matter phenomenology Preskill  _et al._ (1983); Abbott and Sikivie
(1983); Dine and Fischler (1983); Hu  _et al._ (2000); Arias  _et al._ (2012);
Li  _et al._ (2014); Hui  _et al._ (2017); Menci  _et al._ (2017); Diez-
Tejedor and Marsh (2017); Irsic  _et al._ (2017); Visinelli (2017); Kumar
Poddar  _et al._ (2020), which offers possible solutions to three small-scale
cosmological puzzles, namely, cusp-vs-core, missing satellite and too-big-to-
fail problems Moore (1994); Flores and Primack (1994); Navarro  _et al._
(1997); Klypin  _et al._ (1999). All these puzzles are based on comparisons
between simulations and observed cosmological data. Compared to data,
simulations predict high rotation curves towards galaxy centers due to the
large DM density therein, too many satellite galaxies, and too many visible
dwarf galaxies due to the presence of several DM subhalos. Fuzzy dark matter,
being delocalized due to its large De Broglie wavelength, would be more
susceptible to tidal disruptions and lead to broader halos possibly solving
the three aforementioned problems.

Regardless of the axion solution to the strong $CP$ problem or the Fuzzy DM
regime, nonzero abundance of ultralight bosonic DM could lead to a large
occupation number for this field, which would behave as a vacuum expectation
value (VEV). The idea is that the DM field undergoes a phase transition in the
early universe and starts oscillating around its new VEV. Hubble friction
damps these oscillations, but nevertheless, they may be still lead to
observable effects in late times. This oscillating VEV may lead to a curious
phenomenon: time variation of what we define as fundamental constants. Such
effects have been searched for in many different setups, from resonant
microwave cavities to atomic clocks (for a review, see e.g. Ref. Battaglieri
_et al._ (2017)). Here we are interested in the possibility that neutrinos
provide a portal to such ultralight fields. This is motivated by two facts.
First, neutrinos are much lighter than other fermions in the standard model,
and thus even a small time-dependent VEV could lead to a relatively large
impact on neutrino masses and mixings. Second, as the mechanism of neutrino
masses remains unknown, this sector may offer a more natural connection to
ultralight physics beyond the standard model. Although we use Fuzzy DM and
neutrino masses as a motivation, we are indeed more concerned with general
experimental signatures, and thus we adopt a _bottom-up_ approach in this
work. Some of these signatures were explored in previous literature Berlin
(2016); Capozzi  _et al._ (2017); Krnjaic  _et al._ (2018); Brdar  _et al._
(2018); Farzan and Palomares-Ruiz (2019); Liao  _et al._ (2018); Ge and Parke
(2019); Huang and Nath (2018); Ge and Murayama (2019); Cline (2020); Choi  _et
al._ (2020); Farzan (2019). Here we will provide a more detailed analysis,
with special emphasis on aspects directly connected to experimental searches,
and highlighting the different regimes relevant to neutrino oscillation
phenomenology. To be concrete, we will focus on the case of ultralight scalar
fields and defer ultralight vector fields to future analysis. We will use the
DUNE experiment as a case study, even if our conclusions will be valid for any
neutrino oscillation experiment.

##  II Theoretical guidance to experimental searches

First, let us review the theoretical aspects of ultralight scalar fields
coupled to neutrinos. We choose to work in an effective theory framework below
electroweak symmetry breaking. The effective Lagrangian that describes the
system is

$\mathcal{L}_{\rm
eff}=-m_{\nu}\left(1+y\frac{\phi}{\Lambda}\right)\bar{\nu}\nu+{\rm h.c.},$ (1)

where flavor indices are implicit in both $m_{\nu}$ and $y$, $\phi$ is the
ultralight scalar field, $m_{\nu}$ is the neutrino mass matrix and $\Lambda$
is a heavy mass scale 111We do not distinguish Dirac and Majorana neutrinos
here, as the oscillation phenomenology of these cases are identical.. We will
not worry where this effective Lagrangian arrives from and leave
considerations regarding ultraviolet completions, perhaps, to future work.

The local field value can be expressed as

$\phi(x,t)\simeq\frac{\sqrt{2\rho_{\phi}}}{m_{\phi}}\sin\left[m_{\phi}(t-\vec{v}\cdot\vec{x})\right],$
(2)

where the local density $\rho_{\phi}$ should not exceed the local DM density
$\rho_{\rm DM}=0.3~{}{\rm GeV}/{\rm cm}^{3}$, $m_{\phi}$ is the mass of $\phi$
and $v\sim 10^{-3}c$ is the virial velocity. Note that the space-dependent
phase $m_{\phi}\vec{v}\cdot\vec{x}$ of the local field value is much smaller
than $m_{\phi}t$, and thus will be neglected henceforth. The neutrino mass
matrix will therefore receive a contribution from $\phi$ given by

$\delta m_{\nu}=m_{\nu}y\frac{\sqrt{2\rho_{\phi}}}{\Lambda
m_{\phi}}\sin(m_{\phi}t).$ (3)

Without a flavor model (see e.g. Ref. Ding and Feruglio (2020)), the Yukawa
couplings $y$ can have any structure in flavor space, and thus the modulations
of neutrino masses and mixing can bear any correlation. Nevertheless, it is
useful to focus on two simple scenarios: modulation of mass splittings and
modulation of mixing angles.

Besides, the mass of $\phi$ defines the modulation period via

$\tau_{\phi}\equiv\frac{2\pi\hbar}{m_{\phi}}=0.41\left(\frac{10^{-14}~{}{\rm
eV}}{m_{\phi}}\right)~{}{\rm seconds}.$ (4)

In a given experimental setup, there are three characteristic time scales: the
neutrino time of flight $\tau_{\nu}=L/c=3.4(L/1000~{}{\rm km})~{}{\rm msec}$,
the time between two detected events $\tau_{\rm evt}$ (the inverse of the
ratio of events) and the lifetime of the experiment $\tau_{\rm exp}$. From
these, we can identify the following three different regimes for a given
experimental setup.

  * •

Time modulation ($\tau_{\rm evt}\lesssim\tau_{\phi}\ll\tau_{\rm exp}$). When
the period of modulation of $\phi$ is of the same order as the experiment
total run time, a temporal variation of the neutrino signal may be observed.
This is true for the modulation of angles and mass splittings. Experiments
with large statistics and high event rates will be sensitive to time
modulation periods much smaller than the lifetime of the experiments (for
instance, searches for modulations in solar neutrino fluxes are sensitive to
periods ranging from 10 minutes to 10 years Yoo  _et al._ (2003); Aharmim  _et
al._ (2005, 2010)).

  * •

Averaged distorted neutrino oscillations
($\tau_{\nu}\ll\tau_{\phi}\ll\tau_{\rm exp}$). Even when the rate of change of
neutrino oscillation parameters is too fast to be observed as a modulating
signal, the time average oscillation probability may be distorted by such
effects and deviate from the standard model scenario. While the averaging of a
modulating mixing angle can be mapped onto standard oscillations (with some
inferred value of this mixing angle), the averaging of a mass splitting can
lead to distorted neutrino oscillations (DiNOs), smearing out the probability
similarly to an energy resolution smearing Krnjaic  _et al._ (2018). This
regime covers a large range of scalar masses and can be easily searched for in
oscillation experiments, as it boils down to a simple novel oscillation
effect.

  * •

Dynamical distorted neutrino oscillations ($\tau_{\phi}\sim\tau_{\nu}$). As
the modulating period of $\phi$ gets closer to the neutrino time of flight,
the changes in oscillation parameters need to be treated at the Hamiltonian
level and can be modeled by a modified matter effect Brdar  _et al._ (2018).
This matter potential is time-dependent, and thus it changes as the neutrino
propagates towards the detector. When the variation of the matter potential is
too slow compared to the neutrino time of flight, dynamical DiNOs map onto
average DiNOs. In the opposite case, when the variations are too fast compared
to the neutrino time of flight, they cannot be observed and thus one recovers
standard oscillations.

In the next section, we will discuss the phenomenology of each regime in
general terms and highlight their phenomenology in oscillation experiments in
more detail.

##  III General considerations on ultralight scalars and neutrino oscillations

Before going into a more technical discussion, we provide a few insights on
the effects of ultralight scalar fields in neutrino oscillations. To do so, we
will analyze the changes in neutrino oscillation probability in each of the
regimes discussed above, using the DUNE experimental setup as a case study.
For simplicity, we will consider a single parameter modulating at a time.
Modulation of mixing angles are assumed to be of the form

$\theta_{ij}(t)=\theta_{ij}+\eta\sin(m_{\phi}t),$ (5)

where $\theta_{ij}$ represents the undistorted value of the mixing angle and
$\eta$ is the amplitude of modulation of $\phi$. Mass splittings, on the other
hand, are assumed to modulate like (for $\eta\ll 1$)

$\Delta m^{2}_{ij}(t)\equiv m_{i}^{2}(t)-m_{j}^{2}(t)\simeq\Delta
m^{2}_{ij}\left[1+2\eta\sin(m_{\phi}t)\right],$ (6)

where, similarly to above, $\Delta m^{2}_{ij}$ represents the undistorted
value of the mass splitting.

###  III.1 Time modulation

The phenomenology in the time modulation regime is very intuitive: mixing
angles or mass splittings are modulating on time. The oscillation probability
now depends on time via

$P_{\alpha\beta}\equiv
P(\nu_{\alpha}\to\nu_{\beta},t)=P\left[\nu_{\alpha}\to\nu_{\beta};\\{\theta_{ij}(t),\Delta
m^{2}_{ij}(t)\\}\right].$ (7)

For example, in the case of modulating angles, the oscillation probability for
$\nu_{\mu}\to\nu_{\mu}$ disappearance in vacuum, in a simplified two neutrino
framework, would read

$\displaystyle P_{\mu\mu}^{\rm angle}\simeq
1-\sin^{2}\left[2\theta(t)\right]\sin^{2}\left(\frac{\Delta
m^{2}L}{4E}\right)=1-\sin^{2}\left[2\theta+2\eta\sin(m_{\phi}t)\right]\sin^{2}\left(\frac{\Delta
m^{2}L}{4E}\right),$ (8)

where $L$ is the baseline (1300 km for DUNE) and $E$ is the neutrino energy.
Notice that the oscillation probability displays a time modulation via the
$\sin(m_{\phi}t)$ term.

Figure 1: The difference in $\nu_{\mu}$ disappearance oscillation probability
$\Delta P(t)\equiv P_{\mu\mu}(t)-P_{\mu\mu}(0)$ at DUNE for a modulating
$\theta_{23}$ assuming $\eta=0.05$ for several energy bins. We have assumed
normal mass ordering and the best fit values of oscillation parameters from
Ref. de Salas  _et al._ (2018), namely $\Delta m^{2}_{31}=2.5\times 10^{-3}$
eV2, and $\sin^{2}\theta_{23}=0.55$.

In Fig. 1 we show change in $\nu_{\mu}$ disappearance oscillation probability
$\Delta P(t)\equiv P_{\mu\mu}(t)-P_{\mu\mu}(0)$ at DUNE for a modulating
$\theta_{23}$ assuming $\eta=0.05$ for several energy bins. We have assumed
normal mass ordering and the best fit values of oscillation parameters from de
Salas  _et al._ (2018), namely $\Delta m^{2}_{31}=2.5\times 10^{-3}$ eV2, and
$\sin^{2}\theta_{23}=0.55$. There are three important effects to be noted
here. To understand those, we expand $\Delta P(t)$ for small $\eta$ as

$\Delta
P(t)\simeq-2\left[\eta\sin(4\theta)\sin(m_{\phi}t)+2\eta^{2}\cos(4\theta)\sin^{2}(m_{\phi}t)\right]\sin^{2}\left(\frac{\Delta
m^{2}L}{4E}\right).$ (9)

First, the modulation phases across different energies are fully correlated,
as the oscillation probability in all energy bins goes up and down at the same
time. Second, $\Delta P(t)$ depends on $\eta$ and the mixing angle itself. And
finally, when $\sin(4\theta)$ is near zero (e.g. for $\theta_{23}$ near
maximal mixing), we observe a Jacobian effect that shrinks the oscillation
amplitude, as the leading $\eta$ term shrinks and the $\eta^{2}$ term becomes
dominant. In Fig. 1, as the modulation develops, $\sin^{2}[2\theta_{23}(t)]$
goes further away from maximal, making $\Delta P(t)$ larger, and we observe
the larger amplitude. As the phase evolves, $\sin^{2}[2\theta_{23}(t)]$ gets
maximal, then slightly below maximal, maximal again, and back to its original
value. This corresponds to the region to the right in the subplots
($t/\tau_{\phi}>0.5$). If we had chosen to show a modulation of $\theta_{12}$
instead (e.g. in the JUNO experiment An  _et al._ (2016)), the change in
oscillation probability would have been more symmetric.

If instead, the mass splittings modulate on time, the simplified time-
dependent oscillation probability would read

$\displaystyle P_{\mu\mu}^{\rm mass}\simeq
1-\sin^{2}(2\theta)\sin^{2}\left[\frac{\Delta
m^{2}(t)L}{4E}\right]=1-\sin^{2}(2\theta)\sin^{2}\left\\{\left(\frac{\Delta
m^{2}L}{4E}\right)\left[1+2\eta\sin(m_{\phi}t)\right]\right\\}.$

We can see that what modulates is not the amplitude of the oscillation
probability, but rather the position of the minimum. In Fig. 2, we show again
the change in $\nu_{\mu}$ disappearance oscillation probability $\Delta
P(t)\equiv P_{\mu\mu}(t)-P_{\mu\mu}(0)$ at DUNE, but now for a modulating
$\Delta m^{2}_{31}$ assuming $\eta=0.05$ for several energy bins. Again, we
have assumed normal mass ordering and the best fit values of oscillation
parameters from de Salas  _et al._ (2018), namely $\Delta m^{2}_{31}=2.5\times
10^{-3}$ eV2, and $\sin^{2}\theta_{23}=0.55$. The most distinctive features,
in this case, are the correlations and anti-correlations across multiple
energies. This is easy to understand, as when the minimum or other regions of
small probability enter/leave a given energy bin, the average probability
lowers/raises in that bin. Quantitatively, we can see these correlations in
energy by expanding $\Delta P(t)$ around $\eta\to 0$,

$\Delta
P(t)\simeq-2\sin^{2}(2\theta)\sin(m_{\phi}t)\left[\eta\xi\sin(2\xi)+2\eta^{2}\xi^{2}\cos(2\xi)\sin(m_{\phi}t)\right],$
(10)

where we have defined $\xi\equiv\Delta m^{2}L/4E$. Note that this expansion is
valid for $\eta\xi\ll 1$. These correlations are very specific for a given
experimental baseline and energy binning, and one would expect it to be useful
in distinguishing ultralight scalar phenomenology from other time-dependent
effects like Lorentz symmetry violation (see e.g. Ref. Kostelecky and Mewes
(2012)). Besides, due to the moving of the oscillation minimum, the change in
oscillation probability at energies near the minimum can be enhanced. In this
example, $\eta=0.05$ can lead to $\Delta P(t)\simeq 0.2$.

Figure 2: The difference in $\nu_{\mu}$ disappearance oscillation probability
$\Delta P(t)\equiv P_{\mu\mu}(t)-P_{\mu\mu}(0)$ at DUNE for a modulating
$\Delta m^{2}_{31}$ assuming $\eta=0.05$ for several energy bins. We have
assumed normal mass ordering and the best fit values of oscillation parameters
from Ref. de Salas  _et al._ (2018), namely $\Delta m^{2}_{31}=2.5\times
10^{-3}$ eV2, and $\sin^{2}\theta_{23}=0.55$.

A search for time-dependent frequencies in the data would seem suitable to
probe this scenario, for modulations of either mixing angles or mass
splittings. Such studies have already been performed by e.g. SNO, Super-
Kamiokande, and Daya Bay collaborations Yoo  _et al._ (2003); Aharmim  _et
al._ (2005, 2010); Adey  _et al._ (2018). In the next section, we will provide
a detailed time-dependent analysis of mock experimental data (we use the DUNE
set up as an example), employing the Lomb-Scargle method, to estimate
sensitivity to ultralight scalars.

###  III.2 Average Distorted Neutrino Oscillations

The second regime we study here is dubbed average distorted neutrino
oscillations or average DiNOs for short. In this regime, the modulation of
mixing angles or mass splittings is too fast to be observed, but an averaging
effect on oscillation probabilities remains. To see how this comes about,
imagine that the modulating period of the ultralight scalar field is much
shorter than the experimental data collecting time, but still much longer than
the neutrino time-of-flight (such that the oscillation parameters are
essentially constant for each neutrino event). The average oscillation
probability would be given by Krnjaic  _et al._ (2018)

$\displaystyle\langle
P_{\alpha\beta}\rangle=\frac{1}{\tau_{\phi}}\int_{0}^{\tau_{\phi}}dt\,P_{\alpha\beta}(t),$
(11)

where $\tau_{\phi}$ is the period of the ultralight scalar field. First, we
focus on the averaging of mass splittings, which will prove to be more
interesting than averaging of angles. For the vacuum, 2-flavor, muon neutrino
disappearance oscillation, the averaging is

$\displaystyle\langle P^{\rm mass}_{\alpha\beta}\rangle$
$\displaystyle=\frac{1}{\tau_{\phi}}\int_{0}^{\tau_{\phi}}dt\left\\{1-\sin^{2}(2\theta)\sin^{2}\left[\left(\frac{\Delta
m^{2}L}{4E}\right)\left[1+2\eta\sin(m_{\phi}t)\right]\right]\right\\}$ (12)
$\displaystyle\simeq 1-\sin^{2}(2\theta)\left\\{\sin^{2}\left(\frac{\Delta
m^{2}L}{4E}\right)+2\eta^{2}\left(\frac{\Delta
m^{2}L}{4E}\right)^{2}\cos\left(\frac{\Delta m^{2}L}{2E}\right)\right\\},$
(13)

Figure 3: Averaged $\nu_{\mu}$ disappearance probability at DUNE for a
modulating $\Delta m^{2}_{31}$ for three different values of $\eta$ and in the
absence of DiNOs ($\eta$= 0).

where the last expression was expanded to order $(\eta\Delta m^{2}L/4E)^{2}$,
which is typically fine for the first oscillation minimum and $\eta<0.05$. As
shown in Fig. 3 for oscillation parameters from Ref. de Salas  _et al._ (2018)
and different values of $\eta$, the effect of mass splitting averaging is a
smearing in the oscillation probability, similar to the effect of a finite
energy resolution. Therefore, one would expect experiments like DUNE, KamLAND
and JUNO to be ideal to probe such scenarios, as all of them have good energy
resolution and are endowed with broad band beams that allow for the
observation of the full shape of the oscillation probability.

Regarding the averaging of modulating angles, it is easy to see that this
effect simply maps into standard oscillation probability with different
parameters. The average mixing would be Krnjaic  _et al._ (2018)

$\displaystyle\frac{1}{\tau_{\phi}}\int_{0}^{\tau_{\phi}}dt\,\sin^{2}\left[2\theta+2\eta\sin(m_{\phi}t)\right]=\frac{1}{2}\left[1-J_{0}(4\eta)\cos(4\theta)\right]\simeq\sin^{2}(2\theta)(1-4\eta^{2})+2\eta^{2},$
(14)

where $J_{0}$ is a Bessel function of the first kind and the righthand side
was expanded to second order in $\eta$. As one expects, the effect of
averaging is pushing apparent mixing angles away from zero or maximal mixing.
Experimental sensitivity to average DiNO effects on mixing angles depends not
only on the precision with which the experiment can determine the mixing
angles but also on the value of the measured angle itself. Because of that, we
will focus on the modulation of mass splittings in the case of average DiNO.

###  III.3 Dynamical Distorted Neutrino Oscillations

The last regime is the dynamical DiNO, in which the effect of modulation of
the neutrino mass matrix in oscillations needs to be taken into account at the
Hamiltonian level. The matter potential induced by the ultralight scalar field
can be written as Brdar  _et al._ (2018)

$V_{\phi}(t)=\frac{1}{2E}\left[(Ym_{\nu}+m_{\nu}Y)\phi(t)+Y^{2}\phi^{2}(t)\right],$
(15)

where we have defined the coupling matrix $Y\equiv m_{\nu}y/\Lambda$. Then the
Hamiltonian that drives the evolution of the system is given by

$H(t)=H_{\rm vac}+V_{\rm matter}+V_{\phi}(t),$ (16)

where $H_{\rm vac}$ is the vacuum Hamiltonian, that is,
$\textrm{diag}(0,\,\Delta m^{2}_{21},\,\Delta m^{2}_{31})/2E$ in the mass
basis, $V_{\rm matter}$ is the usual matter potential, $V_{\rm
matter}=U^{\dagger}\textrm{diag}(V_{\rm cc},\,0,\,0)U$ in the mass basis, and
$U$ is the PMNS matrix. Since $\phi$ evolves in time, the full Hamiltonian
also depends on time.

Here we propose a simplification of this treatment which is valid when
considering a modulation effect either solely in the mixing angles or solely
in the mass splittings. If the modulation of $\phi$ affects solely the mass
splittings, we can work in the mass basis (denoted by a “0” subscript) and
write

$H_{0}(t)\equiv\frac{1}{2E}\left(\begin{array}[]{ccc}0&0&0\\\ 0&\Delta
m^{2}_{21}(t)&0\\\ 0&0&\Delta
m^{2}_{31}(t)\end{array}\right)+U^{\dagger}\left(\begin{array}[]{ccc}V_{\rm
cc}&0&0\\\ 0&0&0\\\ 0&0&0\end{array}\right)U$ (17)

where the time dependent mass splittings are given in Eq. (6). We write the
neutrino state in the mass basis as $\nu=U^{\dagger}\nu_{\rm fl}$ (the “fl”
subscript denotes flavor).

The evolution of this state is given by $\dot{\nu}=iH_{0}(t)\nu$, which has to
be solved numerically and leads to an “instantaneous” oscillation probability
from flavor $\alpha$ to $\beta$, namely, $P_{\alpha\beta}(t_{0},L)$.
Numerically, the propagation of the neutrino from the source to the detector
can be implemented by dividing the path into $N$ layers of thickness $\Delta
L=L/N$ where the oscillation parameters do not change substantially in each
layer. The oscillation probability from flavor $\alpha$ to flavor $\beta$ is
given by

$P_{\alpha\beta}(t_{0},L)=|\bra{\nu_{\alpha}}U\left\\{\prod_{n=1}^{N}\exp\left[iH_{0}(t_{n})\Delta
L\right]\right\\}U^{\dagger}\ket{\nu_{\beta}}|^{2}~{}~{}~{}~{}~{}~{}~{}~{}~{}(\textrm{for
}\Delta m^{2}\textrm{ modulation}),$ (18)

where $U^{\dagger}\ket{\nu_{\alpha}}$ is a flavor state $\ket{\nu_{\alpha}}$
written in the mass basis, and $t_{n}\equiv t_{0}+n\Delta L$. Note that due to
the presence of the MSW term, the Hamiltonian in different layers do not
commute with each other. Throughout the duration of the experiment, all
possible initial phases will be scanned randomly. Therefore, the observed
oscillation probability is the time average of $P_{\alpha\beta}(t_{0},L)$,
namely

$\langle
P_{\alpha\beta}(L)\rangle=\frac{1}{\tau_{\phi}}\int_{0}^{\tau_{\phi}}dt_{0}P_{\alpha\beta}(t_{0},L),$
(19)

where again, $\tau_{\phi}=2\pi/m_{\phi}$ is the period of oscillation of
$\phi$.

For the case of modulating mixing angle, it is more convenient to write the
Hamiltonian in the flavor basis (again denoted by the subscript “fl”) and
write the full Hamiltonian as

$H_{\rm fl}(t)\equiv\frac{1}{2E}U(t)\left(\begin{array}[]{ccc}0&0&0\\\
0&\Delta m^{2}_{21}&0\\\ 0&0&\Delta
m^{2}_{31}\end{array}\right)U^{\dagger}(t)+\left(\begin{array}[]{ccc}V_{\rm
cc}&0&0\\\ 0&0&0\\\ 0&0&0\end{array}\right),$ (20)

where the PMNS matrix is now changing in time. The evolution of a neutrino of
definite flavor $\nu_{\rm fl}$ is $\dot{\nu}_{\rm fl}=iH_{\rm fl}(t)\nu_{\rm
fl}$ and the observed oscillation probability is given by Eq. (19). Following
the previous recipe, the oscillation probability can be calculated with

$P_{\alpha\beta}(t_{0},L)=|\bra{\nu_{\alpha}}\left\\{\prod_{n=1}^{N}\exp\left[iH_{\rm
fl}(t_{n})\Delta
L\right]\right\\}\ket{\nu_{\beta}}|^{2}~{}~{}~{}~{}~{}~{}~{}~{}~{}(\textrm{for
}\theta_{ij}\textrm{ modulation}).$ (21)

Note that we do not need to rotate the initial flavor neutrino by the PMNS
matrix as we decided to work here in the more appropriate flavor basis. The
observed oscillation probability is again given by the average over $t_{0}$,
as in Eq. (19).

We can see the effect of dynamical DiNOs in Fig. 4 for modulations of $\Delta
m^{2}_{31}$ (left) and $\theta_{23}$ (right) for oscillation parameters from
Ref. de Salas  _et al._ (2018) and a modulation amplitude of $\eta=0.1$ and
$L=1300$ km. In the left panel, it is clear to see that as
$\tau_{\phi}/\tau_{\nu}$ gets larger, the effect shrinks, as the changes in
the matter effect become too fast to affect neutrino oscillations. Besides, as
$\tau_{\phi}/\tau_{\nu}$ gets smaller, the effect of dynamical DiNOs
asymptotes to that of average DiNOs, as one would expect. For mixing angle
modulations (right panel), the Jacobian effect discussed in Sec. III.1
suppresses the impact of modulating $\theta_{23}$ in DUNE. Notice also that
there is a small displacement of the minima and maxima of oscillations. This
is simply due to the fact that the effective mass splitting measured in long
baseline $\nu_{\mu}$ disappearance is not exactly $\Delta m^{2}_{31}$, but
rather a function of atmospheric splittings and mixings, typically dubbed
$\Delta m^{2}_{\mu\mu}$ Nunokawa  _et al._ (2005). In this plot, we have
chosen a fixed value of $\Delta m^{2}_{31}$ to obtain the curves. This
displacement could simply be mapped into a different value of the atmospheric
mass splitting.

Figure 4: The $\nu_{\mu}$ disappearance probability at DUNE for a modulating
$\Delta m^{2}_{31}$ (left) and $\theta_{23}$ (right) in the dynamical regime,
for different values of the ratio between the period of the modulation
$\tau_{\phi}$ and the neutrino time of flight $\tau_{\nu}$.

##  IV Case study: DUNE phenomenology

In this section, we will perform a case study of how the oscillation
experiments can probe the ultralight scalars that we introduced in the
previous sections. We will analyze the DUNE sensitivity to the aforementioned
three regimes. For all analyses performed here, we followed the simulation of
Ref. Alion  _et al._ (2016) using the GLoBES software Huber  _et al._ (2005,
2007) We have assumed the DUNE experiment to have a 1.07 MW beam, four far
detectors with a total mass of 40 kton, and a total run time of 7 years,
equally divided between neutrino and antineutrino mode. We took the matter
density to be constant $\rho=2.8~{}{\rm g}/{\rm cm}^{3}$. The neutrino
spectrum at DUNE spans energies roughly from 1-5 GeV, with a peak at around 3
GeV. From Ref. Alion  _et al._ (2016), the main systematic uncertainties are
related to the beam and cross-section modeling. As it is important for some of
the results that will follow, we call the attention that the reconstructed
neutrino energy resolution from Ref. Alion  _et al._ (2016) is approximately
16% at 3 GeV. Finally, in general terms, the ultralight phenomenology is
better probed by the $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ disappearance channels,
as they bring in the largest statistics.

In our numerical simulations, the best fit oscillation parameters were taken
from Ref. de Salas  _et al._ (2018) ($\Delta m^{2}_{21}=7.5\times 10^{-5}$
eV2, $\Delta m^{2}_{31}=2.5\times 10^{-3}$ eV2, $\sin^{2}\theta_{12}=0.32$,
$\sin^{2}\theta_{13}=0.0216$, $\sin^{2}\theta_{23}=0.55$, $\delta=1.2\pi$) and
normal mass ordering was assumed for concreteness. We have used current priors
for the oscillation parameters that were marginalized over in the analysis
performed 222We took the largest one-sided $1\sigma$ error for each
oscillation parameter from Ref. de Salas  _et al._ (2018))..

###  IV.1 Time modulation phenomenology at DUNE

The signature of the time modulation regime is the presence of a periodic
signal in the oscillated neutrino spectrum at the far detector. As can be seen
in Eqs. (5) and (6), the period of modulation is given by the mass of the
scalar field, and thus a positive signal of time modulation at DUNE would
provide a measurement of the mass of this field.

In general, the search for periodic signals in data sets can be performed
using the Lomb-Scargle periodogram Lomb (1976); Scargle (1982), which is an
extension of the classical periodogram for unevenly separated data (see Ref.
VanderPlas (2018) for a pedagogical review). Such searches are not new in the
context of neutrino physics. For instance, many analyses of periodicities in
the solar neutrino flux have been performed in SNO Aharmim  _et al._ (2005);
Ranucci and Rovere (2007) and SuperKamiokande Ranucci (2006), as well as for
other time-varying signals in the context of Lorentz and CPT violation in Daya
Bay Adey  _et al._ (2018).

To evaluate the statistical significance of a modulation in a data set, we
first define the Lomb-Scargle (LS) power for a frequency $\omega$ as

$\begin{split}P_{LS}(\omega)=\frac{1}{2}\Bigg{\\{}\bigg{(}\sum_{n}g_{n}\cos[2\pi\omega(t_{n}-\tau)]\bigg{)}^{2}\bigg{/}\sum_{n}\cos^{2}[2\pi\omega(t_{n}-\tau)]+\\\
\bigg{(}\sum_{n}g_{n}\sin[2\pi\omega(t_{n}-\tau)]\bigg{)}^{2}\bigg{/}\sum_{n}\sin^{2}[2\pi\omega(t_{n}-\tau)]\Bigg{\\}}\end{split}$
(22)

where $g_{n}=g(t_{n})$ is the signal at the time of the measurement $t_{n}$
and $\tau$ is defined by solving $\tan(4\pi\omega\tau)=\sum_{n}\sin(4\pi\omega
t_{n})/\cos(4\pi\omega t_{n})$.

The significance of a given LS power can be quantified with the False Alarm
Probability Test (FAP), which is a measure of how likely it is that a data set
with no signal would give rise to a peak of the same magnitude as a
consequence of spurious background noise. To estimate the FAP we followed the
Baluev approach Baluev (2008), which provides an upper limit to the value
obtained through more sophisticated approaches based on a Bootstrap Method. A
more detailed discussion can be found in Ref. VanderPlas (2018).

In our analysis, we considered larger than nominal energy bins, which allows
us to extend the sensitivity to smaller periods. In order to cover all the
parameter space presented, different time binning of the events has been
explored. In an experimental setup, different choices of the time bins can be
explored a posteriori, ensuring that the whole accessible parameter space is
covered. As discussed in Fig. 2, time modulation effects from ultralight
scalars can lead to correlated (and anti-correlated) modulations in different
energy bins, depending on which oscillation parameter is modulating. In the
case of a positive modulation signal, such correlations can be exploited to
further constrain the model. Different data sets generated for different
values of modulation amplitude and frequency have been tested with the LS
periodogram and Baluev’s approach to the FAP test.

Figure 5: Expected number of events in different time bins at three different
energy bins for a particular realisation of DUNE. Gray dashed lines represent
the number of events expected with $\eta=0$ and the shadowed region
corresponds to the $1\sigma$ statistical uncertainty. Light blue lines are the
expected signal with non-zero $\eta$ and $m_{\phi}$, while dark blue lines
represent a random realisation of the signal which includes statistical
fluctuations. Left panel corresponds to a set of parameters $(\eta,m_{\phi})$
in the 90$\%$ CL region above shown. Right panel corresponds to a point in the
$\eta-m_{\phi}$ plane which LS method would not identify the frequency
correctly.

A typical realization of the signal in DUNE, for modulation of mass
splittings, is presented in Fig. 5 for illustrative purposes. Two modulation
amplitudes are shown, $\eta=0.063$ and $\eta=0.012$. As we will see later,
DUNE is expected to be sensitive to the larger but not to the smaller. The
expected number of events is presented for three different energy bins
assuming 3 years running time in the neutrino mode 333Although we show only 3
energy bins, we use the entire energy range available in our simulations.. The
original period of the data ($\tau_{\phi}$) and the one determined with the LS
method ($\tau_{LS}$) are presented together with the FAP score.

In Fig. 6, we show the DUNE sensitivity to modulations of the mass splitting
$\Delta m^{2}_{31}(t)$ as discussed in Eq. (6). We present the parameter space
for which the Lomb-Scargle method would identify a frequency in the data set
such that one would be able to state that there is a 90$\%$ probability that
the periodic signal found in data is not due to random noise. It is important
to point out that this is different from stating that there is a periodic
signal with a given frequency in the data set.

Figure 6: Region in the modulation amplitude $\eta$ versus period
$\tau_{\phi}$ plane for which the analysis performed using the Lomb-Scargle
periodogram would be capable to state at a 90$\%$ C.L. that the periodic
signal found in 7 years of DUNE data would not be spurious. Upper x-axis
presents the corresponding range for the mass of the ultralight scalar,
$m_{\phi}$.

The smallest period (largest frequency) to which DUNE is expected to be
sensitive corresponds to the case in which the scalar oscillation period is
about a few days. Sensitivity at large frequencies requires smaller time bins,
leading to fewer events per bin which are more prone to statistical
fluctuation. Consequently, the sensitivity to the amplitude decreases when
looking for shorter periods. The largest period is determined by the running
time of the experiment and is expected to be a couple of years. This means
that DUNE would be sensitive to the range of scalar masses $5\times
10^{-23}\lesssim m_{\phi}\lesssim 10^{-20}$ eV. Moreover, the maximum
sensitivity on the amplitude of the modulation $\eta$ is about 4%.

We do not show the DUNE sensitivity to modulations of the mixing angle
$\theta_{23}$. Due to the Jacobian effect discussed in Sec. III.1, the values
of $\eta$ necessary to yield an observable effect at DUNE would need to be
$\mathcal{O}(1)$. Therefore we see the mass splitting modulation search as a
more promising way of detecting the presence of ultralight bosonic dark
matter. Nevertheless, the same Lomb-Scargle technique described above could be
used to probe the modulations of mixing angles.

This analysis is not free from systematic uncertainties. Deadtime intervals of
any neutrino detector may span from months to milliseconds. These would
constitute an important source of systematic errors when analyzing time
modulations. In particular, quasiperiodic deadtimes such as scheduled breaks
of runs or maintenance and calibration operations could lead to peaks in the
Lomb-Scargle power spectrum. It is possible to deal with them by estimating
the window function, see Ref. VanderPlas (2018). This method allows us to
identify the main features of the structure of the window. The existence of a
certain periodicity in data taking would induce peaks both in the Lomb-Scargle
power spectrum and in the window power spectrum and consequently, one can
identify the corresponding frequencies as related to the experiment and not to
the physical phenomenon under study.

Beam unrelated backgrounds, like cosmic rays and atmospheric neutrino, can
also exhibit time modulation, but these are typically negligible in beam
neutrino experiments. Moreover, beam performance can change over time. One
would expect near-to-far ratios to be less sensitive to these variations.
Nevertheless, the observed neutrino beam is not the same in the near and far
detectors (due to different geometry), and oscillations may further enhance
this difference. Thus, experiments should take those systematics into
consideration when performing a search for the time dependence of the neutrino
signal. In our results, we have not considered those systematics, though we
believe they would not degrade the overall sensitivity by much.

###  IV.2 Average distorted neutrino oscillations at DUNE

The second regime we will study in detail in DUNE is the case of average
distorted neutrino oscillations (average DiNOs). In this case, modulations of
mass splittings are too fast to be observed as a time modulation signal, but
the averaging of the modulation still imprint observable effects in the
oscillation probability. As discussed in Sec. III.2, the modulation of mass
splittings leads to more interesting phenomenology than the modulation of
mixing angles. If $\Delta m^{2}_{31}$ varies in time, the maxima and minima
are displaced periodically. For a very fast modulation, such displacement
manifests as a non-trivial averaging and has to be carefully studied in order
to disentangle it from the distortion caused by the finite energy resolution
Krnjaic  _et al._ (2018). Consequently, the searches here presented would
benefit from improvements in the energy resolution, as the one proposed in
Ref. Friedland and Li (2019). As a general rule of thumb, the new physics
effect here is just a modification of the oscillation probability, so all
systematic uncertainties associated to an oscillation search would be relevant
in this case. We include those systematics using the simulation from Ref.
Alion  _et al._ (2016).

The average oscillation probability in Eq. (12), necessary to estimate DUNE
sensitivity to this scenario, was calculated by numerically averaging the
analytic approximations for neutrinos oscillations in the presence of matter
from Ref. Denton  _et al._ (2016). The DUNE sensitivity to fast modulations of
the mass splitting is shown in Fig. 7. The range of periods that DUNE is
sensitive to, in this regime, is approximately from a year
($\tau_{\phi}\ll\tau_{\rm exp}\sim$ 10 years) to tens of milliseconds
($\tau_{\phi}\gg\tau_{\nu}\simeq$ 4.3 msec). This translates into a very large
range of masses $2\times 10^{-23}\lesssim m_{\phi}\lesssim 3\times 10^{-14}$
eV. On the flip side, if a signal of average DiNOs is observed at DUNE, it
would be very challenging to pinpoint the exact mass of $\phi$, as its effects
has been averaged out. As can be seen from the left panel of Fig. 7, DUNE is
expected to be sensitive to roughly 4% modulation amplitudes of the mass
splitting at 90% C.L. In estimating DUNE’s sensitivity, we have marginalized
over all oscillation parameters. Note that the energy resolution is crucial
here, and we expect improvements in energy resolution to translate into better
sensitivity to average DiNOs.

Besides, one may ask what is the impact on the determination of the standard
mixing parameters if an average DiNO effect is present. We show in the right
panel of Fig. 7 the allowed region by DUNE in the plane
$(\eta,\,\sin^{2}2\theta_{23})$. There is some mild degeneracy between
$\theta_{23}$ and the modulation amplitude $\eta$ for small values of $\eta$.
This is because, for small $\eta$, the minimum of the $\nu_{\mu}\to\nu_{\mu}$
oscillation probability lifts slightly thus allowing values of $\theta_{23}$
closer to $\pi/4$. We have checked that there is no degeneracy between
$\delta_{CP}$ and $\eta$.

Figure 7: Left: $\chi^{2}$ profile as a function of $\eta$ and confidence
levels to discriminate with respect to the standard 3$\nu$ picture. Right:
allowed region in the plane $(\eta,\,\sin^{2}\theta_{23})$ in DUNE. Best fit
values chosen for mock data are marked by a black circle (standard
oscillations assumed as true hypothesis).

###  IV.3 Dynamical distorted neutrino oscillations

Now, we proceed to the last regime of ultralight scalar field phenomenology in
neutrino experiments, that of dynamical distorted neutrino oscillations
(dynamical DiNOs). In this regime, the time of flight of the neutrino is
comparable to the magnitude of the scalar period and we have to take full
account of the variation of the oscillation parameters during the journey from
the source to the detector. As before, we consider modulations in
$\theta_{23}$ and $\Delta m^{2}_{31}$ separately, which leads to the
oscillation probabilities discussed in Sec. III.3. As discussed in the
previous section, the new physics effect here is again just a modification of
the oscillation probability, so all systematic uncertainties associated with
an oscillation search would be relevant in this case.

We evaluate the DUNE sensitivity to the dynamical DiNO regime for the case of
mass splitting modulations and present the results in Figs. 8. As expected,
when the modulation period is much smaller than the neutrino time-of-flight
(of about 4.3 msec), the experimental sensitivity degrades as the oscillation
probability tends to the standard one. Besides, for mass splitting modulation,
when the modulation period is sufficiently large, we recover the average DiNO
sensitivity, see Fig. 4. In the case of mixing angle modulation, the averaging
simply maps the modulating mixing angle onto another value of the mixing, and
thus we do not analyze this case here. In principle, a richer flavor structure
in the coupling between ultralight scalar and neutrinos would allow
concomitant modulation of mass splittings and angles with arbitrary
correlations. For simplicity, we do not pursue this possibility in this
manuscript.

Figure 8: DUNE sensitivity in a fast modulation of the $\Delta m^{2}_{31}$
mass splitting with amplitude $\eta$ versus the mass of the ultralight scalar
$m_{\phi}$ plane at 90$\%$ C.L (green) and $3\sigma$ (blue). Upper x-axis
presents the corresponding range for the dimensionless quantity
$m_{\phi}L_{DUNE}/2\pi$.

Finally, we present in Fig. 9 DUNE’s sensitivity to ultralight scalar dark
matter incorporating the three searches proposed in this paper. It is
remarkable that a neutrino oscillation experiment can probe scalar masses
ranging about 10 orders of magnitude. The transition between dynamical DiNOs
and average DiNOs can be seen to be smooth. When time modulation can be seen
(LS labeled region), it provides an even better probe of ultralight scalars.
While its sensitivity spans a very large region in parameter space, DUNE would
only be able to measure the mass of $\phi$ by the determination of the
modulation frequency, which is only possible in the “LS” region. Elsewhere,
even if distorted neutrino oscillations are observed, it is not clear on how
to determine the ultralight scalar mass.

Figure 9: DUNE sensitivity (90% C.L.) to ultralight scalars via modulation of
the atmospheric mass splitting $\Delta m^{2}_{31}$. The amplitude of
modulation is $\eta$ (see Eqs. (5) and (6)), and the modulation period and
mass of the scalar field is denoted by $\tau_{\phi}$ and $m_{\phi}$,
respectively.

##  V Other constraints on ultralight scalars

Besides neutrino oscillations, ultralight bosonic dark matter may also be
probed by other experiments and cosmological observations. Several other
constraints can be found in previous literature. Here we simply summarize and
comment on them.

  * •

CMB: As pointed out in previous work, the cosmological $\phi$ density
redshifts as nonrelativistic matter, and thus the amplitude of modulation of
$\phi$ is expected to be much larger at early times. This could lead to an
increase in the sum of neutrino masses in the early universe, which would be
constrained by observations of the Planck satellite Ade  _et al._ (2016). The
constraint on the amplitude of $\phi$ is found to be $\eta\lesssim 9\times
10^{-3}$ if the atmospheric mass splitting modulates or $\eta\lesssim 0.1$ if
only the solar splitting modulates Berlin (2016); Krnjaic  _et al._ (2018);
Brdar  _et al._ (2018). Note however that this constraint is model dependent.
For example, if neutrino masses arrive from a seesaw mechanism and $\phi$
couples to, e.g., right-handed neutrinos, a large mass of the latter would
actually imply a smaller mass of active neutrinos in the early universe.

  * •

BBN: Following a similar reasoning, if $\phi$ couples universally to all
standard model fermions, the changes in fermion masses could lead to
observable effects in big bang nucleosynthesis, particularly on the abundance
of 4He Sibiryakov  _et al._ (2020). In our scenario, we are only coupling
$\phi$ to neutrinos, and therefore this bound would not apply. In a UV
complete realization, this bound should be taken into account carefully.

  * •

Astrophysical neutrinos and supernova 1987A: The coupling between neutrinos
and the ultralight scalar could lead to scattering of astrophysical neutrinos
on the dark matter background, making the universe opaque to certain
astrophysical neutrinos. In particular, the observation of neutrinos from the
supernova 1987A requires the universe to be transparent to MeV neutrinos.
These constraints are several orders of magnitude weaker than the ones derived
from neutrino oscillation measurements, except for large masses
$m_{\phi}\gtrsim 10^{-8}$ eV where neutrino experiments lose sensitivity
Krnjaic  _et al._ (2018); Brdar  _et al._ (2018).

  * •

Electron mass modulation: Even if we postulate that $\phi$ only couples to
neutrinos at tree level, a loop level coupling to electrons cannot be avoided,
which could potentially lead to stronger constraints. Nevertheless, the
magnitude of the loop-induced coupling is suppressed by $G_{F}m_{e}m_{\nu}$,
too small to be observed Krnjaic  _et al._ (2018).

  * •

Black hole spin measurements: Ultralight bosons can be emitted copiously by
rotating black holes, thereby suppressing their energy and angular momentum.
The constraints from black hole spin measurements are typically in the
parameter space within the region $10^{-21}<m_{\phi}<10^{-16}~{}$eV. Note that
these bounds present some dependence on the systematic uncertainties related
to black hole mass and spin observations, as well as on the details of the UV
model (for instance, scalar self-interactions may weaken those bounds)
Baryakhtar  _et al._ (2017); Davoudiasl and Denton (2019); Zu  _et al._
(2020).

##  VI Conclusion

In this paper, we have investigated in detail the phenomenology of ultralight
bosonic dark matter fields in neutrino oscillation experiments. We describe
the three phenomenological regimes associated with different masses of the
bosonic field, together with their respective signatures at neutrino
oscillation experiments. As a case study, we estimate the sensitivity of the
DUNE experiment to ultralight scalars via detailed simulations and we provide
a description of how to implement an experimental search for all regimes.
Finally, we show that the DUNE experiment is sensitive to about 10 orders of
magnitude in ultralight scalar mass. The searches presented in this paper are
general and can be easily applied to any other oscillation experiment. Current
experiments like MINOS, NOvA, T2K, Daya Bay, and RENO may perform the very
first experimental searches for these scenarios. In the future, besides DUNE,
we also expect the JUNO to be particularly sensitive to modulations of mass
splittings due to ultralight scalars. Our results exhibit an excellent, unique
opportunity to expand the physics program of neutrino oscillation experiments.

##  VII Acknowledgements

We thank Vedran Brdar, Zackaria Chacko, Gordan Krnjaic, Joachim Kopp, Alex
Sousa and Jaehoon Yu for discussions. PMM is grateful for the kind hospitality
of Fermilab theory group during the development of this work. This manuscript
has been authored by Fermi Research Alliance, LLC under Contract No. DE-
AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office
of High Energy Physics. PMM is supported by the FPU grant FPU18/04571, the
Spanish grants FPA2017-90566-REDC (Red Consolider MultiDark), FPA2017-85216-P
and SEV-2014-0398(MINECO/AEI/FEDER, UE), as well as PROMETEO/2018/165
(Generalitat Valenciana). PMM received partial support from the EU Horizon
2020 project InvisiblesPlus (690575-InvisiblesPlus-H2020-MSCA-RISE-2015). AD
is supported in part by the National Science Foundation under Grant Number
PHY-1914731.

## References

  * Tanabashi  _et al._ (2018) M. Tanabashi _et al._ (Particle Data Group), Phys. Rev. D98, 030001 (2018).
  * Tremaine and Gunn (1979) S. Tremaine and J. E. Gunn, Phys. Rev. Lett. 42, 407 (1979).
  * Madsen (1990) J. Madsen, Phys. Rev. Lett. 64, 2744 (1990).
  * Madsen (1991) J. Madsen, Phys. Rev. D44, 999 (1991).
  * Boyarsky  _et al._ (2009) A. Boyarsky, O. Ruchayskiy, and D. Iakubovskyi, JCAP 0903, 005 (2009), arXiv:0808.3902 [hep-ph] .
  * Peccei and Quinn (1977a) R. Peccei and H. R. Quinn, Nuovo Cim. A 41, 309 (1977a).
  * Peccei and Quinn (1977b) R. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977b).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Vafa and Witten (1984) C. Vafa and E. Witten, Phys. Rev. Lett. 53, 535 (1984).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie (1983) L. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Hu  _et al._ (2000) W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astro-ph/0003365 [astro-ph] .
  * Arias  _et al._ (2012) P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, JCAP 1206, 013 (2012), arXiv:1201.5902 [hep-ph] .
  * Li  _et al._ (2014) B. Li, T. Rindler-Daller, and P. R. Shapiro, Phys. Rev. D89, 083536 (2014), arXiv:1310.6061 [astro-ph.CO] .
  * Hui  _et al._ (2017) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Menci  _et al._ (2017) N. Menci, A. Merle, M. Totzauer, A. Schneider, A. Grazian, M. Castellano, and N. G. Sanchez, Astrophys. J. 836, 61 (2017), arXiv:1701.01339 [astro-ph.CO] .
  * Diez-Tejedor and Marsh (2017) A. Diez-Tejedor and D. J. E. Marsh, (2017), arXiv:1702.02116 [hep-ph] .
  * Irsic  _et al._ (2017) V. Irsic, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, Phys. Rev. Lett. 119, 031302 (2017), arXiv:1703.04683 [astro-ph.CO] .
  * Visinelli (2017) L. Visinelli, Phys. Rev. D96, 023013 (2017), arXiv:1703.08798 [astro-ph.CO] .
  * Kumar Poddar  _et al._ (2020) T. Kumar Poddar, S. Mohanty, and S. Jana, Phys. Rev. D 101, 083007 (2020), arXiv:1906.00666 [hep-ph] .
  * Moore (1994) B. Moore, Nature 370, 629 (1994).
  * Flores and Primack (1994) R. A. Flores and J. R. Primack, Astrophys. J. 427, L1 (1994), arXiv:astro-ph/9402004 [astro-ph] .
  * Navarro  _et al._ (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), arXiv:astro-ph/9611107 [astro-ph] .
  * Klypin  _et al._ (1999) A. A. Klypin, A. V. Kravtsov, O. Valenzuela, and F. Prada, Astrophys. J. 522, 82 (1999), arXiv:astro-ph/9901240 [astro-ph] .
  * Battaglieri  _et al._ (2017) M. Battaglieri _et al._ , in  _U.S. Cosmic Visions: New Ideas in Dark Matter College Park, MD, USA, March 23-25, 2017_ (2017) arXiv:1707.04591 [hep-ph] .
  * Berlin (2016) A. Berlin, Phys. Rev. Lett. 117, 231801 (2016), arXiv:1608.01307 [hep-ph] .
  * Capozzi  _et al._ (2017) F. Capozzi, I. M. Shoemaker, and L. Vecchi, JCAP 1707, 021 (2017), arXiv:1702.08464 [hep-ph] .
  * Krnjaic  _et al._ (2018) G. Krnjaic, P. A. N. Machado, and L. Necib, Phys. Rev. D97, 075017 (2018), arXiv:1705.06740 [hep-ph] .
  * Brdar  _et al._ (2018) V. Brdar, J. Kopp, J. Liu, P. Prass, and X.-P. Wang, Phys. Rev. D97, 043001 (2018), arXiv:1705.09455 [hep-ph] .
  * Farzan and Palomares-Ruiz (2019) Y. Farzan and S. Palomares-Ruiz, Phys. Rev. D99, 051702 (2019), arXiv:1810.00892 [hep-ph] .
  * Liao  _et al._ (2018) J. Liao, D. Marfatia, and K. Whisnant, JHEP 04, 136 (2018), arXiv:1803.01773 [hep-ph] .
  * Ge and Parke (2019) S.-F. Ge and S. J. Parke, Phys. Rev. Lett. 122, 211801 (2019), arXiv:1812.08376 [hep-ph] .
  * Huang and Nath (2018) G.-Y. Huang and N. Nath, Eur. Phys. J. C 78, 922 (2018), arXiv:1809.01111 [hep-ph] .
  * Ge and Murayama (2019) S.-F. Ge and H. Murayama, (2019), arXiv:1904.02518 [hep-ph] .
  * Cline (2020) J. M. Cline, Phys. Lett. B 802, 135182 (2020), arXiv:1908.02278 [hep-ph] .
  * Choi  _et al._ (2020) K.-Y. Choi, E. J. Chun, and J. Kim, Phys. Dark Univ. 30, 100606 (2020), arXiv:1909.10478 [hep-ph] .
  * Farzan (2019) Y. Farzan, Phys. Lett. B 797, 134911 (2019), arXiv:1907.04271 [hep-ph] .
  * Ding and Feruglio (2020) G.-J. Ding and F. Feruglio, (2020), arXiv:2003.13448 [hep-ph] .
  * Yoo  _et al._ (2003) J. Yoo _et al._ (Super-Kamiokande), Phys. Rev. D 68, 092002 (2003), arXiv:hep-ex/0307070 .
  * Aharmim  _et al._ (2005) B. Aharmim _et al._ (SNO), Phys. Rev. D72, 052010 (2005), arXiv:hep-ex/0507079 [hep-ex] .
  * Aharmim  _et al._ (2010) B. Aharmim _et al._ (SNO), Astrophys. J. 710, 540 (2010), arXiv:0910.2433 [astro-ph.SR] .
  * de Salas  _et al._ (2018) P. F. de Salas, D. V. Forero, C. A. Ternes, M. Tortola, and J. W. F. Valle, Phys. Lett. B782, 633 (2018), arXiv:1708.01186 [hep-ph] .
  * An  _et al._ (2016) F. An _et al._ (JUNO), J. Phys. G 43, 030401 (2016), arXiv:1507.05613 [physics.ins-det] .
  * Kostelecky and Mewes (2012) A. Kostelecky and M. Mewes, Phys. Rev. D 85, 096005 (2012), arXiv:1112.6395 [hep-ph] .
  * Adey  _et al._ (2018) D. Adey _et al._ (Daya Bay), Phys. Rev. D98, 092013 (2018), arXiv:1809.04660 [hep-ex] .
  * Nunokawa  _et al._ (2005) H. Nunokawa, S. J. Parke, and R. Zukanovich Funchal, Phys. Rev. D 72, 013009 (2005), arXiv:hep-ph/0503283 .
  * Alion  _et al._ (2016) T. Alion _et al._ (DUNE), (2016), arXiv:1606.09550 [physics.ins-det] .
  * Huber  _et al._ (2005) P. Huber, M. Lindner, and W. Winter, Comput. Phys. Commun. 167, 195 (2005), arXiv:hep-ph/0407333 [hep-ph] .
  * Huber  _et al._ (2007) P. Huber, J. Kopp, M. Lindner, M. Rolinec, and W. Winter, Comput. Phys. Commun. 177, 432 (2007), arXiv:hep-ph/0701187 [hep-ph] .
  * Lomb (1976) N. R. Lomb, Astrophys. Space Sci. 39, 447 (1976).
  * Scargle (1982) J. D. Scargle, Astrophys. J. 263, 835 (1982).
  * VanderPlas (2018) J. VanderPlas, The Astrophysical Journal Supplement Series 236, 1 (2018), arXiv:1703.09824 [astro-ph] .
  * Ranucci and Rovere (2007) G. Ranucci and M. Rovere, Phys. Rev. D75, 013010 (2007), arXiv:hep-ph/0605212 [hep-ph] .
  * Ranucci (2006) G. Ranucci, Phys. Rev. D73, 103003 (2006), arXiv:hep-ph/0511026 [hep-ph] .
  * Baluev (2008) R. V. Baluev, Mon. Not. Roy. Astron. Soc. 385, 1279 (2008), arXiv:0711.0330 [astro-ph] .
  * Friedland and Li (2019) A. Friedland and S. W. Li, Phys. Rev. D99, 036009 (2019), arXiv:1811.06159 [hep-ph] .
  * Denton  _et al._ (2016) P. B. Denton, H. Minakata, and S. J. Parke, JHEP 06, 051 (2016), arXiv:1604.08167 [hep-ph] .
  * Ade  _et al._ (2016) P. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Sibiryakov  _et al._ (2020) S. Sibiryakov, P. Sorensen, and T.-T. Yu, (2020), arXiv:2006.04820 [hep-ph] .
  * Baryakhtar  _et al._ (2017) M. Baryakhtar, R. Lasenby, and M. Teo, Phys. Rev. D 96, 035019 (2017), arXiv:1704.05081 [hep-ph] .
  * Davoudiasl and Denton (2019) H. Davoudiasl and P. B. Denton, Phys. Rev. Lett. 123, 021102 (2019), arXiv:1904.09242 [astro-ph.CO] .
  * Zu  _et al._ (2020) L. Zu, L. Feng, Q. Yuan, and Y.-Z. Fan, (2020), arXiv:2007.03222 [astro-ph.HE] .