Relaxing Cosmological Neutrino Mass Bounds with Unstable Neutrinos

At present, cosmological observations set the most stringent bound on the neutrino mass scale. Within the standard cosmological model ($\Lambda$CDM), the Planck collaboration reports $\sum m_\nu<0.12\,\text{eV}$ at 95% CL. This bound, taken at face value, excludes many neutrino mass models. However, unstable neutrinos, with lifetimes shorter than the age of the universe $\tau_\nu \lesssim t_U$, represent a particle physics avenue to relax this constraint. Motivated by this fact, we present a taxonomy of neutrino decay modes, categorizing them in terms of particle content and final decay products. Taking into account the relevant phenomenological bounds, our analysis shows that 2-body decaying neutrinos into BSM particles are a promising option to relax cosmological neutrino mass bounds. We then build a simple extension of the type I seesaw scenario by adding one sterile state $\nu_4$ and a Goldstone boson $\phi$, in which $\nu_i \to \nu_4 \, \phi$ decays can loosen the neutrino mass bounds up to $\sum m_\nu \sim 1\,\text{eV}$, without spoiling the light neutrino mass generation mechanism. Remarkably, this is possible for a large range of the right-handed neutrino masses, from the electroweak up to the GUT scale. We successfully implement this idea in the context of minimal neutrino mass models based on a $U(1)_{\mu-\tau}$ flavor symmetry, which are otherwise in tension with the current bound on $\sum m_\nu$.


Introduction
Neutrino oscillation experiments show that neutrinos are massive and there is flavor mixing not only in the hadronic sector but also in the leptonic sector. At present, three out of the four parameters of the three family leptonic mixing matrix, U , are experimentally determined (θ 12 , θ 23 and θ 13 ) and there is a ∼ 3 σ hint pointing to a CP-violating value of the Dirac CP phase δ [1,2]. Regarding neutrino masses, neutrino oscillation measurements can provide information only about the neutrino mass squared differences ∆m 2 ij = m 2 νi − m 2 νj . In particular, we currently know that ∆m 2 21 0.0086 eV and | ∆m 2 32 | | ∆m 2 31 | 0.05 eV [3][4][5]. Thus, concerning light neutrino masses, there are still two open fundamental questions: i) what is the neutrino mass ordering [6], i.e. m ν3 > m ν1 (Normal Ordering, NO) or m ν3 < m ν1 (Inverted Ordering, IO), and ii) what is the absolute neutrino mass scale or, equivalently, m ν ≡ m ν1 + m ν2 + m ν3 . Near future neutrino oscillation facilities will be able to definitely answer i) but question ii) can not be addressed in this type of experiments.
The best laboratory constraint on the absolute neutrino mass scale comes from the KATRIN experiment that reports the following Feldman-Cousins upper limit: mν e ≡ i |U ei | 2 m 2 νi < 0.9 eV at 95% CL [7]. Using the information from neutrino oscillation data on U , this bound can be translated into the subsequent upper limit on the sum of the light neutrino masses [8]: m ν 2.7 eV at 95% CL. This should be compared with the substantially stronger constraints arising from Cosmology [9][10][11]. In particular, by using Cosmic Microwave Background (CMB) measurements in conjunction with Baryon Acoustic Oscillations (BAO) data, and within the framework of the standard cosmological model ΛCDM, the Planck collaboration reports [12]: m ν < 0.12 eV at 95% CL. Thus, at face value, the cosmological bound from Planck is more than one order of magnitude stronger than the one reported by KATRIN. From the experimental viewpoint, the CMB bound is very robust since systematic effects in Planck CMB measurements have been shown to be very small [12][13][14][15] and because, given a cosmological model, the theoretical error in predicting CMB spectra is negligible [16][17][18][19][20]. However, all cosmological bounds on m ν are cosmological model dependent and m ν < 0.12 eV (at 95% CL) only applies if ΛCDM is the actual model describing our Universe. Observationally, ΛCDM is an extremely successful cosmological model [12] but, of course, when comparing cosmological and laboratory bounds on m ν , one wonders to what extent the cosmological ones depend upon the data set and cosmological model under consideration. In order for the reader to have a feeling of this matter, we summarize in Table 1 a suite of cosmological constraints on m ν arising from analyzing various data sets, and by using the same data set but within different cosmological models. From Table 1 we can draw some important conclusions: 1. Within the framework of ΛCDM current cosmological bounds on m ν are driven by Planck. This is relevant since, as argued before, Planck constraints are very robust.

Cosmological bounds on
m ν are not substantially altered in standard extensions of ΛCDM. In particular, when the modifications of ΛCDM entail non-standard dark energy, non-standard inflationary perturbation spectra, curvature, extra dark radiation, or some of these together. 3. Neutrino mass bounds can considerably be alleviated if neutrinos possess non-standard properties such as a time dependent m ν or if neutrinos decay on cosmological timescales.
In particular, the analysis performed in [21] which includes Planck+BAO+SNIa+SDSS large scale structure (LSS) data, yields m ν < 0.9 eV at 95% CL, provided that neutrinos decay invisibly to massless BSM states with lifetimes in the range 10 −4 τ ν /t U 0.1 1 , where t U = 13.8 Gyr = 4.35 × 10 17 s [12] is the age of the Universe. This bound is to be compared with the one obtained within the ΛCDM model in the stable neutrino framework, m ν < 0.12 eV at 95% CL [12]. This remarkable relaxation has important implications for neutrino mass models since many of them predict neutrinos with masses m ν 0.12 eV. For instance, most of the two-zero neutrino mass textures, and therefore the flavor models realizing these textures, lead to m ν > 0.12 eV [23][24][25]. Along these lines, we refer to [26] for a very recent discussion of the implications of cosmological neutrino mass bounds for models leading to quasi-degenerate neutrinos.
Moreover, upcoming galaxy surveys such as DESI [35] and Euclid [36] will have a 1σ-sensitivity of σ( m ν ) 0.02 eV and are expected to detect neutrino masses within the next ∼5-10 years if the cosmological model describing our Universe is ΛCDM, see e.g. [37][38][39][40][41]. A neutrino mass detection by DESI/Euclid would provide us with extremely valuable information about the neutrino lifetime [42,43], implying τ ν t U and, thus, excluding cosmologically relevant neutrino decay scenarios. Perhaps even more interestingly, if neutrino masses are not detected by DESI/Euclid and only upper bounds on m ν are reported, neutrino decays will become a prime scenario to explain the non-detection of neutrino masses in cosmological observations. Notice that this would have also a crucial impact in the context of neutrinoless double beta decay experiments, since the theoretical prediction for the decay rate strongly depends on both the absolute neutrino mass scale and the neutrino ordering (for two recent reviews on the topic see for instance [44,45] /meV at 95% CL from various cosmological data sets and within various cosmological models. In the m lightest ν column, the value quoted on the left (right) corresponds to NO (IO). The cosmological bounds listed here assume that neutrinos are degenerate, but relaxing this assumption has a negligible impact for these bounds, see e.g. [30,34]. In the last row we show the constraint on mν from the KATRIN experiment, obtained from precision measurements of the electron spectrum from Tritium β decay [7]. *IO excluded at 95% CL.
In this work, motivated by the potential relaxation of cosmological neutrino mass bounds, we set up a comprehensive study of the models that can lead to neutrino decays with lifetimes τ ν t U . We focus on invisible decay channels because the ones with final states that can interact electromagnetically are bounded to have rates τ ν > (10 2 − 10 10 ) t U [46][47][48][49][50][51][52][53]. The possibility of neutrino decays was considered already in the 70's [54], and in fact, two neutrino mass eigenstates decay within the Standard Model [55][56][57] albeit at a very slow rate, τ SM We note that several invisible neutrino decay models have been put forward, see e.g. [58][59][60][61][62][63]. However, to our knowledge, the majority of theoretical investigations of invisible neutrino decays in cosmology were developed prior to the discovery of neutrino oscillations. Therefore, we find timely to develop new models for invisible neutrino decay and to perform an exhaustive identification of the required particle content, coupling strengths, and mass scales of the particles needed to trigger invisible neutrino decays at a rate τ ν t U .
The structure of this work is the following: In Section 2, at the phenomenological level, we categorize the different decay modes of invisible neutrino decays. In particular, in terms of 2-body and 3-body neutrino decay final states as mediated by neutrinophilic scalars and vector bosons. In Section 3, we identify the relevant neutrino decay channels and lifetimes that can potentially lead to the alleviation of the m ν bounds and briefly review the most stringent constraints on invisible neutrino decays. We then map this information to the properties of the particles involved in the decay, identifying the region of the coupling and mass parameter space compatible with τ ν t U and all the relevant present constraints. In Section 4, we propose a simple extension of the seesaw mechanism in which neutrinos naturally can decay loosening the neutrino mass bounds up to m ν ∼ 1 eV, without spoiling the light neutrino mass generation mechanism. Then, we embed this extension within a minimal model based on a spontaneously broken U (1) µ−τ flavor symmetry [64,65], which predicts m ν > 0.12 eV. Finally, in Section 5 we summarize the main results obtained in this study and draw our conclusions.

Taxonomy of Invisible Neutrino Decays
In this section we categorize the various possibilities for invisible neutrino decays, namely, those in which the decay products do not interact electromagnetically. We then write the most general effective Lagrangians parametrizing them.

Categories
We classify invisible neutrino decays according to two criteria: i) nature of the decay products and ii) number of particles in the final state. The categorization for each case goes as follows: i) Nature of the decay products a) At least another active neutrino mass eigenstate ν j .
Since energy and angular momentum conservation ensure the lightest neutrino state to be stable, the cosmological constraint on m ν cannot be relaxed by more than 0.06 eV and 0.1 eV for NO and IO, respectively. We refer the reader to Section 3 for more details.

