Constraining light sterile neutrino mass with the BICEP2/Keck array 2014 B-mode polarization data

We explore the thermal light sterile neutrino situation from cosmological perspective in the ΛCDM+r0.05+Neff+mseff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda \text {CDM} + r_{0.05} + N_{\text {eff}} + m^{\text {eff}}_{\text {s}}$$\end{document} model using combinations of latest data sets available. Among CMB datasets, we use Planck 2015 temperature and low-l polarization data and the latest data release on the B-mode polarization from the BICEP2/Keck collaboration (BK14). We also use the latest Baryon Acoustic Oscillations (BAO) data from SDSS-III BOSS DR12, MGS, and 6dFS; and a Gaussian prior (HST) on the Hubble constant (H0=73.24±1.74\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_0 = 73.24 \pm 1.74$$\end{document} km/s/Mpc) from direct measurements by Hubble Space Telescope. We find that inclusion of BK14 data makes the constraints on the effective mass of sterile neutrino (mseff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m^{\text {eff}}_{\text {s}}$$\end{document}) slightly stronger by preferring higher σ8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _8$$\end{document} values. The bound of mseff<\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m^{\text {eff}}_{\text {s}}<$$\end{document} 0.46 eV (95% C.L.) is found for the combination of Planck 2015, BAO and BK14 datasets, whereas the bound is mseff<\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m^{\text {eff}}_{\text {s}}<$$\end{document} 0.53 eV (95% C.L.) without the BK14 data. Our most aggressive bound of mseff<\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m^{\text {eff}}_{\text {s}}<$$\end{document} 0.28 eV (95% C.L.) is obtained with Planck 2015, HST and BK14. Our analysis indicates that fully thermalized sterile neutrinos with mass ∼1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 1$$\end{document} eV are slightly more disfavoured with the inclusion of BK14 data. It also seems to make the agreement between Planck 2015 and CFHTLenS (weak gravitational lensing data) worse due to the higher σ8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _8$$\end{document} values.


