A tale of two modes: Neutrino free-streaming in the early universe

# A tale of two modes: Neutrino free-streaming in the early universe

## Abstract

We present updated constraints on the free-streaming nature of cosmological neutrinos from cosmic microwave background (CMB) temperature and polarization power spectra, baryonic acoustic oscillation data, and distance ladder measurements of the Hubble constant. Specifically, we consider a Fermi-like four-fermion interaction between massless neutrinos, characterized by an effective coupling constant , and resulting in a neutrino opacity . Using a conservative flat prior on the parameter , we find a bimodal posterior distribution with two clearly separated regions of high probability. The first of these modes is consistent with the standard CDM cosmology and corresponds to neutrinos decoupling at redshift , that is before the Fourier modes probed by the CMB damping tail enter the causal horizon. The other mode of the posterior, dubbed the “interacting neutrino mode”, corresponds to neutrino decoupling occurring within a narrow redshift window centered around . This mode is characterized by a high value of the effective neutrino coupling constant, (68% C.L.), together with a lower value of the scalar spectral index and amplitude of fluctuations, and a higher value of the Hubble parameter. Using both a maximum likelihood analysis and the ratio of the two mode’s Bayesian evidence, we find the interacting neutrino mode to be statistically disfavored compared to the standard CDM cosmology, and determine this result to be largely driven by the low- CMB temperature data. Interestingly, the addition of CMB polarization and direct Hubble constant measurements significantly raises the statistical significance of this secondary mode, indicating that new physics in the neutrino sector could help explain the difference between local measurements of , and those inferred from CMB data. A robust consequence of our results is that neutrinos must be free streaming long before the epoch of matter-radiation equality in order to fit current cosmological data.

\affiliation

McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, 5000 Forbes Ave., Pittsburgh, PA 15213, USA \affiliation Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, United Kingdom \affiliation Department of Physics, Harvard University, 17 Oxford St., Cambridge, MA 02138, USA \affiliation Physics Department, University of California, Davis, CA 95616, USA \emailAddlachlanl@princeton.edu \emailAddfcyrraci@physics.harvard.edu \emailAddlknox@ucdavis.edu \emailAddzhpan@ucdavis.edu \keywordsNeutrino, Cosmology, Cosmic Microwave Background

## 1 Introduction

Due to their extremely weak interactions with other known particles of the Standard Model (SM), deepening our understanding of neutrinos has become a game of large-scale experiments and statistics. The cosmic microwave background (CMB) and other probes of the large-scale structure of the Universe have proven themselves to be excellent probes of certain aspects of neutrino physics. This is best exemplified by the limits that have been placed on the sum of neutrino masses eV [1, 2, 3, 4, 5, 6, 7], or by the constraints on the effective number of neutrino species [8].

Within the SM, neutrinos interact exclusively through the electroweak force mediated by the and vector bosons [9]. These interactions, due to their weak nature, fall out of equilibrium early on in the history of the Universe when the SM plasma temperature was [10]. After decoupling from other SM particles, neutrinos are thought to free stream throughout the Universe, only interacting with other species through their gravitational interactions. Within the standard cosmological paradigm, free-streaming neutrinos are the primary source of anisotropic stress in the early Universe. This causes the neutrino background to have a significant impact on the growth of cosmological perturbations in the radiation dominated universe. Most notably, a nonvanishing anisotropic stress results in an offset between the two Bardeen potentials in the conformal Newtonian gauge (see e.g. Ref. [11]), which in turn affects the evolution of CMB fluctuations leading to measurable amplitude shifts in the CMB temperature and polarization spectra \citep1972ApJ…176..323S,1973ApJ…180….1P,Chacko:2003dt,Bashinsky:2003tk,Friedland:2007vv,Hou:2011ec. The supersonic (relative to the photon-baryon plasma) perturbations of free-streaming neutrinos additionally cause a phase shift in the oscillations of the power spectra which would not otherwise be observed \citepBashinsky:2007yc,Bashinsky:2003tk, Follin:2015hya, Baumann15, Pan16,2017arXiv170310732B. Moreover, the absence of neutrino clustering on scales entering the causal horizon before neutrinos become nonrelativistic leads to a suppression of the matter power spectrum which could in principle be observed with large-scale structure and CMB lensing observations [23].

So far, the standard evolution of cosmological neutrinos outlined above appears in good agreement with observations [24, 25, 26, 8]. Still, the existence of neutrino mass [27, 28, 29] implies that some sort of physics beyond the SM must couple to neutrinos. Minimal modifications to the SM have been proposed (see e.g. refs. [30, 31, 32, 33]) to accommodate neutrino masses in ways that leave the standard evolution of cosmological neutrinos essentially unchanged. However, it is possible that the physics giving rise to neutrino mass does change the cosmological evolution of the neutrino bath, leading to detectable effects on the CMB and large-scale structures. For instance, neutrino interactions mediated by a not-yet- observed particle [34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 14, 44, 45, 46, 47, 48, 49, 50, 51, 16, 52, 25, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68] could dramatically alter the free-streaming property of neutrinos on scales probed by the CMB. The key questions to ask then are whether current and upcoming data can pinpoint this new physics, and whether the latter can bias current and future constraints on the sum of neutrino masses.

Several works [14, 44, 48, 50, 58, 56, 69, 70] have shown that current cosmological data allow for some flexibility in the underlying theoretical assumption of neutrinos interacting solely via the weak interaction. Previous studies have usually taken a very phenomenological approach to constrain models of neutrino properties, relying either on the effective fluid parameters and (the rest-frame sound-speed and viscosity parameter, respectively) [24, 71, 72, 73, 74, 75, 76, 77, 78, 8, 79], or on an effective redshift [25, 56] meant to characterize the kinetic decoupling or re-coupling epoch of neutrinos, depending on the model considered. While the former approach based on the effective fluid parameters has shown consistency with the standard neutrino cosmology, it is in general difficult to interpret these results in terms of new interactions in the neutrino sector [58, 63]. The analyses using an effective re-coupling or decoupling redshift can capture more general neutrino interaction scenarios, but still leave out the details of the neutrino visibility function which could contain important clues on the nature of the yet-unseen neutrino interaction.

In this work, we compute updated cosmological constraints on the free-streaming nature of neutrinos in the early Universe. We focus on models where a new Fermi-like four-fermion interaction leads to a neutrino self-interaction rate which suppresses free streaming at early times (here is a Fermi-like constant controlling the strength of the new interaction). Such models could occur for instance if SM neutrinos couple to a new scalar or vector mediator with mass above keV [35, 40, 51, 52, 54, 57], or, perhaps more plausibly, if SM neutrinos mix with a sterile neutrino which is itself coupled to a new vector or scalar mediator [62, 65, 68]. We do not consider here models where neutrinos (either active or sterile) couple to a massless or nearly massless particle such as a majoron [36, 37, 40, 39, 14, 44, 45, 47, 16, 64, 66] since these typically lead to a neutrino self-interaction rate scaling as at energies relevant to the CMB. For simplicity, we also focus exclusively on massless neutrinos since the presence of neutrino mass significantly complicates the computation of the self-interaction collision integral. In an upcoming work [80], we show that the results presented here are robust to the inclusion of neutrino masses.

Cosmological constraints on neutrino interaction models with were placed in refs. [58, 56]. There, it was found that the posterior distribution for the interaction strength controlling neutrino self interaction was bimodal. While one of the modes shows broad consistency with the standard evolution of cosmological neutrinos, the second mode of the posterior corresponds to a radically different evolution of cosmological neutrinos in which neutrinos are strongly self interacting until close to the epoch of matter-radiation equality. Assessing the statistical significance of this somewhat surprising “interacting neutrino mode” is complicated by the fact that it depends on the choice of the prior probability assigned to each mode [58]. Here, we adopt a conservative choice of prior and use the latest CMB, baryon acoustic oscillation (BAO), and Hubble parameter data to test for the presence of new physics in the neutrino sector.