s) Only BSM species.
In this case, the constraint on m ν can be significantly relaxed to the level of m ν 1 eV at 95% CL, provided that the BSM particles are massless, see Section 3.
ii) Number of particles in the final state 2) 2-body decays.
Angular momentum conservation requires the decay products to be one fermion and one boson with spin 0 (scalar φ) or 1 (vector Z ). Regarding the fermion, we will consider two possibilities: a light neutrino mass eigenstate ν i or a sterile neutrino ν 4 .

3) 3-body decays.
For sufficiently massive bosons, i.e. m φ, Z > m νi , any 2-body decay is kinematically closed. However, such a boson can mediate off-shell a 3-body decay, which becomes the dominant channel.

4) 4-body decays and beyond.
We will not consider this possibility here since, as we will show in Section 3.2, already the 3-body decay channels are only capable of rendering τ ν < t U across a narrow window of the parameter space. (3a1)

2-body Decays:
3-body Decays: Here the indices i, j, k take values from 1 to 3. ν4 refers to a mainly sterile mass eigenstate with small mixing with the active neutrinos να. φ and Z are neutrinophilic scalars and vector bosons, respectively.
Therefore, according to this classification, there are six different possible topologies for the neutrino decays that we show in Figure 1. The diagrams are labelled by their number of final state particles N = 2, 3 and by whether there are active neutrinos in the final state a or only sterile species s.

Effective Lagrangians
Now that we have clarified the particles that can participate in the decays, in this section we will introduce the most general effective Lagrangians that parametrize their interactions. The details about the computation of the decay rates required for the phenomenological analysis are given in Appendix A.1. Approximated expressions of the decay rates for all the channels under consideration in the relevant limits are presented in Table 2.
On the one hand, the most general renormalizable Lagrangian that describes the interactions between the light neutrinos, a sterile state, and a scalar field φ in the neutrino mass basis is given by: where ν i and ν 4 are the mainly active and mainly sterile neutrino mass eigenstates, respectively. Here h and λ parametrize scalar and pseudo-scalar Yukawa couplings between active/sterile neutrinos and φ. For the sake of concreteness, we shall assume that neutrinos are Majorana particles. The formulae for the decay rates shown in Appendix A.1 and Table 2 can be applied to the Dirac neutrino case doing the following mapping: On the other hand, the interaction Lagrangian parametrizing active/sterile neutrino interactions with a Z is: where P L(R) is the left (right) handed projection operator and g L(R) represent left (right) handed couplings.

Invisible Neutrino Decay Rates
Decay Channel Case Decay Rate 2-Body Decays Table 2. Rates for invisible neutrino decays. For neutrinos coupled to a scalar we have considered pseudoscalar couplings only for concreteness. We have simplified the 3-body decay rates by integrating out the mediating boson and φ/Z * means off-shell mediators. See Appendix A.1 for detailed formulae. The reference numbers for each decay are chosen so as to match relevant scales/parameters, see Section 3.2.
Notice that, in order to translate Equations (2.1) and (2.2) to the neutrino flavor basis, active (ν α ) and sterile (ν s ), we just need to perform the correspondent rotation: ν α = U αi ν i + θ α4 ν 4 , ν s = θ si ν i + θ s4 ν 4 . Since the mixing between the sterile and active neutrinos θ α4 and θ si should be small, U αi with α = e, µ, τ and i = 1, 2, 3 is given by the PMNS matrix, up to subleading corrections driven by the active-sterile neutrino mixing.
Given the phenomenological approach of this section, we will assume that the coupling constants present in the effective Lagrangians introduced in Equations (2.1) and (2.2) are independent parameters. However, in UV complete models some couplings may be absent and, if all are present, they are expected to present correlations among them. In Section 4 we will map the effective Lagrangians presented here into UV complete models able to realize invisible neutrino decays in the region of the parameter space that can relax cosmological bounds on m ν . In any case, from the UV perspective, one would typically expect the following possible relationships between λ ij : λ i4 : λ 44 2 : • λ i4 = 0, λ 44 = 0: Scenario with no light sterile neutrinos. This is the simplest possibility and we will consider neutrino mass models realizing invisible neutrino decay of this sort in Section 4.2.
neutrino are charged under a gauge lepton flavor symmetry. This case is expected to generically suffer from stringent constraints since the bosons mediating this interaction will couple to both charged leptons and neutrinos with couplings of the same order.
This should not be considered as a complete list of all possible scenarios but rather as a sample of typical cases. Indeed, in Section 4.1 we will propose a neutrino decay model with λ i4 λ ij and λ 44 = 0, that does not belong to any of the above items in the list.

Relaxing Neutrino Mass Bounds with Invisible Neutrino Decays
The aim of this section is twofold. Firstly, in 3.1, we discuss the present bounds on invisible neutrino decays and their impact on the extraction of cosmological neutrino mass bounds. We also show the viable region of the τ ν − m ν parameter space where the cosmological constraint on m ν can be relaxed. Secondly, in 3.2, we highlight the range of values for the coupling constants and masses involved in the possible invisible decays (classified in Figure 1) in which this goal can be achieved.

