Signatures of Ultralight Dark Matter in Neutrino Oscillation Experiments

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.


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 [1]. 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) [2][3][4][5].
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. [23]). 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 [24][25][26][27][28][29][30][31]. 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 where flavor indices are implicit in both m ν and y, φ is the ultralight scalar field, m ν is the neutrino mass matrix and Λ is a heavy mass scale 1 . 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 where the local density ρ φ should not exceed the local DM density ρ DM = 0.3 GeV/cm 3 , m φ is the mass of φ and v ∼ 10 −3 c is the virial velocity. Note that the space-dependent phase m φ v · x of the local field value is much smaller than m φ t, and thus will be neglected henceforth. The neutrino mass matrix will therefore receive a contribution from φ given by Without a flavor model (see e.g. Ref. [32]), 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 φ defines the modulation period via In a given experimental setup, there are three characteristic time scales: the neutrino time of flight τ ν = L/c = 3.4(L/1000 km) msec, the time between two detected events τ evt (the inverse of the ratio of events) and the lifetime of the experiment τ exp . From these, we can identify the following three different regimes for a given experimental setup.
• Time modulation (τ evt τ φ τ exp ). When the period of modulation of φ 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 [33][34][35]).
• Averaged distorted neutrino oscillations (τ ν τ φ τ 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 [26]. 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 (τ φ ∼ τ ν ). As the modulating period of φ 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 [27]. 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.

LATIONS
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 where θ ij represents the undistorted value of the mixing angle and η is the amplitude of modulation of φ. Mass splittings, on the other hand, are assumed to modulate like (for η 1) where, similarly to above, ∆m 2 ij represents the undistorted value of the mass splitting.

A. 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 For example, in the case of modulating angles, the oscillation probability for ν µ → ν µ disappearance in vacuum, in a simplified two neutrino framework, would read 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 φ t) term.
In Fig. 1 we show change in ν µ disappearance oscillation probability ∆P (t) ≡ P µµ (t)−P µµ (0) at DUNE for a modulating θ 23 assuming η = 0.05 for several energy bins. We have assumed normal mass ordering and the best fit values of oscillation parameters from [36], namely ∆m 2 31 = 2.5 × 10 −3 eV 2 , and sin 2 θ 23 = 0.55. There are three important effects to be noted here. To understand those, we expand ∆P (t) for small η as 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, ∆P (t) depends on η and the mixing angle itself. And finally, when sin(4θ) is near zero (e.g. for θ 23 near maximal mixing), we observe a Jacobian effect that shrinks the oscillation amplitude, as the leading η term shrinks and the η 2 term becomes dominant. In Fig. 1, as the modulation develops, sin 2 [2θ 23 (t)] goes further away from maximal, making ∆P (t) larger, and we observe the larger amplitude. As the phase evolves, sin 2 [2θ 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/τ φ > 0.5). If we had chosen to show a modulation of θ 12 instead (e.g. in the JUNO experiment [37]), 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 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 ν µ disappearance oscillation probability ∆P (t) ≡ P µµ (t) − P µµ (0) at DUNE, but now for a modulating ∆m 2 31 assuming η = 0.05 for several energy bins. Again, we have assumed normal mass ordering and the best fit values of oscillation parameters from [36], namely ∆m 2 31 = 2.5 × 10 −3 eV 2 , and sin 2 θ 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 2 31 assuming η = 0.05 for several energy bins. We have assumed normal mass ordering and the best fit values of oscillation parameters from Ref. [36], namely ∆m 2 31 = 2.5 × 10 −3 eV 2 , and sin 2 θ 23 = 0.55.
can see these correlations in energy by expanding ∆P (t) around η → 0, where we have defined ξ ≡ ∆m 2 L/4E. Note that this expansion is valid for ηξ 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. [38]). 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, η = 0.05 can lead to ∆P (t) 0.2.
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 [33][34][35]39]. 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.

B. 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 [26] where τ φ 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 where the last expression was expanded to order (η∆m 2 L/4E) 2 , which is typically fine for the first oscillation minimum and η < 0.05. As shown in Fig. 3 for oscillation parameters from Ref. [36] and different values of η, 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 where J 0 is a Bessel function of the first kind and the righthand side was expanded to second order in η. 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.

