Direct comparison of sterile neutrino constraints from cosmological data, $\nu_{e}$ disappearance data and $\nu_{\mu}\rightarrow\nu_{e}$ appearance data in a $3+1$ model

We present a quantitative, direct comparison of constraints on sterile neutrinos derived from neutrino oscillation experiments and from Planck data, interpreted assuming standard cosmological evolution. We extend a $1+1$ model, which is used to compare exclusions contours at the 95% CL derived from Planck data to those from $\nu_{e}$-disappearance measurements, to a $3+1$ model. This allows us to compare the Planck constraints with those obtained through $\nu_{\mu}\rightarrow\nu_{e}$ appearance searches, which are sensitive to more than one active-sterile mixing angle. We find that the cosmological data fully exclude the allowed regions published by the LSND, MiniBooNE and Neutrino-4 collaborations, and those from the gallium and rector anomalies, at the 95% CL. Compared to the exclusion regions from the Daya Bay $\nu_{e}$-disappearance search, the Planck data are more strongly excluding above $|\Delta m^{2}_{41}|\approx 0.1\, \mathrm{eV}^{2}$ and $m_\mathrm{eff}^\mathrm{sterile}\approx 0.2\, \mathrm{eV}$, with the Daya Bay exclusion being stronger below these values. Compared to the combined Daya Bay/Bugey/MINOS exclusion region on $\nu_{\mu}\rightarrow\nu_{e}$ appearance, the Planck data is more strongly excluding above $\Delta m^{2}_{41}\approx 5\times 10^{-2}\,\mathrm{eV}^{2}$, with the exclusion strengths of the Planck data and the Daya Bay/Bugey/MINOS combination becoming comparable below this value.


Introduction
The LSND [1], MiniBooNE [2], and Neutrino-4 [3] collaborations have made observations consistent with anomalous neutrino flavour oscillations. Other, related anomalies have been measured with gallium detectors [4] and reactor neutrinos [5]. These observations suggest that additional neutrino flavours may exist at a mass scale of O(1 eV), beyond the three flavours of the Standard Model.
Measurements of the decay width of the Z boson [6] conclusively show that only three neutrino flavours with m ν < m Z /2 couple through the weak interaction; these three flavours are termed "active", and any additional flavours are therefore referred to as "sterile". The existence of a sterile neutrino can have observable effects since neutrino oscillations allow the sterile flavour states to mix with the active flavour states. Such mixing occurs as the neutrino mass eigenstates are related to the flavour eigenstates through a mixing matrix, the PMNS matrix [7][8][9]. The minimal phenomenological 3 + 1 model of sterile neutrinos adds a single sterile flavour state and a fourth mass eigenstate.
Limits on the existence of sterile neutrinos have been set by observations of the cosmic microwave background (CMB) [10] and by numerous neutrino oscillation experiments [11][12][13][14][15][16][17]. In a commonly used model, cosmological measurements set limits on the parameter ∆N eff , the additional number of relativistic degrees of freedom in the universe arising from the additional neutrino states, and m sterile eff , the effective mass of the sterile neutrino. In a 3+1 model, neutrino oscillation experiments set limits on the mass splitting ∆m 2 41 = m 2 4 − m 2 1 , the difference between the squared masses of the additional, fourth mass eigenstate and the lightest neutrino eigenstate, along with the elements of the 4 × 4 PMNS matrix.
Quantitative comparisons of cosmological and neutrinooscillation limits on sterile neutrinos are complicated due to this difference in parameterization. In a previous article [18], a comparison using a phenomenological model in which only the muon-neutrino flavour mixes into the fourth mass eigenstate was presented. This 1 + 1 model allows only comparisons of ν µ disappearance measurements to the cosmological limits. Other studies [19] have investigated the situation in which only the electron-neutrino flavour is assumed to mix into the fourth mass eigenstate. In this article, we extend our treatment to the full 3 + 1 model, allowing us to compare the cosmological limits to both ν e and ν e disappearance measurements, including the gallium and reactor anomalies, and to the LSND and MiniBooNE ν µ → ν e and ν µ → ν e observations. In addition, we develop a treatment that allows us to extend our comparisons into the so-called degenerate region in which the sterile mass-splitting ∆m 2 41 becomes equal to the mass splitting ∆m 2 31 .

