Constraining Feeble Neutrino Interactions with Ultralight Dark Matter

If ultralight $(\ll$ eV), bosonic dark matter couples to right handed neutrinos, active neutrino masses and mixing angles depend on the ambient dark matter density. When the neutrino Majorana mass, induced by the dark matter background, is small compared to the Dirac mass, neutrinos are"pseudo-Dirac"fermions that undergo oscillations between nearly degenerate active and sterile states. We present a complete cosmological history for such a scenario and find severe limits from a variety of terrestrial and cosmological observables. For scalar masses in the"fuzzy"dark matter regime ($\sim 10^{-20}$ eV), these limits exclude couplings of order $10^{-30}$, corresponding to Yukawa interactions comparable to the gravitational force between neutrinos and surpassing equivalent limits on time variation in scalar-induced electron and proton couplings.


I. INTRODUCTION
Ultralight ( eV) bosonic dark matter (DM) φ is characterized by a macroscopic de-Broglie wavelength which exceeds the inter-particle separation, where v φ is the field veloctiy. If φ is misaligned from the minimum of quadratic potential, it oscillates as a classical field about this minimum according to φ( r, t) = 2ρ φ (t) m φ cos[m φ (t + v φ · r ) + ϕ( r )] , (2) and the corresponding energy density redshifts like nonrelativistic matter ρ φ ∝ a −3 , where a is the cosmic scale factor and ϕ is a possible phase. This phase may encode additional information about spatial variation -e.g. different φ domains arising from cosmological initial conditions 1 or the incoherent virialization in the Galaxy leading to variation on the scale of λ φ . If φ couples to Standard Model (SM) particles, their masses, spins, and coupling constants may inherit time dependence from Eq. (2). In the context of charged SM particles, there are many searches for such phenomena, which typically place very strong limits on the φ-SM interaction strength (see Ref. [1] for a review). By contrast, there are relatively few bounds on DM induced time dependence in the neutrino sector [2][3][4][5][6][7][8][9][10][11][12][13][14][15] and the * abhish@fnal.gov † krnjaicg@fnal.gov ‡ pmachado@fnal.gov § hramani@stanford.edu 1 e.g. due to oscillation starting at slightly different times in different Hubble patches when the field becomes dynamical at H ∼ m φ .
corresponding limits constrain comparatively large interaction strengths primarily via flavor oscillations.
In this paper we introduce the possibility that an ultralight DM candidate φ induces a time dependent Majorana mass for right-handed neutrinos where y φ is a coupling constant and the time dependence arises from Eq. (2). When this mass is small compared to the neutrino Dirac mass m D , the mass eigenstates form a pair of pseudo-Dirac fermions; one "active" ν a and one "sterile" ν s (per generation). These states oscillate into each other with a characteristic probability governed by their squared mass difference δm 2 [16] P (ν a → ν s ) = sin 2 (2θ) sin 2 δm 2 L 4E ν , where L is the baseline, E ν is the energy of the propagating neutrino, and θ ≈ π/4 is the mixing angle, which is near maximal in the pseudo-Dirac limit where m M m D . The Majorana mass governing δm 2 in Eq. (4) is time dependent, so the oscillation rate becomes sensitive to the dark matter density and to its cosmic evolution. This dependence can impact various terrestrial and cosmological observables. In this work we extract resultant bounds and impose extremely strong limits on the induced Majorana mass; depending on the value of m φ we find some limits on the coupling y φ corresponding to a φ mediated Yukawa force comparable to that of gravity.
This letter is organized as follows: in section II we present our theoretical framework, in section III we delineate the qualitatively different neutrino oscillation regimes that φ can induce, in section IV we compute the terrestrial bounds, in section V we determine the cosmological bounds on this scenario, and in VI we make some concluding remarks.

