Bounds on light sterile neutrino mass and mixing from cosmology and laboratory searches

We provide a consistent framework to set limits on properties of light sterile neutrinos coupled to all three active neutrinos using a combination of the latest cosmological data and terrestrial measurements from oscillations, $\beta$-decay and neutrinoless double-$\beta$ decay ($0\nu\beta\beta$) experiments. We directly constrain the full $3+1$ active-sterile mixing matrix elements $|U_{\alpha4}|^2$, with $\alpha \in ( e,\mu ,\tau )$, and the mass-squared splitting $\Delta m^2_{41} \equiv m_4^2-m_1^2$. We find that results for a $3+1$ case differ from previously studied $1+1$ scenarios where the sterile is only coupled to one of the neutrinos, which is largely explained by parameter space volume effects. Limits on the mass splitting and the mixing matrix elements are currently dominated by the cosmological data sets. The exact results are slightly prior dependent, but we reliably find all matrix elements to be constrained below $|U_{\alpha4}|^2 \lesssim 10^{-3}$. Short-baseline neutrino oscillation hints in favor of eV-scale sterile neutrinos are in serious tension with these bounds, irrespective of prior assumptions. We also translate the bounds from the cosmological analysis into constraints on the parameters probed by laboratory searches, such as $m_\beta$ or $m_{\beta \beta}$, the effective mass parameters probed by $\beta$-decay and $0\nu\beta\beta$ searches, respectively. When allowing for mixing with a light sterile neutrino, cosmology leads to upper bounds of $m_\beta<0.09$ eV and $m_{\beta \beta}<0.07$ eV at 95\% C.L, more stringent than the limits from current laboratory experiments.


Introduction
The observation of oscillations between different neutrino flavors firmly establishes that at least two out of three neutrino mass eigenstates m i (i = 1, 2, 3) are non-zero [1][2][3][4][5]. Since in the Standard Model of particle physics neutrinos only come with left-handed chirality, it is not possible to generate a mass term for them by the usual Higgs mechanism as for the other leptons. The discovery of neutrino masses then requires new physics beyond the Standard Model, which may involve the existence of additional right-handed states (see e.g., [6] for a review). These states are sterile, which means that they do not take part in weak, strong, or electromagnetic interactions, and their number and mass scale depend on the specific model. For phenomenological purposes, in this paper we consider only one light sterile state with a mass at the eV scale [7].
One way to search for sterile neutrinos is in laboratory experiments. If produced in sufficient numbers in the early universe via oscillations with active states, sterile neutrinos can lead to observable effects in cosmology. As we will discuss in this paper, direct laboratory measurements and cosmological observations are complementary, and both avenues have led to tremendous progress in our understanding of the neutrino sector over the past few years.
Neutrino oscillation experiments have provided high-precision measurements of two mass-squared splittings [8][9][10]: These mass splittings are responsible for the oscillations between the three neutrino flavor eigenstates associated with the charged leptons of the Standard Model. Given the fact that the sign of ∆m 2 31 is not determined, two mass orderings can be possibly realized in nature: the normal ordering, in which the lightest mass state is m 1 and ∆m 2 31 > 0, or the inverted ordering, in which m 3 is the lightest state and ∆m 2 31 < 0. Current neutrino oscillation measurements show a mild statistical preference for the normal ordering [8][9][10][11].
For many years, however, a number of experimental results that cannot be explained within the context of the three-neutrino oscillation paradigm have been reported. Such anomalies are measured by various so-called short-baseline (SBL) experiments. The first result was found by LSND [12] and later partly confirmed by MiniBooNE [13], but anomalies were also detected by the Gallium experiments GALLEX and SAGE 1 [14][15][16] and by a number of observations of reactor antineutrino fluxes at short distances [17][18][19]. Although the situation is not completely clear, (see e.g. [7,20]), the anomalies point to an additional mass splitting, much larger than the other two: The existence of a new mass splitting implies the presence of a new neutrino mass eigenstate, which corresponds to a fourth flavor eigenstate. Since the hypothetical new neutrino has not been found in other weak interaction measurements [21], it cannot take part in Standard Model interactions and is therefore denoted "sterile" (see e.g. [7,22]). Due to the large gap O(eV 2 ) |∆m 2 31 |, ∆m 2 21 , oscillations with the sterile are not expected to affect standard atmospheric and solar neutrino measurements. One typically assumes that the mixing between the three active neutrinos and the fourth mass eigenstate is small, so that the scenario is commonly referred to as the 3+1 model. For recent reviews on experimental searches for eV-scale sterile neutrinos, see also [7,20,23].
In case such a new neutrino exists, its presence can affect the evolution of the Universe as well. Even if the fourth neutrino does not interact with the other standard model particles, it would be produced in the early Universe via oscillations with the active flavors. Excluding direct detection of relic neutrinos [24], the cosmological presence of the sterile state can only be deduced indirectly, for instance from its contribution to the energy density of the Universe [25]. In particular, considering masses around the ∼ eV scale as motivated by the SBL anomalies, the sterile neutrino is produced as a relativistic particle. Therefore, it contributes at early times to the number of relativistic degrees of freedom N eff , which quantifies the amount of energy density from relativistic species beyond photons in the early Universe. Moreover, the free-streaming length associated to light sterile neutrinos is potentially large, so that they would behave as hot dark matter and leave a distinct imprint on structure formation, similar to the one of active neutrinos [26][27][28][29]. Cosmological observations are sensitive to the hot dark matter energy density Ω hot = (Ω ν + Ω s ), neglecting the contribution of other nonstandard particles. Recent measurements of the temperature and polarisation anisotropies in the cosmic microwave background (CMB) and of large-scale structures (LSS) are able to constrain both N eff and Ω hot h 2 to high precision, see e.g. [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44]. Predictions of big-bang nucleosynthesis (BBN) combined with measurements of primordial element abundances also provide independent, albeit less stringent, constraints on N eff in agreement with CMB+LSS findings [45].
There is another reason to consider additional relativistic particles such as light sterile neutrinos in cosmology, independent of the anomalous results from neutrino experiments. Contributions to the energy density in the early Universe reduce the size of the sound horizon and result in a higher value of the Hubble constant H 0 inferred from baryon acoustic oscillations (see e.g. discussions in [46][47][48][49][50]). Currently measurements based on observing the sound horizon in the CMB or in the galaxy distribution at late times are in tension with local distance ladder measurements from supernovae Ia [51][52][53] and lensed quasar time delay measurements [54][55][56], see e.g. [57][58][59] for further discussions. However, a sterile neutrino alone cannot account for the full discrepancy, as discussed in [41].
In this work, we study the phenomenology of a light sterile neutrino from various points of view. While the main goal is to compute, for the first time, cosmological bounds on all the mixing angles and the mass splitting between the light sterile state and active neutrinos, we also compare the cosmological results with other probes, such as constraints from direct mass measurements and neutrinoless double-β decay, and limits from neutrino oscillation experiments. We will focus on light sterile neutrinos (i.e., ∆m 2 41 < 100 eV 2 ), since we are motivated by anomalies in oscillation experiments that could be explained by the existence of light sterile states.
The rest of this paper is structured as follows. First, we give a brief overview on the theoretical aspects and the status of sterile neutrino probes in section 2. The phenomenology of such neutrinos in the early and late Universe is discussed in section 3. In section 4 we discuss the datasets we consider and the method used to obtain the bounds on the sterile neutrino mixing parameters, for which we present upper limits in section 5. We summarize and conclude in section 6.

