Cosmology-friendly time-varying neutrino masses via the sterile neutrino portal

We investigate a consistent scenario of time-varying neutrino masses, and discuss its impact on cosmology, beta decay, and neutrino oscillation experiments. Such time-varying masses are assumed to be generated by the coupling between a sterile neutrino and an ultralight scalar field, which in turn affects the light neutrinos by mixing. Besides, the scalar could act as an ultralight dark matter candidate. We demonstrate how various cosmological bounds, such as those coming from Big Bang nucleosynthesis, the cosmic microwave background, as well as large scale structures, can be evaded in this model. This scenario can be further constrained using multiple terrestrial experiments. In particular, for beta-decay experiments like KATRIN, non-trivial distortions to the electron spectrum can be induced, even when time-variation is fast and gets averaged out. Furthermore, the presence of time-varying masses of sterile neutrinos will alter the interpretation of light sterile neutrino parameter space in the context of the reactor and gallium anomalies. In addition, we also study the impact of such time-varying neutrino masses on results from the BEST collaboration, which have recently strengthened the gallium anomaly. If confirmed, we find that the time-varying neutrino mass hypothesis could give a better fit to the recent BEST data.

It remains a mystery as to why the absolute scale of neutrino masses is more than six orders of magnitude smaller than its charged lepton partner. This is often ascribed to some new physics at very high energy scales in the spirit of the seesaw mechanism [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36], whose experimental test is extremely challenging. Alternatively, an interesting possibility is to attribute the origin of neutrino masses to the dark sectors, for example, dark energy (DE) [37,38] and dark matter (DM) [39]. As the persistent direct detection searches of weakly interacting massive particles (WIMPs) come out with null signals so far, there has been increasing attention to other DM candidates, such as the ultralight DM [40,41]. The ultralight DM candidate, due to its tiny mass m < O(eV), acts as a delocalized classical-number field. Such a scenario could be responsible for addressing generic studies of varying physical constants. One particular class of ultralight DM, the fuzzy dark matter with a mass m φ 10 −22 eV, can also help to alleviate the small-scale structure issues existing between cosmological observations and simulations [42][43][44][45].
Such an ultralight scalar field, coupled to neutrinos, can give rise to intriguing phenomenological consequences. Neutrinos coupled to the scalar will naturally get a contribution to their masses, in analogy to the Higgs mechanism. As the classical DM field oscillates, the generated mass term for the neutrinos becomes time-varying. There have been various attempts of studying the phenomenology of a coupling between neutrinos and the ultralight field in the literature . A simple realization of this possibility is to couple a classical scalar field (Φ) to neutrinos via the operator ννΦ, which generates a time-varying neutrino mass term. In a UV-complete model, this can be achieved by directly coupling Φ to the lepton doublet L ≡ (ν, e) L , as in the type-II seesaw [32][33][34][35]. However, stringent constraints arise because of the accurate determination of electron properties.
A more feasible way to couple light neutrinos to the scalar field is by mixing with gauge singlets, such as a right-handed sterile neutrino N . It is also natural to have a considerable Yukawa coupling between the sterile component and the scalar field. Motivated by this, our working Lagrangian consistent with gauge symmetries is described as, where Φ represents the ultralight DM field, h is the SM Higgs, the active Majorana neutrino mass term m ν ≡ κv 2 /2 and the Dirac mass m D = y D v/ √ 2 will be generated after Higgs takes the vacuum expectation value v/ √ 2 ≡ h = 174 GeV, and m N is the Majorana mass term of the sterile neutrino. On top of that, the Yukawa interaction gΦN c N with a coupling constant g generates a time-varying mass term to the sterile neutrino. In addition, there might be a dimension-five effective interaction with a coupling constant y and the cutoff scale Λ [49,72,73], similar to the Weinberg operator generating light neutrino masses.
In general, the Majorana mass terms m ν and m N can be vanishing by imposing additional symmetries such as lepton number conservation (then the remaining Lagrangian will look similar to the singlet Majoron model [74][75][76][77][78][79][80]). In the most economical case, one can even generate small neutrino masses, with a minimal interaction form y D L hN + gΦN c N , but the number of sterile neutrinos should be extended to at least two in order to explain the oscillation data. Note that even though we assume Majorana neutrinos in this work, the generalization to Dirac neutrinos is not difficult. An alternative scenario is the one in which right-handed neutrinos have no vacuum mass (m N = 0) but get a tiny Majorana mass from their coupling to the scalar field. This can give rise to ultralight scalar-induced Pseudo-Dirac neutrinos [71]. The framework explored in this work is different from the earlier discussions about the mass-varying neutrinos, where the tiny neutrino mass is hypothetically generated from the coupling to the dark energy (e.g., acceleron [38], quintessence [81], etc.). In that case, neutrino masses are constant over any observable laboratory time scales and one cannot distinguish the induced mass from the vacuum mass by looking for the time-varying signals. Late-time neutrino mass generation can also arise in models having a late-time cosmic phase transition [82], or models having an anomaly due to a gravitational−θ term [83]. In fact, Ref. [82] analyzed redshift dependent neutrino masses in the light of recent Planck data, and reported a preference for models of late-time neutrino masses. This was further analyzed in a follow-up work [84], where additional data from Type-Ia supernovae, as well as structure formation was used. More recently, it was demonstrated that a detection of the diffuse supernova neutrino background could shed light on whether neutrino masses turned on at later redshifts [85].
With the renormalizable terms in Eq. (1), the equation of motion of the scalar in our local galaxy is where the addition of Hubble dilution term, 3HΦ with H being the Hubble expansion rate, is necessary if we consider the scalar evolution over cosmological time scales in the expanding Universe. In the absence of the source term in the right-hand side, the scalar field evolves freely in the Universe after production. Because the scalar is assumed to be produced coherently 1 and its occupation number is very high, it is more appropriate to consider Φ as a classical field instead of a quantum state. Neglecting possible spatial variations (arising from structure formation), which are usually suppressed by the DM velocity v φ ∼ 10 −3 in our Milky Way, the field evolution simply follows, Here the time t is calibrated such that Φ(0) = 0, and φ = √ 2ρ/m φ denotes the field strength, which is approximately φ = 2.15 × 10 15 eV · (10 −18 eV/m φ ) with the dark matter energy density ρ ≈ 0.3 GeV · cm −3 in our local galaxy. It is worthwhile to setup the magnitude of the coupling g by noting that gφ ≈ 2.15 eV ·(g/10 −15 )·(10 −18 eV/m φ ). Unless otherwise specified, we use the capital Φ to denote the complete time-varying field and φ for its amplitude.
Through the Yukawa coupling to Φ, the sterile neutrino develops an effective time-varying mass term. One may worry about energy-momentum conservation within such a setup. Let us consider a closed system formed only by the scalar and the sterile neutrino. The total Lagrangian explicitly possesses a temporal translation symmetry, so the energy of the entire system must be conserved according to Noether's theorem. Furthermore, because the Lagrangian with only the neutrino part preserves spatial translation symmetry, the neutrino momentum must be invariant under the time-varying scalar potential. But the neutrino energy might be perturbed by the scalar field. When the neutrino energy evolves with E 2 N = p 2 N + (gφ sin m φ t) 2 , the energy of the scalar system E φ will also change accordingly due to the feedback effect from N source term as in Eq. (2), such that E N + E φ = const.
In the absence of the scalar potential, diagonalization of Eq. (1) leads to light mass eigenstates ν i (for i = 1, 2, 3) and a heavy one ν 4 . If the vacuum mass m 4 is too large compared to the potential gφ, only the effective coupling by mixing g ij ν Φν i ν j will be relevant at low energies, with i and j being the indices of light neutrino mass eigenstates. However, as has been realized, such a scenario faces inevitable constraints from the observation of cosmic microwave background (CMB) [46][47][48]68]. The growth of φ with the redshift will render cosmic neutrinos non-relativistic in the early Universe, which is in contradiction with the free-streaming property of relativistic neutrinos from the CMB observation. Sensitivities of most of the laboratory searches are not even comparable to the CMB constraint in the order of magnitude. The situation is different in the other regime, when the sterile neutrino mass m 4 is smaller than gφ. In this case, as we go back to the dense early Universe (gφ m 4 ), the potential induced for light neutrinos ν i will be suppressed by m 2 D /(gφ). This is very similar to the seesaw mechanism which naturally generates the small neutrino masses with the suppression of heavy degrees of freedom. The details for this observation are given in Appendix A. The sterile neutrino can be as light as O(eV), which is relevant for short baseline anomalies [86][87][88][89][90][91]. A general issue with the light sterile neutrino is the possibility of its thermalization in the early Universe, leading to an increase in ∆N eff [92][93][94][95]. In fact, the parameter space suggested by the short-baseline anomalies is almost ruled out by the ∆N eff constraints from Big Bang nucleosynthesis (BBN) and CMB [96,97]. As has also been noticed previously [49,56], the framework with a time-varying DM field coupled to sterile neutrino can provide a bonus to suppress the production of light sterile neutrinos in the early Universe.
In this work, we systematically explore the consequences of a time-varying mass of sterile neutrinos in the early Universe, outlining the impact on BBN and CMB. Furthermore, we discuss the possibility of current generation beta decay experiments like KATRIN to probe the time-varying sterile neutrino hypothesis. Finally, we also highlight how short-baseline neutrino experiments fare in light of such a hypothesis. The mass range of the scalar relevant for various probes is illustrated in Fig. 1, including the neutrino oscillation experiment DUNE, the beta-decay experiment KATRIN, the short-baseline experiment BEST, as well as the equivalence time scales in the early Universe. The corresponding time scales are indicated on the top axis. It is worth mentioning that black hole superradiance constraints exists in several mass ranges for the ultralight scalar mass [98,99].
The structure of the rest of the work is as follows. In Sec. II, we separate the time-varying masses via the sterile neutrino portal into two different scenarios. In Sec. III, we investigate the cosmological consequences if the ultralight DM interacts with the sterile neutrino, including both heavy and light sterile neutrino scenarios. In Sec. IV, we explore the distortion effect induced by the time-varying potential in beta-decay spectrum, taking KATRIN experiment as an example. We proceed in Sec. V to discuss the impact of a time-varying light sterile neutrino on the short-baseline experiments. Finally, we discuss our results, and conclude in Sec. VI.