Of course, CMB and other cosmological data are not the only probe of neutrino physics beyond the SM. A summary of constraints on “secret” neutrino interaction has been presented in ref. [61]. For a fairly low mass mediator ( MeV), SN 1987A [40], Big Bang nucleosynthesis [81], the CMB, and the detection of PeV neutrinos at IceCube [61, 65] provide the strongest constraints, with the latter bound having the potential of being the most stringent. There are also flavor-dependent bounds based on measurements of lepton and meson decays [82], although these could potentially be evaded by an appropriate choice of the flavor structure of the matrix coupling the new interaction to SM neutrinos. Other limits [83, 35, 84] coming from Z-boson decay do not apply at the energy scale probed by the CMB. We note that elastic collisions caused by the new interaction do not change the time it takes for neutrinos to escape supernovae [85], although they could lead to new interesting phenomena (see e.g. refs. [86, 87, 88, 89, 90]). Finally, supernova cooling puts bounds on the coupling of majorons to SM neutrinos [91, 92, 93], but the applicability of these to the models considered here is uncertain.

We begin in section 2 by discussing the generic properties of the neutrino interaction model that we consider here and by presenting the key equations describing the cosmological evolution of self-interacting neutrinos. The data used in our analysis is introduced in section 3, while our results are given in section 4. We discuss the implications of our findings in section 5. We finally perform a CMB stage-IV Fisher Matrix forecast in section 6 before concluding in section 7.

## 2 Cosmological evolution with self-interacting neutrinos

### 2.1 Interaction model

We consider a phenomenological model in which SM neutrinos are coupled to a new massive scalar field with mass via an interaction term of the form , where is a dimensionless coupling constant. As discussed in ref. [61], this interaction model should not be taken at face value but instead serves as a rough guide to delimit the parameter space relevant to “secret” neutrino interactions. Within this phenomenological model, we focus on the low-energy effective theory where the neutrino temperature is significantly less than , allowing us to integrate out the massive mediator and treat the interaction as a four-neutrino vertex with a dimensional coupling constant given by1

 Gν=y2νM2ϕ. (1)

As we focus on the regime where the neutrino temperature is significantly less than the mediator mass, we can safely ignore the presence of the mediator in the thermal bath. Of course, the annihilation of the particles when would generally heat the neutrino bath above its standard value, leading to a non-zero value of the number of extra relativistic degrees of freedom, . However, since the specifics of this heating strongly depend on the details of the neutrino interaction model, we do not include it here and thus fix for the remainder of this paper. In an upcoming work [80], we show that our results are robust to this particular choice.

For our Fermi-like phenomenological model, the thermally-averaged neutrino self-interac-tion cross section scales as . Since the neutrino number density itself scales as after weak neutrino decoupling, we get that the neutrino self-interaction rate goes as . Any realistic model of neutrino interaction admitting an equivalent low-energy effective theory as our phenomenological model will lead to a similar interaction rate, up to a pre-factor depending on the details of the interaction. We thus define our neutrino self-interaction opacity as , where is the scale factor describing the expansion of the universe, and is a constant of order unity that is determined by the details of the neutrino interaction model [58]. We note that the overhead dot denotes a derivative with respect to conformal time (which is why there is a scale factor appearing in our definition of opacity), and that the opacity is by convention negative since it is the derivative of the neutrino optical depth. Since the neutrino opacity only depends on the product of and , we define a rescaled coupling constant

 Geff≡√ξGν, (2)

which allows us to write down the neutrino opacity as

 ˙τν=−aG2effT5ν. (3)

As we show below, the CMB is sensitive only to relatively large values of , and it is thus reasonable to assume that , where MeV is the standard Fermi constant. For the remainder of this work, we thus justifiably ignore the electroweak contribution to the neutrino opacity, as is traditionally done in standard CMB analyses.

### 2.2 Cosmological perturbations

We now turn our attention to the evolution of cosmological perturbations in the presence of self-interacting neutrinos. We adopt the notation of ref. [11] in the synchronous gauge. The details of the derivation of the neutrino Boltzmann hierarchy in the presence of self-interaction are presented in ref. \citepoldengott15. The presence of significant momentum transfer in a typical neutrino interaction renders the computation of the collision integrals rather tedious. However, since the neutrino opacity given in eq. \eqrefnu_opacity is a steep function of temperature, there is only a narrow window in time where the details of the neutrino interactions play a role. Indeed, at early times when , the neutrinos form a tightly-coupled fluid such that only the two lowest moments (corresponding to energy density and heat flux) of the Boltzmann hierarchy are nonzero2. On the other hand, after neutrino decoupling when , the collision term plays little role in the evolution of neutrino perturbations. Only near the peak of the neutrino visibility function (see e.g. figure 5 below) can the significant momentum transferred in a neutrino collision modify the neutrino phase-space distribution function. Even in this case, since we are only sensitive to the neutrino perturbations through their gravitational impact on CMB photons, it is unclear whether these distortions of the neutrino phase-space distribution would be observables. Moreover, the absence of an energy sink or source near the epoch of neutrino last-scattering all but guarantees that distortions of the neutrino distribution function would be small3.

While it would be interesting to compute whether the distortions to the neutrino phase-space distribution function can lead to observable effects in the CMB, we adopt here a simplified picture in which we assume that the neutrino distribution function exactly maintains its equilibrium thermal shape throughout the history of the Universe. This allows us to simplify the collision integrals, and write a momentum-integrated hierarchy of neutrino multipole moments. This Boltzmann hierarchy takes the form

 ˙δν=−43θν−23˙h, (4)
 ˙θν=k2(14δν−12Fν2), (5)
 ˙Fν,2=815θν−35kFν,3+415˙h+85˙η+α2˙τνFν,2, (6)
 ˙Fν,l=k2l+1[lFν,(l−1)−(l+1)Fν,(l+1)]+αl˙τνFν,l,   l≥3 (7)

where is the neutrino energy density perturbation, is the divergence of the neutrino heat flux, and are the metric perturbation in the synchronous gauge [11], is the th multipole moment of the neutrino Boltzmann hierarchy, is the magnitude of the comoving wave vector in Fourier space, and the are dimensionless coefficients of order unity that depend on the specifics of the neutrino interaction model. Noting that a change to the value of can be absorbed into the definition of , and that the higher order coefficients have a small impact on the CMB temperature and polarization spectra, we simply set for all in this work.

We solve equations (4) through (7) numerically along with perturbation equations for other standard cosmological components such as baryonic matter, photons, and cold dark matter using a version of the code CAMB modified for the purposes of this study \citepcamb12,camb00. Due to the tight coupling of the self-interacting neutrinos whenever , equations \eqrefl2 and \eqrefhiL are difficult to solve numerically at early times. In this regime, we adopt a tight-coupling approximation scheme [97] where we set and for . We switch off this approximation whenever the condition is satisfied. We have checked that this switch occurs early enough as to not bias our results. Finally, we truncate our neutrino Boltzmann hierarchy using the procedure outlined in ref. [11].

## 3 Data analysis methods

As the analysis performed in ref. [58] found evidence for a multimodal posterior distribution for the cosmological and neutrino parameters in the presence of self-interactions, we do not use the standard Markov Chain Monte Carlo method based on the Metropolis-Hastings (MH) algorithm to explore the parameter space. Indeed, the MH algorithm displays poor convergence properties for posteriors containing distinct islands of high probability separated by wide valleys of low probability, which is exactly the type of posterior that our self-interacting neutrino model produces, as we will discuss in section 4. To tackle this problem, we use the more robust nested-sampling algorithm [98] as implemented in the software package MultiNest \citepferoz08,feroz09,feroz13 to scan the parameter space and generate random samples of the posterior probability distribution function (PDF). As a by-product, nested sampling can also yield an estimate of the Bayesian evidence, which is useful for model comparison (see section 4.2 below).