Experimental searches for light sterile neutrinos
In this section we briefly summarize the theory of active-sterile neutrino oscillations and review the status of light sterile neutrino searches at various experiments. For comprehensive reviews focused on experimental searches for eV-scale sterile neutrinos, see [7,20,23].

Active-sterile neutrino mixing
The mixing between active and sterile neutrinos can be parameterized by means of the unitary 4 × 4 lepton mixing matrix U . There are several ways to write the mixing matrix, so we will briefly explain our choice of parameterization. In the case of mixing between four neutrino states, U is fully characterized by six mixing angles and three CP-violating Dirac phases, but for the purposes of this work we set the phases to zero and therefore assume that oscillations among neutrinos are identical to those among antineutrinos. Three of the angles characterize the oscillations between the three active neutrinos: θ 13 , θ 23 , and θ 23 . The remaining three mixing angles, namely, θ 14 , θ 24 , and θ 34 , describe the mixing with the sterile state. We choose the following parameterization of the mixing matrix [22]: where each R ij is a real matrix describing a rotation by an angle θ ij . Bounds on the mixing matrix can be provided in a way that is independent of the parametrization if instead of quoting limits on the individual mixing angles θ ij we consider the matrix elements directly. Concerning active-sterile neutrino mixing, the most important elements are those of the fourth column:

4)
|U s4 | 2 = cos 2 θ 14 cos 2 θ 24 cos 2 θ 34 . (2.5) We expect the mixing angles θ i4 to be small, in order not to substantially alter the phenomenology of three-neutrino mixing beyond what is allowed by current limits. Therefore, the matrix elements |U e4 | 2 , |U µ4 | 2 and |U τ 4 | 2 are expected to be small, while |U s4 | 2 should in principle be of order unity. For this reason, it is possible to refer to the fourth neutrino mass eigenstate as the "sterile" neutrino, even though strictly speaking a sterile neutrino is a flavor eigenstate and has no definite mass. In the rest of the paper, we will often refer to the mass m s m 4 as the mass of the sterile neutrino.
In addition to the mixing angles, we need to specify the three mass-squared splittings ∆m 2 21 , ∆m 2 31 , and ∆m 2 41 . For the active neutrinos, we always assume normal ordering 2 (∆m 2 31 > 0) and we fix all the standard mixing parameters to the best-fit values obtained in [8]. When considering the full oscillation pattern, the complete expressions for the oscillation probabilities are rather involved and depend on all the mass splittings and mixing angles, see e.g. [65]. When considering SBL oscillations, however, the terms due to the much smaller solar and atmospheric mass-squared differences are suppressed because they correspond to slower oscillations, and only the effect of ∆m 2 41 is relevant. As a consequence, the oscillation probabilities can be well approximated by a two-neutrino mixing formula with appropriate mixing matrix elements. This is the appearance probability to populate a flavor state: whereas the flux of the initial flavor is modulated by the disappearance probability: where L is the distance traveled by the neutrino, E is its energy, and we use Greek letters to refer to the four flavor states α, β ∈ [e, µ, τ, s]. The effective angles ϑ αα and ϑ αβ depend on the elements of the fourth column of the mixing matrix: The last expression is implied by the unitarity of the mixing matrix. From the probabilities in Eqs. (2.6) and (2.7) it is clear why the oscillatory behavior manifests when 4E ∼ ∆m 2 L. If 4E ∆m 2 L, the oscillatory term vanishes; if 4E ∆m 2 L, the oscillation frequency becomes so high that oscillations cannot be resolved anymore and they are averaged out. The eV scale we focus on arises since various experiments find anomalies pointing towards (L/E) −1 ∼ O(eV 2 ). We emphasize that the absolute mass scale of the neutrino states does not enter in the equations, so oscillation experiments depend on the differences of the masssquared values alone. In other words, oscillations are independent of the lightest neutrino mass m 1 . Both appearance and disappearance channels can be used to measure the effective mixing angles and to constrain the mixing matrix elements.