C. 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 [27] where we have defined the coupling matrix Y ≡ m ν y/Λ. Then the Hamiltonian that drives the evolution of the system is given by where H vac is the vacuum Hamiltonian, that is, diag(0, ∆m 2 21 , ∆m 2 31 )/2E in the mass basis, V matter is the usual matter potential, V matter = U † diag(V cc , 0, 0)U in the mass basis, and U is the PMNS matrix. Since φ 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 φ affects solely the mass splittings, we can work in the mass basis (denoted by a "0" subscript) and where the time dependent mass splittings are given in Eq. (6). We write the neutrino state in the mass basis as ν = U † ν fl (the "fl" subscript denotes flavor).
The evolution of this state is given byν = iH 0 (t)ν, which has to be solved numerically and leads to an "instantaneous" oscillation probability from flavor α to β, namely, P αβ (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 ∆L = L/N where the oscillation parameters do not change substantially in each layer. The oscillation probability from flavor α to flavor β is given by where U † |ν α is a flavor state |ν α written in the mass basis, and t n ≡ t 0 + n∆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 αβ (t 0 , L), namely where again, τ φ = 2π/m φ is the period of oscillation of φ.
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 where the PMNS matrix is now changing in time. The evolution of a neutrino of definite flavor ν fl isν fl = iH fl (t)ν fl and the observed oscillation probability is given by Eq. (19). Following the previous recipe, the oscillation probability can be calculated with 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 ∆m 2 31 (left) and θ 23 (right) for oscillation parameters from Ref. [36] and a modulation amplitude of η = 0.1 and L = 1300 km. In the left panel, it is clear to see that as τ φ /τ ν gets larger, the effect shrinks, as the changes in the matter effect become too fast to affect neutrino oscillations. Besides, as τ φ /τ ν 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 A suppresses the impact of modulating θ 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 ν µ disappearance is not exactly ∆m 2 31 , but rather a function of atmospheric splittings and mixings, typically dubbed ∆m 2 µµ [40]. In this plot, we have chosen a fixed value of ∆m 2 31 to obtain the curves. This displacement could simply be mapped into a different value of the atmospheric mass splitting.

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. [41] using the GLoBES software [42,43] 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 ρ = 2.8 g/cm 3 . The neutrino spectrum at DUNE spans energies roughly from 1-5 GeV, with a peak at around 3 GeV. From Ref. [41], 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. [41] is approximately 16% at 3 GeV. Finally, in general terms, the ultralight phenomenology is better probed by the ν µ andν µ disappearance channels, as they bring in the largest statistics.
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 [44,45], which is an extension of the classical periodogram for unevenly separated data (see Ref. [46] 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 [34,47] and SuperKamiokande [48], as well as for other time-varying signals in the context of Lorentz and CPT violation in Daya Bay [39].
To evaluate the statistical significance of a modulation in a data set, we first define the Lomb-Scargle (LS) power for a frequency ω as where g n = g(t n ) is the signal at the time of the measurement t n and τ is defined by solving tan(4πωτ ) = n sin(4πωt n )/ cos(4πω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 [49], 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. [46].
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 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, η = 0.063 and η = 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 3 . The original period of the data (τ φ ) and the one determined with the LS method (τ LS ) are presented together with the FAP score.
In Fig. 6, we show the DUNE sensitivity to modulations of the mass splitting ∆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 FIG. 6. Region in the modulation amplitude η versus period τ φ 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 φ . that this is different from stating that there is a periodic signal with a given frequency in the data set.
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 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. [46]. 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.

B. 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 B, the modulation of mass splittings leads to more interesting phenomenology than the modulation of mixing angles. If ∆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 [26]. Consequently, the searches here presented would benefit from improvements in the energy resolution, as the one proposed in Ref. [50]. 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. [41].
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. [51]. 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 (τ φ τ exp ∼ 10 years) to tens of milliseconds (τ φ τ ν 4.3 msec).
This translates into a very large range of masses 2 × 10 −23 m φ 3 × 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 φ, 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 (η, sin 2 2θ 23 ). There is some mild degeneracy between θ 23 and the modulation amplitude η for small values of η. This is because, for small η, the minimum of the ν µ → ν µ oscillation probability lifts slightly thus allowing values of θ 23 closer to π/4. We have checked that there is no degeneracy between δ CP and η.

C. 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 θ 23 and ∆m 2 31 separately, which leads to the oscillation probabilities discussed in Sec. III C. 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 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 φ 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.

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 φ density redshifts as nonrelativistic matter, and thus the amplitude of modulation of φ 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 [52]. The constraint on the amplitude of φ is found to be η 9×10 −3 if the atmospheric mass splitting modulates or η 0.1 if only the solar splitting modulates [24,26,27]. Note however that this constraint is model dependent. For example, if neutrino masses arrive from a seesaw mechanism and φ 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 φ couples universally to all standard model fermions, FIG. 9. DUNE sensitivity (90% C.L.) to ultralight scalars via modulation of the atmospheric mass splitting ∆m 2 31 . The amplitude of modulation is η (see Eqs. (5) and (6)), and the modulation period and mass of the scalar field is denoted by τ φ and m φ , respectively. the changes in fermion masses could lead to observable effects in big bang nucleosynthesis, particularly on the abundance of 4 He [53]. In our scenario, we are only coupling φ 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 φ 10 −8 eV where neutrino experiments lose sensitivity [26,27].
• Electron mass modulation: Even if we postulate that φ 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 ν , too small to be observed [26].

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.