Effects on cosmological neutrino mass bounds
In order to understand how invisible neutrino decays can potentially relax cosmological neutrino mass bounds, we will first briefly review the thermal history of neutrinos in the Standard Model and how it is modified if neutrinos decay. This is illustrated in Figure 2 where we show the evolution of neutrino energy density in the Standard Model framework and within several invisible neutrino decay scenarios. The energy density evolution has been computed by solving the relevant Boltzmann equations as in [21]. Note that, for simplicity, in Figure 2 we only show the contribution from the neutrino energy density. However, this does not change our conclusions since even though the contribution from the neutrino decay products is not negligible for t ∼ τ ν , it is diluted as (1 + z) 4 . For example, for m ν = 0.5 eV and τ ν = 10 −4 t U the energy density of the decay products today is similar to that of the photons and represents only 0.5% of the energy density that stable neutrinos with m ν = 0.5 eV would encode. The energy density of the decay products becomes relevant for lifetimes that are somewhat similar to t U and for m ν 0.2 eV. However, as we discuss below, such regions of parameter space are excluded by current cosmological observations.
Within the standard cosmological picture, neutrinos decouple from the SM plasma at T ∼ 2 MeV (t ∼ 0.2 s) and ever since they simply free stream without interacting with any species in the Universe. For temperatures T m e , after electrons and positrons annihilate, the ratio between the neutrino and photon energy densities is fixed and keeps being fixed as long as neutrinos are relativistic, as it can be seen in Figure 2. In particular, T γ /T ν 1.4, which corresponds to N SM eff ≡ 8/7 (11/4) 4/3 ρ ν /ρ γ 3. Eventually, when the average momentum of relativistic neutrinos, given by the Fermi-Dirac distribution, drops below its rest mass, T ν < m ν /3.15, neutrinos become non-relativistic and enhance their contribution to the energy density with respect to the photons. This corresponds to a redshift z 1100 m ν /0.58 eV (see Figure 2). Therefore, neutrinos become non-relativistic after recombination (z 1090) provided that m ν < 0.58 eV. As it can be observed in Figure 2, while neutrinos are still relativistic their contribution to the energy density is independent of their mass. However, after Figure 2. Energy density evolution for several neutrino decay scenarios. The period in which recombination takes place and the ones relevant for the Lyman-α forest and Baryon Acoustic Oscillations data are indicated by the grey regions. We show the contribution to the energy density from photons (red), neutrinos (magenta), dark matter (black), baryons (blue), and dark energy (cyan). Note that for the sake of simplicity we do not show the energy density of the neutrino decay products.
becoming non-relativistic, the total matter energy density gets an additional non-negligible contribution arising from stable neutrinos, which today is given by Ω ν h 2 = m ν /93.14 eV. Essentially, CMB observations are sensitive to m ν via the contribution of massive neutrinos to the total energy density after recombination.
The impact of massive neutrinos on the CMB has been reviewed with great level of detail in the literature, see e.g. [8][9][10][11] 3 . The position of the acoustic peaks in the CMB spectra is extremely well measured and it depends both upon i) the expansion history of the Universe prior to recombination (which is essentially independent of m ν provided m ν 1.5 eV as can be observed in Figure 2) and ii) the angular diameter distance to recombination, that can be modified via the contribution of massive neutrinos to the matter energy density between recombination and today. Thus, perhaps the most direct effect of massive neutrinos on CMB observations, encoded in Ω ν h 2 = m ν /93.14 eV, is due to their contribution to the angular diameter distance to recombination. This effect is, however, strongly degenerate with the value of the Hubble constant. Using data from BAO or local measurements of the Hubble constant this parameter degeneracy can be broken. In addition, when CMB observations are considered in isolation, the most relevant effect of neutrinos that become non-relativistic after recombination (m ν < 0.58 eV) is to affect the lensing of the CMB power spectra via their contribution to the gravitational matter potential between recombination and today. We note that there are other additional effects of m ν = 0 [8][9][10][11] but they are sub-leading, given current data, as compared to the effects mentioned above.
If neutrinos decay while or soon after becoming non-relativistic, their energy contribution to matter is reduced with respect to the stable scenario and can even become completely negligible as it is shown in Figure 2. This would thus relax their main impact on the CMB spectra 4 . In such a case, the cosmological bounds on m ν can be relaxed. Nonetheless, CMB and galaxy survey observations can still be used to constrain the neutrino mass and lifetime in the regime in which neutrinos decay  [22]. The magenta region is excluded by the analysis of late time matter clustering and CMB lensing from [21]. In orange we show the line separating relativistically and non-relativistically neutrino decays given by [21]: 1/τν H0 √ Ωm ( mν /(9Tν,0)) 3/2 . Left panel: bounds on decays of the type νi → νj φ/Z with m φ, Z = 0. Right panel:. bounds on decays of the type νi → ν4 φ/Z with mν 4 = m φ, Z = 0 (the same bounds apply to νi → ν4ν4ν4 decays). Note that the impact of νi → νj φ/Z in late time cosmological observables has not yet been studied in the literature. Such study may yield additional constraints. Therefore, the region of parameter space excluded by the Planck bound on m cosmo ν (see Figure 4) shown in the left panel constitutes an estimated constraint.
non-relativistically [69][70][71][72][73]. In this regime, the mass and lifetime of neutrinos is mainly constrained by the shape of the matter power spectrum as relevant for galaxy surveys and also, although to a less extent, by the CMB lensing spectrum. Very recently, the alleviation of the cosmological m ν bound via neutrino decay has been analyzed in [21]. There, it has been shown that, if neutrinos decay while non-relativistic into invisible massless dark radiation with 10 −4 τ ν /t U 0.1, the cosmological m ν bound can be relaxed up to m ν 0.9 eV. This is illustrated in the right panel of Figure 3 where the bounds from [21] are labeled as "CMB+LSS". We note that there is no study of the impact of relativistically decaying neutrinos on late time cosmological observables such as the matter power spectrum. However, since the energy density evolution in the relativistic decay scenario renders always a smaller energy density than in the non-relativistic decay case, we expect that neutrinos decaying into massless dark radiation would lead to a similar maximum alleviation of the m ν bound in both cases.
Neutrinos with lifetimes τ ν 10 12 s eV/m ν decay while relativistic prior to recombination. In this regime, neutrino decays alter the free streaming property of neutrinos in the early Universe. Given the large energy density encoded in neutrinos prior to recombination, this effect has a strong impact on neutrino perturbations and thus on the CMB spectra which is broadly independent of the cosmological model [74]. Therefore, CMB data can also be used to constrain relativistic invisible neutrino decays [66,67]. The strongest current bound comes from the analysis of Planck legacy CMB observations performed in [22]. Provided all neutrinos decay, and irrespectively of the neutrino decay final state, Ref. [22] reports the following bound on the neutrino lifetime: τ ν > 1.3 × 10 9 s (m ν /0.05 eV) 3 . This constraint is shown in Figure 3 labeled as "CMB ν free streaming". Complementary bounds on the neutrino lifetime from laboratory experiments and astrophysics are summarized in Appendix A.2.
Finally, we note that a detailed analysis of the implications for late time observables of ν i → ν j φ/Z decays is absent from the literature. Such an analysis is beyond the scope of this paper. However,  based on fermion number conservation and the previous discussion on the neutrino contribution to matter energy density, we estimate that the effective sum of neutrino masses inferred from the neutrino energy density today, Ω ν h 2 ≡ m cosmo ν /93.14 eV, will be given by m cosmo ν = 3 m lightest provided that all neutrinos except the lightest one decay at a rate τ ν t U /10. In other words, we estimate that the potential maximum relaxation of the bound on the actual sum of neutrino masses, . This is illustrated in Figure 4 where we show the estimated value of m cosmo ν as a function of the actual sum of neutrino masses m ν . On the one hand, given the Planck bound m cosmo ν < 0.12 eV, if two neutrinos decay via the ν i → ν j φ/Z channel with τ ν t U /10, the region of the actual neutrino masses parameter space that would be excluded is m ν > 0.15 eV ( m ν > 0.17 eV) for NO (IO). In such a case, neutrino decays can only slightly relax the current bound. We note, however, that the constrain we report here should be taken as an estimate since a detailed analysis could lead to a slightly different result 5 . On the other hand, notice that the relaxation of the neutrino mass bound is larger for small values of m ν with a maximum alleviation given by a 0.06 eV (0.1 eV) shift for NO (IO) in the massless lightest neutrino limit. This might be very relevant for upcoming galaxy surveys as DESI or Euclid since if, lets say, a bound m cosmo ν < 0.04 eV is eventually reported it would not necessarily exclude NO or IO, being compatible with both NO and IO in the neutrino decay scenario, with m ν < 0.08 eV ( m ν < 0.12 eV) for NO (IO).
We can conclude from this section that decays into purely massless dark radiation represent the best avenue to relax current cosmological constraint on m ν . Neutrino decays including active neutrinos in the final state can potentially play a relevant role in the near future with upcoming data from DESI/Euclid. This can be clearly observed already in Figure 3, where we show (in white color) the viable region of the τ ν − m ν parameter space where the cosmological bound on m ν can be relaxed. In the next section we will explore the range of values for the coupling constants and masses of the particles involved in the decays for which this can be realized.

Allowed parameter space
Using the effective Lagrangians considered in Subsection 2.2, in Appendix A.1 we have computed the invisible neutrino decay rates for all the relevant channels summarized in Figure 1. In this section, we will determine the region of couplings and masses parameter space that can lead to τ ν < t U , including the most constraining model independent bounds on neutrinophilic bosons and invisible neutrino decays. According to the categorization presented in Section 2.1, and following a model independent approach, this will be done separately for each decay channel using the decay rates computed in Appendix A.1. In order to have a feeling of the main dependence of these channels on the different parameters we display in Table 2 a summary of approximated expressions for the neutrino decay rates in some relevant phase space limits. Table 2 has been built so as to highlight relevant scalings with masses and couplings for interesting cases, i.e., for neutrino masses O(0.3 eV) when the decay channel allows m ν ∼ 1 eV and O(0.05 eV) in the case of decay to one active neutrino, which could be relevant for future surveys. We note that there is no interference between the different decay diagrams and the rates only depend upon the modulus of the individual coupling constants. Therefore, to simplify the notation, throughout this section we shall use λ, h and g to denote their absolute value.