Current status of oscillation data and global fits
There is not a single SBL anomaly, but several different experiments with short baselines find anomalous neutrino fluxes of varying degree. We will briefly summarize the current situation in this section.
From the historical point of view, the first anomaly was found by LSND [12], which reported the unexpected appearance of electron antineutrinos in a beam of muon antineutrinos produced from π + decays. Years later, also MiniBooNE [13] confirmed the excess ofν e events, with a similar experimental setup, in partial agreement with LSND even though the MiniBooNE excess is too large to be explained by a sterile alone.
Another anomaly with comparable L/E was found by ν e disappearance measurements by the Gallium neutrino detectors GALLEX and SAGE [14,15,66,67]. Both experiments measured the electron neutrino flux in proximity to a radioactive source, and found a lack of events at significance of ∼ 3σ. A similar effect is observed in measurements ofν e disappearance in close proximity to nuclear reactors [17]. The lack of electron antineutrinos also reaches a significance of ∼ 3σ and was noticed after new calculations found a higher expected initial flux [18,19].
The LSND and MiniBooNEν e appearance data require a non-zero value of ϑ eµ , while in order to explain the reactor and Gallium disappearance measurements one needs ϑ ee > 0. Taken together, this also implies a non-zero mixing ϑ µµ and |U µ4 | 2 . However, this matrix element can be measured independently by muon (anti) neutrino disappearance experiments, and no corresponding anomaly for the relevant L/E values is found by either measurements of the atmospheric muon neutrino flux by IceCube [68,69] or by the MINOS+ [70] collaboration using an accelerator beam. Therefore, a combination ofν e appearance data by LSND and MiniBooNE and the disappearance results from electron and muon (anti) neutrinos in a global fit is currently problematic [71,72].

Laboratory searches for the absolute neutrino mass scale
In addition to neutrino flavor oscillation experiments, laboratory searches for massive neutrinos also include kinematic measurements of β-decay and searches for neutrinoless double-β (0νββ) decay events. In this section, we briefly summarize the consequences that a light sterile neutrino mixing with the active flavors can have for these searches.
One way to probe the absolute neutrino mass scale directly is to observe the cutoff of the electron energy spectrum emitted from β-decay. It occurs at the effective electron neutrino mass m β , given by the incoherent sum of squared masses and mixing parameters: where j = 4 is the contribution from the additional sterile neutrino. So far, β-decay measurements have only been able to set upper limits on the mass scale, with the latest bound m β < 1.1 eV at 90% confidence level (C.L.) recently published by the KATRIN collaboration [73]. Neutrinoless double-β decay (0νββ) searches constrain the half-life T 1/2 of the isotope involved in the decay (see e.g. [74] for a recent review). Assuming that the mechanism responsible for lepton number violation manifested in 0νββ events is the mass mechanism, constraints on the half life can be translated into constraints on the effective Majorana mass m ββ : where m e is the electron mass, G 0ν is the phase space factor and M 0ν is the nuclear transition matrix element for the decay. The effective mass parameter m ββ can be expressed as a coherent sum of mass eigenstates and mixing matrix parameters: where α j are Majorana phases and j = 4 again contains the contribution from the sterile neutrino. So far, no such event has been detected and only upper limits on the lifetime T 1/2 of various isotopes are available, as described in section 4.3. These bounds are usually converted into a range of upper limits on the Majorana mass m ββ depending on theoretical uncertainty in the calculation of the nuclear matrix elements.