Sterile neutrinos in oscillation experiments
In the 3 + 1 model, four neutrino flavour eigenstates, ν l (l = e, µ, τ, s), are related to four neutrino mass eigenstates, ν i (i = 1, 2, 3, 4), with masses m i , by a 4 × 4 extension of the PMNS matrix, U: Throughout this paper, we assume all neutrino and antineutrino oscillation probabilities are equal and therefore use the symbol ν to also refer to ν. If a neutrino of energy E is produced in a flavour eigenstate ν l , the probability that it is detected in flavour eigenstate ν l after traveling a distance L is An experiment searching for ν e or ν µ disappearance thus measures where ∆m 2 ji = m 2 j − m 2 i are the mass splittings. Each mass splitting therefore defines an observable oscillation wavelength, with the elements of the PMNS matrix governing the amplitudes of those oscillations.
Over the majority of the parameter space relevant to sterile-neutrino searches, |∆m 2 41 | |∆m 2 31 | > |∆m 2 21 |. Thus, we can choose L and E to probe only the oscillations at the ∆m 2 41 wavelength, allowing us to approximate the disappearance probabilities in Eq. 3 to Here, we have introduced the mixing angles θ i j that are used to parameterize the PMNS matrix. We refer to this approximation of the oscillation probabilities as a 1 + 1 model since it assumes only one mass splitting, neglecting the effects of ∆m 2 31 and ∆m 2 21 , and assuming only one flavour state at a time (either electron or muon) mixes into the fourth mass eigenstate. The mixing angle θ 14 quantifies how much electron flavour mixes into the fourth mass eigenstate, and the angle θ 24 quantifies this mixing for the muon flavour. In this paper, we use the 1 + 1 model for an analysis of ν e disappearance.

Data from oscillation experiments
We use data from collaborations that report allowed regions consistent with sterile neutrino oscillations. Such regions have been reported by the LSND, MiniBooNE, and Neutrino-4 collaborations, in addition to the regions allowed by the reactor and gallium anomalies. We then compare to the exclusion region from Daya Bay, combined with Bugey-3 and MINOS data, which provides stronger exclusion at lower values of the mass of the fourth mass eigenstate, where the sensitivity of the Planck results decreases.

LSND
The Liquid Scintillator Neutrino Detector (LSND) took data from 1993-1998 at the Los Alamos Meson Physics Facility. A 167 t liquid scintillator detector was placed 30 m away from a stopped-pion source that produced ν µ with energies up to 52.8 MeV [21]. Appearance of ν e was observed in the detector with a total excess of 87.9 ± 22.4(stat.) ± 6.0(syst.) ν e events above the expected background [1]. To explain this excess through oscillations, a mass splitting ∆m 2 41 0.03 eV 2 is required.
We determine the 90% Confidence Level (CL) allowed region by requiring χ 2 − χ 2 min = 4.605 between the observed positron energy spectrum and an estimated spectrum. The appearance spectrum is simulated with pseudo-experiments, producing a reconstructed neutrino energy from a reconstructed positron energy and angle, and integrating the reconstructed neutrino energy over the same binning as in Ref. [1]. The true positron energy, E e + , is the difference between the true neutrino energy and the threshold energy of 1.806 MeV. The ν e → e + cross section is estimated to be linear in E e + . The reconstructed positron energy is smeared by a Gaussian function of the form 7%/ E e + /52.8 MeV, and its angle is Gaussian-smeared by 12 • . The distance L that the neutrino has travelled is uniformly spread in the range [25.85, 34.15] m, and a 14 cm Gaussian smearing is applied to produce a reconstructed distance. The flux is determined for pions decaying at rest to ν µ , with an L −2 weighting applied. The true neutrino energy and distance is used to calculate the oscillated ν e flux with Eq. 6.

MiniBooNE
The MiniBooNE experiment was an 818 t mineral oil Cherenkov detector [22] 541 m away from the neutrino-production target of the Booster Neutrino Beam [23]. The beam could be configured to produce either ν µ or ν µ with mean energy of ≈ 800 MeV. By searching for the appearance of either ν e or ν e , the experiment was sensitive to oscillations driven by a similar range of ∆m 2 41 as LSND. An excess of activity consistent with ν e and ν e was observed. We use the CL contours from the Collaboration's public data release [24].

