Discovering leptonic forces using non-conserved currents

Differences in lepton number (i.e., $ L_e - L_\mu $, $ L_e - L_\tau $, $ L_\mu - L_\tau $, or combinations thereof) are not conserved charges in the Standard Model due to the observation of neutrino oscillations. We compute the divergence of the corresponding currents in the case of Majorana or Dirac-type neutrinos and show that, in the high energy limit, the vector interactions map onto those of a light scalar coupled to neutrinos with its coupling fixed by the observed neutrino masses and mixing. This leads to amplitudes with external light vectors that scale inversely with the vector mass. By studying these processes, we set new constraints on $ L_i - L_j $ through a combination of semi-leptonic meson decays, invisible neutrino decays, neutrinoless double beta decays, and observations of Big Bang Nucleosynthesis/supernova, which can be much stronger than previous limits for vector masses below an eV. These bounds have important implications on the experimental prospects of detecting $ L_i - L_j $ long-range forces.


I. INTRODUCTION
Light vector bosons (X µ ) are simple, technically natural, extensions to the Standard Model of particle physics, which can play a key role in explaining experimental anomalies [1][2][3][4][5][6][7], act as a mediator to the dark sector [8][9][10], or make up dark matter themselves [11][12][13][14][15][16]. If a vector mass (m X ) is well below the weak scale then, depending on the type of current involved, experimental signals have drastically different behavior as a function of m X . For a vector coupled to the electromagnetic current (i.e., a kinetically mixed dark photon), then as m X → 0 it becomes an unobservable correction to electromagnetism resulting in all constraints deteriorating at low masses. On the other hand, if it couples to any other conserved current, j µ , (i.e., a current which satisfies ∂ µ j µ = 0 at the quantum level) then the constraints become massindependent as m X → 0. Lastly, if a vector is coupled to a non-conserved current, there exist processes whose amplitudes are proportional to 1/m X , leading to constraints that get stronger at low masses (see [17][18][19][20][21][22][23][24] for cases where this property was previously used to set bounds on light vectors).
The Standard Model with massless neutrinos has a set of five linearly-independent conserved currents: where j EM is the electromagnetic current, B and L are baryon and lepton number respectively, and L i denotes lepton number of generation, i. Once neutrino masses are introduced the only currents which remain conserved are electromagnetism and B − L (assuming neutrinos are of Dirac-type); The "flavored" leptonic currents, j Li−Lj , cannot explain the observed neutrino masses and mixing without introducing a small breaking, typically assumed to be from integrating out some scalar(s) and * jdror@lbl.gov right-handed neutrinos. [25][26][27][28] 1 While the inevitable breaking has negligible impact when m X > ∼ eV, it can have dramatic effects at lower masses, in the form of enhanced production of the vector's longitudinal mode. Our goal is to explore these effects.
Since the symmetry is broken by the neutrino masses it will lead to amplitudes proportional to (g X m ν /m X ), and hence most prominent at low m X . As we will show, a combination of constraints from semi-leptonic meson decays (K ± → ± νX, π ± → ± νX), lepton-flavor violating processes (ν i → ν j X, X-emitting neutrinoless double beta decays), and neutrino annihilations (changing ∆N eff /supernova), result in important bounds on the L i − L j parameter space. To keep the discussion general we do not specify a particular mass mechanism for the neutrinos. While this has the drawback that processes may develop sensitivity to the UV, most do not, and we focus on bounds that are robust.
Aside from the constraints we present here, the best constraints on gauged L i −L j depend strongly on whether the difference in lepton number includes electron number. For vectors coupled to j Le−Lµ and j Le−Lτ with m X < ∼ O(0.1 eV) the most stringent current limits are from fifth force/EP-violation searches [29][30][31][32] and modifications to neutrino oscillations due to a long-range potential induced by the Earth, Sun, or whole galaxies [29,[33][34][35][36][37][38], with the gauge coupling (g X ) being constrained to be weaker then gravity for m X < ∼ O(10 −3 eV). For L µ − L τ , where common matter does not induce a long-ranged force, the dominant low mass constraints are from dynamics of neutron star inspirals [39] as well 1 These mechanisms are generically easier to build when m X /g X > ∼ 100 GeV since then the vacuum expectation values of any U(1) L i −L j -breaking scalar can then easily contribute evenly to all the neutrinos, thereby preventing any of the lepton mixing angles from being too small. However, we do not take this as a bound as one can always evade this theoretical challenge by introducing multiple scalar fields.
as from pulsar binaries [39,40]. For m X > ∼ 10 −11 eV (the inverse size of a neutron star) these bounds do not apply and the most stringent bounds are from (non-enhanced) neutrino annihilations to vectors increasing ∆N eff . If the universe reheated above the muon mass after inflation and there were no additional entropy dumps after muons froze-out, then the best constraints are from muon annihilations to vectors µ + µ → γX constraining g X < ∼ 4 × 10 −9 [41]. Otherwise, the strongest limits are due to neutrino annihilations constraining g X < ∼ 5 × 10 −6 [42]. We will consider the same process but focusing on the enhanced contributions in section VI. For light bosons there are additional constraints from the observation of near-extremal black holes and the null observation of superradiance [43,44]. While we include these, they could change dramatically if the scalar field used to provide to break U(1) Li−Lj is light enough to induce selfinteractions within the superradiance cloud [44] The constraints presented here, dramatically shape the low mass L i − L j parameter space. In particular, for L e − L µ and L e − L τ the enhanced processes are the strongest bounds for m X < ∼ 10 −22 − 10 −17 eV, depending the theory assumptions. For L µ − L τ the new constraints can be the strongest for all m X < ∼ 100 eV or for eV < ∼ m X < ∼ 10 −10 eV and m X < ∼ 10 −15 eV, depending on the theory assumptions.
The paper is organized as follows. In section II we compute the divergence of L i − L j current and study the correspondence between L i − L j vectors and light scalars. We use these to compute constraints from semileptonic meson decays in section III, neutrino decays in section IV, neutrinoless double beta decays in section V, and neutrino annihilations prior to Big Bang Nucleosynthesis (BBN) and inside supernova in section VI. We conclude by discussing the implications of these new bounds with a focus on the prospects of seeing deviations to neutrino oscillations from long-range forces in section VII.

II. Li − Lj AND NON-CONSERVED CURRENTS
We begin by summarizing gauged L i − L j models, with an emphasis on the parts relevant for non-conserved currents. The coupling of X µ to the Standard Model fermions is given in terms of 2-component Weyl spinors by, 2 where we denote right-handed fields using a c and write all quantities in the flavor basis withˆ. The U(1) Li−Lj charge matrix,Q ij depends on the particular flavor symmetry and is = diag(1, −1, 0), diag(1, 0, −1), or diag(0, 1, −1) for L e −L µ , L e −L τ , or L µ −L τ respectively. 2 We follow the conventions of 2-component spinors laid out in [45].
In addition, there may be couplings of right-handed neutrinos to X, X µQc iiν c † iσ µν c i , with an a priori arbitrary charge matrix,Q c i . We discuss this in some detail below. Imposing one of the lepton flavor symmetries, in combination with electromagnetism forbids any mixing in the charged lepton sector (i.e.,ˆ i = i ) and hence the observed lepton mixing which makes up the Pontecorvo-MakiNakagawaSakata (PMNS) matrix, U , arises entirely from the neutrino mixing matrix (i.e.,ν i = U ij ν j ). As such, performing a rotation to the mass basis, changes the form of the coupling to neutrinos but leaves the charged lepton interaction diagonal. The neutrino mass term depends on whether neutrinos are Majorana (M) or Dirac (D), 3 In either case imposing L i − L j cannot reproduce the PMNS matrix. This is simple to see if neutrinos are Majorana since L i − L j would restrictm ij from having any mixing between one of the flavors from the other two and hence cannot possibly reproduce nonzero values for all three 3 mixing angles. If neutrinos are Dirac then the situation is more subtle. Upon diagonalization the mass term is given by, where V is the right-hand neutrino rotation matrix necessary to diagonalizem ij such that V †m U = m ≡ diag(m 1 , m 2 , m 3 ). V is undetermined until a choice is made for the charges of the right handed neutrinos,Q c i . However, even allowing for arbitrary right-handed neutrino charges, the neutrino mass term violates L i − L j symmetry. To see this explicitly consider the effect of a L i − L j rotation on the neutrino mass masses. In order for the term to be L i − L j invariant we requirê For an arbitrary matrix, VQ c T V † , this equation has no solution since the form ofQ requires the left hand side to have a vanishing column while mU has been measured to be entirely non-vanishing. We conclude that the neutrino masses and their mixing always break U(1) Li−Lj .
The consequence of this breaking can be seen by computing the divergence of the current and we carry this out in detail in appendix A. The result is, or in terms of the mass eigenstates, where we have introduced some convenient notation to eliminate extra U s when possible: Q ≡ U TQ U (Majorana) or Q ≡ U †Q U (Dirac) and m 1,2,3 are the neutrino masses in ascending order. 4 We will use Q orQ depending on which is more convenient.
In the Majorana case, the divergence is manifestly nonzero if Q ij and m i are nonzero, while for the Dirac case one can show that there is no choice of Q c ij that allows the Dirac term to vanish, which is simply a restatement of the above observation that the presence of neutrino masses inevitably breaks the L i − L j symmetry. For all processes studied in this work, a right-handed neutrino coupling to X would only add to the rates at leading order in the neutrino mass (which will be an expansion parameter). To be conservative we simply set Q c = 0 throughout but note that if the right-handed neutrinos are charged under U(1) Li−Lj , it would strengthen our prospective bounds.
The impact of this term becomes clear employing Goldstone boson equivalence theorem (GBET), which states that, in the high energy limit (energy of X much larger than its mass), we can compute amplitudes of X µ employing the replacement, X µ → g X i∂ µ ϕ X /m X , and computing the amplitudes for a scalar ϕ X instead (see appendix A for more details). The ϕ X interactions, in the mass basis, are then governed by, To obtain this expression we have integrated by parts to replace ∂ µ j µ Li−Lj using (9), setting Q c → 0. For boosted external X µ , this results in amplitudes proportional to g X m ν /m X . These will always be parametrically the most important processes as m X → 0. We see that, at high energies L i − L j vectors are effectively Majorons [47] for Majorana neutrinos and neutrino-philic scalars for Dirac neutrinos. Of course one is always free to either use the Goldstone boson expression, (10) or the vector Lagrangian, (3) (though this correspondence breaks down when X µ is virtual). Any UV complete model will contribute new sources to the right-hand side of 9 through the field(s) that spontaneously break the symmetry and adding these could change our results. In particular since we can be in the regime where m X /g X ∼ eV, it is likely other degrees of freedom can be produced on-shell in processes closely related to the ones studied here. While this would add to the rates, these are more model-dependent and we do not include them. Furthermore, using these additional degrees of freedom to eliminate the processes computed here would likely introduce an enormous fine-tuning and would be challenging given the multitude of possible signals.
Lastly, we comment on the minimal requirement for theories with a non-conserved current -ensuring unitarity; If g X /m X is too large, then in any diagram with neutrinos, we can add an additional external X and increase the amplitude, leading to a breakdown of unitarity. Since a g X /m X -enhancement will always cost a power of a neutrino mass, the unitarity condition is, The factor of 4π arises from the phase space cost of adding in an extra particle. The unitarity bound applies equally to both Majorana and Dirac neutrinos and is robust against the details of the UV completion (up to O(1) corrections). We will impose this throughout as an orange region in our constraint plots.
As our first enhanced process we study semi-leptonic meson decays, P ± → ± b ν a X (P = B, D, K, π). This process is depicted in the Goldstone picture in Fig. 1  (top-left), where the scalar is radiated off the neutrino leg. In the vector picture there are two diagrams that one must add together, and a delicate cancellation takes place for the longitudinal mode leaving only the piece proportional to m ν g X /m X . This process was previously considered in [48] to probe L i − L j however they did not include the important enhanced contribution for small m X . Energy-enhanced semi-leptonic decays to νX were discussed for parity violating muon couplings [18] and for vectors coupled only to neutrinos [19,49,50]. Such vectors are coupled to currents whose divergence is proportional to the lepton masses and hence have amplitudes that grow as ∝ m /m X .
The decay rate is different in the case of Majorana and Dirac neutrinos and we compute both. The final result is, where phase space integral is given in term of the energy fractions x 1 ≡ 2E X /M P , x 2 ≡ 2E νa /M P , and we defined α ≡ M /M P , the ratio of the lepton to meson mass, and neglected the X µ and external neutrino masses. The spin-averaged square of the matrix element is, V CKM is the CKM matrix element on the hadronic side ( = 1 for pion decay, the sine of the Cabibo angle for Kaon decay, etc.), and f P + is the meson factor (f π 130.2 ± 0.8 MeV [51], f K 155.7 ± 0.3 MeV [51], etc.), g 0.6 is the weak coupling, m W is the W boson mass. The dimensionless function, I(x 1 , x 2 ), is given by, where the momenta are defined in Fig. 1 (top-left). The phase space integral has an IR divergence when the internal neutrino goes on-shell, a feature emphasized and treated carefully in [52]. This arises from loop corrections of X µ to the neutrino mass. The sensitivity to this divergence is mild and can be regulated by introducing a small mass for the internal neutrino and we employ this strategy here.
There are many experiments that looked for semileptonic charged meson decays. Due to the number of events, the best constraints arise from pion and kaon decays and we focus on those. Decays to electrons have a significant advantage as an additional vector can lift the usual 2-body helicity suppression, reducing the background rates. For rare pion decays a dedicated search for π ± → e ± ν(X → inv) was performed in [53] constraining the branching ratio < ∼ 4 × 10 −6 , which we approximate as our limit. In addition, we note that the recent PIENU experiment has about ∼ 10 3 times as many pions and could likely significantly improve this limit as was recently done for heavy neutrino searches in pion decays [54,55]. No analogous search exists for π ± → µ ± νX though, given its additional experimental challenges, it is not likely to significantly improve the limits.
There are no dedicated searches for K ± → e ± νX or K ± → µ ± νX however E949 looked for K ± → µ ± ννν [56] constraining the branching ratio to < ∼ 2.4 × 10 −6 (improving on a limit set in [57]). Since the neutrino and X channels share a similar final state, we use this limit to approximate our bound. We note that the constraints may be able to be improved by searching for the electron channel due to reduced Standard Model backgrounds.
The limits are depicted in the case of Majorana neutrinos for L µ − L τ and L e − L µ in Figs. 2 and 3. Given their qualitative similarity we reserve the rest of the cases to appendix B. We see that both these constraints are a powerful bound on the parameter space. Furthermore, they are insensitive to the UV completion, do not make any assumptions about cosmology, and apply to both Majorana and Dirac neutrinos. This robustness against the modifications of the model makes this channel important in searching for light L i − L j vectors.
Gauging L i − L j induces (neutrino mass-suppressed) lepton-number violating processes. The simplest possibilities are invisible neutrino decays, as depicted in Fig. 1  (bottom-left), where a heavier neutrino decays into a lighter neutrino, emitting off a L i −L j vector. In principle its possible to look for visible neutrino decays instead of invisible decays through the neutrino magnetic moment, . While the bounds on such decays are much more stringent, with searches looking for spectral distortions in the Cosmic Microwave Background (CMB) constraining lifetimes reaching O(10 20 sec) [58], the insertion of a magnetic moment suppresses the rate relative to the invisible decay by a parametric factor of ∼ (m ν /m W ) 4 which is too large a suppression to be relevant for us here and we focus on invisible decays.
The spin-averaged matrix element squared has a mild dependence on whether neutrinos are Majorana (top line) or Dirac (bottom line): where ReQ ab and ImQ ab denote the real and imaginary parts of Q ab . In the limit the PMNS matrix is real, the rates differ by a factor of 2. The strongest constraint on invisible neutrino decays come from cosmology. If neutrinos are present during recombination then their decay, in combination with subsequent coalescence with X, can prevent free-streaming from efficiently wiping out structure at small enough-in Cosmic Microwave Background [59,60]. A recent study using Planck 2018 data [61] These are the strongest limits on L i − L j through enhanced processes. In addition, the cosmological bound could be drastically improved with the observation of a nonzero neutrino mass sum [62][63][64]. A null result would improve the robustness of the cosmological bound and resultantly significantly change the allowed parameter space making cosmological neutrino lifetime measurements the most promising avenue to discover ultra light L i − L j vectors. While significant, the present constraints assume neutrinos are around during recombination, which need not be the case if additional decay or annihilation channels were active after Big Bang Nucleosynthesis producing other purely free-streaming sources of radiation [65]. Alternatively, one could consider cases where the neutrino masses vary in the late universe [66], rendering them stable during recombination. While perverse, these options reflect a clear model-dependency in the constraints and hence we also include constraints on the neutrino lifetime from terrestrial experiments. The strongest such bound is from studying the flavor composition of solar neutrinos [67]. Observations of the Sudbury Neutrino Observatory in combination of other solar neutrino experiments [68], constrain the lifetime of the second neutrino mass eigenstate, While there is a huge gap between the ability of cosmological and terrestrial bounds, future studies of the flavor content of astrophysical neutrinos at IceCube could improve this sensitivity by roughly four orders of magnitude if astrophysical uncertainties can be sufficiently well understood [69][70][71] (with even some hints at decays in current data [72]). We show the constraining power of these processes in red in Fig. 2 for L µ − L τ and in Fig. 3 for L e − L µ with the cosmological bound as a dashed line and the, more robust, terrestrial bound as a solid line.

V. 0ν2βX DECAYS
If neutrinos are Majorana-type, then L i − L j vectors can induce in neutrinoless double β decays. For L e − L µ and L e − L τ the dominant production mode is single X emission as shown in Fig. 1 (top-right). For L µ − L τ this process is forbidden. This is easiest to see when working in the flavor basis, where the W and X interactions are both flavor diagonal and all the mixing is in the mass terms. Since then the X interactions are diagonal, and it does not couple to electrons for L µ − L τ , an additional neutrino mass insertion is needed to induce neutrinoless double beta process making it highly suppressed. As such, the dominant mode for L µ − L τ is double-X emission.
In either case, X-induced neutrinoless double beta decays are only allowed for Majorana neutrinos, where the high energy limit matches precisely to that of a Majoron with appropriate identification of the couplings, allowing us to make use of dedicated Majoron studies. Searches  [39], black hole superradiance [44], and ∆N eff measured through Big Bang Nucleosynthesis, with the latter depending on whether the universe reheated above the muon mass [41] (dashed) or below [42] (solid). Enhanced constraints come from unitarity (orange) meson decays (blue) with K → νX being most stringent, ∆N eff measured with Big Bang Nucleosynthesis due to enhanced neutrino annihilations (green), and neutrino decays (red). The neutrino decay bounds arise from both terrestrial (solid) and cosmological (dashed) searches, the latter which assumes the neutrinos to be present during recombination.
for neutrinoless double beta decays in association with a Majoron have been carried out by KamLAND-Zen [73] and NEMO [74,75] with the strongest limit from the former, setting the bound on the Majoron coupling to electron neutrinos, |g ee | < ∼ (0.8 − 1.6) × 10 −5 , depending on the nuclear matrix element. The corresponding coupling in terms of L i − L j parameters is, (recall that the hats denote the matrices within the flavor basis). The prospective sensitivity to neutrinoless double beta decays with multiple X emission was considered for Majorons [76] and is comparable in overall rate sensitivity, however, since the rate is proportional to (g X m ν /m X ) 4 , we do not expect these to significantly improve the bounds given other constraints and as such we do not include them here.
The sensitivity of 0ν2βX searches is shown in pink in Fig. 3. As discussed above, these do not apply for Dirac neutrinos, nor for L µ − L τ .

VI. NEUTRINO ANNIHILATIONS
X bosons can also be produced through neutrino annihilations in a thermal bath in the early universe or during a supernova, as depicted in Fig. 1 (bottom-right). If both Xs are transverse then the amplitude is ∝ g 2 X , if one X is transverse and the other is longitudinal, then the amplitude is ∝ g 2 X m ν /m X while if both Xs are longitudinal then the diagram is ∝ g 2 X m 2 ν /m 2 X . Since we focus on the case where g X 1, the dominant contribution will be from double-longitudinal emission and we neglect the rest.
Unlike the rest of the processes discussed so far, which can be computed in either the Goldstone boson equivalent Lagrangian or using the vector Lagrangian, double-X emission processes must be computed using the Goldstone boson picture to get a reasonable estimate of the rate (without providing a UV completion for the source of neutrino masses). The reason for this can be traced back to additional contribution from the UV completion to the divergence of the current in (9). In a U(1) Li−Ljinvariant model there must be a condensate coupling to X µ (e.g., a L i − L j -breaking fundamental scalar). Since the condensate breaks the symmetry, it must also have linear couplings to X µ X µ and will always contribute to processes with double-X emission. These are crucial since otherwise one may mistakenly conclude that double-longitudinal emission amplitudes for νν → XX scattering can scale be ∝ 1/m X , as opposed to the expected ∝ 1/m 2 X . This confusion does not arise if one works in the Goldstone boson equivalent theory since then the condensate coupling is not needed to ensure this cancellation.
Note that a U (1) Li−Lj -breaking condensate will not have a significant impact on single-X emission considered above since single-X emission using additional scalars In gray are previous constraints from fifth forces [29], deviations to neutrino oscillation data [38], and black hole superradiance [44]. The enhanced constraints come from unitarity (orange), meson decays (blue), dominated by K → νX, neutrinoless double beta decay searches (pink), Big Bang Nucleosynthesis/supernova (green), and neutrino decays (red). Neutrino decay bounds depend on whether or not neutrinos are present during recombination (dashed vs solid) must arise from operators with derivatives on the scalar field(s) which can only interfere with the contributions we computed at loop level. Aside from this somewhat technical point this highlights an important point: bounds from neutrino annihilations have additional O(1) uncertainties due to the presence of additional diagrams that contribution at the same order. In this sense multiple-X emission bounds are less robust, prior to choosing a UV completion, then single-X emission bounds. Keeping this in mind, we can still roughly estimate the bounds using places where neutrino scattering to invisible states can be actively measured, namely, Big Bang Nucleosynthesis and supernova.
Since we are interested in the high energy limit, we can again make use of the literature on Majorons and neutrino-philic scalars. Bounds from Big Bang Nucleosynthesis (BBN) on Majorons were computed in [42]. Assuming flavor diagonal couplings they find the Majoron coupling to must be < ∼ 2 × 10 −5 at low masses to avoid disturbing the measured light element abundances. The constraints for L i − L j vectors will be similar with appropriate re-weighting of the coupling. To estimate the bounds we set this equal to ∼ g X m 3 / √ 3m X , where the factor of 1/ √ 3 attempts to account for the fact that only one neutrino is interacting to first approximation (the most massive one). We require the same bound for Dirac neutrinos. The bounds are then shown in Fig. 2 and Fig. 3 for L µ − L τ and L e − L µ respectively. Supernova cooling bounds using supernova 1987A were computed for Majorons in the massless case [77], constraining Majoron couplings to the level of O(10 −6 ).
These could slightly strengthen the bounds from BBN, however are sensitive to the assumption on the inner temperature and could potentially be irrelevant if there was no proto-neutron star present after the supernova at all [78]. Since these do not improve the bounds significantly, we do not include these but note that a supernova in our galaxy or improvements in supernova simulations would improve upon the annihilation bounds.
There are a few other scattering processes one may consider. Charged lepton annihilations can proceed through → XX and → Xγ however these will not be enhanced at tree level. Alternatively neutrinos scattering off leptons or quarks through virtual electroweak boson exchange and radiating off a single X can result in enhanced processes. These will be efficient in a thermal bath at high enough temperatures however any limits set to this way will be inherently sensitive to the physics prior to BBN and will not be a robust bound on the parameter space. Alternatively one could try to use neutrinonucleus or neutrino-electron scattering experiments however these experiments do not have the sensitivity necessary to probe the parameter space considered here.

VII. DISCUSSION
In this work we presented new constraints on L i − L j using the non-conserved nature of the currents. These bounds get stronger for smaller vector masses, making them particularly relevant for searches for long-range forces associated with gauging differences in lepton num-ber. In particular, with current data we find that for L e − L µ or L e − L τ , the new limits are now the strongest for m X < ∼ 10 −18 eV with a traditional cosmology or m X < ∼ 10 −24 eV with a modified cosmology to evade CMB limits on the neutrino lifetime. For L µ − L τ , where searches for long range forces are less constraining, we find the new limits are the most important bounds on the parameter space for m X < ∼ 10 −2 eV with a traditional cosmology or for a modified cosmology: 10 −10 eV < ∼ m X < ∼ 10 −2 eV in conjunction with m X < ∼ 10 −18 eV. In either case our results have profound implications on the prospects of detecting L i − L j vectors using neutrino oscillation searches. For gauged L e − L µ or L e − L τ , the bounds are no longer the strongest to arbitrarily small masses and observing such long-range forces on galactic scales [38] is not possible. Interestingly, if a measurement of the neutrino lifetime improves to cosmological scales (e.g., if a non-vanishing sum of neutrino mass is detected), the new bound could make it impossible to see deviations to neutrino oscillations even from the Earth/Sun. This would require a measurement of τ νa→ν b X > ∼ 10 8 year, within the reach possible with the upcoming Euclid Satellite in combination with Planck data [64].
In addition, there have been proposals to probe gauged L µ − L τ with a mass-mixing of the vector with the Z boson (i.e., introducing a term ε Z m 2 Z X µ Z µ ), and again looking for deviations to neutrino oscillations patterns constraining the product, ε Z g X < ∼ O(10 −51 ) for m X < ∼ 10 −15 eV [79,80]. Combining previously studied bounds on the mass-mixing parameter, ε Z < ∼ 10 −30 (m X /10 −15 eV) [23] with CMB bounds on neutrino lifetime from this work, g X < ∼ 10 −24 (m X /10 −15 eV), we find searching for mass-mixed L µ − L τ is no longer viable. Relaxing the CMB bounds allows a small window for masses above the neutron star binary bounds.
We conclude by commenting on other prospective enhanced signals, which do not turn out to be effective at looking for L i − L j vectors. Firstly, one can look for lepton-number violating processes such as charged current decays at one-loop through neutrino mass insertions. For Dirac neutrinos, the chiral structure of the electroweak interactions forbid any contributions to the amplitude at order O(m a ) leading to suppression by an additional power of m a /m W , rendering the decays unobservable. This suppression turns out to be present also for Majorana neutrinos. This is easiest to see in the Goldstone picture for X. In this case the only diagram that can contribute to a → b X decays are from radiating off ϕ X from an internal neutrino leg. The chirality structure of the neutrinos and the weak interaction then requires an additional neutrino mass or ϕ X insertion for the process to take place. This suppression is significant enough to make this process uninteresting for our purposes. One can also look for enhanced W → νX decays as studied in [20]. The W branching ratios are known to O(0.1%). When combined with the phase space penalty for producing a 3 body final state, this channel is not effective for looking for L i − L j vectors relative to other signals.

ACKNOWLEDGMENTS
We thank Doug Bryman, Simon Knapen, Robert Lasenby, Toby Opferkuch, Maxim Pospelov, Nicholas Rodd, and Ofri Telem for useful discussions. This work was supported primarily by the DOE under contract DE-AC02-05CH11231. The study was initiated and performed in part at the Kavli Institute for Theoretical Physics, which is supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

Appendix A: Goldstone boson equivalence
In this section we prove that the following replacement rule for external X µ : (A1) The switch from X µ to its Goldstone mode ϕ X follows from Goldstone boson equivalence theorem and we take it as a given. Our goal is then to prove that the divergence of the current can be replaced with the mass term, even when the fermions are off-shell. The derivation is in close analogy with that of the Schwinger-Dyson equations.
Consider a general correlation function, Z = D · · · (ψ 1 ψ 2 · · · )e iS(X,ν,ˆ ) (A2) where S contains the kinetic and mass terms for the fermionic and X fields as well as the source term d 4 yj µ (y)X µ (y). We use the shorthands, D · · · ≡ DXDνD and ψ i =ν i (x i ),ˆ i (x i ), or X(x i ). Consider an infinitesimal field redefinition corresponding to a local U(1) Li−Lj rotation, (leaving the vectors untouched). This field redefinition leaves the measure is unchanged since L i − L j transformations are anomaly free. Working in the flavor basis, its clear that the only terms in the action that are changed are the mass and interaction terms, i.e., δS = d 4 y − ∂ µ α(y)j µ (y) − im ijQiνi (y)ν j (y) + h.c.
(A5) We take the neutrinos to be Majorana form here, though a similar equation will hold for Dirac neutrinos.

Appendix B: Additional results
In the main text we only showed the constraints on L e − L µ and L µ − L τ and only for Majorana neutrinos. For those interested in quantitative limits we provide the limits for the other cases here. The summary of limits for L e − L µ , L e − L τ , and L µ − L τ are shown in Fig. 4 in the top, center, and bottom panels respectively. On the left we present the limits for Majorana neutrinos and on the right we present the limits for Dirac neutrinos. The non-enhanced bounds are all shown in gray and come a combination of fifth force/equivalence principle violation searches, cosmology, and long-range forces affecting neutrino oscillation patterns. In orange we show the unitarity bound, in blue the limits from meson decays (dominated by Kaon decays for all cases except L e − L τ for Dirac neutrinos, where K → µνX is suppressed by additional powers of neutrino masses), in pink neutrinoless double beta decay search limits, in green constraints from neutrino annihilations through observations of Big Bang Nucleosynthesis/supernova, and in red the bounds from neutrino decays (dashed are the cosmological bounds).