We employ the cosmological parameter estimation package CosmoSIS \citepzuntz15 to wrap together the MultiNest inference code, the CAMB Boltzmann solver, and the different data likelihoods used in this work (see section 3.1). Our parameter space consists of the 6 standard CDM cosmological parameters, together with the base-10 logarithm of the neutrino self-interaction strength, . We use flat priors for the standard cosmological parameters as described in [103], as well as a flat prior on . We justify this choice of prior by noting that models with MeV have no impact on scales currently probed by the CMB, while models with MeV have their neutrino decoupling epoch occurring in the matter-dominated era when neutrinos play a subdominant role in the evolution of CMB and matter fluctuations. We perform nested sampling with 2000 live points and terminate the chains when the accuracy reaches a tolerance of in the logarithm of the Bayesian evidence. We have verified the stability of our results by both increasing the number of live points and lowering the lower bound of the prior on and subsequently checking that we obtain the same parameter constraints.

### 3.1 Data sets

As our main dataset, we use CMB data from the Planck 2015 release [103]. In particular, we include both the low- and high- temperature auto-correlation data, which we shall collectively refer to as the “TT” dataset4. To reduce the computational burden, we use the “lite” version of the high- likelihood which internally marginalizes over the nuisance parameters of the data reduction \citepplanck15nuisance. We also combine the CMB temperature data with the high- E-mode polarization auto-correlation and temperature cross-correlation data, a combination that we shall collectively refer to as “TT + Pol”5.

We also use data from the Sloan Digital Sky Survey (SDSS-III) Baryon Oscillation Spectroscopic Survey (BOSS) data release 12 [105]. In particular, we include their constraints on the volume averaged distance scale , obtained from the measurement of the BAO feature at three different redshifts , , and . We shall refer to this dataset as “BAO” in the rest of the paper.

Finally, we also use the direct measurement of the local Hubble Parameter from ref. [106], km/s/Mpc, which is obtained from a traditional distance ladder approach using Cepheids and supernovae. We refer to this likelihood as “” in what follows.

## 4 Results

### 4.1 Posterior distributions and parameter limits

We show in figure 1 the marginalized posterior distribution of the parameter for the different combinations of likelihoods considered in this work. As found in ref. [58], we find that the posterior is bimodal, with a first mode characterized by small values of the effective neutrino coupling constant and consistent with the standard CDM paradigm, and a second mode characterized by large values of the effective coupling constant and corresponding to a cosmology in which the onset of neutrino free-streaming is severely delayed compared to the standard scenario. We see that this “interacting neutrino mode” is present for every combination of data sets we used here, although its relative statistical weight compared to the standard CDM mode varies from one data set to the next. We discuss the statistical significance of each mode further in section 4.2 below. We note that since the two modes have very different widths, it is difficult to choose a single smoothing kernel for the whole range of the posterior shown, leading to an over-smoothing of the interacting neutrino mode. To illustrate this, we show in the inset of figure 1 the posterior of the interacting neutrino mode plotted using a smoothing kernel appropriate for the width of that mode, leading to changes of the smoothed PDF at the to level. We however emphasize that the choice of smoothing kernel is only relevant when plotting the posterior, and the quantitative parameter limits presented in this work do not depend on the choice of kernel.

In figures 2, 3, and 4 we present one and two-dimensional marginalized PDFs for the four cosmological parameters that are affected most by the introduction of neutrino self-interactions at early times: , , , and . In each of these plots we compare the effects of using Planck temperature and polarization data with simply using the temperature data (along with the other likelihoods that we list in the legend). The quantitative parameter confidence limits are given in table 1 for the CDM mode, and in table 2 for the interacting neutrino mode. For the CDM mode of the different posteriors we computed, the parameter constraints on , , , and are broadly consistent with those of the standard 6-parameter CDM paradigm [8], with parameter shifts of less than one sigma. As visible in figures 2, 3, and 4, larger values of within the CDM mode are correlated with smaller values of the scalar spectral index. This reflects the impact of neutrino self-interaction on the amplitude of the CMB temperature spectrum: the absence of neutrino free-streaming damping [15] is being compensated by a lower spectral index. However, the fact that the contours of the CDM mode close at finite values of implies that it is not always possible to compensate for the absence of neutrino free-streaming by adjusting the primordial spectrum of fluctuations. We return to this point in section 5.

For the CDM mode, the values of the optical depth to CMB last scattering for all data set combinations are slightly larger than the values presented in ref. [8], but still fully consistent given the error bars. These larger values of are compensated by somewhat higher amplitudes of the scalar fluctuations (compared to ref. [8]) in order to obtain the correct amplitude of the CMB temperature spectrum, which itself depends on the product . Within the CDM mode, all dataset combinations considered here lead to an upper bound on the strength of the neutrino self-interaction approximately given by (95% C.L., see table 1 for the exact limits). This can be translated to a constraint on the neutrino decoupling redshift . This limit is nearly identical to that presented in ref. [58], which is based on Planck 2013 CMB temperature data [107]. It may appear surprising that the limit on effective neutrino self-coupling does not significantly improve with the addition of the 2015 Planck temperature and polarization data. As we will discuss further in section 5, the strength of this limit is largely a function of the maximum CMB multipole included in the dataset. Since both Planck data releases extend to similar (), it is sensible that their constraining power of is similar.

Table 2 summarizes the parameter limits within the interacting neutrino mode. This mode is characterized by a higher value of the Hubble parameter km/s/Mpc, a lower value of the scalar spectral index , and a somewhat lower value of the scalar amplitude of fluctuation . The mean value of within this mode can be translated to a neutrino decoupling redshift , that is, still significantly before matter-radiation equality. The values of , , and are broadly consistent with those of the CDM mode. These results are similar to those presented in ref. [58]. The main difference with this latter work is that the preferred value of we obtain is slightly higher than theirs, and the 1-sigma uncertainty on the mean has shrunk by approximately a factor of 2. This is largely driven by the different calibration of the 2013 and 2015 Planck data [8] and by the addition of CMB polarization data, as can be seen in the inset of figure 1. Given that the -mode polarization data used in our analysis might contain systematics due to temperature-to-polarization leakage [8], we caution that this result should be interpreted with care. We expand on this point in section 5. While CMB polarization data tend to raise the amplitude of the interacting neutrino mode, the addition of the BOSS DR12 BAO data [105] tends to reduce its significance for all cases we considered. To understand why, note that the absence of a free-streaming phase shift [15, 17, 19, 20] for models in this corner of parameter space is compensated by adjusting the angular diameter distance to the last scattering surface (via a change in ). When a low redshift probe such as BAO is added, it provides an extra constraint in this angular-to-physical scale conversion, hence helping to break possible degeneracies. Our results suggest that improved low-redshift measurements of the BAO scale and of its phase [108] could help rule out the interacting neutrino mode as a viable possibility.

Interestingly, the higher value of the Hubble parameter preferred by the interacting neutrino mode helps alleviate the tension between CMB and local [106, 109] measurements of the present-day expansion rate. This is clearly illustrated in figure 4 where we show that the addition of the local Hubble parameter measurement to our analysis raises the statistical significance of the interacting neutrino mode. Comparing the values of given in table 2 with the latest measurement from ref. [109], km s Mpc, indicates that suppressing neutrino free-streaming reduces the tension but does not entirely remove it. Given that we have kept the effective number of neutrinos fixed to its standard value of and have not included neutrino masses in our analysis, it is possible that nonstandard neutrino properties could play a key role in reconciling the CMB and local measurements of the Hubble parameter once the full parameter space of neutrino properties is explored. We leave such exploration to future work [80].

### 4.2 Relative statistical weight of the two modes

We now turn our attention to the relative statistical importance of the interacting neutrino mode compared to the standard CDM mode. For the purpose of this analysis, we define the interacting neutrino mode as occupying the range , while the CDM mode consists of all models with . Upon comparison with figure 1, one can see that this is a reasonable splitting of the parameter space since the plane given by corresponds to a local minimum of the marginalized posterior. As a first statistical comparison, we compute the maximum-likelihood ratio of the interacting neutrino mode to the CDM mode, that is,

 RSIν≡max(LSIν)max(LΛCDM), (8)