Light sterile neutrinos in cosmology
As explained in section 1, sterile neutrinos do not take part in weak, strong or electromagnetic interactions. As a consequence, they will not be produced by Standard Model scattering or annihilations in the very early Universe. In this work, we consider a production mechanism via non-resonant oscillations with active states, the so-called Dodelson-Widrow production mechanism [75,76]. In the very early universe, densities are high and weak interactions frequent. This generates an effective matter potential that suppresses neutrino oscillations. Therefore the sterile state is only populated once densities are low enough for oscillations with the active neutrino eigenstates to occur (see e.g. [77] for a detailed discussion). The mass splitting sets the timescale for oscillations with the fourth neutrino and determines the time when flavor oscillations can start to arise. Here we focus on eV-scale sterile neutrinos and consequently consider mass splittings 10 −2 ≤ ∆m 2 41 /eV 2 ≤ 10 2 . This range corresponds to plasma temperatures between O(100) and O(1) MeV when the sterile starts to be populated.
The thermalization process is described by a set of differential equations that encode the evolution of the momentum distribution functions f α (p) of the various neutrino flavors α ∈ [e, µ, τ, s], and an additional one for the evolution of the photon temperature, as described in detail in [78]. As far as cosmological effects of neutrinos are concerned, what matters is the momentum distribution function of all the neutrino states after their decoupling from the thermal plasma and at the end of electron-positron annihilations into photons, which happens at temperatures around 0.1 MeV.
The final momentum distribution function for the active neutrinos is very close, but not exactly equal, to a Fermi-Dirac shape with the neutrino temperature T ν (see e.g. [78][79][80][81][82]).
The deviation from thermal equilibrium is mostly due to the fact that neutrino decoupling does not occur instantaneously, so there are small distortions at high momenta that come from the energy transferred during electron-positron annihilations to the few neutrinos still coupled to the plasma. For the sterile neutrino, the momentum distribution function f s depends on the degree of thermalization that it reaches. Initially the sterile is absent, and if the mass splitting and the mixing angles are not large enough, oscillations either start too late or are not efficient enough to bring the fourth neutrino into equilibrium with the active flavors. The degree of thermalization of the sterile can be expressed in terms of the effective number of relativistic species (N eff ), which can be constrained by big-bang nucleosynthesis (BBN) [45] and more tightly by CMB observations [83]. After electron-positron annihilations, this parameter can be expressed as with the energy density in photons ρ γ and active plus sterile neutrinos (ρ ν + ρ s ). They can be computed from the integrated distribution functions, and for negligible contributions from the sterile we recover the standard value N 3ν eff = 3.043−3.045 [78,[80][81][82]84]. On the other hand, if the mixing parameters are sufficiently large, there is time for neutrino oscillations to fully bring the fourth neutrino to equilibrium with the active flavors. In such case, the final distribution function of sterile neutrinos f s will also be very close to a Fermi-Dirac spectrum, with the same temperature as the active neutrinos, and N eff 4.05. Intermediate cases correspond to an incomplete thermalization, the sterile neutrino contributes with ∆N eff = N eff −N 3ν eff < 1.01 and its distribution function is significantly non-thermal. In such cases, however, it turns out that since the sterile is populated through oscillations from the thermal active states, the distribution function f s is still proportional to a Fermi-Dirac shape [78] with temperature T ν . Then, we generally have [75,76] where p is the neutrino momentum and T ν is the temperature of active neutrinos. Note that for all cases considered here, the sterile neutrino distribution function reaches its asymptotic values well before primordial BBN at T ∼ 0.1 MeV [78]. This means that, under our assumptions, the value of N eff does not change between the time of BBN and CMB decoupling, which would have consequences for the abundance of light elements [45] important for the calculation of cosmological perturbations (see also [85]). The problem of active-sterile oscillations in the early Universe has been studied in many previous papers, some of them published more than 30 years ago (early references include e.g. [86][87][88][89][90], see the review [79] for an extensive list). Solving the Boltzmann kinetic equations for different neutrino energies is a complex issue, due to the simultaneous presence of neutrino interactions via weak processes and flavor oscillations in an expanding medium. Thus, past analyses [91][92][93][94][95][96][97][98][99][100][101][102][103] have considered various approaches that approximated the multi-momentum calculations to the evolution of an average momentum and/or reduced the number of active neutrino states. In particular, the LASAGNA code [95] solves the quantum kinetic equations in the 1+1 case (assuming one active coupled via oscillations to one sterile neutrino state) with full collision integrals. This code has been used in previous works [94,[100][101][102] to map the active-sterile mixing parameters onto two other quantities relevant for cosmology (N eff and the effective sterile neutrino mass).
When considering a simplified model with only one active neutrino ν a and one sterile ν s , it is possible to relate the degree of thermalization ∆N eff directly to the mixing parameters. For a mass splitting δm 2 as and mixing angle ϑ as , this results in [92] δm 2 where the numerical coefficient is slightly different for electron, muon, or tau flavor neutrinos. If one applies this relation to the 3+1 case, with δm 2 as → ∆m 2 41 and sin 2 (2ϑ as ) 4|U α4 | 2 |U s4 | 2 , one gets that in order to have a fully thermalized sterile neutrino with a mass splitting around 1 eV 2 , a mixing matrix element |U α4 | 2 10 −3 is required. From this relation we can also see that a larger mixing matrix element generally increases N eff towards 4. For larger mass splittings, a smaller mixing is sufficient to generate the same level of thermalization since the oscillations are faster, see Eq. (2.6)-(2.7). While Eq. (3.3) is a rough estimate, we will see that it is a quite good approximation of the full calculation, as long as one mixing angle is varied at each time.
Including also the mixing among the active neutrino states has also been considered before, from the early simplified analyses [92,93] to more recent multi-angle studies [97,98,103] performed within the averaged-momentum approximation. Although some authors have studied the multi-angle (2+1 scenario) and multi-momentum problem [99], only very recently the evolution in the early Universe of the momentum-dependent kinetic equations for the full 4 × 4 density matrix of neutrinos was calculated with a dedicated numerical code, FortEPiaNO 3 , as described in [78].
At late times, the sterile neutrino becomes non-relativistic and it contributes to the matter energy density. At such point in the evolution, its contribution to the energy density is [76,77] where we have introduced the effective mass m s, eff ≡ m s ∆N eff . As mentioned in section 1, light sterile neutrinos might behave as a hot dark matter component and affect the evolution of matter perturbations in a similar way as active neutrinos [25]. The form of the distribution function (3.2) implies that the sterile neutrinos have the same average momentum as the active species. Consequently, the maximum free-streaming length of the sterile is equal to the one of an active neutrino with mass m ν = m s , corresponding to a comoving wavenumber k fs [77]: if neutrinos become nonrelativistic during the matter-dominated era, or if they become non-relativistic during the radiation-dominated era. Here m can be either the mass of an active neutrino, or the mass of sterile neutrino produced through non-resonant oscillations. Their large velocities prevent neutrinos from clustering at scales smaller than the free-streaming length, so the collective effect of active and sterile neutrinos is a step-like suppression of the amplitude of matter perturbations below the free-streaming scale [76,77].
where P ν and P refer to the matter power spectrum with or without neutrinos respectively. The size of the suppression at small scales depends on the fraction of matter density provided by neutrinos, f hot = (Ω ν + Ω s )/Ω m . For sterile neutrinos the onset of the suppression is determined by the free-streaming scale and therefore m s , while the size of the suppression is proportional to m s, eff = m s ∆N eff . The cosmological effects of the sterile neutrino are therefore completely determined once m s (or equivalently Ω s ) and ∆N eff are specified.
In figure 1 we show the predicted contribution of the sterile neutrino to N eff (upper row) and the fraction of total dark matter energy density it represents (lower row, assuming a total Ω c h 2 = 0.119 from the latest Planck measurements [41]), for two different choices for the mixing matrix elements. From the top row of Fig. 1, it is clear that if at least one of the mixing matrix elements is much larger than 10 −3 , one has ∆N eff 1 for m s ∼ 1 eV in agreement with the expectation from Eq. (3.3) and the results of previous analyses of activesterile neutrino oscillations in the early Universe (see e.g. [78,94]). Such a value of ∆N eff is at odds with what is inferred from Planck measurements of CMB anisotropies, or from astrophysical determinations of the abundances of light elements. In general, we note that in some parts of the parameter space probed by our analysis, corresponding to large mixing angles and/or large masses, the contribution of sterile neutrinos to both N eff and the dark matter energy density is large. The latter is also a problem, since cosmological observations also constrain the fraction of hot dark matter to be small. As we shall see in the next section, these regions of parameter space will be indeed excluded by data, in agreement with our expectations.
As we will discuss in detail in section 5, large contributions to N eff are in tension with cosmological constraints, so in order to allow a sterile neutrino to form all the dark matter one needs to move to higher masses m s ∼ keV and very small mixing angles. In this mass range, a sterile neutrino produced by oscillations would behave as warm (or even cold) dark matter, so limits on the hot component would not apply. Then even modest contributions to N eff are sufficient to provide the required dark matter density, although other astrophysical limits apply, as reviewed for instance in [25,104,105]. The analysis of the mass range ∆m 2 41 > 10 2 eV 2 , however, is beyond the scope of this paper, as we are interested in light sterile neutrinos motivated by anomalies in neutrino oscillation experiments. The main question we want to address here is whether such a sterile state is compatible with cosmological bounds.

Datasets and methods
In this section, we present the data sets employed in this analysis, and we discuss the method adopted in this work, including the set up adopted for the FortEPianO code. Since we perform our analysis in a Bayesian framework, we also report the prior choices for the parameters varied in the analysis.

Cosmological data
We consider measurements of the CMB anisotropies in temperature and polarization [106], as well as determinations of the CMB lensing power spectrum [83], as provided by the latest Planck 2018 data release [41,107]. The data sets are incorporated by using the publicly available clik likelihood code described in [106].
To the CMB dataset, we add late-time distance and expansion rate baryon acoustic oscillation (BAO) measurements from the 6dFGS [108], SDSS-MGS [109], and BOSS DR12 [110]  surveys. These measurements are fully consistent with the Planck results and greatly help to break parameter degeneracies.

Oscillation experiments
As described in section 2.2, various oscillation experiments find conflicting results for the active-sterile mixing parameters and it is not possible to combine the measurements in a selfconsistent way. Therefore we group the non-conflicting oscillation measurements and present their results separately.
Theν e flux from nuclear reactors is measured by Bugey-3 [111], DANSS [112,113] NEOS [114], and PROSPECT [115]. Note that these experiments use differential measurements at a varied distance from the source, so a calibration of the intrinsic neutrino flux is not necessary. A combined frequentist analysis results in a preference for non-zero mixing with a sterile at about ∼ 2.5σ with two nearly degenerate best-fit points, located at (∆m 2 41 0.4 eV 2 , |U e4 | 2 0.01) and (∆m 2 41 1.3 eV 2 , |U e4 | 2 0.01) [72] (see also [71,116,117] for previous analyses).
We also include bounds on the mixing with muon neutrinos from measurements of the atmospheric flux by IceCube [68,69] and from the MINOS+ experiment [70], which uses an accelerator beam with two detectors, one close (∼ 500 m) and one far (∼ 800 km) from the source, to put strong constraints covering a wide range of ∆m 2 41 and |U µ4 | 2 . Both are combined in a frequentist analysis [72] and we refer to the combination of both data sets  [64,123]. The upper limits on m ββ depend on the uncertainty of the nuclear matrix element and are not used directly in our analysis. Details on the individual limits and the effect of the nuclear matrix elements are summarized in appendix A.
as "ν µ disappearance". Neither IceCube nor MINOS+ find anomalous events and therefore provide upper bounds on the mixing parameters. IceCube also has a limited sensitivity on |U τ 4 | 2 thanks to the low-energy data from DeepCore [68].

Decay experiments
In this work, we consider the latest measurements of the tritium β decay spectrum released by the KATRIN collaboration [73] that sets a limit m β < 1.1 eV at 90% C.L. We include the constraint by using the approximated analytical likelihood function presented in Eq. (B.3) of Ref. [118]. Since the modification to the decay energy spectrum induced by an extra light sterile species is below the resolution of KATRIN, this likelihood is also valid for the additional sterile neutrino. Concerning 0νββ searches, we consider the results from the KamLAND-ZEN collaboration phase-I [119], and phase-II [120], from the GERDA collaboration [121], and from the EXO collaboration [122]. The constraints are included by using the approximated likelihoods given in [123]. A summary of all data sets used, the isotopes employed by each collaboration and the respective bounds on the lifetime T 1/2 can be found in table 1, and we show the limits on m ββ from the individual experiments in appendix A.
While this paper was in preparation, new results from 0νββ collaborations have been published from the GERDA collaboration [124], the EXO200 collaboration [125], the CUORE collaboration [126]. We note that the tightest constraints on m ββ are still those by KamLAND-ZEN phase-II. As a result, the inclusion of other released data from 0νββ searches would not change the conclusions of this work.

Method
In order to use cosmological data to constrain the sterile neutrino properties, we have to map the fundamental parameters ∆m 2 41 , |U e4 | 2 , |U µ4 | 2 , |U τ 4 | 2 onto N eff and m s , which determine its impact on cosmological observables. As described in section 3 we use FortEPiaNO to calculate the evolution of the sterile distribution function including mixing with all three active neutrinos. We pre-compute a N eff table on a four-dimensional grid spanning the parameter space we are interested in. For the mass splitting, we consider nine logarithmically spaced samples over the range ∆m 2 41 ∈ [10 −2 , 10 2 ] eV 2 , while for each mixing matrix element we take eleven logarithmically spaced samples in the range |U α4 | 2 ∈ [10 −6 , 10 − six parameters 4 , complemented with the lightest neutrino mass m 1 , the sterile mass splitting ∆m 2 41 and the three mixing matrix elements |U α4 | 2 , α = e, µ, τ . We calculate the evolution of cosmological perturbations using the numerical Einstein-Boltzmann solver CLASS [127,128] 5 . Since we map the four active-sterile mixing parameters into the two cosmological parameters N eff and m s , a very large number of samples are necessary to obtain a well-converged posterior in the full sterile parameter space. To speed the calculations up, we therefore proceed in two steps. First we use MontePython [129,130] to run a standard Markov Chain Monte Carlo (MCMC) with CMB+BAO data. For this run, we vary all ΛCDM parameters and include three massless neutrinos with N eff = 3.045 and one additional massive state with varying mass and ∆N eff . Note that with the current sensitivities the data cannot distinguish between various combinations of neutrino parameters as long as they result in the same ∆N eff and Σm ν , so distributing the mass in different ways over the four available states does not affect the resulting constraints in a significant way (this reflects the fact that cosmology is not yet able to constrain the mass hierarchy, see e.g. [33,34,60,62,[131][132][133]). Given the prior N eff > 3.045, the cosmological data sets a 95%-limit of Σm ν < 0.15 eV and N eff < 3.44. We then use a Gaussian kernel density estimate (KDE) to interpolate the marginalized N eff − Σm ν likelihood surface. Evaluating this estimated likelihood is ∼ 10 5 faster than running the full CMB likelihood, and we explicitly make sure that sampling from it results in the same constraints as obtained from the MCMC.
In a second step, we use emcee [134] together with the estimated likelihood to sample the neutrino parameter space as follows: we vary the lightest neutrino mass m 1 , the mass splitting ∆m 2 41 and the three mixing matrix elements |U α4 | 2 . From the last four, we interpolate ∆N eff from the pre-computed FortEPiaNO grid. Since the mass splittings ∆m 2 21 and ∆m 2 31 are much smaller than the sensitivity of current constraints on neutrino masses, we assume three degenerate active neutrinos m 1 = m 2 = m 3 so we get a total Σm ν = 3m 1 +∆N eff m 4 . We then evaluate the KDE of the cosmological likelihood for N eff = 3.045 + ∆N eff and Σm ν . In order to compare the constraints to direct experiments, we also vary the nuclear matrix elements for 136 Xe and 76 Ge defined in Eq. (2.11) and the three Majorana phases α i in Eq. (2.12) to derive limits on m β and m ββ . A summary over all non-ΛCDM parameters varied in the analysis, their prior bounds and shapes is given in table 2.
For constraints involving the β-decay and 0νββ likelihoods described in section 4.3, we marginalize over all parameters listed in table 2 as well. The constraints from reactors and ν µ appearance measurements from MINOS+ and IceCube (see section 4.2), on the other hand, are derived in a frequentist framework and do not depend on the prior choices made here.

Results
In this section we present the results of our analysis. First we discuss the complementarity of the various data sets either in terms of the fundamental parameters of the sterile mixing matrix, or in terms of the derived parameters m β and m ββ directly constrained by terrestrial experiments. Then we explore the effect of assuming other priors on the sterile mass parameter for the cosmological datasets. In the same context we compare our constraints with the ones previously obtained in a simplified scenario where the sterile is only coupled to one neutrino species.

Constraints from cosmology and direct experiments
Cosmology provides tight limits on the sum of neutrino masses Σm ν and the effective number of relativistic species N eff , that can be translated to tight bounds on the parameter space of the mixing matrix with significant relativistic energy contributions from the sterile neutrino. After marginalising over all other cosmological and nuisance parameters, we obtain the constraints on the mixing matrix elements in figure 2 from CMB+BAO. As any mixing reaching a value of |U α4 | 2 ≈ 10 −3 starts to populate the sterile state and leads to detectable N eff contributions as seen in figure 1, this is where cosmological data becomes constraining. There are only small differences between mixing for the various flavors, so the limits on all of the matrix elements are very similar. Since the constraints on both the mass m 4 derived from the total sum Σm ν = 3m 1 + ∆N eff m 4 and N eff vanish for small amplitudes of the sterile distribution function set by ∆N eff , the cosmological data in principle allow large sterile masses as long as the mixing is small and the sterile state is not thermalized or populated to a significant level. As mentioned in section 3, much higher mass ranges in the keV range and larger are in principle possible if the mixing angles are small enough. In such a case the free-streaming length k fs of the sterile neutrino would be pushed to smaller scales, and it would form a warm or even a cold dark matter component [25,104,105].
In figure 3 (left panel), we compare the cosmological results with limits obtained from neutrinoless double-β decay and tritium decay measurements by KATRIN [73] in the ∆m 2 41 − |U e4 | 2 plane, since the decay experiments are not sensitive to the other matrix elements. While sensitivity of the neutrinoless double-β decay experiments and β-decay measurements from KATRIN on the sterile parameters are comparable, cosmological bounds are orders of magnitude stronger. To address the sterile neutrino interpretation of the reactor anomaly, we also show the parameter space preferred by a combined fit of short-baseline measurements of the reactor antineutrino flux. As explained in section 1, such experiments observe a preference in favor of a sterile with ∆m 2 41 ∼ O(eV 2 ) and |U e4 | 2 ≈ 10 −2 . While these parameter values are compatible with current measurements from 0νββ and KATRIN, they are in severe tension with the CMB+BAO data. A mixing of the size needed to explain the reactor data would be more than sufficient to bring the sterile in thermal equilibrium in the early universe. On the right of figure 3 we show a similar comparison for the mixing angle |U µ4 | 2 and present cosmological constraints together with bounds from ν µ disappearance measurements from IceCube and MINOS+ described in section 4.2. While the experimental sensitivity is stronger than on |U e4 | 2 , we still find the CMB+BAO data set to be more constraining. While the cosmological limits rely on model assumptions and can be slightly relaxed in extended parameter spaces, accommodating N eff ≈ 4 is very challenging [135,136].
We also map the cosmological bounds onto the parameter space m β or m ββ directly probed by decay experiments in figure 4. On the left side we show limits from CMB+BAO together with the latest results from KATRIN [73]. Since m β defined in Eq. (2.10) receives contributions from the sterile, the resulting limit from cosmology in the extended 3 + 1 parameter space is slightly higher compared to the expectation m β ≈ m 1 ≈ Σm ν /3 from standard neutrinos. The CMB+BAO data together with the prior assumption N eff > 3.045 Parameter experimental upper limit (95%) cosmological upper limit (95%) P(log ∆m 2 14 ) P(m 4 ) P(∆m 2 14 ) m 4 [eV] -1.6 4.4 6.8 0.08 (0νββ) 0.07 Table 3. Upper limits (95%) for the sterile neutrino mass, the parameters of the mixing matrix, and m β and m ββ . For the limits from direct experiments we take a conservative approach and only consider probes that are not in tension with cosmology. For those, we quote the strongest bound on each parameter. We present cosmological limits for the different prior choices for the mass splitting either on log ∆m 2 41 used in section 5.1, or on either m 4 or ∆m 2 41 discussed in section 5.2, but note that the prior choice barely affects the resulting constraints on m β and m ββ .  Marginalized 68% and 95% constraints on the mass splitting ∆m 2 41 and mixing matrix element |U e4 | 2 from cosmology (blue), from the tritium β-decay measurements by KATRIN (green) and from neutrinoless double-β-decay experiments (0νββ, red), compared with the preferred frequentist regions by the joint fit [72] of reactor experiments discussed in section 4.2, which are in strong tension with cosmological bounds. Right: Cosmological 68% and 95% marginalized constraints on the mixing matrix element |U µ4 | 2 compared to (frequentist) ν µ disappearance results [72] from IceCube and MINOS+ (grey).
lead to an upper bound on the neutrino mass sum of Σm ν < 0.15 eV, corresponding to a limit m cosmo β < 0.05 eV assuming only three active neutrinos, which is slightly degraded to m cosmo β < 0.09 eV if the sterile is included. However, this constraint is still a factor of ∼ 10 stronger than current experimental limits from KATRIN.
On the right-hand side of figure 4 we show a similar comparison of the derived bound on m ββ from cosmology and direct 0νββ searches. Note that the posterior from the combined 0νββ data has a maximum at slightly non-zero values due to an excess of events observed by EXO200 compared to the background expectation [125]. In table 3, we present a summary of all 95% bounds on the sterile mass m 4 , the mixing matrix elements |U α4 | 2 for each flavor, m β and m ββ , comparing the cosmological bound to the respective strongest bound from direct searches. Cosmology provides the tightest limits on all parameters, for most of them by orders of magnitude.

Priors and parameter space volume effects
Since all cosmological bounds presented in figure 2 and discussed in the previous section 5. 1 provide only upper bounds on the sterile parameters in a Bayesian framework, the choice of priors affects the limit. In this section we focus on various priors for the mass splitting ∆m 2 41 . While we adopt a logarithmic prior as a standard case, since cosmological data is ignorant of the order of magnitude of the mass splitting, one can also argue that the parameter of interest is either the mass-squared difference itself -or the sterile mass scale, since cosmological data is approximately sensitive to the energy density parameter Ω s ∝ m s ≈ m 4 . We therefore repeat the previous cosmological analysis using a flat prior either on log 10 ∆m 2 41 , on ∆m 2 41 , or on m 4 , always considering the same range specified before in table 2. Note that while sampling over m 4 we enforce the additional constraint ∆m 2 41 > 10 −2 eV 2 to make sure that the neutrino mass hierarchy is unchanged.
We present the resulting limits on the mixing matrix elements in figure 5 and in table 3. The different choices indeed affects the results, and bounds on the mixing matrix elements are stronger for the flat priors on m 4 or ∆m 2 41 . This is a consequence of the mapping between parameter spaces: the most important constraint from the data is on N eff (which also contributes to Σm ν through ∆N eff m 4 ), and there is a strong degeneracy between the mass splitting and the mixing matrix elements as explained in section 2 and seen from the upper left plot in figure 1. For a fixed higher sterile mass splitting, the mixing angles have to be smaller to stay within the allowed N eff region. Therefore the more weight the priors give to large masses m 4 or ∆m 2 41 , the smaller the matrix elements |U α4 | 2 have to be to fulfill the ∆N eff constraint. We want to emphasize that the choice of priors does not affect the Cosmological marginalized constraints on the mixing matrix elements |U α4 | 2 for flat priors on either log ∆m 2 41 (blue), m 4 (orange) or ∆m 2 41 (green), from CMB+BAO data. Off-diagonal panels show 68% and 95% confidence level probability contours. Panels along the diagonal show one-dimensional probability distributions. The more weight the prior gives to higher sterile masses, the lower the allowed |U i4 | 2 values in order to stay within the N eff range allowed by the cosmological data.
degree of tension with reactor andν e appearance measurements. A different choice of priors changes the weight of the mapping (∆m 2 41 , |U α4 | 2 ) → N eff , but the region in parameter space necessary to explain the anomalies leads to an almost thermal sterile with ∆N eff ≈ 1. This is excluded by cosmological data independent of the assumptions made in the analysis.
The resulting limits on the mixing matrix elements summarized in table 3 are robustly constrained to be |U α4 | 2 < 10 −3 for all flavors. While mixing with electron neutrinos is the most important channel and the resulting bounds on |U e4 | 2 are slightly more stringent, 1+1 ν e 1+1 ν µ 1+1 ν τ Figure 6. Left: Cosmological 68% and 95% marginalized constraints assuming either a 1 + 1 case with the sterile only coupled to ν e (purple), compared to the full 3 + 1 scenario including the full mixing matrix (blue) from figure 2. The difference is mostly a parameter space volume effect, since the limits obtained with a modified coupling only to ν e with an effective mixing |U eff | 2 = Σ α |U α4 | 2 (α ∈ [e, µ, τ ], red) are almost the same as for the 3 + 1 case. Right: Cosmological 68% and 95% marginalized constraints on |U eff | 2 for a single sterile coupled only to one neutrino species ν e (red, same as left), ν µ (green) or ν τ (grey). The differences between coupling to the different flavors are very small.
there is overall only a minor difference between mixing with the different active neutrinos. It is therefore interesting to understand the difference between a simplified case where the sterile is only coupled to one active neutrino (often assumed to be ν e ) and the full mixing with all flavors explored in this paper. We present the comparison between the cosmological constraints on ∆m 2 41 and |U e4 | 2 assuming either a 1 + 1 or a 3 + 1 scenario on the left side of figure 6. The results clearly differ and the 1 + 1 mixing allows for higher mass splittings of the sterile. However, due to parameter space volume effects it is expected that the limits look different. For every point in the (∆m 2 41 , |U e4 | 2 ) plane, in the 3 + 1 scenario there are more ways to achieve a higher ∆N eff and Σm ν by increasing the other mixing matrix elements. Higher mass splittings with small |U e4 | 2 that lead to negligible cosmological effects in the 1+1 model are then still unlikely in the 3 + 1 case since the other mixing angles have to be small as well. We explicitly test this explanation by adding another comparison case where the sterile is still only coupled to electron neutrinos, but we consider an effective mixing matrix element and we sample over three distinct contributions |U i4 | 2 to make up for the larger parameter space volume in the full 3 + 1 mixing scenario. As can be seen on the left side of figure 6, this parameter space volume effect almost completely accounts for the difference between the scenarios. On the right-hand side of figure 6 we test this effective mixing scenario by accounting for the parameter space effect in the same way, but coupling the sterile to either one of electron-, muon-or tau-neutrinos respectively. Again, differences between the individual flavors are very small, and an effective coupling to ν e provides a very good approximation to the full dynamics.
As a consequence we find that the production of light sterile neutrinos via mixing with to three active states can be modelled within a 1+1 model with a single sterile mixing with ν e as long as the effective mixing matrix element in Eq. (5.1) is used. Therefore, the abundance of a single sterile computed in the extended 3+1 model is, to an excellent approximation, very similar to that found in the effective 1+1 scenario, provided that the 1+1 squared mixing matrix element is the sum of the individual squared mixing matrix elements in the 3+1 case.

Conclusion
Several anomalies observed in short-baseline oscillation experiments hint at a new neutrino mass state at the eV scale. In this paper, we have provided a consistent framework to constrain the additional neutrino mass splitting and the active-sterile mixing angles with cosmological data. This also allows us to compare the bounds from cosmological data with results from oscillations, β-decay and 0νββ measurements in a common parameter space.
For the first time, we performed this analysis in a full 3 + 1 scenario where the sterile state is mixed with all three active neutrinos. In order to map the sterile neutrino masssquared splitting and mixing matrix parameters for each flavor ∆m 2 41 and |U α4 | 2 onto the cosmological observables, we have solved the evolution of the 4 × 4 neutrino density matrix in the early universe with the FortEPiaNO [78] code and find that the resulting distribution function of the fourth neutrino state is well approximated by a Dodelson-Widrow form [75].
A combination of CMB and BAO data currently provides the strongest available bounds on the sterile mass splitting ∆m 2 41 and the mixing matrix elements |U α4 | 2 . While the limit on the mass scale depends on prior assumptions, we robustly find that once any mixing matrix element reaches a level of |U α4 | 2 ≈ 10 −3 , the new state would give rise to a detectable relativistic energy contribution in the early universe not seen in cosmological data. The parameter space needed to explain short baseline anomalies with a sterile neutrino leads to a fully thermalized relativistic species with N eff ≈ 4 and is in strong tension with the CMB bounds. We also derive limits on the effective electron neutrino mass m β and the Majorana mass parameter m ββ , measurable in βand 0νββ-decay experiments, from cosmological data, finding m β < 0.09 eV and m ββ < 0.07 eV at 95% C.L. These constraints are tighter than the ones obtained from the latest direct laboratory measurements. Our main results in terms of limits on the sterile mass scale m 4 , the mixing matrix parameters for each flavor |U α4 | 2 , m β , and m ββ , are summarized in table 3.
We have explored the effect of prior choices on the cosmological bounds repeating the analysis with different priors on the mass splitting. Since a sizeable sterile contribution to N eff can be produced either by a large mass splittings or large mixing angles, the main effect of priors is to shift the weight between the two quantities. The maximal contribution to N eff allowed by the data is fixed, so the more weight the prior gives to high mass splittings ∆m 2 41 , the lower the mixing angles have to be to stay within the allowed region.
Almost all limits on sterile neutrinos from cosmology previously reported in the literature were derived in a simplified 1 + 1 scenario where the sterile is only coupled to one active neutrino. Even though the mixing with the different active flavors is almost equivalent, we show that the results differ significantly from the ones obtained with a 3 + 1 mixing scheme assumed for this work. This can be largely attributed to parameter space volume effects and can be accounted for by coupling the sterile with one active neutrino with effective mixing as given by Eq. (5.1). In any case, our results show that allowing for more than one active-sterile mixing does not relax the tension between the active-sterile neutrino parameters favored by the oscillation anomalies and present cosmological observations. Since cosmological data currently dominates the constraints, we do not perform a joint analysis with laboratory experiments at this time. However, the framework provided here can be the basis for new global constraints on light sterile neutrino properties once new data from laboratory searches becomes available.

A Individual 0νββ likelihoods
For completeness we also present the individual constraints on m ββ from GERDA, KamLAND ZEN I & II and EXO200 in figure 7. The width of each curve represents the uncertainty in the nuclear matrix element given in table 2. A modest excess of events observed at EXO200 makes the corresponding probability distribution for m ββ to peak at slightly non-zero values of m ββ . The combination of the individual constraints presented in Fig. 7 gives the probability distribution reported as a red solid line labeled 0ν2β in the right panel of Fig. 4.