Neutrino-4
The Neutrino-4 experiment [3] searches for the disappearance of ν e from the SM3 reactor in Russia. A gadoliniumdoped liquid scintillator detector is divided into 50 sections that can be placed at various distances, from 6 to 12 m, from the reactor core. The data analysis yields an oscillatory pattern to the ν e detection rate as a function of L/E that is interpreted in terms of a sterile neutrino with best-fit oscillation parameters ∆m 2 41 = 7.34 eV 2 , sin 2 (2θ 14 ) = 0.44. We take the 95% CL allowed region directly from Ref. [3].

Reactor anomaly
The reactor anomaly, first described in Ref. [5], is the observation that, with more modern flux calculations, many short-baseline reactor-ν e searches show a deficit compared to the expected flux. This observation can be interpreted as ν e disappearance due to oscillations involving a sterile neutrino. We use the 95% CL allowed region calculated in Ref. [25].

Gallium anomaly
The gallium anomaly, first described in Ref. [4], measured the ν e rate from radioactive calibration sources in the SAGE and GALLEX solar-neutrino detectors. A deficit in the measured rate compared to the expectation can be interpreted as ν e disappearance due to oscillations involving a sterile neutrino. We use the 95% CL allowed region calculated in Ref. [25].

Daya Bay
The Daya Bay experiment consists of eight gadoliniumdoped liquid scintillator detectors that measure the disappearance of electron antineutrinos from the Daya Bay and Ling Ao nuclear power plants in China [26]. The arrangement of eight detectors and six reactor cores provides a range of baselines between 358 m and 1925 m. The Daya Bay experiment was designed to be sensitive to oscillations driven by ∆m 2 31 and θ 13 [27]; however, by looking for non-standard ν e disappearance, Daya Bay can also search for oscillations driven by ∆m 2 41 and θ 14 in the range 10 −4 ∆m 2 41 0.1 eV 2 [13]. We use the Daya Bay data release [28] to recreate the χ 2 surface, and follow the prescribed approach [29], based on the CL s method [30,31], to produce the 95% CL exclusion contour.

Bugey-3
The Bugey-3 experiment took data in the early 1990s. The experiment used two lithium-doped liquid-scintillator detectors [32] to search for the disappearance of ν e at distances of 15 m, 40 m and 95 m from the Bugey nuclear power plant in France [33]. The shorter baseline provides sensitivity to sterile neutrinos at a higher range of ∆m 2 41 compared to Daya Bay.

MINOS
The MINOS experiment used two steel-scintillator calorimeters [34] to search for the disappearance of muon neutrinos and antineutrinos from the NuMI beam at Fermilab [35] at baselines of 1.04 km and 735 km. MINOS was designed to be sensitive to oscillations driven by ∆m 2 31 and θ 23 [36]. By searching for non-standard ν µ and ν µ disappearance at higher energies, it is also sensitive to oscillations driven by the sterile-neutrino parameters ∆m 2 41 and θ 24 [11].
3.9 Combination of Daya Bay, Bugey-3, and MINOS Data The Daya Bay limit was combined with that of Bugey-3 and MINOS to produce limits on the parameters ∆m 2 41 and sin 2 (2θ µe ) that govern ν µ → ν e appearance [14]. In performing this combination, the analysis of the Bugey-3 data was updated to use a more recent calculation of the neutron lifetime in the cross-section of the inverse-β decay process that is used for ν e detection. In addition, the ILL+Vogel flux model [37,38] was replaced with the Huber-Mueller model [39,40]. We use the combined CL s surface of Ref.

