Distorted neutrino oscillations from time varying cosmic fields

Cold, ultralight ( ≪ eV) bosonic fields can induce fast temporal variation in neutrino couplings, thereby distorting neutrino oscillations. In this paper, we exploit this effect to introduce a novel probe of neutrino time variation and dark matter at long-baseline experiments. We study several representative observables and find that current and future experiments, including DUNE and JUNO, are sensitive to a wide range of model parameters over many decades in mass reach and time-variation periodicity.


I. INTRODUCTION
Many popular extensions to the Standard Model (SM) feature ultralight bosonic fields with a present-day cosmic abundance. In these scenarios, the field becomes displaced in the early Universe and begins to oscillate about its minimum. As the Universe expands, the energy density in these oscillations eventually redshifts like nonrelativistic matter and can account for the dark matter (DM) in our Universe (for reviews, see Refs. [1,2]).
If this boson couples to SM fields, it induces periodic time variation in particle masses and couplings. Although there are strong laboratory bounds on time variation in charged particles [3][4][5][6][7] and many new ideas to improve these efforts [8][9][10][11][12][13], time variation in the neutrino sector is very poorly constrained by comparison.
In this paper, we propose a novel strategy to probe time variation in neutrino parameters at future longbaseline experiments and reinterpret existing experiments to extract bounds on this scenario. Earlier work has shown that slowly varying (∼ min −10 years) periodicity in neutrino parameters is constrained by the temporal stability of astrophysical neutrino fluxes [14]. 1 However, we find that qualitatively novel neutrino oscillation signatures can arise at long-baseline experiments even when the time variation is too fast (∼μs − min) to yield any observable flux modulation from astrophysical sources. Even modest couplings to the new field, over a wide range of masses, can significantly modify neutrino oscillation probabilities leading to distorted neutrino oscillations (DiNOs).
To illustrate this effect, consider an ultralight scalar ϕ with a Yukawa coupling to active neutrinos, where g ≪ 1 is a dimensionless coupling whose flavor indices have been suppressed. Normalizing to the scale of neutrino masses, we define a scale Λ ≡ m ν =g so that all neutrino bilinears can be written as where ϕ=Λ ∼ δm ν =m ν represents a fractional shift in the overall scale of neutrino masses and y encodes the flavor structure of this interaction. In the following discussion we consider the effects of this interaction on neutrino oscillations while remaining agnostic about the origin of this operator (see Appendix A for a discussion). In the Galactic halo, the local field value can be written as where ρ ⊙ ϕ ≤ ρ ⊙ DM ¼ 0.3 GeV cm −3 is the field's energy density, m ϕ is the mass, and v ∼ 10 −3 is the virial velocity.
In the presence of this ϕ background, the values of neutrino masses are modulated by the DiNO amplitude which we define in terms of the local ϕ density. Note that η ϕ can be sizable even if ϕ is only a small fraction of the DM. Note also that there are qualitative differences if ϕ affects solar parameters (θ 12 , Δm 2 21 ) or atmospheric parameters (θ 13 , θ 23 , Δm 2 31 ), as we will see below. Although the impact of such modulation on neutrino mass-squared differences and mixings depends on the matrix y, we analyze the effect of modifying one parameter in isolation assuming y ∼ Oð1Þ. Thus, each mass-squared difference can be written as where Δm 2 ij;0 is the undistorted value and ϕ evolves according to Eq. (3), with a corresponding dependence on η ϕ , where we have dropped the subdominant velocitydependent corrections. Similarly, a shift in the mixing angles can be written as where θ ij;0 is the undistorted mixing angle. Note that if the ϕνν interaction is flavor blind, then the rotations that diagonalize the vacuum mass matrix are unaffected and the mixing angles are ϕ independent. In a two-flavor neutrino formalism, the instantaneous vacuum probability for α → α survival is where L is the experiment baseline, E is the neutrino energy, and both θ and Δm 2 depend on ϕ through Eqs. (5) and (6). If the scalar oscillation period τ ϕ ≡ 2π=m ϕ is longer than the characteristic neutrino time of flight T ν , but shorter than the total experimental run time, then neutrinos emitted at different times will sample different values of ϕ over the course of a given experiment. In this regime, the effective oscillation probability is the ensemble average where for a given experimental baseline L ¼ c=T ν , there is a characteristic m ϕ below which standard oscillation probabilities can be distorted. If τ ϕ ≳ 10 minutes, the scalar oscillation can induce observable time variation in neutrino oscillation measurements (e.g., periodicity in the solar ν e flux) [14]. In this work, we study the opposite, highfrequency regime and find that scalars with τ ϕ ≪ min distort neutrino oscillation probabilities even if this time variation cannot be resolved. The effect of fast averaging is intrinsically different for neutrino mixing angles and mass-squared differences. For mixing angles, the net effect of averaging over ϕ induces a shift in the observed mixing angle relative to its undistorted value. Note that the observed sin 2 2θ after averaging can never be zero or maximal since, from Eqs. (3) and (6), we have where J 0 is a Bessel function of the first kind and, to quadratic order in η ϕ , the correction to the sin 2 2θðtÞ distribution is negative (positive) for maximal (minimal) mixing. Thus the observations of nonzero θ 13 [24] and nearly maximal θ 23 [25,26] already constrain the available parameter space. If the scalar primarily affects mass-squared differences (e.g., through flavor-blind Yukawa couplings), the time averaging has a more complicated functional dependence, FIG. 1. Example neutrino oscillation probabilities for a variety of scenarios at JUNO (top) and KamLAND (bottom). For both plots, the thick red curve is the standard oscillation prediction for each setup including the effect of energy resolution smearing (following the prescription in Appendix B). The green and turquoise curves also include the additional effect of ϕ-induced smearing separately distorting Δm 2 31 and Δm 2 21 , respectively. For KamLAND we have assumed a mean baseline between the nuclear reactors and the detector of hLi ¼ 180 km.
which leads to additional L=E smearing and distorts the functional form of oscillation probabilities, particularly near maxima and minima. Thus, the DiNO effect from Eq. (8) adds an irreducible smearing to the oscillation probability signal, similar to an experimental energy resolution, but at the probability level. This effect is shown in Figs. 1 and 2, which present both instantaneous and ϕ-averagedν e →ν e survival probabilities as a function of neutrino energy for JUNO [27] and KamLAND [28,29] as well as ν μ → ν μ and ν μ → ν e oscillation probabilities at the future experiment DUNE [30] (see Appendix B for a discussion of the signal calculation).