where is the likelihood within the interacting neutrino mode, and is the likelihood within the CDM mode. We list the maximum-likelihood ratios for the different dataset combinations considered in this work in table 3. We observe that the interacting neutrino mode is always significantly subdominant in terms of its maximum likelihood, except when both polarization data and the local measurement are included. In this case, it can reach of the maximum likelihood of the CDM model. This is largely due to the fact that the interacting neutrino mode favors higher values of which are more consistent with distance ladder measurement of from ref. \citepriess11. It is also apparent from the ratios shown in table 3 that Planck polarization data tend to favor the interacting neutrino mode more strongly than the temperature data alone. Additionally, as mentioned above, BAO data tend to suppress the maximum-likelihood ratio of the interacting neutrino mode compared to the standard CDM mode.

A perhaps more robust measure of the relative statistical importance of the interacting neutrino mode can be obtained by computing its Bayesian evidence, and comparing it to that of the CDM mode. The Bayesian evidence is defined in terms of an integral over the entirety of the parameter space as

 Z≡∫ΩML(d|θ,M)Π(θ|M)dθ, (9)

where is the vector of parameters, is the volume of the model’s parameter space, is the probability of the data given the parameters and the model, or the likelihood function, and is the prior distribution placed on the model. Note that the Bayesian evidence is then just as dependent on the priors as it is on the model, . In this way we can define a Bayesian evidence for the CDM mode, as well as separately for the interacting neutrino mode based on the choice of prior given at the beginning of this subsection. We perform separate MultiNest [99, 100, 101] runs for each mode of the posterior, requesting a tolerance on the Bayesian evidence estimate provided by the code. These results are then used to compute the Bayes factor6, which is defined as

 BSIν≡ZSIνZΛCDM, (10)

where refers to the Bayesian evidence of model according to the definition in eq. \eqrefevidence. Essentially, is the ratio of the average likelihood within each mode. We list in table 3 the values of and their uncertainty for the different dataset combinations considered in this work. In all cases, the Bayes factor is much less than unity, implying that the data favor the standard CDM mode over the interacting neutrino cosmology.

We thus conclude, assuming massless neutrinos with and a flat prior on , that current data show no statistical evidence for the presence of strong self-interaction in the neutrino sector. While it is intriguing that the likelihoods within the two distinct modes are somewhat comparable (especially when local measurements are included), the small volume of parameter space where the interacting neutrino mode provides a good fit to the data ultimately results in a Bayes factor favoring the CDM mode. We note that the maximum-likelihood ratio and the Bayes factor show slightly different trends when moving from one dataset to the next, but our conclusions are unaffected by this scatter as the numbers remain small nonetheless. While we do not consider here the impact of an alternate choice of priors, it is clear that using, for instance, a flat prior on itself could dramatically alter our conclusions. Arguably, our current choice of priors could be considered conservative since it does not pick a particular energy scale for the new neutrino physics. As long as no external information is known about new physics in the neutrino sector, we strongly believe that a flat, non-informative prior on is the most reasonable choice. Of course, it would be interesting to determine how the statistical significance of the interacting neutrino mode is affected once neutrino masses are included and is allowed to deviate from its standard value.

## 5 Discussion

### 5.1 Bimodality

One of the key questions that our results raise is why we obtain a bimodal distribution. On the one hand, the existence of a mode consistent with the standard CDM cosmology in which neutrinos begin free-streaming before the onset of Big Bang nucleosynthesis is not too surprising. That this mode spans values of the effective neutrino coupling constant that are more than seven orders of magnitude above the standard Fermi constant is simply a reflection that the scales probed by the current CMB data are insensitive to the onset of neutrino free-streaming if it happens early enough. Indeed, for MeV, neutrino decoupling occurs before the Fourier modes probed by the Planck data enter the causal horizon, implying that they are unaffected by the new neutrino interactions and receive the standard phase and amplitude shift associated with neutrino free-streaming. We illustrate this in figure 5 where we show the neutrino visibility function as a function of conformal time.7 The gray band shows the approximate time interval in which the multipoles enter the causal horizon. This multipole range corresponds to scales encompassing all well-measured CMB temperature peaks except for the first one. We see that the visibility functions of models with MeV have no overlap with the time interval at which the modes probed by the current Planck data are entering the horizon. This explains why the posterior shown in figure 1 flattens out for this range of neutrino self-interacting strength.

The sharp suppression of the CDM mode of the posterior distribution for MeV indicate that these values delay neutrino free-streaming long enough for the length scales probed by the CMB damping tail to enter the horizon. This is supported by figure 5 where the red dashed line shows the models corresponding to the 95% upper limit of the CDM mode (MeV), whose visibility function has significant overlap with the modes probed by the CMB. This indicates that the allowed upper limit on within the CDM mode strongly depends on the highest multipole probed by the data since higher are capable of probing an earlier onset of neutrino free streaming and thus smaller values of . It is thus not surprising that our constraint on the CDM mode is similar to that from ref. [58] since the value of between the Planck 2013 and 2015 data release did not appreciably change.

The deep trough in the posterior distribution for indicate that CMB data strongly disfavor neutrino decoupling occurring while the modes corresponding to the Silk damping tail are entering the horizon. Indeed, models with in this range have a neutrino visibility function peaking within the gray band of figure 5. The dot-dashed blue line shows the neutrino visibility function for the best fit model within the interacting neutrino mode. We observe that this visibility function peaks right as the multipoles corresponding to the first CMB temperature peak (green shaded region) begin to enter the Hubble horizon. In this case, none of the CMB temperature peaks in the range receives the phase and amplitude shift usually associated with neutrino free streaming, hence requiring the other cosmological parameters, notably , , and , to absorb the resulting difference in the temperature spectrum (see figure 7). From the perspective of CMB polarization, the visibility function of the best fit interacting neutrino model has a maximum near the epoch when the second peak of the -mode polarization spectrum at is entering the horizon.

To understand the impact of this late neutrino decoupling, it is instructive to look at the CMB temperature and -mode polarization spectra on the relevant scales (), as shown in figure 6. There, we observe that the interacting neutrino mode predicts a slightly lower amplitude for the first peak of the temperature spectrum compared to the standard CDM cosmology, while displaying more power than the standard paradigm at low multipoles. This indicates that the degeneracy with the cosmological parameters , , and is not exact and that the interacting neutrino mode is a compromise between having a large enough amplitude of the first temperature peak while not overproducing power at very low multipoles. This actually hints at why the interacting neutrino mode is statistically subdominant compared to the CDM mode: it leads to a worse fit of the low- temperature data, hence suppressing its overall likelihood. This also indicates why the interacting neutrino mode is only viable for such a small range of : larger values of require even lower values of and and even higher values of in order to fit the damping tail and first peak, which invariably leads to a CMB temperature spectrum with too much power at low multipoles.

The picture that emerges is that self-interacting neutrinos are only viable if they either decouple before the bulk of the modes probed by the CMB enter the horizon, or if they decouple within a narrow window centered around . In all cases, neutrinos must begin to free stream before matter-radiation equality (even within the interacting neutrino mode, ). For this latter case, the large cosmic variance-dominated error bars at low multipoles allow for an approximate multi-parameter degeneracy to correct for the impact of self-interacting neutrino on the CMB damping tail, at the price of having a slightly lower amplitude of the first temperature peak. It is however remarkable that the parameter changes required to make the interacting neutrino mode fit the Planck temperature data reasonably well also lead to an -mode polarization spectra that is in very good agreement with the data. This is illustrated in the right panel of figure 6 (see also next subsection).

### 5.2 Broad structure of temperature and polarization spectra