II. ULTRALIGHT DARK MATTER AND PSEUDO-DIRAC NEUTRINOS
We consider a scalar DM candidate φ with lepton number 2 and a cosmic abundance due to misalignment. In Weyl fermion notation, the Lagrangian in this scenario contains where y ν is the neutrino Yukawa coupling, H is the SM Higgs doublet, is the SM lepton doublet, and N is a SM neutral fermion, i.e. a right-handed neutrino. As we will see next, the presence of a feeble interaction between the scalar DM and the right-handed neutrino can have dramatic effects in neutrino oscillation phenomenology.
To understand the impact of φ on neutrino oscillations, it is instructive to describe the "1+1" scenario, in which there is only one generation of and N . For simplicity, assume that the active state here is an electron flavor neutrino. In the broken electroweak phase, the first term in Eq. (5) generates a Dirac mass of neutrinos. When the φ field is misaligned according to Eq. (2), the second term in Eq. (5) generates a Majorana mass for N , so we have for the Dirac and Majorana contributions, respectively, where v = 246 GeV is the Higgs vacuum expectation value. When m M m D , we obtain two nearly degenerate neutrino mass-squared eigenstates and we define δm 2 ≡ y φ m D 2ρ φ /m φ , where for the splitting between Weyl fermions as opposed to the usual ∆m 2 ij measured in oscillation experiments; here we have taken the local density to be ρ φ = 0.4 GeV/cm 3 [19]. The active-sterile mixing angle in this case is which is nearly maximal, θ ≈ π/4 in our full parameter space of interest. The diagonalization of the mass terms in Eq. (6) is obtained by defining the flavor fields in terms of the mass eigenstates approximately as The time evolution of a ν e state is given by which yields a ν e → ν e survival probability Using Eqs. (2) and (6) we obtain where we have absorbed the v φ dependence in Eq. (2) into the definition of ϕ for brevity. Thus, for a neutrino emitted at t = 0 and observed at some later time t, the resulting electron-neutrino disappearance probability can be written as where we have treated the phase ϕ as a constant over the propagation time.
Generalization for more neutrino flavors is straightforward and can be derived following similar steps as those taken in Ref. [20]. Moreover, to simplify the discussion on the constraints and because the electron-neutrino admixture in ν 3 is small (|U e3 | 1), when φ couples to ν 1 or ν 2 we will only consider nonstandard ν e disappearance, while when φ couples to ν 3 we will only consider nonstandard ν µ,τ disappearance; in both regimes, we treat the active-sterile oscillation in a two-flavor (active-sterile) framework.
As written in Eq. (2), the phase ϕ need not be constant over the full neutrino trajectory. Indeed, in the Galaxy, virialization will disrupt any constant phase value down to coherence patches of order the de-Broglie wavelength in Eq. (1). Thus, the full oscillation probability will depend crucially on the relative size of the oscillation baseline and this coherence scale.
Finally, we note that our scalar mass is not protected by any symmetry, so it will be sensitive to irreducible one-loop corrections of order from the interactions in Eq. (5). Thus, for small y φ in the pseudo-Dirac limit, this contribution does not destabilize the ultralight scalar mass, assuming no φ couplings to heavier states. 2

FIG. 1: Left:
Parameter space for φ coupled only to ν 1 or ν 2 mass eigenstates, which is predominantly constrained ν e oscillation bounds. Here we show bounds from CMB and BBN from Sec. V, Milky Way satellites from Sec. V B, scalar thermalization with neutrinos from Sec. V C, solar neutrino oscillations from Sec. IV A, and model independent limits on light DM from ultra faint dwarf (UFD) heating [17]. For points below the gray dotted line, the φ mediated force between right handed neutrinos is weaker than gravity, which is theoretically disfavored by the weak gravity conjecture [18] Right: same as the left panel, only φ now couples only to ν 3 , so the limits are driven by ν µ,τ oscillations for which the solar bound is subdominant to the atmospheric bound described in Sec. IV B.