Introduction
Sterile neutrinos still remain nothing short of an enigma in neutrino physics. Standard model predicts 3 massless neutrinos, while neutrino oscillation experiments have confirmed that neutrinos have mass, albeit very small. However, presence of anomalies in some short-baseline oscillation experiments [1][2][3][4][5][6] have been explained with an extra species of neutrino, namely a sterile neutrino, of mass 1 eV, which amply mixes with the active neutrinos but is uncharged under the standard model gauge group. Again, there are analyses [7][8][9][10][11][12][13][14][15] which indicate that all the results cannot be explained comfortably with the sterile neutrino hypothesis. A recent result [16] from the MiniBooNE collaboration finds present electron neutrino and anti-neutrino appearance data still consistent with an extra sterile neutrino.
Apart from terrestrial oscillation experiments, in recent years, cosmology has become a very powerful probe of neutrino physics. In a situation where standard model of particle physics is augmented with only an extra sterile neutrino species, there are two parameters of utmost importance. One is the effective number of relativistic neutrino species, N eff , whose theoretically predicted value, considering only the standard model of particle physics with 3 massless active neutrinos is N SM eff = 3.046 [17,18], but is supposed to increase when contribution from the sterile neutrino is counted. The other is the effective mass of the sterile neutrino, m eff s = N 3/4 eff m ph s , where N eff = N eff − 3.046 and m ph s is the physical mass of the sterile neutrino. Cosmology can provide strong constraints on these two parameters. N eff , in general, can have contribution from any relativistic species which is not a photon, and hence it is not restricted only to the neutrino sector. Also, in certain scenarios like very low-reheating scenarios with sterile neutrinos [19] or self-interacting sterile neutrinos [20], N eff can be negative. However, we do not consider such scenarios in this work, and consider only a non-interacting extra species of sterile neutrino.
Provided we are only considering an extension to standard model with neutrino oscillations in a 3+1 scenario, as long as the sterile neutrino is of similar mass to an active neutrino and amply mixes with the active ones, its cosmological implications are identical to the active neutrino. Sufficient mixing will lead to almost complete thermalization [21,22]. However, even if there is partial thermalization, it will, in general, increase N eff , leading to a delayed matter-radiation equality and a higher value of the Hubble parameter, H (z dec ), at the CMB decoupling (given other parameters are kept fixed). This has two main consequences [23] on the CMB anisotropy power spectrum, first being an increase of the first peak of the spectrum due to early Integrated Sachs Wolfe (ISW) effect, and the second being a horizontal shift of the peaks towards higher multipoles. Along with a horizontal shift, there will also be a vertical shift which will decrease the amplitude of the peaks at high multipoles, a phenomenon related to Silk damping. These effects of an additional relativistic sterile neutrino can be partially compensated if other cosmological parameters are simultaneously varied. For example, if the total matter density ω m is also increased without altering the baryon density, the redshift of matter-radiation equality can be kept fixed. These degeneracies tend to degrade the constraints on N eff . However, the CMB power spectra won't be exactly the same even after such adjustments with other parameters, especially because of the neutrino anisotropic stress arising from the quadrupole moment of the cosmic neutrino background temperature anisotropies which alters the gravitational potentials [24,25]. Hence constraints can be put on N eff from CMB power spectra data.
If a light sterile neutrino has a mass 1 eV, it only starts to become non-relativistic during CMB, and hence the effect of the mass is not strong on CMB power spectra. Sterile neutrinos with masses much smaller than 1 eV will have negligible effect on CMB power sprectra. However, when CMB power spectra data is used with other cosmological observations like constraining the Hubble parameter from direct measurements via a Gaussian prior or using the Baryon Acoustic Oscillation (BAO) data or both, better bounds on the mass of the sterile neutrino can be obtained [26]. A plethora of papers [26][27][28][29][30][31][32] are available on the effects of neutrino masses on cosmology. Current bounds on sterile neutrinos from cosmological data imply that fully thermalized sterile neutrinos of mass 1 eV are disfavoured and can only be accommodated with partial thermalization. See previous analyses on constraining sterile neutrino properties with cosmological data .
In this paper, we have, for the first time, used the BK14 data, the data on the B-mode polarization of CMB from BICEP2/Keck collaboration, to constrain the parameters associated with sterile neutrinos in an extended CDM model, which can be simply denoted with CDM + r 0.05 + N eff + m eff s . BK14 constrains the tensor-to-scalar ratio to r 0.05 < 0.07 at 95% C.L, when combined with Planck 2015 and other datasets [56]; while exclusion of the BK14 data leads to a significantly less strong bound of r 0.05 < 0.12 [57]. BK14 data also contains information on gravitational lensing. Thus we expect this data to affect the constraints on the sterile neutrino parameters. We also provide results with N eff fixed at 4.046 and 3.5 separately, i.e., assuming full and partial thermalization of the sterile neutrinos respectively, and this model is denoted as CDM + r 0.05 + m eff s . This paper has the following structure: in Sect. 2 we provide details about our model parameters and other analysis details and briefly describe the datasets used, in Sect. 3 we provide the results of our analysis, and we conclude in Sect. 4.