We compare in figure 7 the CMB temperature anisotropy spectrum of the best fit models for both the CDM mode and the interacting neutrino mode, along with two reference models with MeV and MeV for which the other cosmological parameters are kept fixed at their Planck 2015 values. The lower panel illustrates the fractional residuals between these four models and the Planck 2015 [8] best fit model from their “TT, TE, EE + lowP” data combination. The first thing to notice is that the residuals of the best fit models of the two modes show a similar oscillatory pattern, but with the interacting neutrino mode displaying more power at large multipoles, where the error bars are relatively big. Over the range , the best fit interacting neutrino cosmology does not deviate from the Planck 2015 model by more than 2%, which is somewhat remarkable given how radically different this model is compared to the standard CDM paradigm. This significant difference can be understood by examining the two other plotted models (MeV and MeV) for which the large amplitude and phase shift are readily noticeable. The existence of a degeneracy between the neutrino self-interaction strength and the parameter combination , , and that nearly compensates for these large changes to the CMB temperature spectrum is remarkable. It also highlights the importance of correctly modeling neutrino physics when performing CMB analyses since this degeneracy is missed when the popular parametrization is used to study deviation from standard free-streaming neutrinos.

Figures 8 and 9 display the and CMB spectra, respectively, for the same four models. Again, we observe that the best fit model within the interacting neutrino mode (green dashed line) has residuals similar to that of the best fit model within the CDM mode (red solid line). This indicates that the multi-parameter degeneracy responsible for compensating the impact of neutrino self-interaction on the CMB temperature spectrum also leads to an -mode polarization spectrum that closely matches the Planck data. As noted in table 2, including polarization data does yield a slightly higher value of the neutrino self-interaction strength, but this shift is entirely consistent with the value inferred from temperature-only data at less than 1-. We caution though that the presence of potential systematics within the Planck polarization data could affect our results, and they should therefore be interpreted with care. However, since our results with and without the Planck polarization data are largely consistent with one another, it is very unlikely that these systematics would dramatically change our conclusions. Indeed, we emphasize that our upper limit on within the CDM mode and the existence of the somewhat subdominant interacting neutrino mode are both robust to the possible presence of systematics in the data.

## 6 CMB Stage-IV forecast

We end our analysis by performing a Fisher matrix forecast for the expected sensitivity of Stage-IV CMB experiments (CMB-S4) [111]. The Fisher formalism assumes a Gaussian PDF for the parameters. Since our posterior is bimodal and thus cannot be accurately approximated by a single multivariate Gaussian, we perform separate forecasts for each mode of the posterior. We emphasize that our Fisher analysis cannot be used to determine whether future data can rule out the interacting neutrino mode. Instead, it is used to provide an estimate of how much future data will tighten the parameter constraints within each mode. Following ref. \citepwufisher14, we assume a Gaussian likelihood for the parameter vector given some data vector

 L(θ|d)∝1√detC(θ)exp(−12d†[C(θ)]−1d), (11)

where is the covariance matrix of the data vector. The Fisher matrix itself is then built by taking second-order partial derivatives of the likelihood, or its curvature, at the fiducial values of the parameters

 Fij=⟨∂2logL∂θi∂θj∣∣∣θ=θ0⟩, (12)

and the uncertainty of parameter is given by

 σi≡σ(θi)=√(F−1)ii. (13)

For CMB experiments, we choose our data vector to be , with being the spherical harmonic coefficients of the temperature field and the E-mode polarization field, respectively; consequently its covariance matrix for each multipole is given by

 Cl(θ)≡(CTTl+NTTlCTElCTElCEEl+NEEl),

where the Gaussian noise is defined as

 NXXl=Δ2Xexp(l(l+1)θ2FWHM8log2), (14)

with being the experiment pixel noise level, and being the full width at half-maximum beam size in radians. Combining the above eqns. \eqreflikelihood and \eqreffisher1, we find

 Fij=lmax∑l=lmin2l+12fskyTr(C−1l(θ)∂Cl∂θiC−1l(θ)∂Cl∂θj), (15)

where is the fraction of covered sky. In our forecast, we take partial derivatives directly over rather than its logarithm as this is better suited to the linear approximation of the Fisher formalism. The fiducial parameters for each mode are chosen as the mean values from our inference analysis using the “ TT + Pol” data combination (see tables 1 and 2). The only exception is the value of in the CDM mode, which we take to be MeV. We assume the following specifications for CMB-S4: , , K-arcmin, , and arcmin. We list the forecasted 1- error bars for each mode in table 4. Our forecast suggests that CMB-S4 could significantly improve the upper bound on within the CDM mode, with a projected limit of at 95% confidence level. We note that this improvement is largely driven by the larger and the smaller error bars of the projected CMB-S4 dataset.

## 7 Conclusion

In this paper, we have presented updated cosmological bounds on the free-streaming nature of neutrinos in the early Universe, focusing on models where neutrino self-scattering is well-described by a Fermi-like four-fermion interaction. Using the Planck 2015 CMB temperature and polarization data, the BAO measurement from BOSS DR12, and the local measurement of a Hubble parameter, we have determined that the onset of neutrino free streaming must occur either before any of the Fourier modes probed by the CMB enter the causal horizon () or within a narrow window centered around redshift . This latter case (i.e. the interacting neutrino mode) is possible due to a multi-parameter degeneracy between the neutrino self-interaction strength and the combination of parameters . This degeneracy is however not exact, and the interacting neutrino mode ends up predicting an excess of power in the temperature spectrum at low multipoles compared to what is observed by Planck, leading to the mode being statistically subdominant compared to the more standard CDM mode. In addition, the BAO measurement considered in this work also tends to disfavor the interacting neutrino mode, while the Planck -mode polarization data has the opposite effect and increase the statistical significance of this secondary mode. Interestingly, the interacting neutrino cosmology can naturally accommodate a higher value of the Hubble parameter, and can thus help reconcile to some extent the tension between the local measurements based on a distance ladder with those derived from CMB data.

The surest way to improve on the constraints presented here is to extend the measurement of the CMB damping tail to higher value of . While such measurements already exist (e.g. [113, 114]) and are bound to improve in the near future, the presence of significant foregrounds at in the temperature spectrum might partially nullify the gain of going to smaller scales. However, it is likely that better CMB polarization measurements on small scales could significantly improve the present constraints, as shown with our Fisher forecast for CMB-S4. In any case, probing larger multipoles within the CMB damping tail and reducing the measurement error bars could help rule out the interacting neutrino mode since the latter tends to display more power on small scales compared to its CDM counterpart (see figures 7 and 9). Including additional low-redshift probes could also in principle help further constraining the interacting neutrino mode. We note however that care must be taken when using these types of probes since the late neutrino decoupling within the interacting neutrino mode modifies the matter power spectrum [80] on scales where BAO, galaxy clustering, and weak lensing all have significant sensitivity. Since the analyses performed on the raw data often assumed a CDM background cosmology, (see e.g. ref. [105]) the constraints on derived parameters such as the BAO scale or extracted from these analyses are difficult to interpret within the interacting neutrino mode context. On the bright side, the continued detection of ultra-high energy neutrinos at IceCube could eventually ruled out this peculiar mode of the posterior while significantly improving the bounds within the CDM mode [61, 65].

In this work, we have focused for simplicity on massless neutrinos with a standard thermal history (that is, we fixed ). In reality, we know that neutrinos have mass and that the introduction of non-standard neutrino interaction is likely to change their thermal history and lead to a different value for the effective number of neutrinos. On the one hand, we expect our constraints on the largest possible value of within the CDM mode to be robust to the introduction of neutrino mass since neutrinos are still ultra-relativistic at the onset of free streaming in this case. On the other hand, we expect neutrino masses to have a larger effect on the interacting neutrino mode since neutrino decoupling within this model occurs close to the epoch at which neutrinos become nonrelativistic. It would be interesting to see how the statistical significance of the interacting neutrino mode changes once the neutrino mass and are allowed to vary, especially when local Hubble parameter measurements are included. We leave such consideration to future work [80].

Note added at Publication: We would like to note that after the submission of this work, the preprint [115] appeared on arXiv, which performs a very similar analysis to this work. The main differences between [115] and this work are their exact treatment of the details of the momentum dependence of the neutrino interaction under study here and how varying changes the derived results. While this work employs a much more detailed analysis of this interaction than our work, the results that are obtained are in good agreement with ours. As is pointed out in [115], this provides a good justification for the thermal approximation of the collision term of the Boltzmann hierarchy that we employ here.