II. GENERATING TIME-VARYING ACTIVE NEUTRINO MASSES FROM STERILE NEUTRINOS
In the limit that the sterile neutrino is very heavy compared to the scalar potential, i.e., m 4 gφ, the light neutrinos will receive an effective mass in addition to the original vacuum one (see Appendix A), The active-sterile mixing angle θ can simply be absorbed by redefining g ν ≡ sin 2 θ · g such that the neutrino mass correction reads as g ν φ sin m φ t. In such a case, it is technically indistinguishable whether active neutrinos couple directly to the scalar field or by mixing with the sterile neutrino, and the discussion will be reduced to the usual scenario explored in the literature [46-48, 50-55, 57, 58, 60-70].
This potential is highly testable in neutrino oscillation experiments, provided the DM oscillation cycle, whose period is given by 2π/m φ , is at least of the same order as the neutrino time of flight (i.e., the DM field is not rapidly oscillating). If the DM field is oscillating with a period of the order of the experimental running time, i.e. few days to years, the neutrino mass eigenstates can develop oscillation phases continuously with constant effective masses during the flight [46,48,62]. This would manifest as a time modulation of the signal. For shorter DM cycles, the modulation effect will be averaged, but not vanishing. This averaged effect imprints non-trivial distortions to the original neutrino oscillation probability as a function of the neutrino energy. Moreover, it was recently pointed out that even in the rapidly oscillating (dynamical and beyond) regime, the neutrino flavor transition induced by the scalar field is still possible but with a decreased impact [62].
Even though the time-varying analysis of neutrino oscillation experiments itself is interesting and very rich in phenomenology, the available model parameters are severely constrained from cosmology. In fact, a reasonable cosmological scenario, without fine-tuning the DM evolution, might rule out all the parameter space to which neutrino oscillations are sensitive. Irrespective of whether the scalar field accommodates all the DM abundance or not, its field strength averaged over space in the early Universe as a function of redshift z will read as up to the time when the Hubble expansion rate becomes comparable to the scalar mass, i.e., H ≈ m φ . To see how an experimentally testable coupling is disfavored by cosmology, we take g ν φ = ∆m 2 21 ≈ 8.7 × 10 −3 eV, where the local DM overdensity is approximately φ 2 ≈ 10 5 φ(0) 2 [48,49]. During the era of matter-radiation equality (z eq ≈ 3000), we obtain g ν φ(z eq ) = 4.5 eV, which exceeds the neutrino temperature at that time T ν ≈ 0.5 eV by almost one order of magnitude. As a result, neutrinos become nonrelativistic, and do not free-stream at the speed of light before recombination, thereby spoiling CMB observations. The above conclusion is actually rather conservative and requires DM particles to be populated just before the matter-radiation equality. On the other hand, if the DM production is not fine-tuned and takes place early before BBN at T ν ≈ 1 MeV (corresponding to z BBN ≈ 6 × 10 9 ) such as the misalignment production of QCD axions, a severe limit g ν φ < 7 × 10 −7 eV can be obtained by conservatively requiring g ν φ(z BBN ) < 1 MeV. The allowed tiny local effective time-varying mass clearly rules out all the possibility of realistic laboratory searches.
The above picture changes if m 4 < gφ(z) during the photon decoupling and/or neutrino decoupling era. In Fig. 2, we show the lighter effective neutrino mass m L ≡ Min{ m 1 , m 4 } within a DM oscillation period for different gφ. It is important to note that when the gφ sin m φ t = −m 4 − m 1 is satisfied, it is possible for m L to swap from m 1 to m 4 due to the resonance encountered. We choose the lighter eigenstate m L , because the lighter neutrino is defined to have a dominant overlap with active neutrino, providing the active-sterile mixing is small. For the case of gφ = 100 m 4 , the lighter mass m L is given by m 1 for gΦ > −m 1 − m 4 and m 4 for gΦ < −m 1 − m 4 . Previous analyses simply assume an ad hoc sinusoidal variation, e.g., Eq. (4), in the active neutrino mass or mixing parameters, which need not necessarily to be the case. As shown in Fig. 2 and is clear from Eq. (A2), the variation of light neutrino parameters can have a more complicated form, especially when gφ is large compared to the sterile neutrino mass. In the limit of 0 eV < m 4 gΦ, the mass-variation is instead given by (see Appendix A for a derivation) In the extreme case of m 4 gΦ and sin θ 1, a constant shift ∼ m 4 sin 2 θ with a time variation with magnitude ∼ m 2 4 sin 2 θ/(gΦ) in the neutrino mass will be induced. In this case, we do not expect a large effective neutrino mass in the very early Universe as long as m 4 sin 2 θ is chosen to be small.