III. NEUTRINO OSCILLATION REGIMES
In what follows, we will consider three distinct regimes for neutrino oscillations in the presence of the ultralight scalar fields. These regimes arrive from the relation between the neutrino oscillation length and the modulation frequency of φ or the coherence length that defines the overall phase ϕ. Instead of performing a detailed fit of experimental data, we will recast existing constraints on pseudo-Dirac neutrinos from Ref. [21] on our parameters of interest, y φ and m φ . As neutrinos are ultra-relativistic, we identify t = L in Eq. (14).
In the low frequency m φ L 1 regime, the neutrino encounters a constant phase ϕ domain over the course of its propagation.
Expanding Eq. (14) around m φ L → 0 yields an oscillation probability We can interpret this oscillation probability as follows.
Since the period of the field φ is too long compared to the neutrino time-of-flight, the pseudo-Dirac mass splitting induced by the field is constant for each neutrino. Nevertheless, as an experiment collects data, the mass splitting will evolve as the field φ displays time modulation. In practice, several neutrino experiments have a high enough rate of events to observe time modulation of oscillation probabilities with periods as short as 10 minutes, which would correspond to m φ ∼ 10 −18 eV [3][4][5]12].
Since any small pseudo-Dirac mass splitting leads to maximal mixing, time modulation of neutrino oscillation probabilities due to φ modulation would lead to large, observable effects on oscillation data. Both constant and time dependent pseudo-Dirac mass splittings would be ruled out by neutrino data if observed, and can be used to set limits on the coupling strength y φ for a given m φ . Since sterile neutrino oscillation constraints are typically reported as bounds on δm 2 , we can define an effective mass-squared δm 2 eff by equating the arguments of Eq. (4) and Eq. (16) to obtain assuming cos ϕ ∼ 1. Recasting pseudo-Dirac neutrino limits on δm 2 in Eq. (17) allows to constrain where we have identified δm 2 eff with the constrained value δm 2 lim . Note that, depending on context, ρ φ can either be the cosmological DM density at a given cosmic era or the present day local density.
When the φ modulating frequency is high, m φ L 1, the accumulated phase due to propagation is sufficient to induce many modulation cycles on φ over the neutrino trajectory. However, as long as m φ v φ L 1, the neutrino time-of-flight is shorter than separation time of φ wave packets. A neutrino propagating in this regime will encounter the same value of ϕ across its trajectory, that is, the modulation of φ throughtout the neutrino trajectory is coherent. Without loss of generality, we can set the initial condition ϕ = 0. The effective oscillation probability in this regime is given by a time-average of Eq. (14) over the duration of propagation where we have assumed that ρ φ does not change appreciably across the baseline. In this intermediate regime we repeat the argument leading up to Eq. (18) and constrain C. Random walk: Finally, in the m φ v φ L 1 regime, the neutrino timeof-flight is longer than the wave packet separation of φ, so the neutrino traverses a random sample of φ field patches, each with a different phase ϕ. Along this trajectory, there are approximately m φ v φ L patches whose contributions add incoherently, so the effective phase can be approximated by ϕ eff ∼ m φ v φ L, assuming random distribution of phases ϕ and the phase averaged probability can be written The corresponding limit on the coupling reads

IV. TERRESTRIAL OBSERVABLES
We now consider various terrestrial bounds on pseudo-Dirac neutrinos in the context of our scenario. Depending on the values of y φ and m φ , a particular constraint can apply in any of the three regimes outlined in Sec. III, so the relationship between y φ and m φ will differ in each case