Sterile neutrinos in cosmological measurements
The presence of one or more sterile neutrinos can affect the power spectrum of the CMB. The effective mass of the sterile neutrino is defined as m sterile eff = 94.1 Ω sterile h 2 eV, where h = H/100 with the Hubble parameter H, and Ω sterile is the contribution of sterile neutrinos to the matter energydensity in the Universe. The neutrino number density, n ν , is expressed as a function of the number of effective neutrino species, N eff , as where n γ is the number density of photons in the CMB. Standard cosmology predicts N eff = 3.046, since the process of neutrino decoupling from the CMB was not instantaneous, and neutrinos still interacted with leptons in the primordial plasma [42]. This allows us to define the effective number of additional radiative degrees of freedom, equivalent to the effective number of additional neutrino species, as ∆N eff = N eff − 3.046. We relate m sterile eff and the mass of the fourth neutrino mass eigenstate, m 4 using the standard relationship [10] m sterile Here, we assume a thermally distributed sterile neutrino with a temperature T s that may differ from the active neutrino thermalisation temperature T ν A fully thermalized sterile neutrino with temperature T s = T ν corresponds to a measured ∆N eff = 1 and m sterile eff = m 4 . An alternative relationship between m sterile eff and m 4 , the Dodelson-Widrow mechanism [43], assumes that ∆N eff acts as a linear scaling factor, m sterile eff = ∆N eff m 4 . The choice of this function does not significantly impact our results.

The Planck experiment
The Planck satellite made detailed observations of anisotropies of the CMB between 2009 and 2013, over a frequency range from 30 to 857 GHz [44,45]. The Planck Collaboration combines data from the TT, TE and EE power spectra, the low-multipole EE power spectrum (LowE), CMB lensing, and baryon acoustic oscillations (BAO) to set limits of N eff < 3.29 and m sterile eff < 0.23 eV [10]. These results arise from the use of a flat prior in the range 0 < m sterile eff < 10 eV. A more restrictive prior results in more constraining limits. A flat prior in the range 0 < ∆N eff < 1 is also used. The Planck analysis assumes a normal neutrino-mass ordering and active states with masses m 1 = m 2 = 0 and m 3 = 0.06 eV.
To obtain these limits on sterile neutrinos, the Planck Collaboration fits the data using a ΛCDM + m sterile eff + ∆N eff model, varying neutrino and nuisance parameters to build a large number of points in the parameter space. The number density of these points is proportional to the probability density. We produce exclusion limits in the (∆N eff , m sterile eff ) space by integrating over 95% of the probability density relative to the maximum value, using the Markov chain Monte Carlo data released by the Planck collaboration [10,46].