## Acknowledgements

The work of this paper could not have been completed without the advice and review provided by members of the Knox research group: Dr. Marius Millea and Brigid Mulroe. We would especially like to thank Dr. Joe Zuntz, leader of the CosmoSIS team, for his excellent and quick feedback on technical issues. We thank Shadab Alam for his help putting together the likelihood updates for the most recent BOSS DR12 data. We also thank the McWilliams Center for Cosmology, for the use of their ‘Coma Cluster’ computer cluster. F.-Y. C.-R. wishes to thank Prateek Agrawal for useful discussions. The work of L. L. was supported in part by NSF grant PHY-1263201 as part of the Research Experience for Undergraduates (REU) program. L. L.wishes to thank the University of California, Davis, for hospitality while part of this work was completed. F.-Y. C.-R. acknowledges the support of the National Aeronautical and Space Administration ATP grant NNX16AI12G at Harvard University.

### Footnotes

1. We note that our definition of is slightly different than that used in ref. [61].
2. Here, is the conformal Hubble parameter.
3. We note for instance that for the SM neutrino decoupling near the epoch of nucleosynthesis, the bulk of the neutrino spectral distortions comes from the energy injected by annihilation [94].
4. For clarity, we use the commander_rc2_v1.1_l2_29_B and plik_lite_v18_TT likelihoods at low- and high-, respectively.
5. In addition to the low- temperature likelihood, the exact likelihood used in this case is plik_lite_v18_TTTEEE.
6. Note that we are not using the Bayes factor for model comparison (as it is commonly used), but rather to compare the statistical weight of two modes of the same model.
7. Much like the better-known CMB visibility function, the neutrino visibility function is a probability density function for the time at which neutrinos last scatter.

### References