A. Solar Neutrinos
For electron neutrinos, the pseudo-Dirac splitting can be constrained by measurements of the solar neutrino flux. In the standard three neutrino oscillations paradigm, 8 B neutrinos undergo an adiabatic evolution due to large matter effects in the Sun [22]. This leads to a survival probability P (ν e → ν e ) c 4 13 s 2 12 + s 4 13 0.3, where s ij and c ij are the sines and cosines of mixing angle θ ij . Low energy solar neutrinos, on the other hand, are not affected by matter effects, and thus P (ν e → ν e ) c 4 13 (c 4 12 + s 4 12 ) + s 4 13 0.55. These probabilities describe well experimental data [23][24][25][26][27][28][29][30][31]. This can be used to extract an order of magnitude bound on the splitting in our scenario by demanding that this prediction is not affected by an order one amount. Here we use ρ φ and v φ ≈ 10 −3 with L = 1.5 × 10 8 km, which requires and can be translated into a bound on our model parameters using the relations in Sec. III, where the appropriate regime is determined by m φ . Since solar neutrinos are essentially almost pure ν 2 or incoherent ν e , and ν e has but a small admixture ν 3 mass eigenstate, the corresponding solar limit on the y φ applies only to the right-handed partners N 1,2 . Applying the solar limit from Eq. (23) to the three regimes from Sec. III A, and assuming that the Dirac mass of ν 1 satisfies m 2 D = ∆m 2 21 = 7.4 × 10 −5 eV 2 [32], we find the following constraints.
For m φ 10 −18 eV, solar neutrinos are in the constant φ regime, so from Eq. (18), we find a limit , for m φ < 10 −18 eV. (24) Note that if m φ 10 −24 eV, the period of φ is larger than 20 years, and the observation of pseudo-Dirac mass splittings become dependent of the initial condition ϕ. These results are plotted in the left panel of in Fig. 1, which shows constraints on φ coupled only to ν 1 or ν 2 , corresponding to ν e oscillations measurements.

B. Atmospheric neutrinos
Measurements of the atmospheric neutrinos can place limits on the φ coupling to ν 3 since muon neutrinos have a large admixture of the ν 3 eigenstate. If ν 3 is split in a pseudo-Dirac pair, a substantial deficit of atmospheric ν µ flux would be observed, contradicting experimental data [33][34][35][36]. The characteristic atmospheric baseline is the Earth's radius L ≈ 6000 km, and the Super-Kamiokande constraint on constant pseudo-Dirac mass splittings is [21] δm 2 lim < 10 −4 eV 2 , which translates into a bound on the φ coupling to which is valid for m φ 3×10 −14 eV. For larger φ masses in the modulating φ regime of Sec. III B, we impose the limit These bounds are presented in the orange shaded region of Fig. 1 (right panel). Note that for m φ 10 −10 eV, atmospheric oscillations are in the long baseline regime of Sec. III C, but the bound in this mass range is subdominant to other constraints in Fig. 1 and is not shown. In principle atmospheric neutrinos also bound the φ coupling to ν 1,2 , but solar constraints are stronger.

A. Scalar Evolution
Throughout our analysis, we assume that the φ potential can be written as where, in principle, the size of the quartic is unconstrained by symmetry arguments and can take on any value. However, there is an irreducible contribution to the quartic interaction generated through a Coleman-Weinberg interaction with the neutrinos which is always present in the absence of fine tuning. In an expanding Friedmann-Robertson-Walker universe, φ satisfies the equation of motion where the prime denotes a derivative with respect to φ. If φ is initially displaced from its minimum, it is frozen by Hubble friction until Hφ ∼ V , so if the mass term dominates the potential, V ∼ m 2 φ φ, the field becomes dynamical when m φ ∼ H and oscillate about φ = 0 while redshifting like non-relativistic matter ρ φ ∝ a −3 .
In this scenario, the initial misalignment amplitude φ i during inflation sets the DM abundance. Since m φ /H i 1, where H i is the Hubble scale during inflation, φ generates isocurvature perturbations, which are constrained by CMB measurements [37]. However, as long as φ i /H i 1, isocurvature perturbations can be parametrically suppressed, so for a given H i , a suitable choice of φ i can account for the DM abundance while being safe from this constraint. Furthermore, since m φ /H i 1, φ evolution is predominantly classical during inflation, so the initial amplitude φ i remains a free parameter throughout our analysis and can be chosen to yield the observed DM density [38].