A. Neutrino decoupling and Big Bang nucleosynthesis
The measurements of primordial helium and deuterium abundances agree very well with the theoretical predictions, assuming standard electroweak interactions with massless neutrinos. The weak interaction rates are expected to be altered if neutrinos are very massive at the redshift of decoupling, e.g., if m i (z BBN ) ∼ T ν . A severe constraint can be imposed on gφ if ultralight DM particles are assumed to be populated before the BBN era. This can translate into a constraint on the effective neutrino mass for increasing values of gφ as described by Eq. (5). Fig. 3 shows the average of the lighter neutrino mass m L as a function of the potential gφ. The vacuum mass and mixing angle have been taken as m 1 = 0 eV, m 4 = 3 eV and θ = π/12. We observe that starting from very small field strength with gφ < m 4 , the average mass m L increases linearly with gφ. At the point around gφ = m 4 , the average mass m L stops growing as expected from Eq. (A11). The turning point around gφ = m 4 is the key to alleviate the tension between testable time-varying signals and BBN.The maximum of m L , which is given by m L max ≈ m 4 sin 2 θ, is under control no matter how gφ changes in the early Universe. Hence, in order not to spoil the BBN observations with a large potential, the sterile neutrino parameter should stay within the range m 4 sin 2 θ 1 MeV. Light sterile neutrinos (lower than the neutrino decoupling temperature, i.e., m 4 < 1 MeV) have an additional risk of thermalization and increasing the amount of extra radiation in the early Universe, measured by ∆N eff [100,101]. It has been noticed that such a risk can be evaded by introducing secret interactions among sterile neutrinos. An effective potential suppressing the active-sterile mixing will hence be induced in the presence of just a small background of sterile neutrinos [102][103][104][105][106]. The spirit is similar in our scenario. In the presence of a DM potential, sterile neutrino production will be strongly suppressed. In the following, we numerically investigate this possibility by directly solving for ∆N eff .
In a simplified two-neutrino setup, the evolution equation of the sterile neutrino phase-space distribution f N is given by [107] df where T ν is the neutrino temperature, f νa is the active neutrino phase-space distribution, and Γ ∝ T 5 ν /m 4 W encodes the net neutrino interaction rate. Here, θ m is the effective mixing angle in matter, where ∆ = ( m 2 4 − m 2 1 )/(2E) gives the vacuum oscillation frequency arising from the mass difference between active and sterile species, and denotes the time-averaging operation. Here V T ∝ (T 5 ν /m 4 W ) is a measure of the forward scattering thermal potential experienced by the neutrinos. The scattering term in the denominator encodes the Quantum Zeno effect, where the flavor conversion is suppressed for a large scattering rate. Furthermore, the self-scattering process 2N → 2N may freeze-in and become relevant post BBN, leading to a scattering-induced decoherent population of sterile neutrinos [104]. However, note that this effect is proportional to the coupling g, and is sub-dominant in our case with an extremely small coupling.
In the limit where f νa can be approximated by a Fermi-Dirac function (or any generic function of p/T ν ), Eq. (8) can be solved approximately to obtain the contribution of sterile neutrinos around the time of BBN [108], where g * is the effective number of relativistic degrees of freedom. In Fig. 4, we illustrate the dependence of ∆N eff on the local DM potential gφ for a given sterile neutrino parameter choice ∆m 2 41 = 9 eV 2 and sin 2 2θ = 0.25. Notably, the presence of just a tiny potential, e.g., gφ 10 −7 eV, is able to reduce ∆N eff to a negligible level, i.e., ∆N eff 0.01. Increasing the potential gφ will further reduce ∆N eff , thereby making this scenario safe from primordial abundance constraints. One might worry if the presence of the light scalar Φ itself can act as extra radiation, and run into trouble with BBN predictions. This is again prevented in our scenario due to tiny coupling g considered, i.e., according to the relation gφ ≈ 2.15 × 10 −7 eV · (g/10 −22 ) · (10 −18 eV/m φ ).
Sterile neutrinos with masses in the keV range and produced through such a freeze-in mechanism can also act as possible dark matter candidates, as was pointed out by Dodelson and Widrow [107]. However, such a mechanism is in tension with the non-observation of X-rays originating from the decay of sterile neutrino dark matter [109]. Nonetheless, in a scenario such as ours, the sterile neutrino need not be a DM candidate, or can only be a tiny fraction of the DM density of the Universe. As a result, all the X-ray bounds will be rescaled by the fractional density of sterile neutrinos in the DM energy budget, and hence, not be relevant for our scenario [110]. Interestingly, introducing new interactions among the sterile neutrinos, as well as the active neutrinos have been a popular way to relax these X-ray bounds on these models [111,112].

B. Cosmic microwave background and large scale structure
After decoupling around T ν = 1 MeV, neutrinos evolve freely without scattering in our scenario. However, the neutrino masses will keep varying with the background scalar field, and can be possibly large during the recombination epoch. Meanwhile, the CMB and large scale structure (LSS) observations are very sensitive to the absolute scale of neutrino masses. In fact, a world-leading constraint on the sum of neutrino masses has been set, i.e., Σ i m i < 0.12 eV using the dataset Planck TT, TE, EE + lowE + lensing + BAO [113]. Recently, stronger limits, m ν < 0.09 eV at 95% C.L. were derived [114,115] One of the key effects of large neutrino masses during recombination is to reduce the free-streaming length λ fs compared to the massless neutrino case. Free-streaming neutrinos can damp all the perturbation modes for distances less than λ fs , which is well consistent with the current cosmological data. Hence, any effect which reduces the neutrino free-streaming length λ fs during recombination can be constrained from CMB data.
This scenario is different from the BBN epoch, where neutrinos scatter very rapidly before decoupling. In this case, it is important to figure out how a free neutrino mass or flavor eigenstate propagates within the oscillating scalar field. To demonstrate this idea, we plot in Fig. 5 the two mass eigenvalues m 1 and m 4 described by Eqs. (A2) and (A3), within one DM oscillation cycle. Other parameters are taken as m 1 = 0.1 eV, m 4 = 3 eV and θ = π/12. The solid curves stand for the mass m 1 , which coincides with the vacuum value when φ is vanishing. On the other hand, the dotted curves are for the m 4 . As Φ(t) oscillates, we note a non-trivial evolution pattern for gφ = 20 eV (left panel, in blue). At t = π/(2m φ ), the neutrino masses reach their local maximal value. The lighter neutrino mass is just around m 1 ≈ 0.3 eV. After t > π/m φ , Φ(t) becomes negative, and m 1 continuously grows to ∼ 10 eV, which is around the original m 4 value. Note that the mixing sin 2 θ (in gray, unitless) also reaches values around one, implying that the corresponding relations of flavor and mass eigenstates are swapped. The other eigenvalue m 4 follows an opposite behavior.
During recombination, if these two mass eigenstates evolve adiabatically as Φ(t) oscillates, we would expect the neutrino flavors, being a linear combination of the mass eigenstates, to be relativistic half of the time, and nonrelativistic the other half. This is equivalent to setting neutrino velocity to O(0.5 c), which will hence reduce the free-streaming length λ fs by a factor of around two. However, we notice that two eigenvalues also critically hit each other at the point of gΦ(t) = −m 4 , enforcing a non-adiabatic transition. As a consequence, each time a neutrino state, say | ν 1 , reaches this point, there will be a certain probability for it to transit to | ν 4 , and vice versa. The net effect is to reduce the neutrino free-streaming length, by at most a factor of two. We do not investigate numerically the details about the energy-momentum conservation in the scenario where the scalar, active neutrinos and sterile neutrinos form a highly entangled system, but instead give our remarks. As in the case with only sterile neutrino and scalar, the momentum of neutrinos must be always conserved because of spatial translation symmetry. When | ν 1 transits to | ν 4 , the energy of each neutrino state is increased by approximately gφ if p gφ. However, we notice that | ν 4 at t π/m φ has more overlap with the sterile state |N , hence increasing the contribution of the source term in Eq. (2). This feedback effect should balance the energy between the neutrino and scalar systems.
Even though such a scenario, to our knowledge, has never been strictly investigated with cosmological simulations, one might get an intuition from the neutrino mass limit. The neutrino mass information one can extract depends on the redshift. Firstly, the quoted limit from Planck i m i < 0.12 eV is derived from the data of all available redshifts (from z = 0 to z = 3000). Because the scalar field strength is heavier in the dense early Universe, we are more interested in the consequence at higher redshifts. There are fits using the CMB and LSS data by allowing neutrino masses to vary with the redshift, which find the mass limit derived from the data at z > 1100 can only be Σ i m i < 0.40 eV at 95% CL [84]. This can be translated it into a constraint on neutrino velocity, v ν > 0.97 c, by using the temperature T ν (z = 1100) ≈ 0.18 eV and the relation p ν ≈ 3T ν . Even without a dedicated analysis, it indicates that the scenario with v ν > 0.5 c is in tension with CMB observations, if gφ is too large in the early Universe.
However, the above estimate is too stringent, because the impact on CMB perturbation due to a smaller freestreaming is only one of the effects of finite neutrino masses. The inclusion of more effects, such as background evolution effect, should lead to a more conservative limit on the neutrino velocity, which can only be obtained with a more detailed analysis. Furthermore, it is important to note that while the scalar field keeps oscillating with the frequency m φ /(2π), the overall field strength decreases with Eq. (5) as the Universe expands. At certain point when gφ < m 4 the two neutrino mass eigenvalues start to separate. This is demonstrated in the right panel of Fig. 5, where we show the eigenstates for a smaller value of gφ = 2 eV (blue curve).
The neutrino velocity is reduced in the large gφ case due to the conversion between the lower and upper mass eigenstates when gΦ(t) = −m 4 . However, such an issue does not arise if the mass-variation is due to the higher dimensional term, governed by yφ 2 /Λ case with y > 0. This is demonstrated via the red curves in Fig. 5. In this case, no matter how large the potential becomes, the eigenstate with m 1 , which has dominant mixing with active neutrinos, always stays below m 4 sin 2 θ. The threats from cosmological observations can thus be removed by this higher-dimensional operator. In such a case, the local potential yφ 2 /Λ can be large without spoiling BBN, CMB and LSS to have observable time-varying effect for active neutrinos at laboratories. Note that our only requirement here is m 4 sin 2 θ < 0.40 eV, where the sterile neutrino mass can actually be heavy, e.g., m 4 = 1 TeV with θ = 10 −7 . This requirement is well compatible with the current collider searches of heavy right-handed neutrinos. For instance, the ATLAS and CMS collaborations at LHC have set constraints |V µN | 2 , |V eN | 2 10 −3 [116][117][118] at m 4 = 100 GeV, corresponding to a very loose result m 4 sin 2 θ < 0.1 GeV.
In this section, we explored the effects of either gφ or yφ 2 /Λ term in the early Universe, by keeping the other subdominant. In principle, we can simultaneously have these two effects, where yφ 2 /Λ is suppressed by some cutoff scale. However, as we go to the early Universe, the effective term, yφ 2 /Λ, growing faster might dominate over the gφ term, and save the model from cosmological bounds. In the remaining part of the work, we will explore the consequences of a time-varying neutrino mass on beta decays and light sterile neutrino phenomenology. For these analyses, we shall ignore the higher dimensional operator, and focus on the renormalizable interaction gΦ for simplicity.

IV. TRITIUM BETA DECAYS
A time-varying mass of the sterile neutrino can also leave potentially observable imprints in beta-decay experiments, depending on the mass of the sterile neutrino, as well as the mass and amplitude of the scalar field Φ. When the sterile neutrino is heavy (e.g., m 4 > MeV), larger than gφ and decoupled from the energy scale of beta-decay experiments, the scalar potential can only affect the beta-decay spectrum by mixing with active ones [67]. In the standard 3ν scenario, beta-decay experiments, such as Mainz [119], Troisk [120,121] and KATRIN [15,17], measure the effective neutrino mass m β ≡ 3 i=1 |U ei | 2 m 2 i , which receives incoherent contributions from three generations of neutrinos. By analyzing the electron spectrum from tritium beta decays, stringent limits have been set, e.g., m β < 0.8 eV from the combination of the first and second campaign of KATRIN (KNM1 + KNM2) [17].
In order to show the effect of a time-varying scalar on beta-decay spectrum, we consider a simplified modification to the effective neutrino mass, in the limit gφ < m 4 , as where g ν ≡ sin 2 θ 14 · g represents the effective neutrino coupling suppressed by the active-sterile mixing angle. For Eq. (11), a uniform mixing of the sterile neutrino to three generations of neutrinos has been assumed, such that U ei (φ) = U ei and m i (φ) = m i + g ν Φ (for i = 1, 2, 3) hold, thereby leading to the simplified relation in Eq. (11). In general, the diagonalization of the sum of vacuum mass matrix and DM-induced mass matrix will result in complex dependence on the scalar field. For the accuracy of current beta-decay experiments, the major observable is ascribed to a single mass parameter m β . Different coupling patterns are assumed to not affect much the overall magnitude of modifications. The beta spectrum with the effective neutrino mass m β (φ) can be parameterized as with the spectral function Here, G F is the Fermi coupling constant, V ud is the weak mixing matrix element, g V = 1 and g A = 1.247 are the vector and axial-vector weak coupling constants of Tritium and the function F (Z, E e ) is the ordinary Fermi function describing the spectral distortion in the atomic Coulomb potential. Throughout this work, we use E e and K e to distinguish the total and kinematic electron energies and K end,0 is the electron endpoint energy in the massless neutrino limit. The step function, which is necessary to make the square root real, is not explicitly shown. The ultralight scalar may manifest itself as a modulation effect to the beta spectrum, if the DM cycle can be covered by KATRIN runs, and also be resolvable for the duration of KATRIN spectrum scans. If the above condition is not satisfied, one can also look for the distortion effect by averaging over the DM oscillations. In this circumstance, we investigate analytically how the averaged scalar field modifies the beta spectrum. For the current KATRIN sensitivity, it is a good approximation to keep up to the first order of perturbative expansions on m 2 β (unless g ν φ is very large) in the spectral function, namely, The square of effective neutrino mass as in Eq. (11) averaged over one DM cycle reads Hence, in such a case the presence of the scalar field directly adds a constant term to the square of effective neutrino mass. For the sensitivity of first KATRIN campaign, the averaged scalar effect is degenerate with a usual neutrino mass, but it leads to large neutrino mass cosmology [122]. This degeneracy is non-linear and obvious from Fig. 6, where we have fitted the parameter space of m 2 β and gφ using the KATRIN data from the first campaign (KMN1). Following the analysis strategy from the KATRIN Collaboration, we allow m 2 β to become negative during the fit. From Fig. 6 and Eq. (15), one can see that in this scenario, the effective neutrino mass measured, m 2 β is always larger than the true m 2 β and hence, the upper limits derived when assuming gφ = 0 are conservative. Up to this point, we limited the discussion to the case in which the sterile neutrino is heavy. However, when the sterile neutrino is light enough, additional emission channel of beta decays will be open. In the (3 + 1)ν scenario, we can split the 3ν and sterile neutrino contributions as [123] with | U e4 | = sin θ. The φ(t)-dependent mixing angle θ as well as masses m β and m 4 can be calculated from Eqs. (A2), (A3) and (A4), respectively. For simplicity, we assume that the φ(t) dependence in m β is approximately that of m 1 . This is well justified in light of the current KATRIN sensitivity to the effective neutrino mass. For a general light sterile neutrino, one has to integrate the exact spectral function over DM modulations. A generic numerical treatment can be performed in the analysis without making approximations. Aiming to provide a deeper comprehension of the experimental signatures, we show in the upper panel of Fig. 7 the beta-decay spectra for different scenarios, with a configuration similar to the first KATRIN campaign. The lower panel demonstrates the difference of rates of various scenarios with respect to the standard one with m β = 0 eV. The region where KATRIN can obtain most information about neutrino masses is slightly above the beta-decay endpoint. The statistics are very limited close to the endpoint, while far above the endpoint the beta-decay rate is too large and hence, insensitive to very small distortions. A finite neutrino mass will induce not only kink structure, but also provide a constant shift to the beta-decay spectrum away from the endpoint. This is illustrated in Fig. 7 for four different values of m β = {0, 0.4, 0.7, 1} eV, which correspond to the four gray lines, from right to left. However, the kink structure is still unresolvable given the sensitivity of current generation of beta-decay experiments, resulting in the current limit of m 2 β < 0.9 eV from KATRIN's second campaign [17]. A light sterile neutrino will induce a second kink in the decay spectrum. Such spectral feature, which is expected around K e − K end,0 m 4 , can be well separated from a normal neutrino mass term providing the kink is away from the endpoint, as shown by black dashed curve. The size of the distortion is related to the size of the mixing |U e4 | 2 . Consequently, beta-decay experiments set limits to the sterile neutrino mass and mixing [124,125]. Red curves in Fig. 7 illustrate how adding large scalar potentials to the sterile neutrino will smooth these distortions, mainly because the averaged mixing will be suppressed by the potential. Consequently, this scenario allows to open up the sterile neutrino parameter space in the context of short-baseline anomalies, as we will discuss later in Fig. 9. The case in which the sterile neutrino is heavy is given by the the blue curve. We have shown that the effect of the scalar is degenerate with the mass term (see Fig. 6).  [129]. The purple band is the original rate expectation without any deficits, and the gray one is the best-fit scenario with a constant neutrino flux deficit. A sterile neutrino with the time-varying potential generates the blue curves. To be more specific, the parameters are chosen as m 4 = 3 eV, θ = π/5, gφ = 4.2 eV and m φ = 6.6 × 10 −22 eV.