1. S. Vagnozzi, E. Giusarma, O. Mena, K. Freese, M. Gerbino, S. Ho et al., Unveiling secrets with cosmological data: neutrino masses and mass hierarchy, 1701.08172.
2. E. Di Valentino, S. Gariazzo, M. Gerbino, E. Giusarma and O. Mena, Dark Radiation and Inflationary Freedom after Planck 2015, Phys. Rev. D93 (2016) 083523, [1601.07557].
3. E. Giusarma, M. Gerbino, O. Mena, S. Vagnozzi, S. Ho and K. Freese, Improvement of cosmological neutrino mass bounds, Phys. Rev. D94 (2016) 083522, [1605.04320].
4. E. Di Valentino, E. Giusarma, M. Lattanzi, O. Mena, A. Melchiorri and J. Silk, Cosmological Axion and neutrino mass constraints from Planck 2015 temperature and polarization data, Phys. Lett. B752 (2016) 182–185, [1507.08665].
5. E. Di Valentino, E. Giusarma, O. Mena, A. Melchiorri and J. Silk, Cosmological limits on neutrino unknowns versus low redshift priors, Phys. Rev. D93 (2016) 083527, [1511.00975].
6. R. de Putter, O. Mena, E. Giusarma, S. Ho et al., New Neutrino Mass Bounds from Sloan Digital Sky Survey III Data Release 8 Photometric Luminous Galaxies, Astrophys. J. 761 (2012) 12.
7. E. Giusarma, R. de Putter and O. Mena, Testing standard and non-standard neutrino physics with cosmological data, Phys.Rev. D87 (2013) 043515.
8. Planck Collaboration, Ade, P. A. R., Aghanim, N., Arnaud, M., Ashdown, M., Aumont, J. et al., Planck 2015 results xiii. cosmological parameters, AA 594 (2016) A13.
9. Particle Data Group collaboration, K. Olive et al., Review of Particle Physics, Chin. Phys. C38 (2014) 090001.
10. R. H. Cyburt, B. D. Fields, K. A. Olive and T.-H. Yeh, Big Bang Nucleosynthesis: 2015, Rev. Mod. Phys. 88 (2016) 015004, [1505.01076].
11. C.-P. Ma and E. Bertschinger, Cosmological perturbation theory in the synchronous and conformal Newtonian gauges, Astrophys.J. 455 (1995) 7–25.
12. J. M. Stewart, Perturbations in an Expanding Universe of Free Particles, Astrophys. J. 176 (Sept., 1972) 323.
13. P. J. E. Peebles, The Role of Neutrinos in the Evolution of Primeval Adiabatic Perturbations, Astrophys. J. 180 (Feb., 1973) 1–6.
14. Z. Chacko, L. J. Hall, T. Okui and S. J. Oliver, CMB signals of neutrino mass generation, Phys. Rev. D 70 (2004) 085008.
15. S. Bashinsky and U. Seljak, Neutrino perturbations in CMB anisotropy and matter clustering, Phys.Rev. D69 (2004) 083002.
16. A. Friedland, K. M. Zurek and S. Bashinsky, Constraining Models of Neutrino Mass and Neutrino Interactions with the Planck Satellite, 0704.3271.
17. Z. Hou, R. Keisler, L. Knox, M. Millea and C. Reichardt, How massless neutrinos affect the cosmic microwave background damping tail, Phys. Rev. D 87 (Apr., 2013) 083008, [1104.2333].
18. S. Bashinsky, Mapping Cosmological Observables to the Dark Kinetics, Submitted to: Astrophys. J. (2007) , [0707.0692].
19. B. Follin, L. Knox, M. Millea and Z. Pan, First Detection of the Acoustic Oscillation Phase Shift Expected from the Cosmic Neutrino Background, Phys. Rev. Lett. 115 (2015) 091301, [1503.07863].
20. D. Baumann, D. Green, J. Meyers and B. Wallisch, Phases of New Physics in the CMB, JCAP 1601 (2016) 007, [1508.06342].
21. Z. Pan, L. Knox, B. Mulroe and A. Narimani, Cosmic Microwave Background Acoustic Peak Locations, Mon. Not. Roy. Astron. Soc. 459 (2016) 2513–2524, [1603.03091].
22. C. Brust, Y. Cui and K. Sigurdson, Cosmological Constraints on Interacting Light Particles, ArXiv e-prints (Mar., 2017) , [1703.10732].
23. M. Kaplinghat, L. Knox and Y.-S. Song, Determining neutrino mass from the CMB alone, Phys. Rev. Lett. 91 (2003) 241301, [astro-ph/0303344].
24. R. Trotta and A. Melchiorri, Indication for primordial anisotropies in the neutrino background from WMAP and SDSS, Phys. Rev. Lett. 95 (2005) 011305.
25. A. Basboll, O. E. Bjaelde, S. Hannestad and G. G. Raffelt, Are cosmological neutrinos free-streaming?, Phys. Rev. D 79 (2009) 043512.
26. M. Archidiacono, E. Giusarma, S. Hannestad and O. Mena, Cosmic dark radiation and neutrinos, Adv. High Energy Phys. 2013 (2013) 191047, [1307.0637].
27. V. N. Gribov and B. Pontecorvo, Neutrino astronomy and lepton charge, Phys. Lett. B28 (1969) 493.
28. B. Pontecorvo, Neutrino Experiments and the Problem of Conservation of Leptonic Charge, Sov. Phys. JETP 26 (1968) 984–988.
29. SNO Collaboration collaboration, Q. Ahmad et al., Measurement of the rate of interactions produced by solar neutrinos at the Sudbury Neutrino Observatory, Phys. Rev. Lett. 87 (2001) 071301.
30. M. Gell-Mann, P. Ramond and R. Slansky, Complex Spinors and Unified Theories, Conf. Proc. C790927 (1979) 315–321, [1306.4669].
31. T. Yanagida, Horizontal Symmetry and Masses of Neutrinos, Prog. Theor. Phys. 64 (1980) 1103.
32. R. N. Mohapatra and G. Senjanović, Neutrino mass and spontaneous parity nonconservation, Phys. Rev. Lett. 44 (Apr, 1980) 912–915.
33. J. Schechter and J. W. F. Valle, Neutrino masses in su(2) u(1) theories, Phys. Rev. D 22 (Nov, 1980) 2227–2235.
34. Z. Bialynicka-Birula, Do Neutrinos Interact between Themselves?, Nuovo Cim. 33 (1964) 1484–1487.
35. D. Y. Bardin, S. M. Bilenky and B. Pontecorvo, On the nu - nu interaction, Phys.Lett. B32 (1970) 121–124.
36. G. Gelmini and M. Roncadelli, Left-Handed Neutrino Mass Scale and Spontaneously Broken Lepton Number, Phys.Lett. B99 (1981) 411.
37. Y. Chikashige, R. N. Mohapatra and R. Peccei, Spontaneously Broken Lepton Number and Cosmological Constraints on the Neutrino Mass Spectrum, Phys.Rev.Lett. 45 (1980) 1926.
38. V. D. Barger, W.-Y. Keung and S. Pakvasa, Majoron Emission by Neutrinos, Phys.Rev. D25 (1982) 907.
39. G. Raffelt and J. Silk, Light Neutrinos as Cold Dark Matter, Phys.Lett. B192 (1987) 65.
40. E. W. Kolb and M. S. Turner, Supernova SN 1987a and the Secret Interactions of Neutrinos, Phys.Rev. D36 (1987) 2895.
41. R. V. Konoplich and M. Yu. Khlopov, Constraints on triplet Majoron model due to observations of neutrinos from stellar collapse, Sov. J. Nucl. Phys. 47 (1988) 565–566.
42. A. V. Berkov, Yu. P. Nikitin, A. L. Sudarikov and M. Yu. Khlopov, Possible Manifestations of Anomalous 4 Neutrino Interaction., Sov. J. Nucl. Phys. 48 (1988) 497–501.
43. K. M. Belotsky, A. L. Sudarikov and M. Yu. Khlopov, Constraint on anomalous 4nu interaction, Phys. Atom. Nucl. 64 (2001) 1637–1642.
44. J. F. Beacom, N. F. Bell and S. Dodelson, Neutrinoless universe, Phys. Rev. Lett. 93 (2004) 121302.
45. S. Hannestad, Structure formation with strongly interacting neutrinos - Implications for the cosmological neutrino mass bound, Journal of Cosmology and Astroparticle Physics 0502 (2005) 011.
46. Z. Chacko, L. J. Hall, S. J. Oliver and M. Perelstein, Late time neutrino masses, the LSND experiment and the cosmic microwave background, Phys. Rev. Lett. 94 (2005) 111801, [hep-ph/0405067].
47. S. Hannestad and G. Raffelt, Constraining invisible neutrino decays with the cosmic microwave background, Phys.Rev. D72 (2005) 103514, [hep-ph/0509278].
48. N. F. Bell, E. Pierpaoli and K. Sigurdson, Cosmological signatures of interacting neutrinos, Phys.Rev. D73 (2006) 063523.
49. R. F. Sawyer, Bulk viscosity of a gas of neutrinos and coupled scalar particles, in the era of recombination, Phys.Rev. D74 (2006) 043527.
50. M. Cirelli and A. Strumia, Cosmology of neutrinos and extra light particles after WMAP3, JCAP 0612 (2006) 013, [astro-ph/0607086].
51. G. Mangano, A. Melchiorri, P. Serra, A. Cooray and M. Kamionkowski, Cosmological bounds on dark matter-neutrino interactions, Phys. Rev. D 74 (2006) 043517.
52. D. Hooper, Detecting MeV gauge bosons with high-energy neutrino telescopes, Phys. Rev. D 75 (2007) 123001, [hep-ph/0701194].
53. P. Serra, F. Zalamea, A. Cooray, G. Mangano and A. Melchiorri, Constraints on neutrino – dark matter interactions from cosmic microwave background and large scale structure data, Phys.Rev. D81 (2010) 043507.
54. L. G. van den Aarssen, T. Bringmann and C. Pfrommer, Is dark matter with long-range interactions a solution to all small-scale problems of CDM cosmology?, Phys. Rev. Lett. 109 (2012) 231301, [1205.5809].
55. K. S. Jeong and F. Takahashi, Self-interacting Dark Radiation, Phys.  Lett.  B 725, 134 (2013) , [1305.6521].
56. M. Archidiacono and S. Hannestad, Updated constraints on non-standard neutrino interactions from Planck, JCAP 1407 (2014) 046, [1311.3873].
57. R. Laha, B. Dasgupta and J. F. Beacom, Constraints on New Neutrino Interactions via Light Abelian Vector Bosons, Phys.Rev. D89 (2014) 093025, [1304.3460].
58. F. Cyr-Racine and K. Sigurdson, Limits on Neutrino-Neutrino Scattering in the Early Universe, Phys. Rev. D 109 (Dec., 2014) 1071–1085.
59. M. Archidiacono, S. Hannestad, R. S. Hansen and T. Tram, Cosmology with self-interacting sterile neutrinos and dark matter - A pseudoscalar model, Phys. Rev. D 91 (2015) 065021, [1404.5915].
60. K. Ioka and K. Murase, IceCube PeV–EeV neutrinos and secret interactions of neutrinos, PTEP 2014 (2014) 061E01, [1404.2279].
61. K. C. Y. Ng and J. F. Beacom, Cosmic neutrino cascades from secret neutrino interactions, Phys.Rev. D90 (2014) 065035, [1404.2288].
62. J. F. Cherry, A. Friedland and I. M. Shoemaker, Neutrino Portal Dark Matter: From Dwarf Galaxies to IceCube, 1411.1071.
63. I. M. Oldengott, C. Rampf and Y. Y. Y. Wong, Boltzmann hierarchy for interacting neutrinos I: formalism, Journal of Cosmology and Astroparticle Physics 4 (Apr., 2015) 16, [1409.1577].
64. M. Archidiacono, S. Hannestad, R. S. Hansen and T. Tram, Sterile neutrinos with pseudoscalar self-interactions and cosmology, Phys. Rev. D93 (2016) 045004, [1508.02504].
65. J. F. Cherry, A. Friedland and I. M. Shoemaker, Short-baseline neutrino oscillations, Planck, and IceCube, 1605.06506.
66. M. Archidiacono, S. Gariazzo, C. Giunti, S. Hannestad, R. Hansen, M. Laveder et al., Pseudoscalar—sterile neutrino interactions: reconciling the cosmos with neutrino oscillations, JCAP 1608 (2016) 067, [1606.07673].
67. G. Dvali and L. Funcke, Small neutrino masses from gravitational Î¸-term, Phys. Rev. D93 (2016) 113002, [1602.03191].
68. F. Capozzi, I. M. Shoemaker and L. Vecchi, The Dark Side of MSW: Solar Neutrinos as a Probe of Dark Matter-Neutrino Interactions, 1702.08464.
69. F. Forastieri, M. Lattanzi and P. Natoli, Constraints on secret neutrino interactions after Planck, JCAP 1507 (2015) 014, [1504.04999].
70. F. Forastieri, M. Lattanzi, G. Mangano, A. Mirizzi, P. Natoli and N. Saviano, Cosmic microwave background constraints on secret interactions among sterile neutrinos, 1704.00626.
71. A. Melchiorri and P. Serra, Anisotropies in the neutrino background: An Update, Phys. Rev. D 74 (2006) 127301.
72. F. De Bernardis, L. Pagano, P. Serra, A. Melchiorri et al., Anisotropies in the Cosmic Neutrino Background after WMAP 5-year Data, Journal of Cosmology and Astroparticle Physics 0806 (2008) 013.
73. T. L. Smith, S. Das and O. Zahn, Constraints on neutrino and dark radiation interactions using cosmological observations, Phys.Rev. D85 (2012) 023001.
74. M. Archidiacono, E. Calabrese and A. Melchiorri, The Case for Dark Radiation, Phys. Rev. D 84 (2011) 123008.
75. M. Archidiacono, E. Giusarma, A. Melchiorri and O. Mena, Dark Radiation in extended cosmological scenarios, Phys.Rev. D86 (2012) 043509.
76. M. Gerbino, E. Di Valentino and N. Said, Neutrino Anisotropies after Planck, Phys.Rev. D88 (2013) 063538, [1304.7400].
77. A. Melchiorri, M. Archidiacono, E. Calabrese and E. Menegoni, Constraints on neutrino physics from cosmology, Journal of Physics: Conference Series 485 (2014) 012014.
78. E. Sellentin and R. Durrer, Detecting the cosmological neutrino background in the CMB, Phys. Rev. D92 (2015) 063012, [1412.6427].
79. B. Audren, E. Bellini, A. J. Cuesta, S. G. A. Gontcho, J. Lesgourgues, V. Niro et al., Robustness of cosmic neutrino background detection in the cosmic microwave background, Journal of Cosmology and Astroparticle Physics 3 (Mar., 2015) 036, [1412.5948].
80. C. Kreisch, F.-Y. Cyr-Racine, O. Doré and K. Sigurdson, in prep., .
81. B. Ahlgren, T. Ohlsson and S. Zhou, Comment on “Is Dark Matter with Long-Range Interactions a Solution to All Small-Scale Problems of ? Cold Dark Matter Cosmology?”, Phys.Rev.Lett. 111 (2013) 199001, [1309.0991].
82. A. P. Lessa and O. L. G. Peres, Revising limits on neutrino-Majoron couplings, Phys.Rev. D75 (2007) 094001.
83. M. S. Bilenky, S. M. Bilenky and A. Santamaria, Invisible width of the Z boson and ’secret’ neutrino-neutrino interactions, Phys.Lett. B301 (1993) 287–291.
84. M. S. Bilenky and A. Santamaria, ’Secret’ neutrino interactions, hep-ph/9908272.
85. D. A. Dicus, S. Nussinov, P. B. Pal and V. L. Teplitz, IMPLICATIONS OF RELATIVISTIC GAS DYNAMICS FOR NEUTRINO-NEUTRINO CROSS-SECTIONS, Phys.Lett. B218 (1989) 84.
86. H. Davoudiasl and P. Huber, Probing the origins of neutrino mass with supernova data, Phys. Rev. Lett. 95 (Nov, 2005) 191302.
87. M. Sher and C. Triola, Astrophysical consequences of a neutrinophilic two-higgs-doublet model, Phys. Rev. D 83 (Jun, 2011) 117702.
88. P. Fayet, D. Hooper and G. Sigl, Constraints on light dark matter from core-collapse supernovae, Phys. Rev. Lett. 96 (Jun, 2006) 211302.
89. K. Choi and A. Santamaria, Majorons and supernova cooling, Phys. Rev. D 42 (Jul, 1990) 293–306.
90. M. Blennow, A. Mirizzi and P. D. Serpico, Nonstandard neutrino-neutrino refractive effects in dense neutrino gases, Phys.Rev. D78 (2008) 113004, [0810.2297].
91. M. Kachelriess, R. Tomas and J. Valle, Supernova bounds on Majoron emitting decays of light neutrinos, Phys.Rev. D62 (2000) 023004, [hep-ph/0001039].
92. Y. Farzan, Bounds on the coupling of the Majoron to light neutrinos from supernova cooling, Phys.Rev. D67 (2003) 073015, [hep-ph/0211375].
93. S. Zhou, Comment on astrophysical consequences of a neutrinophilic 2HDM, Phys.Rev. D84 (2011) 038701, [1106.3880].
94. G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti and P. D. Serpico, Relic neutrino decoupling including flavor oscillations, Nucl. Phys. B729 (2005) 221–234, [hep-ph/0506164].
95. C. Howlett, A. Lewis, A. Hall and A. Challinor, CMB power spectrum parameter degeneracies in the era of precision cosmology, Journal of Cosmology and Astroparticle Physics 4 (Apr., 2012) 27, [1201.3654].
96. A. Lewis, A. Challinor and A. Lasenby, Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models, The Astrophysical Journal 538 (Aug., 2000) 473–476, [astro-ph/9911177].
97. F.-Y. Cyr-Racine and K. Sigurdson, Photons and Baryons before Atoms: Improving the Tight-Coupling Approximation, Phys.Rev. D83 (2011) 103521.
98. J. Skilling, Nested sampling for general bayesian computation, Bayesian Anal. 1 (12, 2006) 833–859.
99. F. Feroz and M. P. Hobson, Multimodal nested sampling: an efficient and robust alternative to Markov Chain Monte Carlo methods for astronomical data analyses, Monthly Notices of the Royal Astronomical Society 384 (Feb., 2008) 449–463, [0704.3704].
100. F. Feroz, M. P. Hobson and M. Bridges, MULTINEST: an efficient and robust Bayesian inference tool for cosmology and particle physics, Monthly Notices of the Royal Astronomical Society 398 (Oct., 2009) 1601–1614, [0809.3437].
101. F. Feroz, M. P. Hobson, E. Cameron and A. N. Pettitt, Importance Nested Sampling and the MultiNest Algorithm, ArXiv e-prints (June, 2013) , [1306.2144].
102. J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Manzotti, S. Dodelson et al., CosmoSIS: Modular cosmological parameter estimation, Astronomy and Computing 12 (Sept., 2015) 45–59, [1409.3409].
103. Planck Collaboration, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi et al., Planck 2015 results. XI. CMB power spectra, likelihoods, and robustness of parameters, ArXiv e-prints (July, 2015) , [1507.02704].
104. Planck Collaboration, R. Adam, P. A. R. Ade, N. Aghanim, M. I. R. Alves, M. Arnaud et al., Planck 2015 results. X. Diffuse component separation: Foreground maps, ArXiv e-prints (Feb., 2015) , [1502.01588].
105. BOSS collaboration, S. Satpathy et al., BOSS DR12 combined galaxy sample: The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: On the measurement of growth rate using galaxy correlation functions, Submitted to: Mon. Not. Roy. Astron. Soc. (2016) , [1607.03148].
106. A. G. Riess, L. Macri, S. Casertano, H. Lampeitl, H. C. Ferguson, A. V. Filippenko et al., A 3% Solution: Determination of the Hubble Constant with the Hubble Space Telescope and Wide Field Camera 3, Astrophys. J. 730 (Apr., 2011) 119, [1103.2976].
107. Planck Collaboration, P. A. R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown et al., Planck 2013 results. XVI. Cosmological parameters, Astron. Astrophys. 571 (Nov., 2014) A16, [1303.5076].
108. D. Baumann, D. Green and M. Zaldarriaga, Phases of New Physics in the BAO Spectrum, 1703.00894.
109. A. G. Riess et al., A 2.4% Determination of the Local Value of the Hubble Constant, Astrophys. J. 826 (2016) 56, [1604.01424].
111. K. N. Abazajian, P. Adshead, Z. Ahmed, S. W. Allen, D. Alonso, K. S. Arnold et al., CMB-S4 Science Book, First Edition, ArXiv e-prints (Oct., 2016) , [1610.02743].
112. W. L. K. Wu, J. Errard, C. Dvorkin, C. L. Kuo, A. T. Lee, P. McDonald et al., A Guide to Designing Future Ground-based Cosmic Microwave Background Experiments, The Astrophysical Journal 788 (June, 2014) 138, [1402.4108].
113. K. T. Story, C. L. Reichardt, Z. Hou, R. Keisler, K. A. Aird, B. A. Benson et al., A Measurement of the Cosmic Microwave Background Damping Tail from the 2500-Square-Degree SPT-SZ Survey, Astrophys. J. 779 (Dec., 2013) 86, [1210.7231].
114. S. Das, T. Louis, M. R. Nolta, G. E. Addison, E. S. Battistelli et al., The Atacama Cosmology Telescope: temperature and gravitational lensing power spectrum measurements from three seasons of data, JCAP 1404 (2014) 014, [1301.1037].
115. I. M. Oldengott, T. Tram, C. Rampf and Y. Y. Y. Wong, Interacting neutrinos in cosmology: exact description and constraints, JCAP 1711 (2017) 027, [1706.02123].
137388