Model
Below we list the vector of parameters we have varied in this work in two cosmological models.
For CDM + r 0.05 + N eff + m eff s model: with N eff fixed to the value 4.046, which corresponds to full thermalization of the sterile neutrino with active neutrinos and to the value 3.5, which corresponds to partial thermalization.
The first six parameters correspond to the CDM model. Here ω c = c h 2 and ω b = b h 2 are the physical cold dark matter and baryon densities at present, respectively. s is the angular sound horizon, i.e., the ratio between sound horizon and the angular diameter distance at decoupling. τ is the reionization optical depth. n s and A s are the power-law spectral index and power of the inflationary power spectrum, respectively, at the pivot scale of k * = 0.05h Mpc −1 .
r 0.05 is the tensor-to-scalar ratio, also defined at the pivot scale of k * = 0.05h Mpc −1 . N eff , effective number of relativistic species which are not photons, is given by, where T γ is the temperature of the photons and ρ r is the radiation density. In our model, apart from photons, we only have 3 active and one sterile neutrino as relativistic species. The sterile neutrino is assumed not to have any self-interactions, or interactions with other particle species. In our work, we have fixed the active neutrino sector to give a contribution of N SM eff = 3.046 to N eff , with two massless and one massive neutrino with mass of 0.06 eV. Thus the contribution to N eff from the sterile species is simply N eff = N eff − 3.046.
Here it should be mentioned that in certain scenarios like low reheating Note that while we use the value 3.046 which is predominant in literature, a recent study [58] had found N SM eff = 3.045. When the sterile neutrino is relativistic at early times, assuming the only radiation species are photons and neutrinos, contribution of a light sterile neutrino to N eff is given by [59], where T ν is active neutrino temperature, p is the neutrino momentum, and f s ( p) is momentum distribution function of the sterile neutrino. At late times its energy density is parametrized as an effective mass [59,60]: ph s π 2 ρ c dp p 2 f s ( p), (2.5) where ρ c is the critical density, s h 2 is the sterile neutrino energy density. Since sterile neutrinos don't have electroweak interactions and they have mixing with the active neutrinos, they cannot decouple after the decoupling of active neutrinos. Active neutrinos decouple at a temperature T ∼ 1 MeV, when all of them are relativistic. Hence f s ( p) doesn't depend on the physical mass of the sterile neutrino, m ph s . However f s ( p) depends on the production mechanism of the light sterile neutrino. If the production is through a thermal process, one can simply write f s ( p) = (e p/T s +1) −1 , the usual Fermi-Dirac distribution function, where T s is the sterile neutrino temperature. In this case, it can be shown that, Non-thermal production, on the other hand, can lead to various possible scenarios. One of the popular scenarios is the Dodelson-Widrow (DW) mechanism [61], for which where β is a normalization factor. In this case, one gets [59], So, the m eff s parametrization can accommodate two different scenarios of sterile neutrino production. Also notice that in the CDM +r 0.05 + m eff s model, fixing N eff = 4.046 leads to m eff s being same as m ph s . In our work, we conduct a Bayesian analysis to derive constraints on the sterile neutrino parameters. For all the parameters listed in Eqs. (2.1), and (2.2), we impose flat priors. We also limit the physical mass of the sterile neutrino to m ph s ≤ 10 eV. The prior ranges are provided on the Table 1. We run chains using the November 2016 version of the Markov Chain Monte Carlo (MCMC) sampler CosmoMC [62] which incorporates CAMB [63] as the Boltzmann code and the Gelman and Rubin statistics [64] to estimate the convergence of chains.

Datasets
We use separate combinations of the following datasets: Cosmic Microwave Background: Planck 2015: Measurements of the CMB temperature and low-l polarization from Planck 2015 [65] are used. We consider the high-l (30 ≤ l ≤ 2508) TT likelihood, and also the low-l (2 ≤ l ≤ 29) TT likelihood. We refer to this combination as TT. We also include the Planck polarization data in the low-l (2 ≤ l ≤ 29) likelihood, and denote this as lowP. We also use the Planck lensing potential measurements via reconstruction through the four-point correlation functions of the Planck CMB data [66]. We call this simply as lensing. Residual systematics may be present in the the Planck 2015 high-l polarization data [57], so we refrain from using it.
B mode polarization data of CMB: Considering the B-mode polarization of CMB, we incorporate the recent dataset publicly available from BICEP2/ Keck collaboration which includes all data (multipole range:  20 < l < 330) taken up to and including 2014 [56]. This dataset is referred to as BK14.

Baryon acoustic oscillations (BAO) measurements and related galaxy cluster data:
In this analysis, we include measurements of the BAO signal obtained from different galaxy surveys. We make use of the SDSS-III BOSS DR12 [67] LOWZ and CMASS galaxy samples at z eff = 0.38, 0.51 and 0.61, the DR7 Main Galaxy Sample (MGS) at z eff = 0.15 [68], and the 6dFGS survey at z eff = 0.106 [69]. We call this complete combination as BAO. Here z eff is the effective redshift of a survey.
Hubble parameter measurements: We use a Gaussian prior of 73.24 ± 1.74 km/s/Mpc on H 0 , which is a recent 2.4% determination of the local value of the Hubble parameter by [70] which combines the anchor NGC 4258, Milky Way and LMC Cepheids. We denote this prior as HST.