2-Body Decays:
Case2a. In a realistic scenario, once the possibility of a neutrino decay is given, we expect all possible decay channels to be opened, i.e., ν 3 → ν 1,2 φ/Z and ν 2 → ν 1 φ/Z (ν 2 → ν 1,3 φ/Z and ν 1 → ν 3 φ/Z ) for NO (IO). Furthermore, theoretically one would expect that once couplings are turned on they are of similar strength, given the observed mixing among active neutrinos. The resulting allowed parameter space for a pseudo-scalar φ and Z in the final state is shown in Figure 5. In Appendix A.3 we show the result for each individual mode, including also the case of φ being a scalar boson, assuming that the rest of the channels are switched off. Notice that, as we have discussed in the previous section, the maximum relaxation of the m ν bound is obtained for the case in which two neutrinos decay (as it is shown in Figures. 4 and 5). The present Planck constraint can only be very slightly alleviated by 0.03 eV (0.05) for NO (IO).
Case 2a: Pseudo-scalar couplings in ν i → ν j φ decays. From the left panels of Figure 5 (see also the detailed plots for each individual mode in Appendix A.3), we can extract the following two pieces of information: i) CMB bounds on invisible neutrino decays restrict pseudo-scalar couplings to be: λ 31 , λ 32 10 −10 , and λ 21 10 −8 ; ii) the present cosmological neutrino mass bound can only marginally be alleviated for coupling strengths in the range 10 −15 λ 31 , λ 32 10 −10 and 10 −13 λ 21 10 −8 , which give rise to decays occurring on timescales shorter than the age of the universe without substantially altering neutrino free streaming in the early Universe.
Case 2a: Scalar couplings in ν i → ν j φ decays. As far as the mass of the active neutrino in the final state can be neglected, the results are very similar to the pseudo-scalar case, becoming well differentiated when the decaying neutrino mass is close to the final state neutrino mass, due to the extra suppression for degenerate neutrinos in the case of pseudo-scalar couplings, see Equation (A.1). We find that CMB constraints on invisible neutrino decays bound the couplings responsible for neutrino decay to be h 31 , h 32 10 −11 (0.05 eV/m ν ) and h 21 5×10 −11 (0.05 eV/m ν ) (see Figure 9). We also find that couplings in the range 10 −15 h 31 , h 32 10 −11 (0.05 eV/m ν ) allow a slight relaxation of the cosmological m ν bound, in the same way as we discussed in the previous scenario. Case 2a: ν i → ν j Z decays. From the right panels of Figure 5 we can see that CMB bounds on invisible neutrino decays imply m Z /g L 31 , m Z /g L 32 1 GeV and m Z /g L 21 10 MeV. We find that the cosmological bounds on m ν could be slightly relaxed as in the scalar cases, provided 1 GeV m Z /g L 31 , m Z /g L 32 25 TeV. In order to plot our results as a function of m Z /g L ij , we have  Figure 5. Parameter space for case 2a (see Figure 1). Left panel: νi → νj φ decays with φ being a pseudoscalar boson. Right panel: νi → νj Z decays. Here we have assumed that all the couplings are different from zero and all neutrinos decay except the lightest one. The blue region is excluded by the constraint on invisible neutrino decays from Planck CMB observations [22]. In dashed green we highlight the current bound on mν within the ΛCDM [12] and stable neutrino framework. The arrows indicate to what extent it can be relaxed, see also Figure 4. The grey region is excluded by the bound m cosmo ν < 0.12 eV reinterpreted in the neutrino decay scenario.
used the leading order expanded expression for the decay rate as given in Table 2, where m Z /m ν ∼ 0 is the expansion parameter. This will allow us to easily make the connection between the effective coupling m Z /g L ij and a vacuum expectation value within a concrete model realization. The parameter spaces arising when this approximation or the full formula are considered only differ in a tiny region of m Z (where m Z has a significant impact in the phase space). Our approximation is thus in excellent agreement with the full formula given in Equation (A.2).
Case2s: This case corresponds to ν i → ν 4 φ/Z , with φ being a scalar or pseudo-scalar and is shown in Figure 6. Both, scalar and pseudo-scalar, give rise to the same decay rate if m ν4 , m φ,Z m νi and h = λ.
Case 2s. ν i → ν 4 φ decays. We display the relevant parameter space for this scenario in the left panel of Figure 6 for the case m ν4 = m φ = 0. The upper bound is driven by the modification of the CMB due to neutrino free streaming when relativistic neutrino decays are included, and is given by λ i 4 Z Figure 6. Parameter space for case 2s (see Figure 1). Left panel: parameter space for νi → ν4 φ decays. Right panel: parameter space for νi → ν4 Z decays. Both correspond to the 2s case in Figure 1. We have chosen m φ, Z mν i . For νi → ν4 φ we consider only pseudo-scalar couplings but note that the rates are identical for scalar couplings. In orange we show the line separating relativistically and non-relativistically neutrino decays given by [21]: 1/τν H0 √ Ωm ( mν /(9Tν,0)) 3/2 . The blue region is excluded by CMB constraints on relativistically decaying neutrinos [22]. The magenta region is excluded by cosmological bounds on nonrelativistically decaying neutrinos [21]. The Planck bound on mν extracted within ΛCDM is indicated by the dashed green line. The constraint from [21] does not saturate the Planck bound for τν tU simply because different data sets were considered in the analysis of [21].
requiring the decay to happen faster than the age of the universe, we identify a large plateau in the range 10 −15 λ 4i 10 −12 (0.3 eV/m ν ) where the cosmological neutrino mass bound can be alleviated up to m ν ∼ 1 eV (see also Figure 3). Case 2s: ν i → ν 4 Z decays. The region of parameter space m Z /g L i4 < 100 GeV(m ν /0.3 eV) 3 is excluded by the CMB constraints on relativistically decaying neutrinos, as it can be observed in the right panel of Figure 6. Cosmological bounds on non-relativistically decaying neutrinos exclude the region given by m ν > 0.25 eV and m Z /g L i4 10 TeV. In summary, we find that for 100 GeV(m ν /0.3 eV) 3 m Z /g L i4 10 TeV, this decay mode can lead to an amelioration of the cosmological neutrino mass bounds up to m ν 1 eV.

3-Body Decays:
The 3-body decays rates are fairly independent upon the nature of the bosonic mediator. In particular, the case of Z mediated decays and φ mediated decays with scalar couplings are fully analogous to the φ mediated decay with pseudo-scalar couplings, given the following approximate mapping: Therefore, we only discuss decays solely realized via pseudo-scalar couplings (λ = 0 , h = 0). But the same conclusions that we will draw here for the pseudo-scalar case apply to the other scenarios. Several bounds apply to sub-MeV neutrinophilic bosons as relevant for 3-body invisible neutrino decays. The most relevant ones are: i) CMB constraints on neutrino free streaming that exclude ν i − φ/Z couplings as small as 10 −13 for 0.1 eV < m φ/Z < 300 eV [75], ii) BBN constraints, and iii) supernova cooling constraints [76]. There are several studies dealing with BBN and SN1987A bounds on neutrinophilic bosons, see e.g. [22,77,78] and [79][80][81][82][83][84], respectively. However, given the wide range of masses and types of couplings we consider in our analysis, we have rederived the bounds using our effective Lagrangians, see Equations (2.1) and (2.2). The derivation is described in Appendix A.5.  Figure 1). We have assumed that mν 4 = 0 and that active neutrinos have universal couplings to φ. Right panel: parameter space for νi → νjφ → ν4ν4ν4 decays (case 3s in Figure 1). We show laboratory bounds from KamLAND-Zen [85], CMB constraints on ν-φ interactions from Planck CMB observations [75], and SN1987A and BBN constraints on neutrinophilic bosons -see Appendix A.5 for the derivation of the latter two. Note that for decays happening at the same rate, the difference in the final phase space enforces larger couplings for increasing active neutrino mass if one active neutrino is in the final state, whereas the opposite occurs in the case with only sterile final states. Identical conclusions can be drawn for both scenarios but mediated by a Z .
Case 3a0: ν i → ν jνk ν k . In this case, τ ν < t U requires a light and not very weakly coupled neutrinophilic boson with λ ij 0.6 m φ /10 keV which is already excluded by meson decays, double beta decay searches, and various astrophysical and cosmological constraints (see left panel of Figure 7), and thus this channel is not phenomenologically viable.
Case 3a1: ν i → ν jν4 ν 4 . This decay mode is controlled by λ ij and λ 44 , see Table 2. λ ij is constrained by various laboratory experiments, astrophysics and cosmology. However, λ 44 is largely unconstrained since it parametrizes the interaction between a sterile neutrino and a neutrinophilic boson. In the left panel of Figure 7 we show the most favorable region in which τ ν < t U can be realized. As it can be seen in the figure, even saturating λ 44 = 4π, the bounds from the CMB, BBN and SN1987A on sub-MeV neutrinophilic mediators exclude τ ν < t U for this type of topology.
Case 3s ν i → ν 4ν4 ν 4 . The right panel of Figure 7 shows the relevant parameter space for this decay channel. The couplings controlling this decay mode (λ i4 and λ 44 ) are not particularly well constrained. The strongest bounds on λ i4 arise from BBN data and supernova cooling, as shown in Figure 7, while λ 44 is essentially unbounded. In the right panel of Figure 7 we show the most favorable case with λ 44 = 4π. We can clearly appreciate that for ν i → ν 4ν4 ν 4 decays to be cosmologically relevant, mediators with m φ/Z 10 keV and coupling strengths λ i4 > 10 −14 4π λ44 m φ/Z /eV 2 0.3 eV/m ν are required. For λ i4 = λ 44 , however, only a tiny window with m φ/Z 10 eV would be viable. Therefore, at the phenomenological level, ν i → ν 4ν4 ν 4 decays can still satisfy τ ν < t U for m φ/Z 10 keV. It is important to remark that, even tough the decay rate is independent of λ ij at the model independent level, in a UV complete model λ ij is expected to be present and correlated via mixing with the parameters driving the decay. Since this coupling is subject to stringent cosmological and astrophysical constraints (see left panel of Figure 7), ν i → ν 4ν4 ν 4 decays can only render τ ν < t U provided λ 44 , λ i4 λ ij .

Summary:
We have identified the cosmologically relevant regions of parameter space in our effective Lagrangians parametrizing invisible neutrino decay, see Equations (2.1) and (2.2). We have shown that a wide region of parameter space is excluded by current CMB, BBN and SN1987A constraints. Furthermore, we have also highlighted the range of values for the relevant couplings and masses that can lead to τ ν < t U and in which, therefore, invisible neutrino decays can relax cosmological bounds on m ν . We find that 2-body decays into new light BSM species can alleviate such bounds in a substantial region of the parameter space. Finally, our results show that among all possible 3-body decays, only those of the type ν i → ν 4ν4 ν 4 can lead to τ ν < t U and only in a narrow window of the parameter space.
It is worth mentioning that although 2-body decays with one active neutrino in the final state can only very slightly relax the current Planck bound, such decays may become relevant if future galaxy surveys such as DESI or Euclid do not detect neutrino masses at all, or point towards a NO spectrum. From Figure 4, we see that the absence of a neutrino mass signal could be a hint of unstable neutrinos. Moreover, a measurement of m cosmo ν = 0.06 eV, which would be interpreted as a signal of NO in the standard framework, could instead correspond to IO if the active neutrinos decay to the lightest one. This would be particularly relevant in the context of neutrinoless double beta decay experiments since the theoretical prediction for the decay rate crucially depends on both the neutrino ordering and the absolute neutrino mass scale.

Models for Neutrino Masses and Invisible Neutrino Decays
A large number of models generating the neutrino mass and mixing pattern observed in neutrino oscillation experiments are in strong tension with cosmological bounds on m ν . For instance, 5 out of 7 potentially viable two-zero neutrino mass textures predict m ν > 0.12 eV [25]. This tension with cosmology generically occurs also for inverse Majorana neutrino mass matrices with two zeros textures [23,24]. Any model generating these textures would thus be disfavored by cosmology as it is the case, for instance, for minimal U (1) µ−τ neutrino mass models [86][87][88][89]. Furthermore, plenty of generic models do not necessarily strictly predict m ν > 0.12 eV, but accommodate m ν > 0.12 eV across large regions of the parameter space (see [90][91][92][93][94][95][96][97][98] for reviews on neutrino mass models).
In this section, we first propose a simple extension of the seesaw mechanism with an extra scalar field and a sterile neutrino state, both charged under an additional U (1) X symmetry, in which ν i → ν 4 φ decays with m νi m ν4 m φ = 0 can naturally occur, alleviating the neutrino mass bounds up to m ν ∼ 1 eV. Then, we embed this extension within a minimal neutrino mass model based on a spontaneously broken U (1) µ−τ flavor symmetry [64,65], which predicts m ν > 0.12 eV [86][87][88][89] and therefore is in strong tension with the current Planck bound.

A simple model for ν i → ν 4 φ decays within the seesaw framework
Here we will study a simple extension of the seesaw mechanism able to open new active neutrino decay channels. In particular the 2s channel (ν i → ν 4 φ) that can lead to a substantial relaxation of the cosmological neutrino mass bounds, see Figure 6. The presence of this decay mode necessarily requires to introduce new fields in order to account for the extra light sterile neutrino and the light boson in the decay final state. Of course, the new sector should not spoil the light neutrino mass generation mechanism and, in particular, the new sterile state can not present a large mixing with the active neutrinos in order to satisfy bounds on light sterile neutrinos. From the theoretical point of view, this means that a new extremely light (or massless) sterile neutrino state, essentially decoupled from the light neutrino masses, needs to be generated.
Our proposal is to enlarge the usual content of the seesaw model (i.e., three right-handed neutrinos, N R ) by adding an extra fermion singlet S L , and one complex scalar singlet Φ. To avoid large couplings with the active sector and generate a massless sterile state, we also extend the symmetry of the model with a new global U (1) X . Charging S L and Φ with opposite U (1) X charges and leaving the rest of the fields uncharged, the only new terms allowed by the symmetry are of the form yΦN R S L . Notice that, at tree level, the Majorana mass term S c L S L and Yukawa coupling between S L and the active neutrinos are forced to be zero by the symmetry. U (1) X is dynamically broken when Φ takes a vacuum expectation value (vev) v Φ and, since we chose U (1) X to be global, one of the two components of Φ will remain massless (the Majoron) while the other will have a mass of the order of v Φ . After symmetry breaking, the complete 7 × 7 neutrino mass matrix in the flavor basis is given by where m D and M R are the 3 × 3 Dirac and Majorana matrices already present in the seesaw model, while Λ is a 3 × 1 matrix whose elements are given by Λ a = y a v Φ . Diagonalizing the mass matrix under the assumption Λ m D M R , we obtain Notice that for the sake of simplicity, we have ignored the dependence on the PMNS leptonic mixing matrix U . The light neutrino masses and mixing are generated through the seesaw mechanism, and the mixing among the active and heavy states remains the same as in the standard seesaw. There is a zero mode associated to the new sterile state that presents a mixing with ν i given by Λ/m D and with the heavy N i of Λ/M R . The corrections from the new sector to the active neutrino mass matrix are remarkably suppressed, O(Λ 2 /M R ), and thus negligible as long as Λ m D . Obviously, the same requirement leads to a small mixing between the sterile state and the active neutrinos, avoiding any potential tension with neutrino oscillations or BBN data. In particular the sterile-active neutrino mixing, θ α4 ∼ Λ/m D , is best constrained by BBN data that roughly requires |θ α4 | 2 4 × 10 −6 eV/ |∆m 2 4i | (see Appendix A.6) and neutrino oscillation experiments that set |θ α4 | 2 O(0.01) [99][100][101]. Remarkably, since in our model m 2 ν4 − m 2 νi < 0, the sterile neutrino will resonantly be produced via mixing in the early universe, analogously to the MSW-effect and contrary to the widely studied case with m 2 ν4 −m 2 νi > 0. Consequently, the BBN constraint that we have derived in our scenario is roughly two orders of magnitude more stringent than the one derived for the m 2 ν4 − m 2 νi > 0 case [102][103][104] 6 . Given that |∆m 2 4i | = m νi 1 eV, BBN thus gives a stronger bound than neutrino oscillation data: θ α4 ∼ Λ/m D 10 −3 . Therefore, considering a mild hierarchy Λ/m D 10 −3 the model is in agreement with the constraints from the active neutrino sector.
The structure of the mass matrix given in Equation (4.1), also called Minimal Extended Seesaw (MES), has been proposed in order to address the hints of light sterile neutrinos (but still heavier than the active ones) in some oscillation experiments, and may arise as a consequence of some flavour symmetries [105][106][107][108], or a dark U (1) [109]. Notice that the rank of the matrix is 6, and therefore there is always a massless neutrino at tree level. However typically the hierarchy Λ m D is assumed, so the massless state is mainly active, while in our scenario Λ m D leads to a mainly sterile massless neutrino.
In order to do a mapping to the parameters of the phenomenological Lagrangian, we write now the coupling between the N R and the new fields in the mass basis where φ (σ) is the massless (massive) component of Φ and, for the sake of simplicity, we are again ignoring the dependence on the PMNS mixing matrix. The first line of the above equation give us the mapping to the phenomenological parameters from (2.1): Notice that the hierarchy among the λ couplings does not coincide with the naive expectation λ ij ∼ λ i4 θ ∼ λ 44 θ 2 . In our model, the dominant coupling of the light neutrino (active and sterile) sector with the massless scalar field is given by λ i4 , and is controlled by the mixing between the active and heavy neutrinos. Thus, if the standard seesaw mechanism is at work, we have λ i4 y · 10 −13 m ν 0.1eV matching the required range of values for the 2s channel (ν i → ν 4 φ) to relax the m ν bound (see Figure 6), as we will discuss in more detail in a concrete neutrino mass model below. Note that, if the massive scalar σ has a mass m σ m ν it could potentially mediate 3-body decays. In Section 3.2 we concluded that only the 3s channel is phenomenologically viable, and only provided that h ij h 44 , h i4 . However, in our neutrino decay model h 44 0 which precludes the possibility of realizing the 3s channel.
Our proposal provides an extra massless sterile neutrino, avoiding at the same time large corrections to the light neutrino sector, thanks to the addition of an extra U (1) X global symmetry. This prevents S L to couple to the active neutrinos and forbids the Majorana S L self coupling. A pertinent question arises in this context: is the model radiatively stable? The most relevant 1-loop effect is driven by the finite correction to the Majorana mass term S c L S L δM LL ∼ y 2 (4π) 2 which is completely negligible for the parameter space that can alleviate the m ν cosmological bound tension via active neutrino decay. The model is thus stable at the radiative level, and can be easily embedded in any extension of the standard seesaw scenario. Finally, it is worth to point out that although we have considered a global U (1) X symmetry, the ideas presented here can also be applied to a feebly coupled gauged U (1) X . In such a case, at least two extra fermion singlets S L with opposite U (1) X charges are required to cancel anomalies, while the corresponding gauge boson X associated to the U (1) X symmetry should be extremely light to allow for the 2-body decay ν i → ν 4 X.

Neutrino masses and decays within a U (1) µ−τ symmetry
In order to illustrate how the model independent results presented in Section 3.1 can be mapped to UV complete models, in this section we consider neutrino mass models based on a spontaneously broken U (1) µ−τ flavor symmetry [64,65] which are potentially in tension with current Planck data. This type of flavor symmetries are easily anomaly free, given the Standard Model particle content [110,111], and considering the right choice of charges for the right-handed neutrinos introduced in order to account for light neutrino masses. Neutrino mass models based on a U (1) µ−τ flavor symmetry correctly fit the light neutrino mass and mixing pattern observed in neutrino oscillation experiments. However, they generically predict m ν > 0.12 eV [86][87][88][89]. For the sake of concretness and illustration, we shall concentrate on the most minimal U (1) µ−τ SM extension featuring a SM singlet scalar ϕ and 3 heavy sterile neutrinos with U (1) µ−τ charges Q(ϕ) = +1, Q(N e ) = 0, Q(N µ ) = +1, and Q(N τ ) = −1. In this set up, after electroweak and U (1) µ−τ symmetry breaking, active neutrinos get a mass via the seesaw mechanism. The resulting pattern of masses and mixings within this framework has been analysed in detail in [86][87][88][89]. For this particular setting, only NO is consistent with the measured neutrino sector parameters, and the prediction for m ν is mainly dependent upon the value of θ 23 [87]. For values of θ 23 = 48.3 •+1.1 • −1.9 • (corresponding to the best fit ± 1 σ from NuFIT 4.1 without SK data 7 [3]) one finds that 0.17 eV < m ν < 0.47 eV, while for θ 23 = 40.8 • → 51.3 • (corresponding to the 3 σ range) 0.13 < m ν < 2.7 eV, where the upper bound is given by KATRIN. Therefore, the model is in strong tension with the Planck bound on m ν . In what follows, we will consider a weakly coupled realization of this minimal U (1) µ−τ neutrino mass model and highlight the allowed regions of parameter space for which neutrinos decay within cosmological timescales being able to relax the m ν bound from cosmology. We will first study the case in which U (1) µ−τ is promoted to a gauge symmetry where the ν i → ν j Z decays become the most relevant channels. We then consider the case in which U (1) µ−τ is a global symmetry (à la [60]) and the main decay channels are of the type ν i → ν j φ, where φ is the Goldstone mode associated to ϕ upon U (1) µ−τ spontaneous symmetry breaking. Although from the discussion in Section 3.2 we expect that for these scenarios the neutrino mass bound can only be slightly alleviated, since the decay products contain always an active neutrino and they are quite heavy, we use them as an illustration of the matching of a UV completed model to our effective Lagrangian parameters. Moreover, some of our conclusions are easily applicable to other realizations of the U (1) µ−τ symmetry which lead to slightly lighter active neutrinos, where such decays could relax the tension with the Planck bound more significantly [89]. Finally, we will implement the simple extension described in the previous section in the context of the minimal U (1) µ−τ model, showing that in this case the neutrino decay to only BSM species can naturally occur, ameliorating the neutrino mass bounds up to m ν ∼ 1 eV. Gauge Case: The U (1) µ−τ interaction is mediated by a Z boson. The corresponding Z -neutrino interaction arises through the coupling with the SU (2) lepton doublets and is given by where in the second term we have explicitly written the Lagrangian in the mass basis as in (2.2). In order to illustrate the strength of the couplings that play a role in the relevant decay channels, we map the parameters in Equation (4.7) to the phenomenological parameters of Equation (2.2), introducing the bestfit values of the PMNS mixing parameters from NuFIT 4.1 [3] for NO: |g L 12 |/g µ−τ 0.14 , |g L 13 |/g µ−τ 0.54 , |g L 23 |/g µ−τ 0.82 , (4.8) where the g L ij couplings with i = j lead to invisible neutrino decays of the type ν i → ν j Z provided m Z < m νi − m νj . In the model under consideration, neutrino oscillation data can only be fitted for NO and approximately degenerate light neutrinos. Thus, given the values of the couplings above, ν 3 → ν 2 Z and ν 3 → ν 1 Z decays occur at roughly the same rate.
As can be seen from Table 2, the total rate of invisible neutrino decays depends only on the light neutrino masses, g µ−τ and m Z for which we assume m Z m ν (see e.g. [112] for theoretical arguments motivating why this could be the case). Using Equation (4.8), the results shown in Figure 5 (right panels) can be easily mapped to the present model. As also recently pointed out in [113], we find that bounds on relativistically decaying invisible neutrino decays represent the most stringent constraint on the considered parameter space 10 −10 eV m Z m ν 8 . In particular, the vacuum expectation value of the scalar field should fulfil v µ−τ ≡ m Z /g µ−τ > 1 GeV for m Z m ν . We also find that for v µ−τ < 25 TeV neutrinos have a lifetime τ ν < t U . Therefore, for 1 GeV < v µ−τ < 25 TeV and m Z m ν , in this model neutrino decays are able to slightly relax the tension between the model prediction for m ν and the Planck bound. Global Case: If U (1) µ−τ is a global symmetry, the CP-odd component of the ϕ scalar field becomes one of the final states in the 2-body decay channel 2a 9 . Upon spontaneous symmetry breaking of U (1) µ−τ , the CP-odd component of the scalar field, i.e. φ, becomes a Goldstone boson. This scalar φ can be seen as a flavorful massless majoron with m φ = 0. Quantum gravity is expected to break all global symmetries [114] leading to m φ = 0. However, given that the exact way in which gravity breaks global symmetries is highly dependent upon the unknown nature of the gravitational theory at the Planck scale, we shall assume that m φ m ν . In this setting, h ij = 0 in the effective Lagrangangian given in Equation (2.1) since φ is a CP-odd scalar. Furthermore, the effective coupling λ ij , which controls the ν i → ν j φ decay, is generated via mixing between the heavy sterile and active neutrinos through the interaction L ⊃ −Y eµ ϕN c e N µ − Y eτ ϕ † N c e N τ , where N α are the heavy right handed neutrinos and Y eµ and Y eτ are Yukawa couplings. Taking into account the canonical seesaw scaling of the mixing between the heavy states and the active neutrinos θ N ν m ν /M N , and assuming that Y eµ ∼ Y eτ , we can estimate To be in accordance with the observed light neutrino mass spectrum and mixing, the contribution to the right-handed neutrino Majorana masses from the coupling to the scalar should be of the same order of the mass terms allowed by the U (1) µ−τ symmetry, M ee N c e N e and M µτ N c µ N τ , as we have considered in the above equation.
Using Equation (4.9), the allowed parameter space shown in the left panels of Figure 5, maps to v µ−τ > 1 GeV (neutrinos not spoiling Planck CMB observations) and corresponds to v µ−τ 10 TeV (active neutrino lifetime satisfying τ ν < t U ). Interestingly, the relevant range for the energy scale of the model lies around the electroweak scale, v µ−τ v H . Hence, similarly to the gauge case, such energy scales point towards a low-scale realization of the seesaw mechanism at O(TeV). Decay to sterile neutrinos: Now we propose a particular realization of the scenario described in Section 4.1 within the U (1) µ−τ minimal model (with either gauge or global symmetry), that opens the possibility of active neutrinos decaying to a lighter, mainly sterile state. Recall that the new species are one extra fermion singlet S L and one complex scalar singlet Φ, which have opposite charges under a global U (1) X symmetry. We assign zero U (1) µ−τ charge to S L , in order to keep the anomaly cancellation, while Φ could carry U (1) µ−τ charge. There are only three viable choices for the U (1) µ−τ charge of Φ that allow the required interaction terms yΦN R S L : From the mapping given in Equation (4.4), we see that the coupling driving the ν i → ν 4 φ decay is given by λ i4 y m D /M R , which is linked to the light neutrino mass scale via Equation (4.5) since the seesaw mechanism is at work. In Figure 8 we show the relevant parameter space in our scenario, as well as the theoretical prediction for m ν . We see that our extension can evade the bound on m ν from cosmology across a large region of the parameter space. In particular, for y ∼ 1, our model can naturally evade the bound on m ν for 10 10 GeV M R 10 19 GeV, which points to the canonical seesaw. Such a high scale for M R avoids any potential constrain from the early universe on the new Lagrangian terms given in the second line of Equation (4.3). Any sizeable contribution from the heavy N R to low energy particle physics experiments would also be suppressed. Notice that lighter M R scales can also be easily realized by simply considering smaller values of the y coupling.
In summary, our model provides a natural framework to alleviate the m ν cosmological bound via ν i → ν 4 φ decays (2s). In this scenario, light neutrino masses are generated via the canonical seesaw mechanism and, remarkably, are not spoiled by the presence of the massless sterile neutrino that mixes with the active sector. This is possible just invoking a hierarchy between the vev of the new scalar field and the Higgs vev of v Φ /v H 10 −3 .

Summary, Conclusions and Outlook
The cosmological bounds on the absolute neutrino mass scale within ΛCDM are very stringent. In particular, the Planck collaboration reports m ν < 0.12 eV at 95% CL [12]. This bound, taken at face value, is more than one order of magnitude stronger than the one reported by the KATRIN experiment [7], which can be reinterpreted as m ν < 2.7 eV at 95 % CL. Even though cosmological neutrino mass bounds are cosmological model dependent, in typical extensions of ΛCDM the bound is only slightly relaxed to m ν 0.2 eV, see Table 1. Furthermore, the constraints depend on the data set under consideration, but the ones driven by Planck CMB observations are quite robust (see Section 1). This puts several neutrino mass models in trouble since a large number of them predict m ν 0.12 eV.
In this context, invisible neutrino decays -namely, those in which the decay products do not interact electromagnetically -represent a particle physics avenue to relax the potential tension between theoretical model predictions and cosmology. Indeed, the m ν bound can be loosen up to m ν 1 eV for neutrinos decaying on cosmological timescales [21]. Motivated by this, in this paper we have carried out a comprehensive analysis of models that can accomplish neutrinos decaying invisibly within cosmological timescales. For this purpose, in Section 2, we have presented a taxonomy of all possible invisible neutrino decays as highlighted in Figure 1. At the phenomenological level, we have then written down general effective Lagrangians which describe these decay channels in Equations (2.1) and (2.2). In Section 3, after reviewing the most relevant constraints on invisible neutrino decays, we have identified the parameter space compatible with current limits that can lead to a relaxation of the cosmological m ν bound for each possible decay channel. The main findings of our phenomenological study are: • 2s invisible decay channels (see Figure 1) of the type of ν i → ν 4 φ with m φ = m ν4 = 0, can produce a significant relaxation of the neutrino mass bounds up to m ν 1 eV while satisfying all known present constraints, for ν i −ν 4 −φ interaction strengths given by 10 −15 λ 10 −10see Figure 6. An analogous result is obtained for the ν i → ν 4 Z decay scenario, provided that 1 GeV m Z /g L i4 25 TeV and m Z m νi .
• 2a neutrino decays driven by the ν i → ν j φ modes, with m φ m νi , can only be able to slightly relax current cosmological bounds while simultaneously being compatible with CMB constraints on invisible neutrino decays -see Figure 5. This requires ν i − ν j − φ coupling strengths in the 10 −15 λ 10 −10 range. We reach a similar conclusion for ν i → ν 4 Z decays provided that 1 GeV m Z /g L ij < 25 TeV with m Z m νi .
• Due to a combination of terrestrial, astrophysical, and cosmological constraints on sub-MeV neutrinophilic bosons, 3-body neutrino decays involving at least one active neutrino in the final state (cases 3a0, 3a1 and 3a2 in Figure 1) cannot lead to τ ν < t U . The only topology that remains phenomenologically viable involves only sterile neutrinos as final states and is given by ν i → ν 4ν4 ν 4 (3s in Figure 1). However, this requires a particular hierarchy among the active-sterile-mediator couplings, λ 44 , λ i4 λ ij , and the mediating particle to be lighter than ∼ 10 keV.
In light of the above, we believe that near future experiments provide further motivation to study the phenomenology of neutrinos decaying within cosmological timescales. Indeed, the upcoming galaxy surveys DESI [35] and Euclid [36] will have a 1σ-sensitivity of σ( m ν ) 0.02 eV and are therefore expected to detect neutrino masses in cosmology within the next 5-10 years if the cosmological model describing our Universe is ΛCDM. If such a detection is made, it will lead to a bound on the neutrino lifetime τ ν t U [42,43] and our study shows that wide regions of parameter space in decaying neutrino models would be excluded. On the other hand, the possibility that both DESI and Euclid will not measure m ν , reporting only upper bounds, is even more interesting. Were that to be the case, invisible neutrino decays with τ ν < t U will become a prime scenario to understand such measurements. In our analysis we have already identified the minimal particle content and relevant parameter space required to implement such decays. Furthermore, in the next ∼5 years, the KATRIN experiment is expected to reach a 90% CL sensitivity to mν e of 0.2 eV [7,115]. If KATRIN reports a neutrino detection 10 and only upper bounds on m ν arise from cosmology, invisible neutrino decays will become a very well motivated possibility to reconcile both sets of data.
From the model building perspective, in Section 4, we have proposed a simple extension of the seesaw scenario which contains a (mainly) sterile neutrino lighter than the (mainly) active ones, together with a massless Goldstone boson φ from a spontaneously broken U (1) X global symmetry, so that invisible neutrino decays of the type ν i → ν 4 φ can be cosmologically relevant. Remarkably, in this set up the mixing between active neutrinos and ν 4 are naturally small and the new sterile state does not spoil the standard seesaw mechanism for generating active neutrino masses, so it could be easily embedded in more complicated models within this framework.
As an illustration, we have considered minimal neutrino mass models based on a spontaneously broken U (1) µ−τ flavor symmetry, which generically predict m ν > 0.12 eV and are thus in tension with the current cosmological bound on neutrino masses. We have found that in feebly coupled realizations of such models neutrino decays are possible, and are in somewhat better agreement with cosmological observations. Our analysis shows that the cosmological m ν bounds are slightly relaxed for U (1) µ−τ symmetry breaking scales v µ−τ 25 TeV, while v µ−τ > 1 GeV is required in order to satisfy current CMB constraints on invisible neutrino decays.
Furthermore, we have implemented the minimal extension of the seesaw model described above within the context of the U (1) µ−τ flavor symmetry. In this case the tension between the prediction for m ν from U (1) µ−τ models and cosmology can be fully evaded, as it is shown in Figure 8. Within this scenario, τ ν < t U can be associated to energy scales 1 GeV v µ−τ 10 19 GeV. Such energy scales suggest that UV complete models featuring invisible neutrino decay modes could be linked to phenomena such as baryogenesis, leptogenesis, and dark matter. Although beyond the scope of this paper, we think that it would be very interesting to pursue such avenues and we leave it for future work.
Note Added in Proof: When this paper was accepted for publication and in proof, Ref. [116] appeared on the ArXiv, revisiting the process of invisible neutrino decay in the early Universe. Even though the new analysis from Ref. [116] does not affect significantly our results and conclusions, a detailed discussion of the potential implications for our study can be found in Appendix A.7.
For a massive Z the decay rate is given by Case 2s: The corresponding two decay rates can be simply obtained by setting m νj → 0 in case 2a and performing the following replacement in the above equations.  Case 3a2: The ν i → ν 4 φ * → ν 4ν4 ν j decay rate is given by The expression for the ν i → ν 4 Z * → ν 4ν4 ν j decay rate is given by Equation (A.5) performing the following mapping Case 3s: The result for the decay and for the Z mediated case

A.2 Complementary bounds on invisible neutrino decay rates
In this Appendix we will briefly review the bounds on the invisible neutrino decay lifetime that can be extracted from astrophysics and laboratory experiments, which are complementary to the stronger constraints from cosmology presented in Section 3.1.
Astrophysical Bounds. The observation of neutrinos from the supernova SN1987A constrains the neutrino lifetime to be τ ν /m ν > 5.7 × 10 5 s eV −1 [117], provided all neutrinos decay into very light and non-interacting BSM species. Given the low number of observed neutrinos from SN1987A, there are no solid bounds from SN1987A on decay modes that involve only one neutrino decay or that count with active neutrinos in the final state [118,119]. However, given the sensitivity of current neutrino detectors, data from the next galactic supernova is expected to probe decay modes of the type ν i → ν j φ with rates τ ν 10 7 s eV −1 [118][119][120]. Similarly, neutrino decays can be searched for with neutrino telescopes [121][122][123]. A very recent analysis of IceCube data [124] shows that τ ν1, 2 > (1 − 3) × 10 −3 s eV −1 at 90% CL for IO. For NO no bound has been reported yet, but similar lifetimes are expected to be probed in the years to come [124]. Furthermore, neutrino telescopes are expected to reach sensitivities to neutrino decays up to τ ν /m ν 10 3 s eV −1 from astrophysical sources [121][122][123][124].
Even though laboratory and astrophysical bounds on the neutrino lifetime are orders of magnitude weaker than those derived from cosmological observations, the latter can be avoided if neutrinos decay today but not in the early Universe [141]. This can potentially occur if the neutrino masses and/or couplings with additional BSM species are time dependent. We consider that this possibility is not generic and we have not contemplated it in this work.
A.3 Allowed parameter space for individual ν i → ν j φ/Z decays In Figure 9 we show the parameter space for decays of the type ν i → ν j φ/Z when every decay mode is individually considered while the rest are assumed to be switched off. According to the discussion in Section 3.1, we estimate that the maximal relaxation of the cosmological mass bounds is given by the mass difference between the mother and daughter particle. Hence, it is obvious that, if the only opened decay channel is ν 2 → ν 1 φ/Z , cosmological neutrino mass bounds can not substantially be modified because ∆m 2 21 ≈ 7.4 × 10 −5 eV 2 , which leads to a maximal shift of roughly 10 −3 eV for both NO and IO. The situation is different for the decay modes ν 3 → ν 1 φ/Z and ν 3 → ν 2 φ/Z (ν 1 → ν 3 φ/Z and ν 2 → ν 3 φ/Z ) for NO (IO), since the relevant mass difference involved is |∆m 2 31 | ∼ |∆m 2 32 | = 2.5 × 10 −3 eV 2 ∆m 2 21 . Each of these four channels can potentially relax the present Planck neutrino mass bound up to roughly 0.14 eV. It is important to remark that the excluded grey region in the figures can be considered as a conservative bound since: i) we have considered the different decay modes individually; ii) a detailed analysis of all the potential invisible neutrino decay effects in cosmology would likely lead to a stronger constraint. The upper panels correspond to neutrino decays via pseudo-scalar couplings, the middle panels correspond to scalar couplings, and the lower panels to neutrino decays with a vector boson in the final state. For concreteness, we have considered m φ , m Z mν . We highlight in blue the bound on invisible neutrino decays from Planck CMB observations [22] (see Section 3.1). In dashed green we highlight the current bound on mν within ΛCDM [12], and show with an arrow the extent to which it can be relaxed. We shade in grey regions of parameter space that are be excluded by the current Planck bound m cosmo ν < 0.12 eV.
Here we explicitly revaluate the SN1987A bound from supernova cooling, but considering all possible production channels of φ's for m φ = 0. In order to set the constraint we consider a supernova core of temperature T SN = 30 MeV and radius R SN = 10 km. We also take into account that there is a large number of electrons neutrinos in the SN core [142] and choose an associated chemical potential to be µ νe = 200 MeV (while a negligible one for ν µ and ν τ ones). These values for the temperature, radius and neutrino chemical potentials are typically found in supernova simulations [76,142] and, in particular, are the same values considered in Ref. [82]. In order to set a constraint we require that the luminosity of emitted φ species does not exceed the typical binding energy of a neutron star E b = 3 × 10 53 erg. We take into account production viaνν → φ andνν → φφ and account for the fact that if the produced bosons decay or annihilate back into neutrinos within R < R SN there is no cooling. Explicitly, this means where V = 4π 3 R 3 SN , ∆t = 10 s, K 2 (x) is a Bessel function of the second kind, µ is the chemical potential of a given neutrino, and Γ φ→νν is the decay rate at rest of φ into neutrinos. The number density of the neutrinophilic boson can be approximated by its thermal equilibrium value, n φ 1.2 T 3 /π 2 .
The physical interpretation of Equation (A.21) is the following: the first line accounts for emission of φ's viaνν → φφ. It is the very same expression as for the energy density rate for this process, see e.g. Equation 2.15 of [143], but modulated by an exponential that accounts for the fact that the φ's will become trapped in the core if φφ →νν and φ →νν processes have a mean free path shorter than R SN . The second line in (A.21) represents the same as the first line but for production viaνν → φ inverse decays.
Finally, we obtain our constrain on the relevant coupling strengths by taking the annihilation and decay rates from Equations (A.20), (A. 16), and (A.17), using (A.21) and by requiring ∆E SN < E b 3 × 10 53 erg. For masses below m φ ∼ 100 MeV the SN1987A bound for electron neutrinos is given by while for µ and τ neutrinos it leads to In the right panel of Figure 10 we show the bounds from SN1987A for masses up to 1 GeV. The λ ij couplings (mass basis) are related to λ ee (flavor basis) via a rotation with the PMNS matrix and thus the constraint on λ ij is similar to the λ ee . As a result, in the main text, and in particular in the left panel of Figure 7, we choose to show the bound on the λ ee coupling. We can also use SN1987A data to constrain the coupling among the scalar field, the active neutrino and the sterile state. In this case, the main constraint arises due to processes such asν i ν j → φφ as mediated by a sterile neutrino. In this case, we assume λ 44 λ i4 and that φ decays toν 4 ν 4 . Given these processes, SN1987A excludes for m φ 10 MeV and m ν4 MeV , (A.24) Constraints on ν-φ interactions from supernova deleptonization have also been derived restricting λ 10 −4 [80]. These bounds are however subject to the details of the supernova explosion and, thus, have not being considered in the present paper. We also note that recently, bounds on ν-φ interactions have been derived based on neutrino shockwave considerations within a supernova: λ 10 −2 , but extending up to m φ/Z 1 GeV [144]. Big Bang Nucleosynthesis: We follow [143,145] and use NUDEC BSM to model the physics of neutrino decoupling in the presence of φ bosons and ν 4 neutrinos. We write down evolution equations for the temperature of the electromagnetic plasma, active neutrinos, ν 4 , and φ bosons including all relevant interactions among them. For 2 ↔ 2 processes we take m νi = m ν4 = m φ = 0 since these masses are negligible for temperatures that can impact N eff , namely T > T dec ν 2 MeV. We assume that there is no primordial population of φ and ν 4 particles for T 10 MeV and that for T 10 MeV these species are only produced via their interactions with active neutrinos. We note that the bounds we derive below would become more stringent had we considered a primordial population of such species or additional production channels.
We show the resulting values of N eff as a function of the relevant coupling constants in the left panel of Figure 10. The current BBN bound, N eff < 3.4 [146,147] at 95% CL, is indicated by the orange dashed line. We set conservative bounds on individual neutrino-boson interactions by assuming that the rest of the couplings vanish -otherwise constraints will be more stringent. By doing so we find that successful BBN requires λ ij < 8 × 10 −6 , λ i4 < 1.6 × 10 −5 , (BBN) (A. 25) where these couplings are defined in Lagrangian (2.1). Note that these bounds strengthen by a factor 4 √ 3 1.3 if three neutrinos simultaneously interact with φ. Figure 11. Left panel: N eff as a function of |θ4α| for eV-scale sterile neutrinos with ∆m 2 ≡ m 2 ν 4 − m 2 ν i > 0. Solid lines correspond to our calculation using NUDEC BSM and dashed lines correspond to the results of [102]. We appreciate the excellent agreement. We note that a similar agreement is found for |θ 4µ/τ |. Right panel: 95% CL BBN exclusion contour for |θ4α| for sterile neutrinos lighter than the mainly active neutrinos (∆m 2 < 0) as computed using NUDEC BSM.
In the limit in which m φ is not negligible, φ species can also be produced in the early Universe viaν i ν j → φ processes prior to neutrino decoupling and therefore modify the expansion history of the Universe. If the processν i ν j → φ thermalizes prior to T dec ν it would lead to N eff 3.62. Including this type of processes, we find that for N eff < 3.4 the interaction couplings are required to be λ ij < 1.2 × 10 −9 MeV/m φ , (BBN) (A. 26) which roughly corresponds to requiring Γ(ν i ν j → φ) < H(T = 4 MeV) and applies when 2 m ν m φ < 1 MeV.

A.6 BBN bounds on very light sterile neutrinos
In this appendix we place BBN bounds on the mixings of mostly sterile neutrinos with active neutrinos as relevant for invisible neutrino decays. There are very sophisticated studies dealing with cosmological constraints on eV-scale sterile neutrinos [102][103][104]. Previous studies, however, have focused on the case of sterile neutrinos that are heavier than active neutrinos, ∆m 2 ≡ m 2 ν4 − m 2 νi > 0. We, on the other hand, are interested in sterile neutrinos that are lighter than active neutrinos since they could represent the active neutrino decay final state, ∆m 2 < 0. The early Universe phenomenology of both ∆m 2 > 0 and ∆m 2 < 0 is fairly analogous but with a relevant difference. In our case of interest, ∆m 2 < 0, there is resonant sterile neutrino production via a MSW-like effect. This leads to significantly more stringent constraints on |θ α4 | for ∆m 2 < 0 than for ∆m 2 > 0.
In the early Universe, sterile neutrinos are produced via scatterings and oscillations of active neutrinos. The production of eV-scale sterile neutrinos in the thermal plasma roughly peaks at T ∼ 4 MeV ( |∆m 2 |/0.1 eV) 1/3 [148]. Active neutrinos decouple from electrons and positrons at T ∼ 2 MeV and therefore any sterile neutrinos produced at T 2 MeV will impact N eff as relevant for BBN and CMB observations. In this appendix we use NUDEC BSM [143,145] to model neutrino decoupling in the presence of an eV-scale sterile-neutrino. As in [143,145], we assume all relevant species can be described by thermal equilibrium distribution functions with an evolving temperature and negligible chemical potentials. We then solve for the temperatures of photons, active neutrinos and sterile neutrinos including all relevant interactions among them. We account for the active-sterile neutrino interactions by using the collision term described in Equation 1 of [148] and assuming there is no asymmetry between particles and antiparticles in the early Universe. We validate our method for ∆m 2 > 0 by comparing our results against state-of-the-art references in the left panel of Figure 11. There we can appreciate the excellent agreement between our calculation and that of [102] for the case ∆m 2 > 0. In the right panel of Figure 11 we show the bounds on |θ α4 | for ∆m 2 < 0. We consider the latest BBN bound N eff < 3.4 [146,147] at 95% CL and find a limit on |θ α4 | of We notice that this bound is roughly two orders of magnitude more stringent than for ∆m 2 = m 2 ν4 − m 2 νi > 0.

A.7 Implications of a relaxation of the CMB bounds on invisible neutrino decays
Ref. [116] represents a substantial improvement in our understanding of the cosmological implications of neutrino decays and is the first to write down the full inhomogeneous Boltzmann equation for neutrino decays in the early Universe. Armed with this machinery, the authors argue that the bound on relativistically decaying neutrinos from CMB observations considered in our analysis (see Section 3) τ ν > 1.3 × 10 9 s (m ν /0.05 eV) 3 [22] is relaxed by up to 3 orders of magnitude to τ ν > (0.4 − 4) × 10 6 s (m ν /0.05 eV) 5 [116]. In spite of this relaxation, using this new bound in our study does not significantly modify our results and conclusions. This is explicitly shown in Figure 12 where we highlight the impact of considering the bound from Ref. [116] in our analysis (dashed blue lines) to be compared with our results from [22] (solid blue lines). In Ref. [116] it is also argued that the derived bound on relativistically decaying neutrinos a priori only applies for neutrino masses m ν < 0.2 eV. However, for masses m ν > 0.2 eV, decaying neutrinos will still substantially reduce the neutrino anisotropic stress energy tensor prior to recombination. Indeed, for m ν ∼ 0.5 eV neutrinos become non-relativistic close to recombination which means that the actual rate that should suppress neutrino freestreaming is Γ ∼ τ −1 ν at T ∼ m ν /3. For T m ν /3 one needs to account for the suppression arising from freeze-out inverse decays. Even though the exact constraints in this region have not been studied, we expect that the region of parameter space in which neutrinos decay at a rate 1/τ ν > H(T = T CMB ) would be excluded by Planck data. However, for the purpose of our study we extrapolated the result of [22] that we therefore expect to be a conservative estimate.
Finally, in Ref. [116] the CMB anisotropies including these new details of the neutrino decay process in the early Universe were not calculated. Thus, it remains a task to confirm the expectations from [116] with a full analysis of the Planck CMB data. As we discuss in Section 3 and in the previous paragraph, such analysis will finally clarify what are the precise constraints on ultrarelativistically decaying neutrinos from CMB observations.  Figure 3 (upper panels) and Figure 6 (lower panels) but including the results from [116] (dashed blue lines), to be compared with the bounds presented in our analysis and derived from [22] (solid blue lines).