Tritium beta decay with modified neutrino dispersion relations: KATRIN in the dark sea

We explore beta decays in a dark background field, which could be formed by dark matter, dark energy or a fifth force potential. In such scenarios, the neutrino's dispersion relation will be modified by its collective interaction with the dark field, which can have interesting consequences in experiments using tritium beta decays to determine the absolute neutrino mass. Among the most general interaction forms, the (pseudo)scalar and (axial-)vector ones are found to have interesting effects on the spectrum of beta decays. In particular, the vector and axial-vector potentials can induce distinct signatures by shifting the overall electron energy scale, possibly beyond the usually defined endpoint. The scalar and pseudoscalar potentials are able to mimic a neutrino mass beyond the cosmological bounds. We have placed stringent constraints on the dark potentials based on the available experimental data of KATRIN. The sensitivity of future KATRIN runs is also discussed.

There might be various scenarios in which neutrinos feel a background of particles in the "dark sea", which can be formed by fifth force, ultralight DM, as well as dark energy (DE) field. By collectively interacting with the background fields, the neutrino kinematics is expected to be altered similar to the established MSW effect. In this context, neutrino oscillations in the dark potentials have been widely discussed in the literature [39,[66][67][68][69][70][71][72][73][74]. However, the neutrino oscillation experiments are only sensitive to relative values of the potential. If three generations of neutrinos couple identically to the dark field, similar to the neutral Z-exchange in the Standard Model, the neutrino oscillation will be blind to the dark potential. Furthermore, if ultralight dark matter is responsible for the dark potential, and the dark field is fast oscillating over the neutrino baseline, neutrinos in flight will experience a vanishing averaged effect. In this regard, the experiments probing absolute neutrino masses such as β decays can be a better place to search for such dark potentials. In this work, we perform a systematic analysis of the impact of dark potentials on β-decay neutrino mass searches in which kinematics is used to determine m ν . The emphasis will be put on the ongoing KATRIN experiment [75][76][77], which has recently set the world-leading model-independent constraint on the absolute neutrino mass. However, the framework presented in this work should also in principle apply to other promising neutrino mass experiments, such as Project 8 [78][79][80][81], ECHo [82] as well as PTOLEMY [83]. It is worthwhile to mention that KATRIN data have already been used to set constraints on the Lorentz violation parameter by searching for possible periodic signals [84]. Despite various realizations of the dark potential, our results of analysis will be presented in a form as model-independent as possible, such that they can be applied to various scenarios of model realizations.
The rest of this work is organized as follows. In Sec. 2, starting from general interaction forms we derive the modified equation of motion and dispersion relation of neutrinos. In Sec. 3, we give some examples of specific model realizations for completion. In Sec. 4, we compute the beta-decay rate in the dark potential. In Sec. 5, we explore the experimental effect at KATRIN. We make our conclusion in Sec. 6.

Dark MSW effect
Regardless of the nature of the underlying interaction, neutrinos will feel a "dark potential" that can have five different forms [85,86]: Here, δ M = 1/2 (1) if neutrinos are Majorana (Dirac) particles, and φ, ϕ, V µ , a µ and T µν represent real scalar, pseudoscalar, vector, axial-vector and tensor fields, respectively. The neutrino field is composed of ν = ν L + ν c L for Majorana neutrinos, and ν = ν L + ν R for Dirac neutrinos. For Dirac neutrinos, the coupling constants g φ , g V and g a are Hermitian matrices in general, while g ϕ is anti-Hermitian. For the Majorana case, g φ and g a are real symmetric matrices, while g ϕ and g V are purely imaginary symmetric and antisymmetric matrices, respectively. The fields φ, ϕ, V µ , a µ and T µν do not necessarily correspond to real particles, e.g., they may represent contributions from coherent forward scattering or virtual fifth forces. Note that here we assume the field X = φ, ϕ, V µ , a µ , T µν to be real, and only the overall sign combined with "charge" g X X matters for neutrinos.
A well-known example is already given by the standard MSW effect. For instance, via the charged-current interaction with electrons in the Earth, the neutrino will feel a tiny potential of the form g V V 0 + g a a 0 = √ 2G F n e , where G F is the Fermi constant, n e is the electron number density, and the isotropic spatial components g a a are averaged out. Our results are simplified by imposing the assumption that the dark field is purely time-like, i.e., we assume from the cosmological principle that the preference of spatial orientation of the background is not significant or simply averaged out. Hence, the antisymmetric tensor field T µν [49] is vanishing in our context.
The equation of motion (EOM) of the neutrino wave function given the interactions in Eq. (1) is described by where the neutrino is written in the mass eigenstate basis in vacuum, therefore the mass matrix M ν is diagonal. We can assume that three generations of neutrinos share the same coupling, i.e., g X ∝ 1 for all coupling constants, which in particular implies that there is no effect in neutrino oscillations. Note that for Majorana neutrinos, the diagonal vector interactions will be vanishing. The collective effect of the dark sea is to modify the dispersion relation of neutrinos, which can be obtained by multiplying γ µ (i∂ µ − g V V µ + g a a µ γ 5 ) + ( M ν + g φ φ − g ϕ ϕγ 5 ) to the left of Eq. (2). For the plane-wave solution we have where E ν is yet to be fixed. Ignoring the cross terms of the background fields (i.e., we turn on only one field at a time), we end up with the following dispersion relation for neutrinos: where E ν and p ν are the energy and momentum of the neutrino,p denotes the direction of the momentum, and Σ ≡ γ 5 γ 0 γ stands for the spin operator. Eq. (3) has typically two energy solutions: the upper one (positive for most cases) corresponds to the particle ν (+) , and the lower one (mostly negative), which is not bounded from below, should be interpreted as the antiparticle ν (−) . For the antineutrino, the direction of the momentump should be reversed along with the energy accordingly. Eq. (3) implies that the scalar and pseudoscalar potentials g φ φ and |g ϕ |ϕ simply add to the vacuum mass term. The vector potential g V V 0 shifts the overall energy of neutrinos. The axialvector potential g a a 0 will lead to helicity-dependent energies of the states withp · Σ ν (+) = ±ν (+) for neutrino and −p · Σ ν (−) = ±ν (−) for antineutrino, where '±' stands for the right-and lefthelicity states, respectively. This split of energies is discussed in Appendix A. In the massless limit, i.e., M ν = 0, the difference between the vector and axial vector interactions vanishes for the active neutrino ν L , because the left-and right-handed fields are decoupled. More technical details on the derivation of the dispersion relation and how the neutrino should be canonically quantized in the dark background are presented in Appendices A and B.

Model realization
Even though our results derived from KATRIN will be presented in a model-independent form, we want to discuss a specific model realization here for completion. In general, there are several ways to produce dark potentials for neutrinos in the literature, including the following.
• A background of ultralight dark matter, dark radiation or dark energy coupled to neutrinos.
The ultralight field can be treated as a classical one, and the neutrino dispersion relation is affected simply by assigning an expectation value to the dark field in the Lagrangian, e.g., g φ φ νν for an ultralight scalar. Modified neutrino oscillations in such scenarios have been discussed in detail in the literature [45][46][47][48][49][50][51][52][53][54][55][56][58][59][60][61][62][63]. For the scalar ultralight dark matter, the field evolves in the Universe as φ =φ cos m φ t with m φ being the mass of the dark matter field. The field strength readsφ = √ 2ρ/m φ , where ρ is the energy density of dark matter. For the vector dark matter A µ , the spatial component A µ will be averaged out if the polarization of dark matter is randomized, while the temporal component A 0 is simply vanishing for free vector particles nearly at rest. • Coherent forward scattering of neutrinos with massive dark matter particles, the so-called "dark NSI" [64,87]. In this case, it is the elastic scattering of the neutrino wave function off the dark matter grid which changes the dispersion relation. The form of the dark potential depends on the type of interaction between dark matter χ and neutrinos. However, due to the smallness of the dark matter density, a large coupling and a small mediator mass would be required to achieve an observable potential. • A fifth force sourced by heavy dark matter [39,87] or by ordinary matter [65,66,[66][67][68][69][70][71][72][73][74]88].
This can be regarded as a special case of coherent forward scattering, as the mediator mass is extremely small, such that the interaction can be described by a long-range force. The fifth force as a classical virtual field is able to directly modify the neutrino propagation, similar to the motion of electrons in a Coulomb potential. Note that to avoid severe constraints from the charged lepton sector, such force is usually assumed to be generated by mixing with sterile neutrinos.
For illustration, we will elaborate on the last scenario, which can easily realize all forms of dark potentials described in Eq. (1). In order to generate different forms of the dark potential in Eq. (1), we consider a fermionic dark matter χ to source a long-range scalar or vector potential. The Lagrangian is then given as After integrating out dark matter configurations, we can obtain the fifth force potential sourced by the ambient dark matter. For the vector potential, one ends up with an effective potential A µ = (g χ n χ /m 2 A , 0, 0, 0) [69,71] with n χ being the dark matter number density and m A being the mediator mass, and the spatial components are vanishing because χ is non-relativistic. For the scalar potential, one simply has Φ = g χ n χ /m 2 Φ . With the scalar or vector fifth force, one can obtain the dark potential forms in Eq.
The Debye screening effect will reduce the magnitude of dark potentials in astrophysical environments with dense mobile νν pairs [39,89,90]. In analogy to the electromagnetic screening in a metal with free electrons, the free neutrino pairs will shield the vector field. This is equivalent to giving the vector field a Debye screening mass m DS , which causes drastic exponential decrease in the field strength, e.g., V 0 ∝ e −m DS r /r with r being the distance from the source. Taking the vector potential for example, the resultant potential in the solar system with the screening effect reads where n χ denotes the local DM density, and the screening mass induced by neutrino plasma is roughly m DS ∼ g ν T ν . The long-range force parameters can be set to be very small, e.g., m A = 10 −21 eV and g ν = 10 −17 , such that no laboratory bounds other than beta decays can apply to the parameter space of interest. One might be concerned about astrophysical limits from the likes of BBN, CMB, LSS and supernovae, where the dark matter or neutrino density is much higher than in the Earth. However, various arguments attempting to constrain this type of force are mostly invalid due to the screening effect in the dense environment; for relevant details see Refs. [39,89,90]. The Lagrangian in Eq. (4) will also induce the DM self-interaction. Before we discuss beta decays, let us briefly comment on the consequences of this self-interaction. The self-interacting DM is in fact favored by observations of small scale structures [26,27], which itself is a wellmotivated topic. The major constraint should come from the observations of structure formation. To avoid making DM too collisional, one puts g χ 4 × 10 −3 (m χ /GeV) 3/4 [91]. There is also a collective bound on DM long-range interactions from tidal stream of the Sagittarius satellite [92][93][94], which has excluded g χ /m χ > 10 −19 GeV −1 with m Φ 4 × 10 −28 eV corresponding to the Sagittarius dwarf galaxy orbit 16 kpc. However, we have checked that most of the parameter space that can give O(eV) dark potentials remains unconstrained by those considerations.

The beta-decay rate in the dark sea
The microscopic nature of the nuclear beta decay related to electron neutrinos makes it an excellent complementary probe of the dark potential other than neutrino oscillations. The amplitude for the transition of beta decays, e.g. 3 H → 3 He + e − + ν e , is given by whereḡ V andḡ A stand for the vector and axial-vector coupling constants of the charged-current weak interaction of nucleons, respectively, the higher order magnetic and pseudoscalar form factors of the nucleons are neglected, and p, p , p e and p ν are the momenta of tritium, helium, electron and neutrino, respectively. Summing over the spins of particles other than neutrino, we arrive at where for the unsummed spin bilinear of neutrinos v(p ν , s ν )v(p ν , s ν ) under the impact of dark potentials, Eqs. (24) and (33) in the Appendix should be taken. For the vector dark background, we are ready to sum over the final neutrino spin, i.e. v(p, s)v(p, s) = / p − m. But for the axial-vector case, due to the split of energy levels discussed in Appendix A, the integration over phase space for two helicity states should be performed separately. This introduces extra complexity. After the index contraction, the matrix element for the outgoing neutrino with p ν,s in the axial-vector background is With a vanishing a 0 , the result will be reduced to the standard one, which is consistent with Ref. [95]. The matrix element in the vector background has a similar expression, but one can sum over the neutrino helicity, leading to which is close to the standard results but with p ν in vacuum replaced by p ν = ( E ν , p ν ).
The final beta-decay rate without the sum of neutrino helicity reads Note that the neutrino energy in the phase space factor is different from the vacuum case, namely E ν for the vector background and E = E s for the axial-vector one, such that the normalization and completeness relations of spinors can appreciate the simple forms of Eqs. (23), (24), (32) and (33) in the Appendix. In principle, one can take different normalization conventions, but the final result is invariant. The neutrino energy in the delta function should take the form in Eq. (22) or (31). The integration should be done in the rest frame of the tritium, in accordance with the frame picked out by the dark sea considering the Earth is non-relativistic. After the trivial integration over d 3 p and decomposing d 3 p e = |p e | 2 dp e d cos θ eν dφ eν , we have Since d cos θ eν = E /(|p e | · |p ν |)dE and the neutrino favors no specific direction, the decay rate is simplified to We are left with integrating over the neutrino momentum in order to obtain the differential spectrum with respect to the electron energy. The integration limit of |p ν | can be obtained by requiring that has a solution for any −1 cos θ eν 1 and m e E e E max e . Trivial analytical solutions exist for the standard [96] and vector cases. For the axial-vector case, we integrate the rate numerically. It is worthwhile to remark that at the maximal electron energy (i.e., minimal neutrino energy) in the axial-vector case, |p| = 2sg a a 0 , and the phase space factor in Eq. (12) is not vanishing as in the standard case. This gives rise to a finite decay rate at the endpoint of electron spectrum.

Signals at KATRIN
In the absence of dark field, the rate of beta decays, 3 H → 3 He + e − + ν e , reads [95][96][97][98][99] where N T is the total mass of the tritium sample, and σ(E e ) is the reduced cross section (see, e.g., Ref. [100] for the expression). The kinematics of the beta-decay spectrum is contained in where K end,0 = (m 3 H − m e ) 2 − m 2 3 He /(2m 3 H ) is the endpoint energy in the relativistic theory assuming a vanishing neutrino mass. The actual endpoint energy for the neutrino mass m i is approximately given by As mentioned above, the scalar and pseudoscalar potentials merely add an effective mass term to neutrinos, which is kinetically indistinguishable from the vacuum mass. In fact, in some scenarios they are even postulated to be the origin of small neutrino masses [39]. However, we need to emphasize that since the current dark potential is expected to be different from that in the early Universe, the model-dependent cosmological bounds on the absolute neutrino mass, e.g. Σ < 0.12 eV [101], can be evaded or weakened. This will possibly lead to large signals in future KATRIN runs, which expect no visible effect of neutrino masses if the stringent cosmological bounds are adopted [100].
The vector potential has a profound effect on the beta-decay spectrum. To clearly see that, we pick out from Eq. (3) the vector contribution to the energy, which reads E ν = p 2 ν + m 2 i + g V V 0 for the neutrino mode and E ν = p 2 ν + m 2 i − g V V 0 for the antineutrino mode. The fact that neutrino and antineutrino excitations feel opposite vector potentials implies that neutrino mass experiments using electron capture such as ECHo [82], will see an opposite effect compared to KATRIN, providing a way to independently test the effect. The antineutrino energy of beta decays in the vector background can run into the negative (but bounded), when g V V 0 > p 2 ν + m 2 i . The beta-decay spectrum can hence extend beyond the normal kinematic limit K end,0 . This is not a surprise as the process which is not kinematically allowed in vacuum can take place if the medium modifies the dispersion relations [9], a phenomenon familiar in, e.g., plasmon decay. In the presence of the vector potential, the electron endpoint energy in Eq. (15) will be shifted towards Since the axial-vector interaction distinguishes two helicity states, to exactly calculate the spectrum we have to perform the integration over the unsummed helicity amplitudes with Eq. (12). We assume T 2 as the only tritium source, and only the final-state excitations of 3 HeT + need to be considered. The differential spectrum will have to sum over the final-state distributions of the daughter molecule. Given the accuracy of KATRIN, we use the Gaussian-averaged final-state distributions in Ref. [102]. The ultimate integrated event rate at KATRIN is given by convolution of the differential spectrum dΓ β /dK e with the spectrometer response function. The predicted rate is given by where A s is the normalization factor, N T is the target tritium number, qU is the applied retarding potential, and R bkg is the background rate. Here, we have assumed a constant background rate within the energy window of interest, following the KATRIN collaboration [77]. If the vector potential is too large, the reconstruction of the background rate can be affected. However, we

Count rate (cps)
KATRIN spectrum without background STD g ϕ ϕ, g φ φ = 0.5 eV g V V 0 = -5 eV g V V 0 = -2 eV g a a 0 = -5 eV g a a 0 = -2 eV g ϕ ϕ = g V V 0 = 0.5 eV  Figure 1: The integrated beta-decay spectra for the first KATRIN campaign in various types of dark seas: vector potentials with g V V 0 = ±5 eV (solid red curves) and g V V 0 = ±2 eV (dashed red curves) as well as axial-vector potentials with g a a 0 = ±5 eV (solid blue curves) and g a a 0 = ±2 eV (dashed blue curves). For these curves, the neutrino mass has been fixed to m 1 = 0.1 eV. Contributions from φ and ϕ which mimic a neutrino mass exceeding cosmological bounds are shown as the dashed gray curve. The left (right) panel gives the spectra without (with) the background. In the right panel, the KATRIN data points are also given for comparison. Note that for demonstration, the errors are shown as fifty times the standard deviation.
note that the background and signal follow different spectrum shapes and can be statistically distinguished by the fitting procedure. The response function f res (K e − qU ) is dependent on the surplus energy E e − qU , where qU is the applied electric potential. For the first campaign of KATRIN, we use the response function given in Fig. 2 of Ref. [75] with a column density 1.11 × 10 17 molecules · cm −2 . The predicted rate is to be compared with the measured one, for which the ring-averaged event rate with statistical and systematic errors is available in Ref. [77]. The variables A s , R bkg and m 2 ν will be taken as free parameters during the fit. In Fig. 1, we have illustrated the distortions of beta-decay spectrum in various dark seas at KATRIN. The left panel stands for the ideal case without taking into account the background, while the right one gives a more realistic result with the background. The vector potential (red curves) shifts the whole spectrum to lower or higher endpoints, without changing the spectral shape. In comparison, the axial-vector potential (blue curves) induces a non-trivial distortion to the spectrum near the endpoint, but it indeed converges to the vector case away from the endpoint. By measuring this distortion, one can distinguish between the effects of vector and |g ϕ ϕ|, |g φ φ| g V V 0 , g a a 0 < 0 g V V 0 , g a a 0 > 0 With NP (Future)  Figure 2: The constraints on the dark potentials g V V 0 , g a a 0 , g φ φ and g ϕ ϕ using the data of KATRIN's first campaign, shown as the solid curves. We have normalized the likelihood L, i.e., dividing by the likelihood maximum L max in each case. The sensitivities corresponding to the ultimate KATRIN goal m ν < 0.2 eV at 90% level, i.e., σ(m 2 ν ) = 0.025 eV 2 [105], along with a reduced potential fluctuation by a factor of three are given as the dotted lighter curves (in the presence of a potential) and dashed curves (without new physics). The horizontal gray lines represent ∆χ 2 = 1, 4 and 9, respectively. axial-vector potentials. However, we find numerically that current KATRIN runs are not yet sensitive to this distortion, which would require more statistics. We expect that the Project 8 experiment and PTOLEMY proposal can provide better sensitivities to probe such a distortion near the endpoint, but how good the sensitivity is would require a further study. On the other hand, the scalar potentials φ and ϕ mimic the effect of neutrino masses, shown as dashed gray curves. They can induce an effective neutrino mass beyond the cosmological constraint. Even though we turn on only one dark potential at a time, in the left panel of Fig. 1 we show the case of g φ φ = g V V 0 = 0.5 eV as the dotted green curve to imply the existence of degeneracy when multiple dark potentials are turned on. One can notice that the effects of two dark potentials can counterbalance each other to some extent.
It is worthwhile to mention that the previous anomalous signal detected by Troitsk [103] may be explained by the vector potential, which can shift the endpoint even beyond its maximum value and affect the reconstruction of m 2 ν . However, the anomalous signal has disappeared in later measurements [104]. To avoid false signals, experimental systematics must be well controlled in order to probe such new physics effects.
In the limit of massless neutrinos, the difference between vector and axial-vector dark potentials should be vanishing, as expected from the analysis in Sec. 2. However, as in Fig. 1 the vector and axial-vector cases have an apparent difference and cannot smoothly converge by taking the neutrino mass to be vanishing. Only when the neutrino mass is comparable to time-scale of the background field formation (typically ∼ 1 Gyr corresponding to 10 −32 eV), the axial-vector scenario starts to approach the vector case. For further discussions, see Appendix C.
We continue with fitting KATRIN data to our scenarios. The effect of scalar and pseudoscalar potentials is identical to being from m 2 ν . The consequences of vector and axial-vector potentials are the same at energies away from the endpoint, for which KATRIN collects the most events. For the present KATRIN sensitivity, the major effects of the vector and axial-vector potentials are to shift the endpoint energy E 0 of electrons. Hence their effects are entirely ascribed to the fit of E 0 ≡ K end + g V V 0 at KATRIN. The endpoint E 0 and the squared neutrino mass m 2 ν are regarded as free parameters in KATRIN fits [75][76][77]. We perform our own fit with available data of the first KATRIN neutrino mass campaign (KNM1) [75]. For this purpose, we vary freely the vacuum neutrino mass over m 2 ν ≥ 0. In the official fit of KATRIN, the m 2 ν < 0 region is kept to account for data fluctuations, but eventually it is removed by proper statistical interpretations in obtaining the mass limit. We adopt here a simplified approach with the main interest being the minimized χ 2 , whose function for KNM1 is constructed as where the statistical and systematic uncertainties, σ i sta and σ i sys , are taken from Ref. [77] for each retarding potential qU i . For KNM1, there are 27 set points of qU i in total, which have already been shown in the right panel of Fig. 1. We hence assume the vacuum neutrino mass to be vanishing and vary other parameters E 0 , A s and R bkg freely. Marginalizing over parameters other than E 0 1 , our fit result on the endpoint will be compared with the expected value, and this is then used to set a limit on g V V 0 . Note that we do not attempt to fit the results of the second KATRIN neutrino mass campaign (KNM2) here [77], for which the statistical framework should account for every detector ring and there is not sufficient information to perform the fit ourselves.
The fit yields E 0 = 18573.79 ± 0.02 eV. The actual Q-value is obtained by correcting for molecular recoil (1.72 eV) and potential fluctuations of the tritium source and main spectrometer (−0.2 ± 0.5 eV for KNM1), which gives Q = 18575.31 ± 0.5 eV, slightly smaller than the expectation of 18575.72 ± 0.07 eV. For the (pseudo)scalar potential, we assume the vacuum mass to be vanishing.
The likelihood L = exp (−∆χ 2 /2) for our dark potentials is given in Fig. 2. The 2σ limits corresponding to ∆χ 2 < 4 read −1.4 eV < (g V V 0 , g a a 0 ) < 0.6 eV and (|g φ φ|, |g ϕ ϕ|) < 1.1 eV. For potential values smaller than O(1) eV, the χ 2 is dominated by events away from the endpoint, where the effects of vector and axial-vector become almost the same.
A slight preference of the (axial-)vector potential g V V 0 = −0.42 eV can be noted. Keeping the best-fit values so far, future sensitivities to reach the reference KATRIN target m ν < 0.2 eV and by reducing potential fluctuations by a factor of three are shown as dotted curves. In comparison, assuming no new physics contributions, the sensitivities are instead given by the dashed curves.

Conclusion
We have performed a novel and systematic study on the effects of dark neutrino potentials on beta decays, especially focusing on the ongoing KATRIN experiment. By collectively interacting with background fields, neutrinos will have dispersion relations different from the ones in vacuum, which induces distinct distortions to the beta-decay spectrum. Observable consequences include neutrino mass signals beyond the ones bounded from cosmological constraints, events beyond the kinematical endpoint of the decay, and spectral distortions. We find that the current KATRIN data favor the (axial-)vector potential g V V 0 = −0.42 eV, but more statistics should be required to draw a more robust conclusion. The current KATRIN runs are not yet sensitive to the discrimination between vector and axial-vector cases, but future experiments may be able to achieve this by more accurately measuring the distortions near the endpoint.
The ECHo experiment with the electron capture technique will feel an opposite vector potential compared to KATRIN, which can provide a complementary probe if the dark potential is present. The next-generation beta-decay experiment like Project 8 (with molecular or atomic tritium) measuring the differential spectrum can be an excellent further probe of dark potentials, providing a well controlled energy scale of the spectrometer. With a tritium source as large as 100 g, we expect the PTOLEMY proposal to have an overwhelmingly better sensitivity, which is interesting for a future work.

A Massive neutrinos in an axial-vector background
The plane-wave solutions to the Dirac equations in the axial-vector background differ from the vacuum ones. To see that, we collect the left-and right-handed fields into ν = ν L + ν R (ν R = ν c L for Majorana neutrinos), which satisfies the equation of motion For the positive and negative frequency modes, we have First of all, the energy eigenvalues should be derived. This can be done by multiplying a matrix (p µ γ µ + g a a 0 γ 0 γ 5 + m) from the left to Eq. (20), yielding where Σ = γ 5 γ 0 γ is simply the spin operator, andp · Σ/2 represents the helicity with eigenvalue s = ±1/2. Different from the vacuum case, the neutrino energy is split for the two helicities by the temporal component of the background. The energy eigenvalues for both neutrino and antineutrino (the same for Majorana case) read Note that the above equation does not apply to the massless case with m = 0 eV. By imposing the orthogonality conditions 2 , the structure of the spinors u and v is found to be of the form where E 0 = (g a a 0 ) 2 + m 2 is the neutrino energy at rest. Notice that for the spinor v, we have the relationp · Σ v(p, s) = −2s v(p, s), in comparison top · Σ u(p, s) = 2s u(p, s). The helicity completeness relations of the spinors are It is easy to verify these relations with the help of the orthogonality conditions, and one recovers the standard results in the limit of a 0 = 0. Expanding the field operator as ν = d 3 p(b ν p,s ), we arrive at where b p,s and d † p,s should be interpreted as the particle annihilation and antiparticle creation operators, respectively, for Dirac neutrinos. For Majorana neutrinos, the condition ν = ν c will force b p,s = d p,s , i.e., b p,s annihilates simultaneously the positive-and negative-frequency excitations. Using the orthogonality conditions, the canonical quantization rules of ν(x) consistently lead to The neutrino Hamiltonian with normal ordering can then be expanded as (e.g., for the Dirac case) 2 One can check that these orthogonality conditions guarantee the plane wave to be normalized under the field expansion Eq. (25), e.g., d 3 x ν † p ,s ν p,s = δ 3 (p − p)δ s s .

B Dirac neutrinos in a vector background
The results for the vector background are more straightforward. Given the EOM the plane-wave spinor should satisfy This leads to the energy eigenvalues where "±" in E ± corresponds to the neutrino (u) and the antineutrino (v), respectively. Note again that these results do not apply to Majorana neutrinos. Hence, in principle one can distinguish Majorana and Dirac neutrinos by the experimental signature of the vector potential, if one would know that the interaction is diagonal in flavor.
The orthogonality conditions as well as the completeness relations are consistently given by where the effect of the dark background is absorbed into Here, it is equivalent to replacep · Σ with γ 5 / S, where / S ≡ (|p|/m, Ep/m). Note that we have not yet summed over the spin s in the above expressions, and the standard results can be easily obtained by summing over s. These relations should be used along with the expansion or equivalently the form Eqs. (32), (33) and (35) indicate that one may think of all the relations withp = { E, p} similar as those in the vacuum. The net impact of the dark background is adding an overall phase exp(−ig V V 0 t + ig V V · x) to the neutrino field. Since other fields (e.g., n, p and e) do not feel this phase, it will enter into the factor δ 4 (· · · ± g V V ) which imposes energy momentum conservation. Ultimately, the Hamiltonian of neutrino field is found to be The evolution of neutrino energy as the dark potential g V V 0 (red curve) or g a a 0 (blue curves) adiabatically increases. The neutrino momentum is taken as p = 0.2 eV. Right panel: The dispersion relation with respect to p with g V V 0 (red curve) or g a a 0 (blue curves) being 0.2 eV. In both panels, for the solid curves the neutrino mass has been set to m = 0.05 eV, while for the dotted blue one m = 0 eV.

C Some remarks on the ground state in the dark sea
The formation of the background field typically takes place on cosmological time scales, say 1 Gyr corresponding to ∼ 1/(2 × 10 −32 eV), which is significantly larger than the Compton frequency of neutrinos, i.e., the inverse of mass 1/m ν ≈ 6.6 × 10 −15 s for m ν = 0.1 eV. The neutrino modes will therefore always stay in their energy eigenstates during the adiabatic formation of the background field, meaning that the eigenvalues of neutrino energy should change in a continuous and smooth manner without transitions.
For the vector and axial-vector cases, let us investigate in more detail how the energy of a neutrino mode evolves when the background field g V V 0 or g a a 0 gradually changes from zero to a certain value. Their dispersion relations for the right-helicity antineutrino (corresponding to ν L in the massless limit) are recast as follows: where p ≥ 0 represents the magnitude of neutrino momentum. As long as m = 0, these are indeed smooth functions of the potential field. To be specific, we take the neutrino mass as m = 0.05 eV and set the neutrino momentum to be p = 0.2 eV. Then let g V V 0 and g a a 0 adiabatically change from 0 eV to 1 eV. The evolution of energy is shown in upper panel Fig. 3. In the lower panel, we fix g V V 0 and g a a 0 as 0.2 eV and vary p. For comparison, in both panels of Fig. 3 we give the case of axial-vector potential with vanishing neutrino mass as dotted curves. Special care should be taken when the neutrino mass is vanishing, i.e., m = 0 in Eq. (38). By taking the derivative of Eq. (38), we have ∂E ∂(g a a 0 ) = − ∂E ∂p = g a a 0 − p (p − g a a 0 ) 2 + m 2 .
It is clear that as long as m = 0, the energy E is a smooth function of a 0 and p. However, when m = 0, the derivative becomes ill-defined at p = g a a 0 . A smooth solution to the massless case in g a a 0 l, r l, r r l l r m ν Figure 4: An illustration of the evolution of the neutrino energy eigenstates, as the magnitude of dark potential g V V 0 < 0 (left panel) or g a a 0 < 0 (right panel) is adiabatically increasing. The baseline of E ν = 0 is set by the brown line. The state with an empty circle corresponds to the neutrino, while that with the filled circle to the antineutrino. The energy is split for the left-helicity ('l') and right-helicity ('r') states in the axial-vector case. Neutrinos are not allowed to stay in the shaded region on the right.
the axial-vector potential should be which becomes identical to the vector scenario in Eq. (37). This is exactly what we expect when the neutrino mass is vanishing, for which the difference of results between vector and axial-vector scenarios is supposed to vanish. The adiabatic evolution of the eigenstates of neutrinos is schematically shown in Fig. 4 for the vector and axial-vector cases, respectively. We explain the figure quantitatively in what follows and note that we have verified it by numerically solving the Dirac equation. For the axial vector case, there is an energy barrier set by the neutrino mass which keeps the neutrino state above the zero-point energy, as g a a 0 adiabatically increases. In the massless limit, such a barrier does not exist, and the left-handed field shifts smoothly down as in the vector case. Because only the left-handed neutrino field is responsible for the beta decays (left-helicity for neutrino and right-helicity for antineutrino in Fig. 4), the effects of axial-vector and vector dark potentials on the beta-decay spectrum should be the same for m ν = 0.
For the axial-vector potential, the adiabatic approximation will break down in the extremely narrow parameter space 0 < m ν 10 −32 eV (for which the background field changes faster than the neutrino mass), and one may expect the probability of tunneling crossing the mass barrier (for the massless case without barrier, the tunneling probability is equivalently one). This is very similar to the matter effect of neutrino oscillations in varying matter profile [7]. On the resonance, when the adiabaticity parameter is large γ O(1), the neutrino will always stay in one specific mass eigenstate. But for γ O(1), transition occurs from one neutrino mass eigenstate to another. However, we should notice that at least two neutrino mass eigenvalues should be larger than 0.0086 eV according to neutrino oscillation data, and hence the transition for them is always adiabatic. The above discussion is only relevant for cosmological time scales.
For the time scale relevant to KATRIN runs, the dark potential explored here does not change.
In Fig. 1 of the main text, in the limit of m ν → 0, the result of the axial-vector case seems unable to continuously transit to that of the vector one. This is not a surprise considering that the cosmological time scale, one billion years corresponding to 2 × 10 −32 eV, separates the massless and sizable massive cases.