Results
For convenience, we have separated the results in two subsections for the the two different models. The description of models and datasets are given at Sects. 2.1 and 2.2, respectively. We have presented the results, first in the CDM + r 0.05 + N eff +m eff s model, and then in the CDM+r 0.05 +m eff s model. All the marginalized limits quoted in the text or tables are at 68% C.L. whereas upper limits are quoted at 95% C.L., unless otherwise specified.
3.1 Results for CDM + r 0.05 + N eff + m eff s model In this section, we present the results for the CDM+r 0.05 + N eff +m eff s model. In Table 2 we have provided results without BK14 data, whereas, in Table 3, the results are with BK14, to compare. We have presented constraints on the three parameters r 0.05 , N eff , and m eff s . with which we have extended the CDM model, and also two derived parameters H 0 and σ 8 , which are important in constraining the sterile neutrino mass.
With only TT + lowP, we see that the bound on the sterile mass is relaxed at m eff s < 0.78 eV. The bound gets tightened with BAO data, which partially breaks the degeneracy between m eff s and H 0 present in the TT + lowP data, by rejecting lower values of H 0 [71,72] and leads to a bound of m eff s < 0.53 eV. This effect can be seen pictorially in Fig. 1 where addition of BAO data leads to a significantly smaller magnitude of anti-correlation between m eff s and H 0 . The HST prior also breaks the degeneracy partially, as can be seen in Fig. 1. However, the H 0 values preferred by the HST prior are larger than BAO, which leads to a preference to even smaller masses (m eff s < 0.34 eV) to keep the comoving distance to the surface of last scattering fixed [71]. Adding HST and BAO together with CMB however does not provide better bound than CMB+HST. Also, the lensing data degrades the bound on m eff s . We note that CMB and/or BAO data do not allow full thermalization of sterile neutrinos. However, at 95% C.L., with TT + lowP + HST, we obtained a N eff = 3.63 +0. 44 −0.42 . Such high values of N eff disallow the standard model prediction of N SM eff = 3.046 at 95% C.L. but allow N eff = 4.046, i.e., full thermalization. On the other hand, it is also imperative to consider recent constraints on N eff coming from Big Bang Nucleosynthesis (BBN). Planck 2018 results [73] have provided bound of N eff = 2.95 +0. 56 −0.52 (95% C.L.) (which is independent of the details of the CMB spectra at high multipoles) by combining the helium, deuterium, and BAO data with an almost model-independent prior on θ s derived from Planck data. Another recent study on BBN [74] provide a tight bound of N eff = 2.90 ± 0.22 (68% C.L.), which means at at 95% C.L., there will be only a small overlap in the values of N eff provided by [74] and TT + lowP + HST. It is also to be noted that addition of the HST prior leads to a slightly inferior fit to the data, due to the 3.4σ tension present between Planck and HST regarding the value of H 0 [70]. We find that in this CDM + r 0.05 + N eff + m eff s model, compared to TT + lowP, the dataset TT + lowP + HST degrades the χ 2 -fit by an amount of χ 2 = +3.43.
Akaike information criterion (AIC): To understand the improvement/worsening of the quality of fit with addition of sterile neutrino parameters (N eff and m eff s ) we need to compare the fit to data given by CDM +r 0.05 + N eff + m eff s with that of CDM + r 0.05 . Since the number of parameters in the two models are not same, a popular method to compare the fit is the Akaike information criterion (AIC) [75]. For a particular model and data, AIC is defined as, where k is the number of parameters in the model. The model with lower AIC corresponds to the preferred model. Thus, comparison with another model (with the same data) can be done with AIC = χ 2 +2 k. Usually models with extra parameters provide better fit to the data since they have a larger parameter space. The 2 k term penalises models with extra parameters to prevent any over-fitting. Here 2 k = 4.
We find that for the TT + lowP + HST data: i.e., the CDM + r 0.05 + N eff + m eff s model provides a better χ 2 fit compared to CDM + r 0.05 . But due to the 2 extra parameters, AIC = −0.3. Since this difference is small, it implies that the goodness of fits to the TT + lowP + HST data for the two models are similar.
Since the main aim of this paper is to analyze the role of the BK14 data, Table 3 lists the bounds on the cosmological parameters, now with BK14 data included in each combination. The inclusion of the BK14 data seems to have almost no effect on the bounds of N eff and H 0 , as can be seen by comparing the results of Tables 2 and 3. However, bounds on m eff s improve slightly across all data combinations. The 1-D marginalized posteriors for m eff s and r 0.05 for various datasets are shown in Fig. 2. While for TT + lowP, we had m eff s < 0.78 eV, this bound improves to m eff s < 0.68 eV with TT + lowP + BK14. Addition of BAO data further improves this bound to m eff s < 0.46 eV. Our most aggressive bound in this paper comes with TT + lowP + BK14 + HST: m eff s < 0.28 eV. BK14 data significantly constrains the tensor-to-scalar ratio, r 0.05 . TT + lowP provides r 0.05 < 0.127 whereas TT + lowP + BK14 gives a constraint of r 0.05 < 0.068. However, we found only a very small correlation between r 0.05 and m eff s , and that does not explain the decrease in mass. In fact the correlation coefficient (defined as R i j ≡ C i j / C ii C j j , where i and j are the two parameters being considered and C is the covariance matrix of cosmological parameters) between r 0.05 and m eff s to be R m eff s ,r 0.05 = −0.08 with TT + lowP and R m eff s ,r 0.05 = +0.02 with TT + lowP + BK14, i.e., there is no significant correlation before addition of BK14 and also no significant change after. However we also find slightly increased values of σ 8 across all data combinations when BK14 is included. For instance, for TT + lowP, we have σ 8 = 0.802 +0.040 −0.029 , which increases to σ 8 = 0.814 +0.036 −0.027 with TT + lowP + BK14. Since σ 8 is the normalization of matter power spectrum on scales of 8h −1 Mpc, a higher σ 8 prefers lower sterile neutrino mass, as larger neutrino masses create larger suppressions in the matter power spectrum [26]. Thus σ 8 and m eff s , both are strongly anti-correlated. Indeed, we found R σ 8 ,m eff s = −0.84 with TT + lowP and R r 0.05 ,m eff s = −0.81 with TT + lowP + BK14, and hence, even such small changes in σ 8 should also create small changes in m eff s , which we find is the case here. This has been depicted in Fig. 3. Again, notice that the lensing data prefers a lower σ 8 value. As in Table 3, TT + lowP + BK14 + HST + BAO yields σ 8 = 0.841 +0.025 −0.021 , whereas adding the lensing data to this combination yields a lower σ 8 = 0.820 +0.019 −0.015 . Due to the same anti-correlation between σ 8 and m eff s , we see that inclusion of lensing data degrades the m eff s bounds. Overall, we can say that the BK14 data makes the case for fully thermalized eV scale sterile neutrinos slightly worse. The parameter to justify this statement is m eff s . As we have shown that addition of the BK14 data does not affect the N eff bounds, BK14 data does not affect the thermalisation situation, as far as cosmological data is concerned. However, short baseline oscillation experiments predict a fully thermalised sterile neutrino of mass 1 eV. This requires that both N eff = 4.046 and m eff s 1 eV be allowed by the data. Since adding the BK14 data tightens the bounds on m eff s for all of the cosmological dataset combinations, it  also takes the m eff s value further away from the 1 eV value, while N eff bounds almost remain unchanged. Effect of BK14 data on sum of active neutrino masses ( m ν ) was also studied by us recently in [71], in the CDM + r 0.05 + m ν model, where we had also found slightly increased σ 8 . This is also indirectly confirmed by the recent Planck 2018 results, where they provide a bound of m ν < 0.12 eV with Planck TT,TE,EE+lowE+lensing+BAO data in CDM + m ν model [73], whereas the bound is m ν < 0.11 eV with Planck TT,TE,EE+lowE+lensing+BK14+BAO data in the CDM +r + m ν model [76]. This similar effect was seen to persist even in a 12 parameter extended scenario in a recent study with non-phantom dynamical dark energy [77]. In this paper we have shown that such an effect is also present in an extended CDM cosmology with light sterile neutrinos. CMB B-mode polarization has two known sources [78]. The first one is the inflationary gravitational waves (IGW), i.e., tensors (expected to produce a bump peaked around l 80, the so called 'recombination bump' in the BB-mode CMB spectra) as tensors induce quadruple anisotropies in the CMB within the last scattering surface. The tensor signature can-not be reproduced by scalar perturbations, and the amplitude of the recombination bump depends on the tensor-to-scalar ratio. The second source is gravitational lensing by large scale structure. It leads to deflection of CMB photons at late times, which converts a small part of the E mode power into B mode. This lensing BB spectra is expected to have a peak around l 1000. The BICEP2/Keck experiment has a multipole range 20 < l < 330 aiming to constrain the tensorto-scalar ratio. However since r 0.05 and m eff s are only weakly correlated, the slightly stronger constraints on the neutrino masses is possibly coming from gravitational lensing information encoded in the BK14 data, and not from measurement of r 0.05 .
H 0 and σ 8 tensions: It is also worth noting that in CDM model, with TT + lowP, Planck collaboration [57]   allows for N eff > 3.046 and higher N eff values prefer a higher H 0 , to keep the acoustic scale parameter θ s fixed [57], which is very well constrained by Planck data. Thus H 0 and N eff are strongly correlated. The CDM + r 0.05 + N eff + m eff s model also helps in reconciling the σ 8 tension present in the σ 8 − m plane in CDM model between Planck 2015 and weak lensing survey, like CFHTLenS [79] and KiDS-450 [80]. For instance, the KiDS-450 survey constrains the quantity S 8 ≡ σ 8 √ m /0.3 = 0.745 ± 0.039 which has a 2.3σ tension with Planck TT + lowP, which prefers a much higher value of S 8 = 0.851 ± 0.024 [57]. Planck data also prefers higher values of σ 8 compared to CFHTLenS. With TT + lowP in base CDM model, one gets σ 8 = 0.829 ± 0.014 [57].
However, in this CDM + r 0.05 + N eff + m eff s model, with TT + lowP, we get σ 8 = 0.802 +0.040 −0.029 , which is much lower and thereby the conflict is decreased somewhat. We also get S 8 = 0.824 +0.030 −0.027 , which is better agreement with KiDS-450 than CDM. However, the BK14 data prefers slightly higher σ 8 values and thereby increases the tension between Planck and these weak gravitational lensing surveys. This can be visualized in Fig. 4, where we see that the inclusion of BK14 data drives the 2-D contours upwards to a small extent. In Fig. 4, we have used the CFHTLenS data with conservative cuts as described in [57].
Another important point is that while CDM + r 0.05 + N eff +m eff s helps in relieving the H 0 and σ 8 tensions present in the CDM model, they are not both relieved together in any region of the allowed parameter space. In the right panel of Fig. 4, we can see that the regions where σ 8 has lower values, H 0 also has lower values (while we need higher values of H 0 to relieve the H 0 tension), and similarly, where H 0 has higher values, σ 8 also has higher values (while we need lower values of σ 8 to relieve the σ 8 tension). This in turn implies that the two conflicts are not resolved together in this model. And BK14 data worsens the conflicts even more. The HST prior also doesn't help the issue here. As we can see from Tables 2 and 3, the inclusion of this Gaussian prior leads to a preference for much higher N eff values, and higher σ 8 values as well, increasing the conflict.
3.2 Results for CDM + r 0.05 + m eff s model In this section we verify the stability of the results obtained in the previous section, by going to a smaller parameter space. We stop varying N eff and fix its value to 4.046 and 3.5. The first one corresponds to complete thermalization of sterile neutrinos, while the later one corresponds to partial thermalization. We have restricted ourselves to CMB data only. For N eff = 4.046 and N eff = 3.5, the results are given in Tables 4  and 5 respectively.   We see that BK14 does help in obtaining better constraint on the sterile mass also in this reduced parameter space. For N eff = 4.046, with TT + lowP, we get m eff s < 0.66 eV, whereas inclusion of BK14 leads to a tighter bound of m eff s < 0.50 eV. Similar case of strengthening of mass bound is seen with N eff = 3.5, although these bounds are more relaxed compared to the case N eff = 4.046, as a higher N eff prefers a higher H 0 . Again we see that the BK14 data itself does not affect the H 0 constraints much, but heavily constraints the tensor-to-scalar ratio, and also slightly increases the preferred σ 8 values. The main conclusions made in the previous section on the larger parameter space thus remains unchanged in this smaller parameter space.
It is imperative to note that for sterile neutrinos produced by a thermal process and obeying Eq. 2.6, for N eff = 4.046, we have m ph s = m eff s , whereas for N eff = 3.5, we have m ph s = 1.8m eff s . Hence, for N eff = 3.5 and with TT + lowP + BK14, we have a corresponding bound of m ph s < 1.13 eV. This implies that CMB data allows sterile neutrinos with mass 1 eV, but only with partial thermalization with N eff 3.5. When we compare the quality of fit to the TT + lowP + BK14 data between the CDM + r 0.05 + m eff s model (N eff = 3.5 and 4.046) and the CDM + r 0.05 (with N eff = N SM eff ), we find that, for the N eff = 4.046 case: best-fit ( CDM + r 0.05 + m eff s ) − χ 2 best-fit ( CDM + r 0.05 ) = +7.03 (3.3) whereas, for the N eff = 3.5 case: best-fit ( CDM + r 0.05 ) = −0.22 (3.4) These correspond to AIC = +9.03 (for N eff = 4.046) and AIC = +1.78 (for N eff = 3.5). Thus, the model with partial thermalization of N eff = 3.5 provides only a slightly worse fit to the data compared to the CDM + r 0.05 model (with N eff = N SM eff ), and is preferred by the data much more than the full-thermalization case. This is not surprising as in the previous section we had seen that CMB data alone did not allow complete thermalization.