B. Milky Way Satellites
In order for φ to account for the full DM abundance, it must redshift like non-relativistic matter (ρ φ ∝ a −3 ) in the early universe, starting at least at matter-radiation equality at a critical redshift z ∼ 10 6 , corresponding to a temperature T ∼ keV [39]. Since the φN N interaction in Eq. (5) yields an irreducible quartic scalar self-interaction term, we need to ensure that the ρ φ is not dominated by the quartic contribution at T eq ; otherwise it would redshift like radiation ρ φ ∝ a −4 (or faster if even higher polynomial terms dominate instead) [40]. Avoiding this fate requires where φ ≡ φ(T ). Using the scaling in Eq. (2), we find where ρ c is the present day critical density. The inequality in Eq. (34) defines the gray shaded regions in Fig. 1 where this effect would erase the Milky Way satellites already observed. However, note that this bound is modeldependent as it can be evaded if φ is only a small fraction of the total dark matter abundance, in which case it need not redshift like nonrelativistic matter at early (or even later) times.

C. Avoiding Thermalization
The Yukawa interaction in Eq. (5) enables φν → φν scattering which can bring φ particles in the misaligned condensate into equilibrium with neutrinos if the rate ever exceeds Hubble expansion. Since active neutrinos don't couple directly to φ, the cross section for this process requires two Dirac mass insertions and scales as σ ∼ y 4 φ m 2 D /T 4 . Furthermore, since both φ and the neutrinos are ultralight, the scattering rate scales as Γ ∼ nσ ∝ T −1 , up until T ∼ m D corresponding to the maximum rate relative to Hubble. Demanding that less than 3.8% of the φ population in the condensate is upscattered and becomes relativistic [41] at this temperature implies where g = 3.36. This bound is shown in Fig. 1 as the green shaded region.

D. CMB/BBN
In this section, we investigate the effects of the scalar field in the early universe, specifically active to sterile oscillations, which increase the effective number of neutrino species, ∆N eff .

Cosmological Field Density
If the relic density was set by the misalignment mechanism, then the DM density grows as T 3 and remains as non-relativistic DM until the temperature T H when m φ = 3H(T H ), where H is the Hubble parameter. Above this temperature, the field is constant due to Hubble friction and only contributes to the vacuum energy, so we have where T 0 = 2.72 K is the present day CMB temperature and ρ φ (T 0 ) = Ω dm ρ c is the cosmological DM density, which is related to the local overdensity via ρ φ (T 0 ) ≈ 3 × 10 −6 ρ φ . In what follows, we insert Eq. (36) into m M (T ) = y φ φ(T )/2 using Eq.
(2) to model the Majorana mass as function of cosmic temperature.