II. PHENOMENOLOGY
Although a detailed experimental analysis is outside the scope of this paper, 2 we present estimates of the experimental sensitivities of current and future neutrino experiments in terms of the ratio η ϕ taking into account possible new interpretations of oscillation parameters. In Fig. 3 we summarize our main results as bounds and projections in the m ϕ -η ϕ plane assuming separately that ϕ only affects solar (top panel) and atmospheric oscillations (bottom panel).
(1) Bounds from Δm 2 : The presence of a time-varying ϕ background influences the interpretation of existing neutrino oscillation parameters in Eq. (9). The excellent energy resolution of the KamLAND experiment of 6.4% at 1 MeV [31] constrains η ϕ ≲ 3.2% via Δm 2 21 smearing [e.g., induced by nonzero y 11 in Eq. (1)], as shown in Fig. 1, where we illustrate the effect of distorting Δm 2 31 [e.g., induced by nonzero y 33 in Eq. (1)] and Δm 2 21 on the disappearance probability ofν e . Note that the dip in the oscillation minimum for KamLAND is set by θ 12 , in analogy with Eq. (7). However, the heights of the maxima are affected by θ 13 -induced atmospheric oscillations (which are smeared out due to the experimental resolution), and thus are essentially fixed due to the current precision of θ 13 measurements. DiNOs therefore distinctively affect the height of the maxima in KamLAND.
Dedicated analyses are required to estimate the sensitivity to ϕ smearing in Δm 2 31 distortion at current experiments: at MINOS [32] the neutrino flux peaks away from the oscillation minimum; Super-K atmospheric data [33] has a nontrivial neutrino flux; and both NOνA [34] and T2K [35] are narrow-band beams which do not resolve the functional form of the oscillation probability.
(3) Future neutrino oscillation experiments: The future JUNO experiment aims to measure the neutrino hierarchy by observing the small-amplitude, highfrequency oscillations caused by Δm 2 31 in solar baseline oscillations. Its proposed energy resolution of 3% translates into a sensitivity of η ϕ ∼ 1.5%, as shown in Fig. 1, where we illustrate the disappearance probability Pðν e →ν e Þ in JUNO for ϕ-smearing effects on Δm 2 21 and Δm 2 31 . The effect on Δm 2 21 , although non-negligible, can be mimicked by adjusting the undistorted value of θ 12 , as we verified numerically. The DUNE experiment will also provide a constraint on η ϕ of the order of half of its energy resolution (around 15% at 1 GeV) due to Δm 2 31 smearing as shown in Fig. 2, where we illustrate the ϕ-smearing effect on Δm 2 31 on the disappearing probability Pðν μ → ν μ Þ as well as the appearance probability Pðν μ → ν e Þ. determinations of θ 23 by NOνA [26], T2K [36], DUNE [30], and T2HK [37] can bound η ϕ if the observed of θ 23 is maximal. (4) CMB and BBN: Since the average cosmological ϕ density redshifts as nonrelativistic matter, ρ ϕ ðzÞ ∝ ð1 þ zÞ 3 , where z is the redshift, the amplitude of neutrino mass modulation is much larger at earlier times. As observed in Ref. [14], a sharp increase in the overall scale of neutrino masses in the early Universe can exceed the Planck limit on their sum P i m ν i < 0.23 eV [38]. For concreteness, in the following discussion we assume a normal neutrino mass ordering and we assume the coupling in Eq. (1) remains valid in the early Universe; see Appendix D for a discussion of model dependence.
ϕ coupled to a heavy eigenstate: If the ϕ affects the heaviest neutrino, the absolute neutrino mass scale is modified and we demand the average ϕ-induced correction be no greater than an Oð1Þ effect at recombination, where z rec ≃ 1100 is the redshift at recombination, ρ ϕ ðz ¼ 0Þ is the present-day average cosmological ϕ density, and we have taken m ν ∼ 0.1 eV. Translating this into a bound on the local modulation amplitude in the solar neighborhood, we find where and ρ ⊙ ϕ is the local ϕ density and we have assumed the usual halo overdensity relation ρ ⊙ ϕ ∼ 10 5 ρ ϕ ðz ¼ 0Þ.
The ϕνν interaction can also affect neutrino free streaming during the cosmic microwave background (CMB) epoch via νν → νν scattering [39]. In our regime of interest m ϕ ≲ eV, this bound implies which is presented in Fig. 3 as a diagonal region. ϕ coupled to lighter eigenstates: In this regime, a similar argument applies, but ϕ now couples only to light neutrinos with m 1 ∼ m 2 ∼ ffiffiffiffiffiffiffiffiffiffiffi Δm 2 21 p ≈ 0.008 eV. This assumption translates into the requirement We note that the free-streaming bound from Eq. (14) is also reduced if ϕ couples predominantly to lightermass eigenstates. However, the bounds in Eqs. (13) and (15) apply only if ϕ accounts for all of the dark matter at recombination; if it only constitutes a subdominant fraction of the DM density, it need not be dynamical in the early Universe, so the constraint no longer applies. In this work, wherever the modulation effect exceeds these bounds, we will assume that ϕ oscillation begins after recombination.
Similar considerations apply if ϕ is dynamical during big bang nucleosynthesis (BBN) at z ∼ 10 8 ; however, if the field becomes dynamical in between the BBN and CMB epochs, then ϕ can still constitute all the dark matter as long as Eq. (12) or Eq. (15) are satisfied. See Appendix D for a discussion of the model dependence of these bounds. (5) Solar neutrino periodicity: The observed temporal stability of solar neutrino fluxes by Super-K imposes a tight bound on neutrino mass variation over 10 min-10 year time scales [14,40]. The period of ϕ-induced mass variation in our setup is so as long as m ϕ ≳ 10 −17 eV, the period is too short to have been observed at Super-K. If instead the period is within the Super-K time resolution, the η ϕ amplitude is bounded to be below Oð10%Þ [14]. (6) SN 1987A: Light weakly coupled particles can introduce anomalous energy loss in supernovae (SNe), thereby conflicting with observations from SN1987A [41,42]. Since ϕ only couples to neutrinos, the total energy released due to νν → ϕϕ annihilation for a SN temperature of T ¼ 30 MeV is where Δt is the SN blast duration. The observed energy released in SN1987A is approximately 10 51 erg, so to avoid an order-one correction, we demand Λ ≳ 50 keV (see Appendix C for details). However, the neutrino mass modulation amplitude scales as η ϕ , so fixing the magnitude of this effect and saturating the SN bound on Λ ∼ 50 keV implies a maximum scalar mass where the index in m ν;i refers to heaviest neutrino mass eigenstate to which ϕ couples. This bound defines the diagonal shaded region in Fig. 3 where we take Ω ϕ ¼ Ω DM . (7) Charged lepton mass variation: A time-varying neutrino mass can induce charged lepton mass variation through the loop-level diagram depicted in Fig. 4. However, the amplitude of this variation is estimated to be δm e m e ∼ G F m 2 ν 16π 2 η ϕ ∼ 10 −28 where we have taken m ν ¼ 0.1 eV. The existing limit on electron mass time variation is typically written as δm e =m e ≲ 6 × 10 −20 d m e ð10 −11 eV=m ϕ Þ, where d m e is a dimensionless Yukawa coupling constrained by tests of the equivalence principle (see, e.g., Ref. [10]) to be below 10 −2 for m ϕ < 10 −12 eV. Clearly, the effect induced by this loop is negligible. However, if Higgs-portal operators of the form ϕH † H and ϕ 2 H † H are allowed, there will be additional time variation in fermion masses due to Higgs mixing after electroweak symmetry breaking. The effect of these operators can be naturally suppressed, for instance, by extra dimensional locality (e.g., if ϕ and H live on different branes in a higher-dimensional model); however, the model-building details are beyond the scope of this work (see Appendix E for a discussion).