V. IMPACT ON GALLIUM AND REACTOR ANOMALIES
Until now, we have discussed the impact of a consistent scenario of time-varying neutrino masses on cosmology and beta decay experiments. If the sterile neutrino is light today, it will also have important consequences for the long-standing gallium and reactor anomalies. We address these issues in this section. First, we will focus on the case where the DM cycle matches the time scale of the experiment, giving rise to a signal of time modulation. We use the latest BEST experiment as a case study. Second, we move to higher scalar masses in the rapidly oscillating regime and discuss its impact on the parameter space of light sterile neutrinos.

A. Time modulation and the gallium anomaly results
Solar neutrino detectors GALLEX [126] and SAGE [127] have used gallium to measure the neutrino emission rate from radioactive 51 Cr and 37 Ar sources, and found the rates to be lower than the prediction. This is known as the gallium anomaly, which can be, in principle, explained by adding an eV-mass light sterile neutrino to the SM. This anomaly has been recently strengthened by the BEST collaboration, which detects ν e from a 51 Cr source kept in a two-volume detector, filled with gallium [128,129]. BEST provides high statistics real-time data of the event rates in the inner as well as outer detectors. The purpose of this subsection is to illustrate that the type of signal predicted by the time-varying scenario could manifest as anomalous experimental results like the one reported by BEST, and the time modulation probed would correspond to the ultralight scalar mass range of interest. To do so, we performed a simple fit of the available data considering m 4 , θ 14 and gφ as free parameters.
Our results are shown in Fig. 8, which depicts the neutrino reaction rate at inner and outer targets, starting from 14:02 on 05 July 2019. The event rates with error bars, shown in red, are taken directly from Ref. [129]. The purple band is the original prediction without any electron neutrino disappearance. The gray band represents the best-fit scenario with a constant neutrino flux deficit. The model with time-varying sterile neutrino mass generates the blue curve, which clearly shows a periodic behaviour in addition to a reduction in the expected ν e flux. To generate this curve, we use the best-fit values, m 4 = 3 eV, θ = π/5, gφ = 4.2 eV and a DM mass m φ = 6.6 × 10 −22 eV (corresponding to a cycle of 73 days), which is right in the fuzzy DM regime. We find that the original fit with constant flux deficit gives a global χ 2 ≈ 32, while the time-varying scenario can reduce it to χ 2 ≈ 18. Then, it might seem that the time-varying massive sterile neutrino scenario provides a better explanation of the BEST results than the vanilla sterile neutrino case. Note that the procedure followed did not account for possible time-dependent systematics. Besides that, in order to actually claim the existence of a time modulation in data, a more sophisticated analysis is required in order to address if the periodic behaviour in data is spurious (for instance, due to statistical fluctuations) and its statistical significance. Finally, in order to actually conclude whether this scenario provides a better fit to the data, one would need to perform a proper Monte-Carlo statistical analysis which allows to interpret consistently the meaning of the lower χ 2 obtained.
Apart from that, the interpretation of the anomaly in terms of a time-varying sterile neutrino is in conflict with solar neutrino data and other active neutrino oscillation experiments. The active to sterile oscillation will lead to  a deficit in the solar neutrino flux as well [130,131]. The original BEST fit is already in tension with the solar neutrino data, which cannot be saved by simply adding a modulation, unless there is a large overdensity of DM inside the Sun (e.g., ten times of that on the Earth). Furthermore, by mixing, the active neutrino masses will unavoidably receive time-varying corrections. However, no modulation pattern or averaged effect have been found in existing oscillation experiments [132,133] so far. In particular, modulations with an amplitude larger than 10% are excluded by solar data for modulation periods ranging from O(10 minutes) to O(10 years) [134][135][136]. It remains to be seen whether the modulation pattern persists or not in future data collections. Should the modulation pattern be caused by other uncontrolled time-dependent systematics rather than time-varying sterile neutrinos, we may leave it for a future analysis. Nevertheless, this serves as an illustrative example of how an experiment like BEST could constrain the time-varying neutrino mass.

