Statistical tests of sterile neutrinos using cosmology and short-baseline data

In this paper we revisit the question of the information which cosmology provides on the scenarios with sterile neutrinos invoked to describe the SBL anomalies using Bayesian statistical tests. We perform an analysis of the cosmological data in $\Lambda$CDM$+r+\nu_s$ cosmologies for different cosmological data combinations, and obtain the marginalized cosmological likelihood in terms of the two relevant parameters, the sterile neutrino mass $m_s$ and its contribution to the energy density of the early Universe $N_{\rm eff}$. We then present an analysis to quantify at which level a model with one sterile neutrino is (dis)favoured with respect to a model with only three active neutrinos, using results from both short-baseline experiments and cosmology. We study the dependence of the results on the cosmological data considered, in particular on the inclusion of the recent BICEP2 results and the SZ cluster data from the Planck mission. We find that only when the cluster data is included the model with one extra sterile neutrino can become more favoured that the model with only the three active ones provided the sterile neutrino contribution to radiation density is suppressed with respect to the fully thermalized scenario. We have also quantified the level of (in)compatibility between the sterile neutrino masses implied by the cosmological and SBL results.


Introduction
It is now an established fact that neutrinos are massive and leptonic flavors are not symmetries of Nature [1,2]. In the last decade this picture has become fully established thanks to the upcoming of a set of precise experiments. In particular, the results obtained with solar and atmospheric neutrinos have been confirmed in experiments using terrestrial beams [3]. The minimum joint description of all these data requires mixing among all the three known neutrinos (ν e , ν µ , ν τ ), which can be expressed as quantum superposition of three massive states ν i (i = 1, 2, 3) with masses m i leading to the observed oscillation signals with ∆m 2 21 = (7.5 +0.19 −0.17 ) × 10 −5 eV 2 and |∆m 2 31 | = (2.458 ± 0.002) × 10 −3 eV 2 and non-zero values of the three mixing angles [4].
In addition to these well-established results, there remains a set of anomalies in neutrino data at relatively short-baselines (SBL) (see Ref. [5] for a review) including the LSND [6] and MiniBooNE [7,8] observed ν µ → ν e transitions, and theν e disappearance at reactor [9][10][11] and Gallium [12][13][14] experiments. If interpreted in terms of oscillations, each of these anomalies points out towards a ∆m 2 ∼ O(eV 2 ) [15][16][17][18][19][20][21][22][23] and consequently cannot be explained within the context of the 3ν mixing described above. They require instead the addition of one or more neutrino states which must be sterile, i.e. elusive to Standard Model interactions, to account for the constraint of the invisible Z width which limits the number of light weak-interacting neutrinos to be 2.984 ± 0.008 [24].
Several combined analyses have been performed to globally account for these anomalies in addition to all other oscillation results in the context of models with one or two additional sterile neutrinos [20][21][22][23] with somehow different conclusions in what respects to the possibility of a successful joint description of all the data. Generically these global fits reveal a tension between disappearance and appearance results. But while Refs. [21,22] seem to find a possible compromise solution for 3+1 mass schemes, the analysis in Ref. [23] concludes a significantly lower level of compatibility. In particular, Ref. [21] concludes that a joint solution is found with 0.82 ≤ ∆m 2 41 ≤ 2.19 eV 2 at 3σ (1.1) Alternative information on the presence of light sterile neutrinos is provided by Cosmology as they contribute as extra radiation to the energy density of the early Universe which can be expressed as where ρ γ is the photon energy density and the value of N eff in the Standard Model (SM) is equal to N SM eff = 3.046 [25]. The presence of extra radiation is then usually quantified in terms of the parameter ∆N eff ≡ N eff − N SM eff . Sterile neutrinos contribute to ρ R (i.e. to ∆N eff ) in a quantity which, in the absence of other forms of new physics, is a function of their mass and their mixing with the active neutrinos which determines to what degree they are in thermal equilibrium with those. In particular, in 3 + N s scenarios with the values hinted by the SBL results (of the order of those in Eq. (1.2)) the sterile neutrinos are fully thermalized (FT) and each one contributes to ρ R as much as an active one (see for example [26,27]), so ∆N 3+Ns,F T eff = N s . (1.3) Analyses of cosmological data have hinted for the presence of extra radiation, beyond the standard three active neutrinos since several years (see for example Ref. [28] and references therein) and several authors have invoked the presence of eV-scale sterile neutrino as a plausible source of the extra radiation [29][30][31][32][33] 1 . In the last four years the statistical significance of this extra radiation (or its upper bound) as well as the overall constraint on the neutrino mass scale has been changing as data from the Planck satellite [38], the Atacama Cosmology Telescope (ACT) [39], the South Pole Telescope (SPT) [40,41], and most recently of the BICEP2 [42] experiment, has become available [43][44][45][46][47][48][49][50].
In particular, the recent measurement of B-mode signals from the BICEP2 [42] collaboration excludes a zero scalar-to-tensor ratio at 7σ and report a value of r = 0.20 +0.07 −0.05 . This high value of r is compatible with previous results from the Planck data [38] if a running of the scalar spectral index dn s /d ln k is considered well beyond the characteristic value of 10 −4 of slow-roll inflation models. Alternatively the tension might be eased by the presence of sterile neutrinos [44][45][46][47][48] without invoking such a large running of the scalar spectral index.
Since cosmological data have the potential to test regions of parameter space of the sterile neutrino scenarios invoked to account for the SBL anomalies, the question of to what degree they support them has received an increasing attention in the literature [33,[51][52][53][54]. The generic conclusion is that, in order to accommodate the cosmological observations within the 3 + N s scenarios motivated by SBL results, some new form of physics is required to suppress the contribution of the sterile neutrinos to the radiation component of the energy density at the CMB epoch with respect to the FT expectation Eq. (1.3). Among others extended scenarios with a time varying dark energy component [31], entropy production after neutrino decoupling [55], very low reheating temperature [56], large lepton asymmetry [26,57,58], and non-standard neutrino interactions [59][60][61], have been considered. All these mechanisms have the effect of diluting the sterile neutrino abundance or suppressing the production in the early Universe.
In this paper we revisit the question of the information which cosmology provides on the sterile scenarios introduced to explain the SBL anomalies using precise Bayesian statistical tests which we briefly describe in Sec. 2. Section 3 contains the results of our cosmological analysis in ΛCDM+r + ν s cosmologies for three representative sets of cosmological data. With these results at hand we answer the question of how much cosmological data favour or disfavour the scenario with sterile neutrino masses invoked by SBL anomalies, with respect to a model without sterile neutrinos? in Sec. 4, and we do so in terms of the departure from the fully thermalized expectation, Eq. (1.3). We also discuss the consistency of sterile parameter constrains implied by cosmology and SBL. In Sec. 5 we summarize our conclusions.

Model comparison
Bayesian inference is a rigorous framework for inferring which, out of set of models or hypotheses H i , are favoured by a data set D. Bayes' theorem is used to calculate the probabilities of each of the hypotheses after considering the data, the posterior probabilities, Here Pr(D|H i ) is the probability of the data, assuming the model H i to be true, while Pr(H i ) is the prior probability of H i , which is how plausible H i is before considering the data. Considering especially the case of a discrete set of models, one can compare two of them by calculating the ratio of posterior probabilities, the posterior odds, as In words, the posterior odds is given by the prior odds Pr(H i )/ Pr(H j ) multiplied by the Bayes factor B ij = Pr(D|H i )/Pr(D|H j ), which quantifies how much better H i describes that data than H j . The prior odds quantifies how much more plausible one model is than the other a priory, i.e., without considering the data. If there is no reason to favour one of the models over the other, the prior odds equals unity, in which case the posterior odds equals the Bayes factor. For a model H, containing the continuous free parameters Θ, Pr(D|H) also called evidence of the model is given by Here, the likelihood function Pr(D|Θ, H) is the probability (density) of the data as a function of the assumed free parameters, which we often denote by L(Θ) for simplicity. The quantity Pr(Θ|H) is the correctly normalized prior probability (density) of the parameters and is often denoted by π(Θ). The assignment of priors is often far from trivial, but an important part of a Bayesian analysis. From Eq. (2.3), we note that the evidence is the average of the likelihood over the prior, and hence this method automatically implements a form of Occam's razor, since usually a theory with a smaller parameter space will have a larger evidence than a more complicated one, unless the latter can fit the data substantially better.
Bayes factors, or rather posterior odds, are usually interpreted or "translated" into ordinary language using the so-called Jeffreys scale, given in Tab. 1, where "log" is the natural logarithm. This has been used in applications such as Refs. [62][63][64] (and Refs. [65,66] in neutrino physics), although slightly more aggressive scales have been used previously [67,68].

Parameters of interest and the marginal likelihood
In Bayesian statistics if one assumes a particular parametrized model to be correct, the complete inference of the parameters of that model is given by the posterior distribution We see that the evidence here appears as a normalization constant in the denominator. Since the evidence does not depend on the values of the parameters Θ, it is usually ignored in parameter estimation problems and the parameter inference is obtained using the unnormalized posterior. For the case in which we have a only a subset of parameters of interest, λ, so Θ = (λ, η), where η denotes the nuisance parameters, P (λ, η) = Pr(λ, η|D, H) ∝ L(λ, η)π(λ, η) , (2.5) and inference on λ is obtained by marginalizing over the nuisance parameters in the usual way P (λ) = P (λ, η)dη , (2.6) with no need to ever consider a likelihood L(λ) depending only on the parameter of interest. This is typically unproblematic in the case where the data is sufficiently informative to eliminate all practical prior dependence. Often, however, this is not the case, and there can be large dependence on the prior chosen. In this case one can consider as a partial step the marginal likelihood function such that P (λ) ∝ L(λ)π(λ). This likelihood function then encodes the information on λ contained in the data (under H and after taking into account the uncertainty on the nuisance parameters), without needing to specify a prior π(λ). Note that, since the marginal likelihood is not a probability density, it is not normalized to unity, and is not sufficient to perform the full inference. Also, it is generally different from the profile likelihood, so regions defined by −2 log(L(λ)/L max ) < C will not, in general, be the same as those using the profile likelihood, although arguments justifying defining regions in this way for profile likelihoods typically also apply to marginal likelihoods. As shown in Ref. [69], the profile and marginal likelihoods are indeed similar for the cosmological data and models considered there. In any case, the marginal likelihood can still be useful in scientific reporting as a rough guide to what information the data contains, for example, by considering regions defined by −2 log(L(λ)/L max ) < C [70].
Furthermore, if two data sets do not share any common nuisance parameters, the two marginal likelihoods can simply be multiplied to obtain the total marginal likelihood. Notice, however, that the marginal likelihood still depends on the priors on the nuisance parameters. If the nuisance parameters are well-constrained this dependence will be small, but in cosmology this is necessarily not always the case.

Data sets
In our cosmological analysis, we use data on Cosmic Microwave Background, large scale structure (LSS) baryon acoustic oscillations (BAO) measurements, Hubble constant H 0 , and galaxy cluster counts. In particular, we define the following data combinations: • CMB: It includes the current Planck data [38] of the temperature anisotropy up to l = 2479, the high multipole values (highL), coming from ACT [39] and SPT data [40,41], that covers respectively the 500 < l < 3500 and 600 < l < 3000 range, and the EE and TE polarization data from WMPA9 [71] (WP). It also includes the CMB lensing potential (lensing) reconstructed by the Planck collaboration [38] through the measurement of the four-point function.
• BAO: It includes the Data Release 11 (DR11) sample of the recent measurements by the BOSS collaboration [72]. The DR11 sample is the largest region of the Universe ever surveyed, covering roughly 8500 square degrees, with a redshift range 0.2 < z < 0.7. The measure of the sound horizon at the drag epoch has been evaluated at redshift z=0.32 and at z=0.57, finding values in agreement with previous BAO measurements 2 .
• BICEP2: the 9 channels of the CMB BB polarization spectrum recently released by the BICEP2 experiment [42].
• HST: the data from the Hubble Space Telescope [77] on H 0 , obtained through the distance measurements of the Cepheids: • PlaSZ: The counts of rich cluster of galaxies from the sample of Planck thermal Sunyaev-Zel'Dovich catalogue [78]. It constrains the combination In order to test the dependence of our results on the inclusion of the recent BICEP2 data and on the tension with local HST and cluster PlaSZ results we perform the analysis with three different combinations of the data sets above that we label: • DATA SET 1: CMB+BAO, where, as described above, CMB=Planck+WP+highL+lensing data, and BAO=DR11. 2 We have verified that the inclusion of the previous determinations of the BAO from the Sloan Digital Sky Survey (SDSS) Data Release 7 (DR7) [73,74], the 6dF Galaxy Survey (6dFGS) [75], and WiggleZ measurements Ref. [76] does not affect our results and therefore for the sake of simplicity we have not included them in the analysis. 3 For simplicity we do not include the cosmic shear data of weak lensing from the Canada-French-Hawaii Telescope Lensing Survey (CFHTLenS) which constraints σ8(Ωm/0.27) 0.46 = 0.774±0.040, [79][80][81]. Although the combinations constrained by PlaSZ and CFHTLenS are not the same, given the much better precision of the PlaSZ data, the impact of including CFHTLenS in our analysis is very small.
• DATA SET 2: CMB+BAO+BICEP2. We add to the previous data set, the results from BICEP2.
• DATA SET 3: CMB+BAO+BICEP2+HST+PlaSZ. We add to the previous data set the results from HST and Planck SZ counts of galaxy clusters.

Cosmological model
We consider in our analysis a ΛCDM cosmology extended with a free scalar-to-tensor ratio, and three active plus one sterile neutrino species with a hierarchical neutrino spectra of the 3+1 type which we denote as ΛCDM+r + ν s . In this case the three active neutrinos have masses m i=1,3 |∆m 2 31 | while the forth sterile neutrino has a mass m 4 ≡ m s . As mentioned in the introduction, in the absence of other form of new physics the contribution of the sterile neutrino to the energy density is completely determined by its mass and its mixing with the active neutrinos. However in extended scenarios this may not be the case. So generically we will consider that in the 3+1 scenario, irrespective of m s and mixings, the sterile neutrino contributes to ρ R as where F N T is an arbitrary quantity which quantifies the departure from the fully-thermalized active-sterile neutrino scenario and which we will consider to be independent of T in the relevant range of T in the analysis. So in what follows we will label as ∆N eff ≡ F N T this parameter. The effect of m s is included in the analysis via the effective parameter with being h the reduced Hubble constant and Ω s ≡ ρ s /ρ c , with ρ s the sterile neutrino energy density and ρ c the current critical density. This effective mass is not equal to the physical mass m s in general, and their relation depends on the assumed phase-space distribution of the sterile neutrinos. For thermally distributed sterile neutrinos characterized by a temperature T s (in general different from the temperature of the active neutrinos T ν ) while if they are produced by non-resonant oscillations (the so-called Dodelson-Widrow scenario) (DW) [82] the resulting phase-space distribution of the sterile neutrinos is equal to that of the active neutrinos up a constant factor. In this case Altogether our analysis contain nine free parameters where ω b ≡ Ω b h 2 , ω c ≡ Ω c h 2 being the physical baryon and cold dark matter energy densities, Θ s the ratio between the sound horizon and the angular diameter distance at

Parameter
Prior decoupling, τ is the reionization optical depth, A s the amplitude of primordial spectrum, n s the scalar spectral index, and r the scalar-to-tensor ratio. To generate the marginalized likelihoods we use the CosmoMC package [83], implemented with the Boltzmann CAMB code [84]. The parameters in (3.6) are assigned uniform priors with limits as given in Tab. 2. Since we are interested in the constraints on sterile neutrino parameters, we follow Sec. 2.2, and aim to evaluate the marginal likelihoods of (m eff , N eff ), which will then be used for the tests presented in Sec. 4. Thus in this analysis we consider {ω b , ω c , Θ s , τ, log[10 10 A s ], n s , r} as the cosmological nuisance parameters 4 . Note that we have employed uniform priors on m eff and N eff in the numerical analysis but that the marginal likelihoods to be presented below do not depend on the priors, since these are "factorized out". Furthermore, since the nuisance parameters are rather well-constrained the precise vales of the limits in Tab. 2 do not not affect any results, and also employing different shapes on these parameters would have a small impact. The exception is r which is not so well-constrained (especially for certain data sets) and for which the physical lower limit is important.

Marginal cosmological likelihoods
The results of our analysis for the three data set combinations described in Sec. 3.1 are shown in Figs. 1 and 2.
In the the upper left panel of Fig. 1 and in the left panels of Fig. 2 we plot the contours of the marginal likelihood L(N eff , m eff ) normalized to the value of L(∆N eff = 0) (which does not depend on m eff ) 5 . The red contour delimits the regions for which m s in Eq. (3.4) exceeds 10 eV for which the sterile states becomes indistinguishable from cold or warm dark 4 We do not include the lensing amplitude AL as a free nuisance parameter, even when adding the local measurements on σ8 and Ωm, for which a significant deviation from the standard value of unity is preferred. We have checked that adding AL only slightly shifts the preferred region for m eff to higher values. 5 For likelihoods more than eight log-units away from the maximum value, we extrapolate using a con-matter [38] and hence for values of the parameters below the curve, the marginal likelihood is not evaluated. Black dashed contours are those of −2 log(L/L max ) < 2.30, 6.18, 11.83, nominally corresponding to the 1σ, 2σ, and 3σ levels in two dimensions.
In the upper right panel of Fig. 1 and the right panels of Fig. 2 we show the contours of the marginal likelihood of (N eff , m s ) for the thermally distributed sterile neutrinos while in the lower left panel of Fig. 1 we show the corresponding marginal likelihood in the DW scenario. We see that for both DW and thermal ν s scenarios, m s becomes increasingly large for decreasing ∆N eff , hence the distinctive appearance of large "flat" regions and weak constraints on m s for small ∆N eff . Also as seen in Fig. 1  , for ∆N eff ≤ 1 the likelihood contours in the DW scenario are shifted to slightly larger masses in an amount which decreases as ∆N eff increases.
The impact of BICEP2 data can be seen by comparing the upper panels in Fig. 1 and Fig. 2. The addition of BICEP2 data gives a preference for large values of the tensor-toscalar ratio r and, due to its correlation with ∆N eff , and therefore leads to the the shift to slightly larger values of ∆N eff observed in upper panels in Fig. 2. Since r is now better constrained, also the preferred region becomes slightly smaller.
The effects of adding HST+PlaSZ in the analysis are displayed in the lower panels in Fig. 2 where we see a shift to larger values of both ∆N eff and m eff . This is so because the constraints on σ 8 of the PlaSZ measurement, which are in tension with the other experiments within this model, can be somewhat alleviated by an increase in m eff , while the inclusion of HST yields an increase ∆N eff .
All these results are in qualitative agreement with those in the analyses in Refs. [44][45][46][47][48][49][50]. However we notice that by showing the marginal likelihood we are explicitly not assigning any priors to the sterile parameters. This is an advantage in the absence of a physical motivation for them, especially m s , since the data at hand is expected to leave a large prior dependence. For example, if one used a prior on m s which is uniform in log m s instead of in m s , as in general the likelihood is non-negligible for a vanishing sterile mass, the derived Bayesian constraints on m eff and ∆N eff would be very different. In principle one could embark on an extensive prior sensitivity analysis, but in this work we will instead focus on analyses for which the results have little or no prior dependence.
Finally, in order to better illustrate how the constraints depend on the sterile mass we plot in Fig. 3, the slices of the marginal likelihood as a function of m s for fixed ∆N eff . Also shown in the figure is the marginal likelihood for the SBL analysis in the 3+1 scenario (marginalized with respect to the lighter neutrino masses and all mixings) as given in Fig.1 in Ref. [52].

Test on sterile neutrinos models
In this section we perform the statistical tests on the 3+1 scenarios invoked to explain the SBL anomalies using the results of the cosmological analysis presented in the previous stant value. It could be made more accurate by using non-constant functions such as polynomials but no qualitative change is expected.  section. In doing so we are going to assume that the contribution of the sterile neutrino to the cosmological observables is independent of the active-sterile neutrino mixing (see discussion around Eq. (3.2)). Under this assumption our tests will require the marginalized likelihood of the sterile neutrino mass m s obtained from the analysis of the SBL anomalies in the 3+1 scenario (marginalized with respect to all other oscillation parameters in the scenario), and the marginalized likelihoods of (m s , ∆N eff ) from the cosmological analysis previously presented.
Concerning the SBL likelihood, we consider two different functions that can be interpreted as the two limiting cases. In the first, we consider the precise SBL likelihood as given in Fig. 1 in Ref. [52] which we reproduce in Fig. 3 (in this case below roughly seven log-units from the maximum value, we set the SBL likelihood to a constant value). We CMB+BAO+BICEP2 Thermal ν s CMB+BAO+BICEP2+HST+PlaSZ Figure 2. Same two upper panels in Fig. 1, but for CMB+BAO+BICEP2 cosmological data (SET 2, upper panels), and CMB+BAO+BICEP2+HST+PlaSZ cosmological data (SET 3, lower panels). label this case in the following as "full SBL likelihood". In the second case, we approximate the SBL likelihood as a top-hat shaped likelihood which is constant and non-zero between 0.86 eV and 1.57 eV and zero otherwise (which we illustrate by the arrow in Fig. 3). This is what we label in the following as "box SBL likelihood".

Sterile neutrinos vs none
The first question we want to address is whether the current data shows evidence of the existence of sterile neutrinos, and how strong this evidence is. Generically in Bayesian analyses this takes the form of model comparison between a model without sterile neutrinos and a model with sterile neutrinos using some posterior odds, Eq. (2.2). There are several ways to go about answering this question. In here we are interested in testing what cosmology has to say on this comparison for the sterile models invoked to explain the SBL anomalies. This is, in this case the first model In all panels we also include the marginal likelihood for the SBL analysis in the 3+1 scenario (marginalized with respect to the lighter neutrino masses and all mixings) as given in Fig.1 in Ref. [52]. We denote by the red arrow the width and height of the box used to define "box SBL likelihood" (see text for details).
is taken to include one sterile neutrino of mass m s as required to accommodate the SBL anomalies and which contributes to relativistic energy density in the early Universe as some fix ∆N eff .
In this case we can define a posterior odds as: and where in the last line we have used that SBL is not sensitive to any parameters effecting cosmology once we assume that there are no sterile neutrinos and hence D SBL and D c are independent under H 0 . The quantity B upd 10 quantifies how much better the prediction of cosmological data assuming sterile neutrino and SBL data is than the prediction assuming no sterile neutrinos. Now, the model H 1 is inherently more complex than H 0 since it contains additional parameters (sterile mass and mixings), and this is as usual compensated for in a Bayesian analysis. In the first row of Eq. In any case, it is the first factor in Eq. (4.1) the factor by which the cosmological data have updated the SBL-only posterior odds to the final SBL+cosmology odds and which we will be using in our quantification. Furthermore, under the additional assumption that ∆N eff is unconstrained by SBL data, Pr(D SBL |H 1 ) does not depend on ∆N eff , and hence B upd 10 is in fact simply proportional to the combined marginal likelihood of ∆N eff , normalized such that B upd 10 = 1 for ∆N eff = 0. Before discussing the results, let us mention that here, as in any Bayesian analysis, the results are in principle always prior dependent, and we should consider how large this dependence is in practice. First, as discussed in Sec. 3, the marginal likelihood depends on the priors on the cosmological nuisance parameters, but this dependence is expected to be small (except possibly for r). More significantly, there is the dependence of the total Bayes factor and odds on the prior on m s , even when it is well constrained by the combined data set. However the value of B upd 10 does not strongly depend on the the shape of the prior nor its upper limit. This is so because the well-constraining SBL data is used to update the prior to a posterior (which is rather insensitive to the prior) which is then used to analyze the cosmological data. We use a uniform prior between 0 and 10 eV as the nominal upper limit, since, as described in Ref. [38], this roughly defines the region where (for the CMB) the particles are distinct from cold or warm dark matter.
The results are shown in the left panel of Fig. 4. We see that for the CMB+BAO and CMB+BAO+BICEP2 data sets, B upd 10 decreases quite steadily with ∆N eff . These cosmological data hence disfavour models with sterile neutrinos required to explain the SBL anomalies over the model without sterile neutrinos independently of how much the contribution of the sterile neutrino to the energy density is suppressed with respect to the fully thermalized expectation. We also see that the addition of the BICEP2 data has a small impact on these conclusions.
For the CMB+BAO+BICEP2+HST+PlaSZ data set, there is, on the contrary, a significant peak for intermediate values of ∆N eff . The ΛCDM+r model is significantly disfavoured by this combination of cosmological data, while increasing ∆N eff increases the cosmological likelihood for SBL-compatible masses. However, further increasing ∆N eff , cosmology requires a too small mass, so the B upd 10 decrease again. Notice also that in fact the ΛCDM+r is so disfavored that it is in the region where the (too large) constant extrapolation is used. Hence, the exact B upd 10 is expected to be even larger.

Consistency of parameter constraints
In addition to comparing the models with and without sterile neutrinos we now address the question of the consistency of the parameter constraints from the different data sets within the 3+1 model. A Bayesian test was formulated in [85], in which a model where both data sets are fitted by the same physical parameters is compared with a model in which each data set uses their own parameters. However, as is often the case in model comparison, the result can depend crucially on the prior on the parameters, in our case m s in particular. Since we cannot motivate using a specific shape or limits on the prior, we instead use the corresponding χ 2 -test (these were also compared in [85]). Although not rigorous as the Bayesian test, it has the clear advantage that it is prior-independent.
In particular, if we want to test how inconsistent the constraints on the sterile mass from the different data sets are once a certain ∆N eff is assumed, i.e., without considering how favoured or disfavoured that ∆N eff is by the cosmological data, we should evaluate where the hat denotes the value at the best fit i.e., optimized over m s (χ 2 SBL does not depend on ∆N eff ).
The results of this test are shown on the right-hand side of Fig. 4. As expected, the results for CMB+BAO and CMB+BAO+BICEP2 are quite similar and show an steady increase of the inconsistency with ∆N eff . Also comparing the DW and thermal scenarios, we see that in general to obtain the same ∆χ 2 larger values of ∆N eff are required in the DW scenario. This is so because, as explained before, in the DW scenario the preferred masses are shifted to larger values by an amount which increases as ∆N eff < 1 decreases. As for the results for the box and full SBL likelihoods, we notice that, typically, using the full SBL likelihood gives a larger inconsistency. This is because the box likelihood is constant over a wide range of m s , and over this wide range the cosmological likelihood typically varies significantly. The combined χ 2 can then easily be reduced by finding the best fit within this range.
For the CMB+BAO+BICEP2+HST+PlaSZ combination, one observes in each curve a peak for small value of ∆N eff , although it does not reach the level of 2.5σ. These peaks are due to the fact that (cf. Fig. 2) cosmology prefers large physical masses, too large to fit the SBL data. As ∆N eff increases, the mass constraints become compatible, but as ∆N eff continues to decrease, cosmology requires the masses to be smaller than those which can fit the SBL data and the inconsistency increases.
Again, we stress that this consistency test is in principle not affected by how favoured or disfavoured the considered value of ∆N eff is, but instead consider that value to be "true", and then test the compatibility of the constraints on the mass. For example, comparing the left and right panels in the last raw of Fig. 4 we see that for the CMB+BAO+BICEP2+HST+PlaSZ data set, even though from the right panel we read that the mass ranges required for cosmology and SBL become highly incompatible for ∆N eff close to one, from the left hand side we see that these large values of ∆N eff are not particularly disfavoured compared to small ∆N eff , for which the mass constraints are compatible. So what we see is that large ∆N eff is disfavoured because the sterile mass required by SBL and cosmology are incompatible and small ∆N eff is disfavoured because it is so in the cosmological (and consequently in the combined) analysis. So small and large ∆N eff have comparable Bayes factors, but this is only because they are both disfavoured.

Summary
In this paper we have revisited the question of the information which cosmology provides on the scenarios with O(eV) mass sterile neutrinos invoked to explain the SBL anomalies (Eq. (1.2)) using Bayesian statistical tests and study how the results depend on the inclusion of the recently CMB polarization results of BICEP2 and on the inclusion of local measurements which show some tension with the Planck and LSS-BAO results when analyzed in the framework of the ΛCDM scenario.
In order to do so we have first performed an analysis of three characteristic sets of cosmological data in ΛCDM+r + ν s cosmologies as described in Sec. 3. The result of our analysis is presented in Figs. 1 and 2 in the form of marginalized cosmological likelihoods in terms of the two relevant parameters, the sterile neutrino mass m s and its contribution to the energy density of the early Universe ∆N eff . The results clearly indicate that as long as the HST and SZ cluster data from Planck are not included, cosmological data favours the sterile neutrino mass m s clearly well below eV unless its contribution to the energy density is suppressed with respect to the expected from a fully thermalized sterile neutrino. The inclusion of the BICEP2 data does not substantially affect this conclusion. Conversely including these HST and SZ cluster data higher sterile masses become favoured.
With these results, we have performed in Sec. 4 two statistical test on their (in)compatibility with the corresponding likelihood derived from the analysis of the SBL results as given in Ref. [52]. In the first test we have asked ourselves whether cosmology favours or disfavours  Fig. 4 implies that as long as the HST and SZ cluster data from Plank is not included, the cosmological analysis disfavour the 3+1 model with respect to 3+0. The inclusion of these cosmological data however favours the 3+1 model for an intermediate range of ∆N eff . We summarize in Tab. 3 the ranges of ∆N eff for which we find the evidence against or in favour of the 3+1 model compared to the model with only the 3 active neutrinos to be weak, moderate or strong from these analyses.
The second test performed deals with the (in)compatibility of the sterile mass constraints as required to describe SBL and cosmology. For this we have evaluated the ∆χ 2 defined in Eq. (4.4) which we plot in the right panels in Fig. 4. Altogether we read that this test yields inconsistency on the m s required by cosmology and SBL larger than 3σ for In summary, we find that the analysis of cosmological results from temperature and polarization data on the CMB as well as from the BAO measurements from LSS data disfavours the 3+1 sterile models introduced to explain the SBL anomalies over the scenario without sterile neutrinos, and also that their allowed/required ranges of m s are incompatible. This is so even if new physics is involved so that the contribution of the sterile neutrino to the energy density of the Universe (and therefore to the cosmological observables) is sup-pressed with respect to that of the fully thermalized case resulting from its mixing with the active neutrinos. When the local measurement of the H 0 by the Hubble Space Telescope, and the cluster SZ cluster data from the Planck mission is included, compatibility can be found between cosmological and SBL data, but still requires a substantial suppression of the ν s contribution to ρ R .  10 as a function of ∆N eff . Right panels: Consistency of mass constraints. In all panels the results are shown for thermal ν s (solid lines) and for the DW scenario (dashed lines) and for the SBL full (black) and Box likelihood (red). We show the results for the three cosmological data sets considered: CMB+BAO (upper row), CMB+BAO+BICEP2 (middle row), CMB+BAO+BICEP2+HST+PlaSZ (lower row).