III. CONCLUSION
In this paper we have found that ultralight scalar dark matter can significantly distort neutrino oscillation phenomenology by introducing rapid time variation in neutrino masses and mixing angles. As neutrinos travel between source and detector, they traverse a spacetime-dependent scalar background, which must be averaged over when calculating flavor transition probabilities. If this effect primarily modifies mixing angles, it can shift the extracted parameters away from extremal values (maximal or minimal mixing), thereby requiring a reinterpretation of θ 13 and θ 23 measurements. Alternatively, if the scalar modulation primarily affects mass-squared differences, it can distort the functional form of oscillation probabilities (e.g., the locations of maxima/minima) akin to an additional energy resolution effect. Thus, experiments with good energy reconstruction (e.g., KamLAND, JUNO, and DUNE) are particularly sensitive to this scenario.
We found that existing neutrino oscillations experiments -including KamLAND, T2K, and Daya Bay-already exclude a large portion of the scalar parameter space, covering over a dozen orders of magnitude in mass reach; however, the applicability of these limits depends on whether the misaligned scalar modifies the parameters extracted at each experiment. Future experiments including JUNO and DUNE are poised to significantly extend this coverage. The main focus of our investigation is the Yukawa interaction gϕνν in Eq. (1), which introduces a quadratic correction to the scalar mass of order where we have taken m ν ¼ 0.1 eV. For representative values of this coupling, we have fractional mass modulation of order which is compatible with the modulation amplitudes that can be probed at long-baseline experiments (as shown in Fig. 3). However, in the absence of additional new physics the mass range below ∼10 −12 eV requires fine-tuning to avoid large corrections to m ϕ . However, other significant corrections can potentially destabilize the scalar mass in a UV-complete model of neutrino masses, for instance, if the interaction in Eq. (1) arises from a type-I seesaw model [43,44] with righthanded neutrinos N coupled to the scalar field via where L, H are the lepton and Higgs doublets, y ν and y are Yukawa couplings, and ϕ is sequestered from other SM fields. Integrating out the N field gives the familiar Weinberg operator ðLHÞ 2 =Λ with an additional ϕ interaction, which reduces to Eq. (1) after electroweak symmetry breaking. If ϕ enjoys an approximate shift symmetry, it can be naturally light and produced cosmologically through the misalignment mechanism, thereby constituting some fraction of the dark matter abundance.