Cosmological Sterile Neutrino Production
To compute the early universe sterile neutrino yield, it is convenient to define r β as the ratio of active/sterile momentum moments where angular brackets · · · s,a denote a thermal average over the sterile and active distributions, respectively. Generalizing the formalism of Ref. [42], r β satisfies the Boltzmann equation where Γ = 7π 24 G 2 F pT 4 , and the mixing angle is where the effective matter potential for each flavor a = e, µ, τ can be written as where η = (n L − nL)/n γ = 6 × 10 −10 is the lepton asymmetry, C 1 = 0.95, C e 2 ≈ 0.61, C µ,τ 2 ≈ 0.17, and the ± refer to neutrinos and antineutrinos [43]. Here the vacuum mixing angle θ 0 in Eq. (39) is φ dependent where we have used Eqs. (6) and (9). Note that the first two moments of the active neutrino distribution yield the number and energy densities ( p 0 a = n a , p 1 a = ρ a ) In the following subsections, we derive detailed ∆N eff limits from BBN and CMB based on ν a → ν s oscillations around T ≈ MeV; later oscillations do not affect light element yields or the Hubble rate. The oscillation probability is maximized when the argument of Eq. (4) is order one, implying where we have used δm 2 = m D m M and m M ∼ y φ φ from Eq. (6), φ ∼ √ ρ φ /m φ ∝ T 3/2 from Eq. (2), and approximated L ∼ (G 2 F T 5 ) −1 as the neutrino mean-free-path, setting T = MeV throughout. Thus, the blue-shifted DM density at BBN greatly enhances the neutrino Majorana mass and yields on order-one oscillation probability for extremely feeble couplings y φ ∼ O(10 −29 ). Note that in our numerical study below, we use the full temperature dependence from Eq. (36) which also accounts for the Hubble damped regime when T > T H and is relevant for the smallest values of m φ we consider.
However, from Eq. (36), for sufficiently large values of y φ and ρ φ , m M > m D so neutrinos are no longer pseudo-Dirac fermions at high temperatures. In this regime, ν a → ν s oscillations are sharply suppressed as θ 0 → π/2 in Eq. (41), so there is a ceiling to the couplings that can be probed in the early universe; this effect yields concave regions for the BBN/CMB regions in Fig. 1.
In terms of r in Eq. (38), sterile production via a flavor oscillations predicts where T νe dec ≈ 3.2 MeV and T νµ,ντ dec ≈ 5.34 MeV are the temperatures of ν aνa → e + e − chemical decoupling [43]. Assuming the ΛCDM cosmological model, the Planck collaboration constraints ∆N eff 0.28 [37] and we show this constraint in Fig. 1 as the blue shaded region alongside projections from future measurements with CMB-S4 [44].

BBN ∆N eff Limit
A nonzero ∆N eff from sterile production also yields a larger initial neutron/proton fraction at the onset of BBN, which increases the primordial helium fraction. As in Eq. (43), for φ coupled to ν µ,τ , the effect on BBN arises purely from the expansion rate via where r 1 is the solution to Eq. (38) with β = 1 evaluated at decoupling, assuming no initial population of steriles. The blue contour of Fig. 1 (right panel) shows parameter space where ∆N BBN eff > 0.5 [45,46] for φ coupled to the ν 3 mass eigenstate, implying oscillations from ν µ and ν τ flavor states.
However, for ν e → ν s oscillations, there are two distinct effects that impact the n/p ratio: oscillations before ν e chemical decoupling at T νe dec ≈ 3.2 MeV change the expansion rate as above, and oscillations after decoupling deplete the ν e density. Both effects can be captured with a shift in the effective Fermi constant via and a simultaneous shift in g via where g ,SM = 10.75 during BBN and T nuc ≈ 0.8 MeV is the temperature at which nucleon inter-conversion freezes out in the SM. Note that p 2 a ∝ T 5 , which sets the weak scattering rate Γ ∼ G 2 F T 5 , so r 2 = p 2 s / p 2 s yields the fractional departure from this rate.
We can economically capture both effects with an equivalent ∆N BBN eff [47]

VI. CONCLUSIONS
In this letter we have presented the first cosmologically viable model in which neutrino masses acquire time dependence through their coupling to ultralight dark matter. In our scenario, the DM interaction sets the right-handed neutrino Majorana mass and neutrinos are pseudo-Dirac fermions with small mass splittings between active and sterile states. Since in the pseudo-Dirac regime the mixing angle between active and sterile is maximal, we extract limits on ultra feeble Yukawa couplings between DM and right-handed neutrinos, constraining values of order y φ ∼ 10 −30 for m φ ∼ 10 −19 eV in the "fuzzy" DM regime [48]; for such small couplings the φ mediated Yukawa force between right-handed neutrinos is comparable to that of gravity.
Throughout our analysis, we have emphasized bounds from solar and atmospheric neutrino oscillations, large scale structure, and the CMB/BBN eras. However, additional limits on this scenario may also be extracted by studying cosmic ray propagation [20] or diffuse supernova background neutrinos [49,50], which we leave for future work.