Electron neutrino disappearance in a 1 + 1 model
To translate from the parameter space (∆N eff , m sterile eff ) to the parameter space (sin 2 2θ 14 , |∆m 2 41 |), we use LASAGNA [47] for calculating ∆N eff as a function of the mass splitting |∆m 2 41 | and mixing angle sin 2 2θ 14 . LASAGNA solves the quantum kinetic equations describing neutrino thermalization in the early universe by evolving the equations over a temperature range for input values of |∆m 2 41 | and sin 2 (2θ 14 ). Limits from neutrino disappearance experiments can be interpreted in the 1 + 1 model, which assumes that only one active flavour state mixes into the fourth mass state and that the three other mass states form a single, mass-degenerate state, ν d . For ν e disappearance experiments, we allow only the ν e flavour to mix into the fourth mass state. This is equivalent to varying θ 14 whilst fixing θ 24 = θ 34 = 0. In this model, we can write ν e = cos θ 14 ν d − sin θ 14 ν 4 , (9) ν s = sin θ 14 ν d + cos θ 14 ν 4 .
LASAGNA calculates the Bloch vectors for neutrinos and (P 0 , P) for anti-neutrinos using the 1 + 1 model. The resulting vector P + s = (P 0 + P 0 )+(P z + P z ) enters the expression where the momentum distribution, p, of the neutrinos is assumed to obey a Fermi-Dirac distribution at temperature T. A temperature range of T = [40,1] MeV covers the period from the beginning to the end of decoupling. We assume the lepton asymmetry, L = (n l − n l )/n γ , to be zero. It was shown (d) Fig. 1 (a, b) Cosmological parameters ∆N sterile eff and m sterile eff calculated, using LASAGNA, in the oscillation space of the 1 + 1 model that is relevant for ν e and ν e disappearance measurements. We use the thermal sterile neutrino mass (Eq. 8) and assume vanishing lepton asymmetry (L = 0). We also show the Daya Bay exclusion contour; the region to the right of the contour is ruled out at the 95% CL. (c, d) The oscillation parameters of the 1 + 1 electron-neutrino disappearance model, ∆m 2 41 and sin 2 (2θ 14 ), in the cosmological space (m sterile eff , ∆N eff ). The region above the blue line is excluded by the Planck TT, TE, EE and low-multipole EE power spectra at the 95% CL. A prior of m 4 < 10 eV is applied, shown by the hatched region that has not been considered in our probability density estimation. in Ref. [18] that the Planck exclusion region is significantly reduced in a 1+1 model for ν µ disappearance for large lepton asymmetries (L = 10 −2 ).
We use LASAGNA to calculate ∆N eff for a grid in the oscillation parameter space of |∆m 2 41 | ≡ |m 2 4 − m 2 d | and sin 2 (2θ 14 ), as shown in Fig. 1(a). Equation 8 allows us to express this result for all relevant combinations of ∆N eff , m sterile eff , sin 2 (2θ 14 ), and |∆m 2 41 | (Figs. 1(b)-1(d)). In Fig. 2(a) we express the Planck exclusion limit in the parameter space (sin 2 (2θ 14 ), |∆m 2 41 |) and overlay the limit from Daya Bay, and the allowed regions from Neutrino-4 and the gallium and reactor anomalies. The equivalent contours translated into the cosmological parameter space (m sterile eff , ∆N eff ) are shown in Fig. 2(b). In both figures, we show the Planck limit with and without the BAO and CMB lensing data.
The limits obtained using the Planck data with and without the BAO and CMB lensing data are strongly constraining in both parameter spaces in the region above |∆m 2 41 | 2 ≈ 0.1 eV 2 and m sterile eff ≈ 0.2 eV, and exclude the allowed regions from the Neutrino-4 experiment, and from the gallium and reactor anomalies. The Daya Bay experiment is sensitive to the regions of low |∆m 2 41 | and m sterile eff , where the cosmological data are less constraining.   . 2 (a) shows, in the neutrino-oscillation parameter space, limits on the electron flavour mixing with the fourth mass state, using a 1 + 1 model. The exclusion region from the Daya Bay oscillation experiment, and the allowed regions from the Neutrino-4 experiment and the reactor anomaly, come from searches for ν e disappearance. The allowed region from the gallium anomaly comes from a search for ν e disappearance. For the Daya Bay line, everything to the right is ruled out at 95% CL. The solid blue line labeled 'Planck' shows the exclusion using the Planck TT, TE, EE and low-multipole EE power spectra, using Eq. 8 to relate m sterile eff to m 4 , with the region to the right ruled out at 95% CL. The dashed blue line shows the impact of further including CMB and BAO data into the Planck limit, again using Eq. 8. The dashed grey line illustrates the impact on the Planck limit (the solid blue) of using the mean momentum approximation (MMA). Graph (b) shows the same set of limits (minus the MMA line) in the cosmological parameter space. The Neutrino-4 and gallium-anomaly lines are no longer visible as they are compressed up along the ∆N eff = 1 axis. The hatched region corresponds to the prior of m 4 < 10 eV assumed in the Planck analysis.