APPENDIX B: DETAILS OF OSCILLATION AVERAGING
In this appendix we present the full procedure for calculating three-flavor neutrino oscillation observables with energy resolution smearing. For neutrino flavor transitions undistorted by ϕ, the two-flavor probability in Eq. (7) generalizes to where U is the Pontecorvo-Maki-Nakagawa-Sakata matrix (in the presence of matter effects [45,46]), α; β ¼ e, μ, τ, i; j ¼ 1-3, and the neutrino oscillation parameters are given in Ref. [47]. All of our numerical results adopt this ansatz for the standard oscillation probability for the integrand in Eq. (8). Although we include the Mikheyev-Smirnov-Wolfenstein effect [45,46] in our calculation of the oscillation probabilities, it does not differentially modify oscillations in the presence of the ϕ background.
To model the experimental energy resolution we evaluate the effective probability where f is a Gaussian distribution modeling the energy reconstruction, namely, p is the energy resolution at a reference energy E 0 , and E t;r are the true and reconstructed energies. Note that here f has been normalized such that R ∞ 0 dE t fðE t ; E r ; σ exp E Þ ¼ 1. Adding the effect of ϕ averaging to the oscillation probability yields which defines the procedure for obtaining the oscillation curves in Figs. 1 and 2.

APPENDIX C: SUPERNOVA BOUNDS
In this appendix we derive an order-of-magnitude bound on η ϕ from SN1987a by making some conservative assumptions on the emission rate of ϕ from the thermalized core of neutrinos inside the supernova with T ∼ 30 MeV. The cross section for ϕ emission via νν → ϕϕ annihilation is estimated to be where we adopt m ν ¼ 0.1 eV. The ϕ emission rate per neutrino in a thermal bath with number density n ν ¼ For a SN radius of R ∼ 10 km, we have approximately N ν ∼ 4πR 3 n ν =3 ∼ 4 × 10 54 neutrinos in the core, which implies an anomalous cooling rate of so the total energy loss per Δt ¼ 10 sec burst is To avoid adding an order-one correction to the energy loss from SN1987a, we demand ΔE ϕ ≲ 10 51 erg, which imposes the modest limit Λ ≳ 50 keV. Since most of the couplings we consider in this work are vastly greater than this bound, this bound imposes essentially no limit on the relevant parameter space, as shown in the bottom panel of Fig. 3. Note that in the top panel, the SN bound assumes that ϕ only affects solar oscillation parameters and, therefore, only couples to the lighter neutrino mass eigenstates.
Since the cross section scales as σðνν → ϕϕÞ ∝ m 2 ν and we now take m 1 ∼ m 2 ∼ 0.005 eV as a reference mass scale, the bound is correspondingly weaker.