B. Reinterpretation of light sterile neutrino parameter space
The gallium anomaly is further strengthened by the presence of the reactor antineutrino anomaly -a series of reactor experiments which observed a deficit of the electron antineutrino flux as compared to some predictions (see [87,137] for recent reviews). In this subsection, we will discuss how the light sterile neutrino parameter space changes in connection with the gallium as well as the reactor anomalies under the oscillating scalar potential hypothesis. Note that a typical reactor neutrino experiment baseline is small, e.g., baselines are of the order 0.01-1 km, and therefore is sensitive to DM modulations driven by a mass (or equivalently frequency) m φ ≈ 10 −7 − 10 −9 eV. For simplicity, we will focus on the scenario where the DM oscillates very rapidly during the neutrino flight from the source to the detector, such that all the time-varying information is effectively lost. Since the DM potential encountered by the neutrinos changes rapidly, neutrinos in flight do not have enough time to develop any nontrivial phases due to this interaction. The Hamiltonian relevant for active-sterile neutrino oscillations is reduced to a vacuum term and a time-dependent one, namely H = H 0 + H I (t) with In the event of rapid DM oscillations, any terms in the Hamiltonian proportional to single powers of gΦ will be averaged to zero. On the other hand, the quadratic term will be averaged to (gφ) 2 /2. This constant term behaves similar to a matter potential, causing a shift in the dispersion relations; see Appendix B for further details. Fig. 9 shows the impact of a rapidly oscillating DM, coupled to sterile neutrinos, on the different types of constraints existing in the sterile neutrino parameter space for four different benchmark potential values, gφ = 0, gφ = 3 eV, gφ = 10 eV and gφ = 50 eV (from top left to bottom right). The red shaded regions indicate our KATRIN limits. The constraint from the ∆N eff requirement [91] is shown as the gray region in the top-left panel, while there are no such constraints for the other three cases. The black dashed lines show the region favored by gallium and reactor data [138], while the brown shaded region indicate the solar neutrino constraint [139]. For convenience, the active neutrino masses have been fixed to zero in all the cases. We find that as gφ increases, the favored parameter space from the global analysis of gallium and reactor data shrink to the top-right corner. This is due to the fact that for gφ m 4 , the active-sterile mixing will be severely suppressed, and one needs to go to larger vacuum mixings to have the same effective parameters. Similar trends are seen for the solar bounds as well. Regarding KATRIN, for increasing values of gφ, larger values of the mixing become allowed for masses below 10 eV. The bumpy structure present in the vanilla sterile neutrino constraint (top-left panel), which arises as a consequence of statistical fluctuations, is kept approximately in all cases. The analysis has been limited to the case of m β = 0 eV. Including this parameter in the fit should smooth the displayed curves. As a related issue, one can see that for non-zero gφ, maximal mixing is not allowed for sterile neutrino masses of O(50 − 100eV). That region of parameter space corresponds to the limit in which sterile neutrinos are integrated out and only the time modulation affecting active neutrinos can manifest in beta-decay experiments. In that case, there is a degeneracy between gφ and m 2 β (see discussion in Section IV), so it is likely that in a more general analysis by including m 2 β as a free parameter, the constraint in that region will be relaxed. As discussed before, for neutrino oscillation experiments, only (gφ) 2 term exists when we average over the rapidly oscillating scalar field. As a result, the mixing angle suppression goes as θ ∼ θ · m 2 4 /(gφ) 2 . On the other hand, for beta-decay experiments, the suppression to the mixing angle is much milder θ ∼ θ · m 4 /(gφ) because the oscillating scalar field is not averaged at the Hamiltonian level as in the neutrino oscillation experiment.
There are additional constraints in the parameters space of light sterile neutrinos coming from oscillation experiments (see Ref. [140] for a recent review and references therein for a detailed description of existing limits). We leave the study of the rich phenomenology expected in that case to a future work.