Electron neutrino appearance in a 3 + 1 model
When considering sin 2 (2θ µe ), both mixing angles θ 14 and θ 24 must be allowed to be non-zero to allow both ν e and ν µ flavours to mix into the ν 4 state, and so we work in the 3 + 1 model with one sterile and three active neutrino flavours, albeit setting θ 34 = 0. This model can be solved exactly [48] but working with the full momentum dependence of the quantum kinetic equations is computationally very intensive. Hence, we use the mean momentum approximation (MMA) following the prescription of Ref. [49] summarized below.
The neutrino density matrix, ρ(x, y) = ρ ee ρ eµ ρ eτ ρ es ρ µe ρ µµ ρ µτ ρ µs ρ τe ρ τµ ρ ττ ρ τs ρ se ρ sµ ρ sτ ρ ss , depends on the mixing angles and mass splittings. It can be written as a function of reduced time, x ≡ m/T, and reduced momentum, y ≡ p/T, where m is an arbitrary mass scale and T is the initial temperature of the thermal, active neutrinos. This matrix is used to calculate ∆N eff for any required values of θ 14 , θ 24 and ∆m 2 41 as The MMA assumes that the momentum dependence of ρ(x, y) can be factorized out as a Fermi-Dirac distribution, ρ(x, y) → f FD (y)ρ(x). The equations of motion for the neutrino and anti-neutrino density matrices is then written assuming that all neutrinos have the same momentum, y . We solve the resulting differential equations of motion numerically with an implicit Runge-Kutta algorithm of order 5, RADAU5 [50], using a publicly available C++ implementation [51]. To evaluate ∆N eff , we evolve the density matrix from T = 100 MeV to T = 1 MeV. To project the cosmological limits onto the sin 2 (2θ µe ) axis, we minimise the value of ∆N eff as a function of θ 14 and θ 24 along a contour of constant sin 2 (2θ µe ); the derived 95% confidence limits therefore assume the maximum possible thermalisation for a given value of θ µe . The resulting values of ∆N eff as a function of ∆m 2 41 and sin 2 (2θ µe ) are shown in Fig 3. In the region |∆m 2 41 | |∆m 2 31 |, the mass splitting ∆m 2 41 is driving neutrino oscillations at wavelengths similar to those driven by the active-neutrino mass splittings. This is referred to as the degenerate region, and in this region the RADAU5 solver slows down drastically due to the stiffness of the problem when degeneracies are crossed. To mitigate this, we increase the tolerance by a factor of 10 after every 100,000 steps of the algorithm, starting from a default tolerance of 10 −10 , reaching a maximum tolerance of 10 −4 required for certain parameters to converge quickly.
We evaluate the impact of the MMA by repeating the ν e disappearance analysis in the 1 + 1 model using this approximation. The result of this is shown in Fig 2(a), illustrating that, under the MMA, the cosmological exclusion contours expressed in the (∆m 2 41 , sin 2 (2θ 14 )) parameter space become slightly weaker.
In Figure 4, we show the Planck exclusion contours, with and without the BAO and CMB lensing data, in the (∆m 2 41 , sin 2 (2θ µe )) parameter space. We compare this to the limits from the Daya Bay/Bugey/MINOS combination, and the allowed regions from the LSND and MiniBooNE ν e → ν µ searches. The Planck exclusion region strongly excludes the entirety of the LSND and MiniBooNE allowed regions. The Daya Bay/Bugey/MINOS combined exclusion region is comparable in its exclusion power to that from the Planck data for mass splittings below ∆m 2 41 ≈ 5 × 10 −2 eV 2 and becomes more constraining below ∆m 2 41 ≈ 10 −3 eV 2 .

Conclusions
The discovery of a sterile neutrino would have major implications for the field of particle physics. The presence of both possible observations from neutrino oscillation experiments such as LSND and MiniBooNE, negative results from other oscillation experiments, and negative results from cosmological experiments, have left the field in an ambiguous situation. A particular challenge in drawing conclusions is quantitative comparison of limits from neutrino oscillation data with those from cosmology, due to the different parameter spaces in which measurements from these two sets are expressed.
In this article, we discuss a procedure to convert limits on sterile neutrinos between the (|∆m 2 41 |, θ 14 , θ 24 ) parameter space of neutrino oscillation physics and the (m sterile eff , ∆N eff ) parameter space of cosmology. We use the LASAGNA software package to solve the quantum kinetic equations of neutrinos in the early universe in a 1 + 1 model, allowing us to compare the exclusion regions obtained from Planck data with both allowed regions and exclusion regions from ν e and ν e disappearance searches. In a 3 + 1 model, we use a mean momentum approximation to solve the quantum kinetic equations, allowing us to compare the Planck exclusion with allowed regions and exclusion regions corresponding to ν µ → ν e searches. We find that the Planck data strongly excludes the allowed regions from the Neutrino-4, LSND and MiniBooNE experiments, as well as from the gallium and reactor anomalies. Compared to the Daya Bay exclusion region from ν e disappearance, Planck is much more constraining above |∆m 2 41 | ≈ 0.1 eV 2 and m sterile eff ≈ 0.2 eV, whereas at lower values, Daya Bay provides a more stringent exclusion on θ 14 . The Planck data provide the strongest exclusion on the θ µe parameter that describes ν µ → ν e appearance above ∆m 2 41 ≈ 5 × 10 −2 eV 2 ; below this value, the Daya Bay/Bugey/MINOS combination becomes comparable in terms of its exclusion power.