APPENDIX D: MODEL DEPENDENCE OF CMB AND BBN BOUNDS
In Fig. 3, we show the naive bounds on this model from the sum of neutrino masses during the CMB epoch if ϕ couples predominantly to atmospheric (top) and solar (bottom) parameters. These bounds are dotted because they assume that the interaction in Eq. (1) remains valid out to arbitrarily large redshift, which need not be true.
For illustration, the toy model in Eq. (A3) realizes the ϕνν Yukawa interaction by integrating out N which couples directly to ϕ in the UV. When ϕ begins to oscillate, the N mass MðϕÞ ≡ Λ − yϕ acquires a timevarying component. For energies below MðϕÞ, the N can be integrated out and the effect on active neutrinos neutrinos becomes where m ν ¼ y 2 ν v 2 =Λ and we realize Eq. (1). However, the field ϕ grows with redshift ϕðzÞ ∝ ð1 þ zÞ 3=2 , so this approximation is no longer valid at early times where y 2 ν ðHLÞ 2 MðϕÞ ≃ y 2 ν v 2 yϕ νν; Λ ≪ yϕ; ðD2Þ and the physical neutrino mass scale ∝ 1=hjϕji decreases in the early Universe. While this toy model fails to address the fine-tuning issues that arise in models of ultralight scalars, it serves to illustrate the model dependence of the CMB and BBN bounds, which naively exclude sizable neutrino time variation if Eq. (1) is extrapolated back to arbitrarily early times. Another possibility is that the field ϕ is not dynamical at early times and only begins its oscillation after the CMB has formed-perhaps as a result of a hidden-sector phase transition. Such a scenario does not allow ϕ to be a large component of the present-day dark matter because DM is known to have been present during the CMB epoch. However, ϕ can still be a subcomponent of the DM and induce a large modulation effect on neutrino masses through the η ϕ dependence on other parameters in Eq. (4). Importantly, if the oscillation starts after BBN, but before the CMB has formed, ϕ can be all of the dark matter, but the CMB bounds in Fig. 3 apply in this case. It is not known whether a realistic model can elegantly UV complete such a scenario, but we leave this investigation for future work.

APPENDIX E: MODEL DEPENDENCE ϕ-HIGGS MIXING
In principle, ϕ-Higgs mixing θ hϕ could be induced by the loop-generated operator μϕH † H with sin θ hϕ ∼ μv=m 2 h . Although a model-independent estimate yields negligible θ hϕ , given the experimental constraints on the temporal modulation of the electron mass, one may suspect that a significant mixing may arise in full models. To exemplify this point, in the toy model of Eq. (A3) we can estimate μ ∼ y 2 ν yΛ=16π 2 . The fractional variation of the electron mass is given by δm e =m e ¼ sin θ hϕ hϕi=v ∼ μhϕi=m 2 h , which is translated into the DiNO amplitude constraint Here, d m e is a dimensionless Yukawa coupling constrained by tests of the equivalence principle to be below ∼0.01 [10]. In this toy-model example, as long as the neutrino mass generation scale Λ is below a GeV, the effect of ϕ-Higgs mixing is small.