VI. CONCLUSION
The ultralight DM, as an alternative DM candidate, can develop large interactions with the invisible neutrinos. Such interactions have triggered studies at various neutrino experiments like oscillations, beta decays and 0νββ decays. However, active neutrinos and charged leptons form doublets in the SM, and hence caution should be taken when we couple neutrinos to the exotic fields. One of the safest ways is to introduce a SM singlet, e.g., sterile neutrino, which can share the exotic forces with light neutrinos by mixing.
In that spirit, we have systematically explored the time-varying neutrino masses by assuming that the sterile neutrino couples to an ultralight scalar field. The original time-varying neutrino masses, if testable in neutrino oscillation experiments, are actually in severe tension with the neutrino mass constraints from Big Bang nucleosynthesis, cosmic microwave background, and large scale structures. We demonstrate that the coupling with light sterile neutrinos can provide time-varying signals in active neutrinos while being in accordance with cosmological observations. However, to evade the constraints from CMB and LSS, it might be necessary to introduce higherdimensional operators. If the sterile neutrino is light, there are additional concerns about producing extra radiation prior to BBN, as measured by ∆N eff . We study this scenario, and show that the BBN constraint on light sterile neutrinos can be safely avoided in this case. Additional thermal production of the light DM is not of concern due to tiny couplings in our scenarios. This allows us to provide a cosmology-friendly scenario of time-varying masses of active neutrinos.
We show that the coupling of the ultralight scalar to neutrinos can manifest at beta-decay experiments like KATRIN. We perform a dedicated analysis of the impact of such time varying neutrino masses in KATRIN. We find that in the presence of additional light sterile neutrinos, the limits set by KATRIN change considerably.
Furthermore, we also study the impact of such models on light sterile neutrino searches, with emphasis on the results from reactor and gallium experiments, as well as solar data. We find that as the value of the DM mass increases, the region favoured by gallium and reactor data shrink to larger values of sterile neutrino mass and mixing angles. This is due to the large suppression of active-sterile mixing angle, leading to the necessity of large vacuum mixing angles and masses to satisfy the data. Additionally, in the latest gallium data from the BEST collaboration, we find that the time-varying neutrino mass can give rise to signatures in the direction of the observed anomalous results, with the hypothesis, giving a lower χ 2 value, ∆χ 2 = −14, with respect to the original fit. However, such a scenario turns out to be in conflict with other neutrino oscillation experiments, pointing to other possible origin of the modulation pattern.
Such a framework of time-varying neutrino masses can not only be tested against cosmological observations, but also provide observable signatures in beta decay experiments like KATRIN. Furthermore, it also holds potent consequences for short baseline anomalies that plague neutrino physics these days. m 4 cos 2θ 14 , we have the approximation m 1 m 1 + sin 2 θ 14 · gΦ , (A5) m 4 m 4 + cos 2 θ 14 · gΦ , (A6) which usually applies to the case with heavy sterile neutrinos. In the early Universe, we often have |gΦ| m 4 , which will lead to another useful approximation depending on the sign of gΦ. In the case of gΦ m 4 > 0 eV, we have the approximation Whereas, in the case of gΦ −m 4 < 0 eV, we instead have When gΦ switches sign, the expressions for m 1 and m 4 are also swapped, but note that the active neutrino always has dominant overlap with the lighter eigenstate. The above results are analogous to the analytical observations of standard MSW effect. In the absence of light sterile neutrino, active neutrinos will develop large effective masses in the early Universe because φ scales as (1 + z) 3/2 with z being the redshift. This significantly constrains the effect which can be probed in laboratory. The introduction of a light sterile neutrino is very intriguing in avoiding the CMB limit. From Eq. (A11), one can observe that the contribution of large gΦ to the effective active neutrino mass is negligible, contrary to the case in Eq. (A5). This can be understood as the fact that the mixing between active and sterile neutrinos is severely suppressed when gΦ becomes extremely large, whereas gΦ only acts on the sterile component. Such a scenario can also be a key to suppress the production of sterile neutrino in the early Universe [141].
Appendix B: Neutrino oscillations with rapidly oscillating scalar potential The scenario with slowly oscillating scalar field (compared to the neutrino time of flight) is trivial. Here we explore the scenario with rapid potential alternation during the neutrino time of flight. Let us focus on the scenario with two neutrino flavors as before, i.e., ν a and N c . Suppose an active left-handed neutrino component ν a,L is emitted from the source. We need to figure out how the time-varying neutrino mass alters the prediction of standard neutrino oscillation picture. The analysis present below should also be applicable to time-varying masses of active neutrinos only. We start by explaining that neutrino-antineutrino oscillation, e.g., ν a,L → ν a,R , is always negligible in our scenario. Since we are interested in light sterile neutrinos here, all mass terms are tiny in comparison to the neutrino beam energy E ν 1 MeV. We further assume that the scalar potential is also negligible compared to E ν , such that neutrinos are always ultra-relativistic during the flight. This guarantees that a left-helicity state is always dominantly composed of the left-handed field. In the mean time, during the propagation the neutrino helicity is a conserved quantum number, because the addition of the time-varying potential Φ(t) does not affect the spatial translation symmetry (conserving three momentum) and isotropy (conserving angular momentum). To see it more explicitly, we write down the quantum mechanical Hamiltonian driving the two-component spinors in the basis of (ν a,L , ν c a,L , N c , N ): The helicity operator is given byp · Σ withp being the momentum direction vector and Σ = σ ⊗ 1 4×4 . One can easily check that [H,p · Σ] = 0 by noticing [p · σ, σ · ∂] = 0. As a consequence, the neutrino produced from the source, which is dominantly left-helicity, will project mainly into the left-handed field operator at the detector. The right-handed component is always negligible providing that all the mass terms and scalar potential are small compared to the neutrino beam energy. We have also verified the behavior numerically. Hence, in the following we shall confine the discussion to the flavor conversion between fields with the same handness (say, left handness). In the flavor basis, the approximated Hamiltonian for a plane wave reads .
The Hamiltonian relevant for active-sterile neutrino oscillations is reduced to a vacuum term and a time-dependent one, namely H = H 0 + H I with H I (t) = 1 2p 0 m D gΦ m D gΦ 2m N gΦ + (gΦ) 2 . (B4) The evolution of neutrino state |ν t = C a |ν a + C N |N c is governed by the equation Here, |C a | 2 gives the survival probability of ν a . When H I is tiny compared to H 0 , one may consider to perturbatively analyze the effect of gφ by treating H I as a small parameter. However, we note that neutrino oscillation is an accumulative effect over a macroscopic baseline L, and small terms accumulated over L are not always appropriate to be treated as perturbations. In a general circumstance, exact calculation should be invoked to derive the oscillation probability. For a rapidly oscillating scalar field with m φ ∆m 2 /(2E), neutrinos feel an averaged potential during propagation. This does not mean that the scalar effect is vanishing. In fact, the average of H I (t) generates a non-vanishing contribution H I (t) = 1 2p which affects neutrino oscillations by modifying the effective mass-squared difference and mixing angle. The ultimate effective Hamiltonian will be H 0 + H I (t) . The above results should also be valid for decoherent scalar field with random phases, just as the normal matter effect works for an assemble of random microscopic particles.