Discussion
Short Baseline (SBL) Oscillation anomalies have hinted towards a fully thermalized sterile neutrino with mass around 1 eV. In this paper we have studied, for the first time, the light eV scale sterile neutrino situation in cosmology in light of the BICEP2/Keck array 2014 CMB B-mode polarization data. We call this dataset BK14. We first considered an extended− CDM scenario with tensor perturbations and sterile neutrino parameters: CDM + r 0.05 + N eff + m eff s model. Apart from BK14, we have used Planck 2015 temperature and low-l polarization data (TT + lowP), latest BAO data and a Gaussian prior on the Hubble constant (HST) from local measurements. We find that inclusion of the BK14 data has almost no effect on the bounds of N eff and H 0 but it strengthens the bounds on m eff s to a small extent by preferring slightly higher values of σ 8 , with which m eff s is strongly anti-correlated. The BK14 data also tightly constraints the tensor-to-scalar ratio, r 0.05 but we find negligible correlation between r 0.05 and m eff s . This makes us think that the effect on mass is coming from the gravitational lensing information encoded in the B-mode polarization and not from the Inflationary Gravitational Waves. The bound of m eff s < 0.46 eV (95% C.L.) is found for the combination of Planck 2015, BAO and BK14 datasets, whereas the bound is m eff s < 0.53 eV (95% C.L.) without the BK14 data. Our most aggressive bound of m eff s < 0.28 eV (95% C.L.) is obtained with Planck 2015, HST and BK14. The HST prior also leads to high N eff values which allow full thermalization of the sterile neutrino (at 2σ ) but such high values are in conflict with bounds from Big Bang Nucleosynthesis. Also, addition of the HST prior to the TT + lowP data leads to a slightly worse χ 2 fit to the data. On the other hand, it is to be noted that as per the Akaike information criterion (AIC) the CDM+r 0.05 + N eff +m eff s model provides equally good fit to the data as the CDM + r 0.05 model, for the TT + lowP + HST data combination. Previous studies have indicated that fully thermalized sterile neutrinos with mass ∼ 1 eV (as predicted by SBL experiments) are disfavoured by cosmological data. Our analysis indicates that it becomes slightly more disfavoured with the inclusion of BK14 data, due to tighter mass bounds. The BK14 data also seems to make the agreement between Planck 2015 and CFHTLenS (weak gravitational lensing data) worse due to the higher σ 8 values.
We would also like to mention that the Planck 2018 results, released during the preparation of this article, indirectly show tightening of bounds on m ν with BK14. They provide a bound of m ν < 0.12 eV with Planck TT, TE, EE + lowE + lensing + BAO data in CDM + m ν model [73], whereas the bound is m ν < 0.11 eV with Planck TT, TE, EE + lowE + lensing + BK14 + BAO data in the CDM + r + m ν model [76]. Thus we expect our main conclusion regarding BK14 helping in improving the bound on sterile neutrino mass will remain unchanged if used with the recent Planck 2018 data instead of Planck 2015 that we have used in this paper.
While this work was still being completed, a new Bmode polarisation data was released publicly, from the same BICEP2/Keck collaboration. This newly released data includes all the measurements upto and including 2015, and thus we call it BK15 [81]. To understand the effect of the new data, we performed an MCMC analysis with TT + lowP + HST + BK15 in the CDM + r 0.05 + N eff + m eff s model (with all other settings remaining unchanged). We found the following bounds: m eff s < 0.27 eV (95% C.L.), r 0.05 < 0.061, and σ 8 = 0.847 +0.026 −0.021 . In the same model, when we had used BK14 instead of BK15, we had found (see table 3), m eff s < 0.28 eV (95% C.L.), r 0.05 < 0.073, and σ 8 = 0.846 +0.026 −0.020 . As we can see, that while the bound on r 0.05 changes, the bounds on m eff s and σ 8 almost remain unchanged. We also checked that other parameters of interest, like H 0 and N eff change negligibly. As before, since r 0.05 and m eff s have only a very weak correlation, it doesn't affect the mass bound. On the other hand, since TT + lowP + HST + BK15 almost doesn't change the bound on σ 8 , the mass bound almost remains the same. Thus, we find that reanalysis with BK15 instead of BK14 will not change the neutrino mass bounds.
This tension between SBL and cosmological datasets has given rise to a number ideas to reconcile the eV-scale sterile neutrinos with cosmology. These include introduction of new "secret interactions" among sterile neutrinos which modifies the background potential and blocks thermalization [20,37,50,52,[82][83][84][85][86][87], modifications to the cosmic expansion rate at the time where sterile neutrinos are produced [88], large lepton asymmetry [89][90][91], time varying dark energy component [92], very low reheating temperature [93]. The recent results that have come from the MiniBooNE collaboration [16] have rekindled interest in the sterile neutrinos. Thus, there seems to be interesting and exciting times ahead in the realm of light sterile neutrinos in cosmology.
Acknowledgements The authors sincerely thank the anonymous referee for the thoughtful comments and efforts towards improving our manuscript. SRC thanks the cluster computing facility at HRI (http:// cluster.hri.res.in). SRC also thanks Steen Hannestad for useful discussions. The authors would also like to thank the Department of Atomic Energy (DAE) Neutrino Project of HRI. This project has received funding from the European Union's Horizon 2020 research and innovation programme InvisiblesPlus RISE under the Marie Sklodowska-Curie Grant Agreement No. 690575. This project has received funding from the European Union's Horizon 2020 research and innovation programme Elusives ITN under the Marie Sklodowska-Curie Grant Agreement No. 674896.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All the cosmological data we have used in this paper to analyze and obtain results are already publicly available from corresponding experimental collaborations.] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .