The Neutrino Puzzle: Anomalies, Interactions, and Cosmological Tensions

TL;DR

This paper explores how neutrino self-interactions delaying free-streaming can reconcile cosmological tensions in H_0 and σ_8, offering a novel model that fits data well despite Bayesian penalties.

Contribution

It introduces a neutrino self-interaction model with non-zero masses that delays free-streaming, addressing cosmological tensions and anomalies in neutrino physics.

Findings

Delaying neutrino free-streaming increases H_0 and decreases σ_8.

The model prefers N_eff ≈ 4, impacting particle physics.

A new localized feature appears in the matter power spectrum.

Abstract

New physics in the neutrino sector might be necessary to address anomalies between different neutrino oscillation experiments. Intriguingly, it also offers a possible solution to the discrepant cosmological measurements of $H_0$ and $\sigma_8$. We show here that delaying the onset of neutrino free-streaming until close to the epoch of matter-radiation equality can naturally accommodate a larger value for the Hubble constant $H_0=72.3 \pm 1.4$ km/s/Mpc and a lower value of the matter fluctuations $\sigma_8=0.786\pm 0.020$, while not degrading the fit to the cosmic microwave background (CMB) damping tail. We achieve this by introducing neutrino self-interactions in the presence of a non-vanishing sum of neutrino masses. This strongly interacting neutrino cosmology prefers $N_{\rm eff} = 4.02 \pm 0.29$, which has interesting implications for particle model-building and neutrino oscillation anomalies. We show that the absence of the neutrino free-streaming phase shift on the CMB can be compensated by shifting the value of other cosmological parameters, hence providing an important caveat to the detections made in the literature. Due to their impact on the evolution of the gravitational potential at early times, self-interacting neutrinos and their subsequent decoupling leave a rich structure on the matter power spectrum. In particular, we point out the existence of a novel localized feature appearing on scales entering the horizon at the onset of neutrino free-streaming. While the interacting neutrino cosmology provides a better global fit to current cosmological data, we find that traditional Bayesian analyses penalize the model as compared to the standard cosmological. Our analysis shows that it is possible to find radically different cosmological models that nonetheless provide excellent fits to the data, hence providing an impetus to thoroughly explore alternate cosmological scenarios.

Tables

Table 1. Table 1 : Adopted prior ranges

Parameter	Prior
${\log }_{10}\left(G_{\mathrm{eff}}{\mathrm{MeV}}^2\right)$	$[-5.5,-0.000001]$
$\sum m_{\nu }$ [eV]	$[0.0001,1.5]$
$N_{\mathrm{eff}}$	$[2.0,5.0]$
${\Omega }_b{h}^2$	$[0.01,0.04]$
${\Omega }_c{h}^2$	$[0.08,0.16]$
$100{\theta }_{\mathrm{MC}}$	$[1.03,1.05]$
$\tau$	$[0.01,0.25]$
$\ln \left({10}^{10}{A}_s\right)$	$[2,4]$
$n_s$	$[0.85,1.1]$
$y_{\mathrm{cal}}$	$[0.9,1.1]$

Table 2. Table 2 : TT + lens + BAO + H 0 subscript 𝐻 0 H_{0} Constraints: Parameter 68 % percent 68 68\% Confidence Limits

Parameter	Strongly Interacting Neutrino Mode	Moderately Interacting Neutrino Mode
${\Omega }_b{h}^2$	${0.02245}_{-0.00033}^{+0.00029}$	$0.02282\pm 0.00030$
${\Omega }_c{h}^2$	${0.1348}_{-0.0049}^{+0.0056}$	${0.1256}_{-0.0039}^{+0.0035}$
$100{\theta }_{\mathrm{MC}}$	$1.04637\pm 0.00056$	${1.04062}_{-0.00056}^{+0.00049}$
$\tau$	$0.080\pm 0.031$	${0.127}_{-0.029}^{+0.034}$
$\sum m_{\nu }$ [eV]	${0.42}_{-0.20}^{+0.17}$	${0.40}_{-0.23}^{+0.17}$
$N_{\mathrm{eff}}$	$4.02\pm 0.29$	$3.79\pm 0.28$
${\log }_{10}\left(G_{\mathrm{eff}}{\mathrm{MeV}}^2\right)$	$-{1.35}_{-0.066}^{+0.12}$	$-{3.90}_{-0.93}^{+1.0}$
$\ln \left({10}^{10}{A}_s\right)$	$3.035\pm 0.060$	${3.194}_{-0.056}^{+0.068}$
$n_s$	$0.9499\pm 0.0098$	${0.993}_{-0.012}^{+0.013}$
$H_0$ [km/s/Mpc]	$72.3\pm 1.4$	$71.2\pm 1.3$
${\Omega }_m$	$0.3094\pm 0.0083$	$0.3010\pm 0.0080$
${\sigma }_8$	$0.786\pm 0.020$	${0.813}_{-0.020}^{+0.023}$
${10}^9{A}_s$	${2.08}_{-0.13}^{+0.11}$	$2.44\pm 0.15$
${10}^9{A}_s{e}^{-2\tau }$	$1.771\pm 0.016$	${1.892}_{-0.017}^{+0.019}$
$r_{\ast }$ [Mpc]	$136.3\pm 2.4$	$139.1\pm 2.3$
$100{\theta }_{\ast }$	$1.04604\pm 0.00056$	${1.04041}_{-0.00064}^{+0.00058}$
$D_A$  [Gpc]	$13.03\pm 0.23$	$13.37\pm 0.21$
$r_{\mathrm{drag}}$ [Mpc]	$138.8\pm 2.5$	$141.6\pm 2.3$

Table 3. Table 3 : Comparison to Λ ​ CDM Λ CDM \Lambda\mathrm{CDM} for TT + lens + BAO + H 0 subscript 𝐻 0 H_{0}

Parameter	Strongly Interacting Neutrino Mode	Moderately Interacting Neutrino Mode
$\Delta {\chi }_{\mathrm{low}\ell }^2$	$0.66$	$-0.75$
$\Delta {\chi }_{\mathrm{high}\ell }^2$	$-1.15$	$1.08$
$\Delta {\chi }_{\mathrm{lens}}^2$	$0.06$	$-0.24$
$\Delta {\chi }_{H_0}^2$	$-6.68$	$-6.12$
$\Delta {\chi }_{\mathrm{BAO}}^2$	$-0.81$	$-0.36$
$\Delta {\chi }_{\mathrm{Total}}^2$	$-7.91$	$-6.39$
$\Delta \mathrm{AIC}$	$-1.91$	$-0.39$

Table 4. Table 4 : Mode Comparison. Here, ℬ SI ​ ν subscript ℬ SI 𝜈 \mathcal{B}_{\mathrm{SI}\nu} is the Bayes factor between the SI ν 𝜈 \nu and the MI ν 𝜈 \nu mode, ℛ SI ​ ν subscript ℛ SI 𝜈 \mathcal{R}_{\mathrm{SI}\nu} is theie maximum likelihood ratio, and Δ ​ χ 2 = χ SI ​ ν 2 − χ MI ​ ν 2 Δ superscript 𝜒 2 subscript superscript 𝜒 2 SI 𝜈 subscript superscript 𝜒 2 MI 𝜈 \Delta\chi^{2}=\chi^{2}_{\mathrm{SI}\nu}-\chi^{2}_{\mathrm{MI}\nu} . The low- ℓ ℓ \ell dataset refers to low- ℓ ℓ \ell TEB if polarization was included and low- ℓ ℓ \ell TT if only temperature was used. Similarly, the high- ℓ ℓ \ell dataset refers to high- ℓ ℓ \ell TT,TE,EE if polarization was included and high- ℓ ℓ \ell TT if only temperature was used.

Parameter	TT,TE,EE	TT,TE,EE + lens + $H_0$	TT + lens + BAO	TT + lens + BAO + $H_0$
${ℬ}_{\mathrm{SI}\nu }$	$0.03\pm 0.01$	$0.10\pm 0.04$	$0.13\pm 0.04$	$0.37\pm 0.10$
${ℛ}_{\mathrm{SI}\nu }$	$0.26$	$0.63$	$0.81$	$2.14$
$\Delta {\chi }_{\mathrm{low}\ell }^2$	$2.47$	$2.18$	$2.00$	$1.41$
$\Delta {\chi }_{\mathrm{high}\ell }^2$	$0.22$	$-0.16$	$-1.53$	$-2.23$
$\Delta {\chi }_{\mathrm{lens}}^2$	–	$1.34$	$0.16$	$0.30$
$\Delta {\chi }_{H_0}^2$	–	$-2.12$	–	$-0.56$
$\Delta {\chi }_{\mathrm{BAO}}^2$	–	–	$-0.20$	$-0.44$
$\Delta {\chi }_{\mathrm{Total}}^2$	$2.69$	$0.92$	$0.43$	$-1.52$

Table 5. Table 5 : Comparison of the interacting neutrino cosmology to Λ ​ CDM + N eff + ∑ m ν Λ CDM subscript 𝑁 eff subscript 𝑚 𝜈 \Lambda\mathrm{CDM}+N_{\mathrm{eff}}+\sum m_{\nu} for TT + lens + BAO + H 0 subscript 𝐻 0 H_{0}

Parameter	Strongly Interacting Neutrino Mode	Moderately Interacting Neutrino Mode
$\Delta {\chi }_{\mathrm{low}\ell }^2$	$2.40$	$0.99$
$\Delta {\chi }_{\mathrm{high}\ell }^2$	$-3.40$	$-1.17$
$\Delta {\chi }_{\mathrm{lens}}^2$	$-0.20$	$-0.50$
$\Delta {\chi }_{H_0}^2$	$-1.32$	$-0.76$
$\Delta {\chi }_{\mathrm{BAO}}^2$	$-0.81$	$-0.36$
$\Delta {\chi }_{\mathrm{Total}}^2$	$-3.33$	$-1.81$
$\Delta \mathrm{AIC}$	$-1.33$	$0.19$

Table 6. Table 6 : Strongly interacting neutrino cosmology: Parameter 68 % percent 68 68\% Confidence Limits

Parameter	TT,TE,EE	TT,TE,EE + lens + $H_0$	TT + lens + BAO	TT + lens + BAO + $H_0$
${\Omega }_b{h}^2$	${0.02219}_{-0.00022}^{+0.00024}$	$0.02257\pm 0.00018$	${0.02239}_{-0.00036}^{+0.00029}$	${0.02245}_{-0.00033}^{+0.00029}$
${\Omega }_c{h}^2$	${0.1189}_{-0.0038}^{+0.0032}$	$0.1222\pm 0.0032$	${0.1311}_{-0.0065}^{+0.0090}$	${0.1348}_{-0.0049}^{+0.0056}$
$100{\theta }_{\mathrm{MC}}$	${1.04590}_{-0.00045}^{+0.00061}$	${1.04622}_{-0.00042}^{+0.00054}$	${1.04623}_{-0.00044}^{+0.00067}$	$1.04637\pm 0.00056$
$\tau$	${0.0633}_{-0.0066}^{+0.0089}$	${0.0624}_{-0.0074}^{+0.0089}$	${0.082}_{-0.036}^{+0.028}$	$0.080\pm 0.031$
$\sum m_{\nu }$ [eV]	${0.166}_{-0.18}^{+0.064}$	${0.069}_{-0.066}^{+0.027}$	${0.39}_{-0.20}^{+0.16}$	${0.42}_{-0.20}^{+0.17}$
$N_{\mathrm{eff}}$	${2.88}_{-0.22}^{+0.19}$	$3.20\pm 0.18$	$3.80\pm 0.45$	$4.02\pm 0.29$
${\log }_{10}\left(G_{\mathrm{eff}}{\mathrm{MeV}}^2\right)$	$-{1.60}_{-0.089}^{+0.14}$	$-{1.55}_{-0.080}^{+0.12}$	$-{1.41}_{-0.066}^{+0.20}$	$-{1.35}_{-0.066}^{+0.12}$
$\ln \left({10}^{10}{A}_s\right)$	${2.995}_{-0.015}^{+0.019}$	$2.994\pm 0.017$	${3.036}_{-0.071}^{+0.054}$	$3.035\pm 0.060$
$n_s$	$0.9273\pm 0.0080$	$0.9412\pm 0.0061$	$0.947\pm 0.011$	$0.9499\pm 0.0098$
$H_0$ [km/s/Mpc]	${66.2}_{-1.9}^{+2.3}$	$70.1\pm 1.3$	$71.1\pm 2.2$	$72.3\pm 1.4$
${\Omega }_m$	${0.327}_{-0.026}^{+0.013}$	${0.2961}_{-0.011}^{+0.0075}$	$0.3115\pm 0.0090$	$0.3094\pm 0.0083$
${\sigma }_8$	${0.799}_{-0.017}^{+0.041}$	${0.824}_{-0.010}^{+0.015}$	$0.786\pm 0.020$	$0.786\pm 0.020$
${10}^9{A}_s$	${1.998}_{-0.030}^{+0.039}$	$1.998\pm 0.034$	${2.09}_{-0.15}^{+0.10}$	${2.08}_{-0.13}^{+0.11}$
${10}^9{A}_s{e}^{-2\tau }$	$1.760\pm 0.014$	$1.763\pm 0.013$	$1.766\pm 0.016$	$1.771\pm 0.016$
$Y_P$	$0.2430\pm 0.0029$	$0.2476\pm 0.0024$	${0.2549}_{-0.0048}^{+0.0060}$	$0.2577\pm 0.0034$
$r_{\ast }$ [Mpc]	$145.8\pm 2.0$	$143.0\pm 1.7$	${138.2}_{-4.3}^{+3.2}$	$136.3\pm 2.4$
$100{\theta }_{\ast }$	${1.04626}_{-0.00046}^{+0.00060}$	${1.04629}_{-0.00044}^{+0.00053}$	${1.04604}_{-0.00046}^{+0.00060}$	$1.04604\pm 0.00056$
$D_A$ [Gpc]	$13.93\pm 0.19$	$13.67\pm 0.16$	${13.21}_{-0.41}^{+0.30}$	$13.03\pm 0.23$
$r_{\mathrm{drag}}$ [Mpc]	$148.5\pm 2.1$	$145.6\pm 1.8$	${140.8}_{-4.3}^{+3.3}$	$138.8\pm 2.5$

Table 7. Table 7 : Moderately Interacting Neutrino Mode: Parameter 68 % percent 68 68\% Confidence Limits

Parameter	TT,TE,EE	TT,TE,EE + lens + $H_0$	TT + lens + BAO	TT + lens + BAO + $H_0$
${\Omega }_b{h}^2$	$0.02203\pm 0.00023$	$0.02246\pm 0.00018$	${0.02254}_{-0.00035}^{+0.00030}$	$0.02282\pm 0.00030$
${\Omega }_c{h}^2$	$0.1191\pm 0.0031$	$0.1220\pm 0.0027$	${0.1220}_{-0.0046}^{+0.0039}$	${0.1256}_{-0.0039}^{+0.0035}$
$100{\theta }_{\mathrm{MC}}$	$1.04085\pm 0.00044$	$1.04063\pm 0.00040$	$1.04086\pm 0.00058$	${1.04062}_{-0.00056}^{+0.00049}$
$\tau$	${0.0642}_{-0.0082}^{+0.0095}$	${0.0645}_{-0.0073}^{+0.0090}$	$0.108\pm 0.033$	${0.127}_{-0.029}^{+0.034}$
$\sum m_{\nu }$ [eV]	${0.150}_{-0.16}^{+0.054}$	${0.052}_{-0.052}^{+0.020}$	${0.28}_{-0.23}^{+0.12}$	${0.40}_{-0.23}^{+0.17}$
$N_{\mathrm{eff}}$	$2.95\pm 0.19$	$3.29\pm 0.16$	${3.44}_{-0.38}^{+0.30}$	$3.79\pm 0.28$
${\log }_{10}\left(G_{\mathrm{eff}}{\mathrm{MeV}}^2\right)$	$-{4.44}_{-0.77}^{+0.58}$	$-4.26\pm 0.69$	$-4.12\pm 0.77$	$-{3.90}_{-0.93}^{+1.0}$
$\ln \left({10}^{10}{A}_s\right)$	${3.059}_{-0.019}^{+0.022}$	${3.067}_{-0.016}^{+0.019}$	$3.150\pm 0.067$	${3.194}_{-0.056}^{+0.068}$
$n_s$	$0.9548\pm 0.0089$	$0.9718\pm 0.0073$	${0.980}_{-0.015}^{+0.014}$	${0.993}_{-0.012}^{+0.013}$
$H_0$ [km/s/Mpc]	${65.3}_{-1.7}^{+2.2}$	$69.3\pm 1.2$	${69.3}_{-1.9}^{+1.7}$	$71.2\pm 1.3$
${\Omega }_m$	${0.335}_{-0.025}^{+0.012}$	${0.3021}_{-0.010}^{+0.0077}$	$0.3075\pm 0.0092$	$0.3010\pm 0.0080$
${\sigma }_8$	${0.798}_{-0.016}^{+0.038}$	${0.826}_{-0.011}^{+0.014}$	${0.809}_{-0.018}^{+0.021}$	${0.813}_{-0.020}^{+0.023}$
${10}^9{A}_s$	$2.132\pm 0.043$	${2.148}_{-0.035}^{+0.039}$	${2.34}_{-0.18}^{+0.14}$	$2.44\pm 0.15$
${10}^9{A}_s{e}^{-2\tau }$	$1.875\pm 0.018$	$1.888\pm 0.016$	$1.880\pm 0.021$	${1.892}_{-0.017}^{+0.019}$
$Y_P$	$0.2439\pm 0.0027$	$0.2486\pm 0.0022$	${0.2506}_{-0.0048}^{+0.0041}$	$0.2550\pm 0.0035$
$r_{\ast }$ [Mpc]	$145.5\pm 1.9$	$142.8\pm 1.5$	${141.9}_{-2.7}^{+3.0}$	$139.1\pm 2.3$
$100{\theta }_{\ast }$	$1.04117\pm 0.00054$	$1.04066\pm 0.00047$	$1.04086\pm 0.00070$	${1.04041}_{-0.00064}^{+0.00058}$
$D_A$ [Gpc]	$13.97\pm 0.17$	$13.72\pm 0.14$	${13.63}_{-0.25}^{+0.28}$	$13.37\pm 0.21$
$r_{\mathrm{drag}}$ [Mpc]	$148.3\pm 1.9$	$145.4\pm 1.6$	${144.5}_{-2.8}^{+3.1}$	$141.6\pm 2.3$

Equations

$${\mathcal{L}}_{\mathrm{int}}=g_{ij}{\bar{\nu }}_i{\nu }_j\phi ,$$

$$g_{ij}\equiv g_{\nu }{\delta }_{ij},$$

$$|\mathcal{M}{|}_{{\nu }_i}^2$$

$$=2G_{\mathrm{eff}}^2\left(s^2+t^2+u^2\right),$$

$$f_{\nu }(\mathbf{x},\mathbf{p},\tau )=f_{\nu }^{(0)}(\mathbf{p},\tau )[1+{\Theta }_{\nu }(\mathbf{x},\mathbf{p},\tau )],$$

$$f_{\nu }^{(0)}(\mathbf{p},\tau )=\frac{1}{e^{p/T_{\nu }}+1},$$

$${\Theta }_{\nu }(\mathbf{x},\mathbf{p},\tau )=-\frac{d\ln f_{\nu }^{(0)}}{d\ln p}\frac{\delta T_{\nu }(\mathbf{x},\tau )}{{\bar{T}}_{\nu }(\tau )},$$

$${\Xi }_{\nu }(\mathbf{x},\mathbf{p},\tau )$$

$${\tilde{\Xi }}_{\nu }(\mathbf{k},\mathbf{p},\tau )=\sum _{l=0}^{\infty }(-i{)}^l(2l+1){\nu }_l(k,p,\tau )P_l(\mu ),$$

$$C_{\nu }^{(1)}[\mathbf{p}]$$

$$\displaystyle\quad\times\sum_{l=0}^{\infty}(-i)^{l}(2l+1)\nu_{l}P_{l}(\mu)\Bigg{(}A\left(\frac{p}{T_{\nu}}\right)$$

$$\displaystyle\qquad\qquad\qquad\qquad+B_{l}\left(\frac{p}{T_{\nu}}\right)-2D_{l}\left(\frac{p}{T_{\nu}}\right)\Bigg{)},$$

$$\frac{\partial {\nu }_l}{\partial \tau }+k\frac{q}{ϵ}\left(\frac{l+1}{2l+1}{\nu }_{l+1}-\frac{l}{2l+1}{\nu }_{l-1}\right)$$

$$-4\left[\frac{\partial \varphi }{\partial \tau }{\delta }_{l0}+\frac{k}{3}\frac{ϵ}{q}\psi {\delta }_{l1}\right]$$

$$\displaystyle\qquad=-a\frac{G_{\rm eff}^{2}T_{\nu}^{5}\nu_{l}}{f_{\nu}^{(0)}(q)}\left(\frac{T_{\nu,0}}{q}\right)\Bigg{(}A\left(\frac{q}{T_{\nu,0}}\right)$$

$$\displaystyle\qquad\qquad\qquad\qquad+B_{l}\left(\frac{q}{T_{\nu,0}}\right)-2D_{l}\left(\frac{q}{T_{\nu,0}}\right)\Bigg{)},$$

$$\frac{\partial F_l}{\partial \tau }+k\left(\frac{l+1}{2l+1}{F}_{l+1}-\frac{l}{2l+1}{F}_{l-1}\right)$$

$$-4\left[\frac{\partial \varphi }{\partial \tau }{\delta }_{l0}+\frac{k}{3}\psi {\delta }_{l1}\right]=-a{G}_{\mathrm{eff}}^2{T}_{\nu }^5{\alpha }_l{F}_l,$$

$${\alpha }_l$$

$${\rho }_{\mathrm{R}}=\left[1+N_{\mathrm{eff}}\frac{7}{8}{\left(\frac{4}{11}\right)}^{4/3}\right]{\rho }_{\gamma },$$

$${\ddot{d}}_{\mathrm{c}}+\frac{\dot{a}}{a}{\dot{d}}_{\mathrm{c}}=-k^2\psi ,$$

$$d_{\mathrm{c}}\equiv {\delta }_{\mathrm{c}}-3\varphi ,$$

$$d_{\mathrm{c}}(k,\tau )=-\frac{9}{2}{\varphi }_{\mathrm{p}}+k^2{\int }_0^{\tau }d{\tau }^{\prime }{\tau }^{\prime }\psi (k,{\tau }^{\prime })\ln ({\tau }^{\prime }/\tau ),$$

$$r_{\ast }={\int }_0^{a_{\ast }}\frac{c_{\mathrm{s}}(a)}{a^{2}H(a)}da,D_{\mathrm{A}}={\int }_{a_{\ast }}^1\frac{da}{a^{2}H(a)},$$

$$\mathcal{Z}\equiv \mathrm{Pr}\left(\mathbf{d}|\mathcal{M}\right)={\int }_{{\Omega }_{\theta }}\mathrm{Pr}\left(\mathbf{d}|\theta ,\mathcal{M}\right)\mathrm{Pr}\left(\theta |\mathcal{M}\right)d\theta ,$$

$${\mathcal{B}}_{\mathrm{SI}\nu }\equiv \frac{\mathrm{Pr}\left({\mathcal{M}}_{\mathrm{SI}\nu }|\mathbf{d}\right)}{\mathrm{Pr}\left({\mathcal{M}}_{\mathrm{MI}\nu }|\mathbf{d}\right)}=\frac{{\mathcal{Z}}_{\mathrm{SI}\nu }}{{\mathcal{Z}}_{\mathrm{MI}\nu }}\frac{\mathrm{Pr}\left({\mathcal{M}}_{\mathrm{SI}\nu }\right)}{\mathrm{Pr}\left({\mathcal{M}}_{\mathrm{MI}\nu }\right)}.$$

$${\mathcal{R}}_{\mathrm{SI}\nu }=\frac{\max \left[\mathcal{L}\left({\theta }_{\mathrm{SI}\nu }|\mathbf{d}\right)\right]}{\max \left[\mathcal{L}\left({\theta }_{\mathrm{MI}\nu }|\mathbf{d}\right)\right]}.$$

$$\mathrm{AIC}=-2\ln \left(\mathcal{L}\right)+2k={\chi }_{\mathrm{Total}}^2+2k,$$

$$\Delta \mathrm{AIC}={\mathrm{AIC}}_{\mathrm{I}\nu }-{\mathrm{AIC}}_{\Lambda \mathrm{CDM}}=\Delta {\chi }^2+2\Delta k,$$

$$\frac{d{f}_{{\nu }_i}}{d\lambda }=\sum _{j,k,l=1}^3{C}_{{\nu }_i+{\nu }_j\to {\nu }_k+{\nu }_l}[f_{{\nu }_i},f_{{\nu }_j}{f}_{{\nu }_k},f_{{\nu }_l}]$$

$$d{s}^2=a^2(\tau )[-(1+2\psi )d{\tau }^2+(1-2\varphi )d{\vec{x}}^2],$$

$$P^{\mu }\equiv \frac{d{x}^{\mu }}{d\lambda },$$

$$\frac{d}{d\lambda }=\frac{d\tau }{d\lambda }\frac{d}{d\tau }=P^0\frac{d}{d\tau }=\frac{E(1-\psi )}{a}\frac{d}{d\tau },$$

$$\frac{d{f}_{{\nu }_i}}{d\tau }=\frac{\partial f_{{\nu }_i}}{\partial \tau }+\frac{\mathbf{p}}{E}\cdot \nabla f_{{\nu }_i}+p\frac{\partial f_{{\nu }_i}}{\partial p}\left[-\mathcal{H}+\frac{\partial \varphi }{\partial \tau }-\frac{E}{p^2}\mathbf{p}\cdot \nabla \psi \right],$$

Full text

I Introduction

The neutrino sector of the Standard Model (SM) of particle physics is a promising area to search for new phenomena that could help pinpoint the Ultraviolet completion of the SM. Indeed, terrestrial neutrino experiments have identified several anomalies that could potentially indicate the presence of new physics in the neutrino sector (see, e.g., Ref. Dentler et al. (2018) for a recent review). Of particular significance are the $\nu_{\mu}\rightarrow\nu_{\rm e}$ appearance results from the MiniBooNE Aguilar-Arevalo et al. (2018) and LSND Aguilar-Arevalo et al. (2001) collaborations which, if interpreted within a neutrino oscillation framework that includes an extra sterile neutrino, would indicate the presence of such a sterile neutrino at very high statistical significance. Within this “3+1” neutrino oscillation framework, these results are, however, very difficult to reconcile with the absence of anomalies in the $\nu_{\mu}\rightarrow\nu_{\mu}$ disappearance as probed by recent atmospheric Aartsen et al. (2017, 2016) and short-baseline Adamson et al. (2017a, b) experiments. If these results are confirmed by future analyses, it is likely that new physics beyond the sterile+active oscillation models would be necessary to resolve the tension between neutrino appearance and disappearance data.

Astrophysical and cosmological observations provide complementary means of probing the properties of neutrinos. This is perhaps best illustrated by the cosmological constraints on the sum of neutrino masses $\sum m_{\nu}<0.12$ eV Aghanim et al. (2018) obtained by combining cosmic microwave background (CMB) data from the Planck satellite with baryon acoustic oscillation (BAO) measurements. Cosmological observables such as the CMB and large-scale structure (LSS) are also sensitive to the presence of new interactions (see e.g. Refs. Bialynicka-Birula (1964); Bardin et al. (1970); Gelmini and Roncadelli (1981); Chikashige et al. (1980); Barger et al. (1982); Raffelt and Silk (1987); Kolb and Turner (1987); Konoplich and Khlopov (1988); Berkov et al. (1988); Belotsky et al. (2001); Hannestad (2005); Chacko et al. (2005); Hannestad and Raffelt (2005); Sawyer (2006); Mangano et al. (2006); Friedland et al. (2007); Hooper (2007); Serra et al. (2010); van den Aarssen et al. (2012); Jeong and Takahashi (2013); Laha et al. (2014); Archidiacono et al. (2015); Ng and Beacom (2014); Cherry et al. (2014); Archidiacono et al. (2016a); Cherry et al. (2016); Archidiacono et al. (2016b); Dvali and Funcke (2016); Capozzi et al. (2017); Brust et al. (2017); Forastieri (2017); Balducci et al. (2018); Lorenz et al. (2019); Choi et al. (2018)) in the neutrino sector that would modify their standard free-streaming behavior during the radiation-dominated epoch following their weak decoupling. In the literature, a phenomenological description based on the $c_{\rm eff}$ and $c_{\rm vis}$ parametrization Hu (1998) has often been used to test the free-streaming nature of neutrinos in the early Universe Trotta and Melchiorri (2005); Melchiorri and Serra (2006); De Bernardis et al. (2008); Smith et al. (2012); Archidiacono et al. (2011, 2012); Gerbino et al. (2013); Archidiacono et al. (2013); Melchiorri et al. (2014); Sellentin and Durrer (2015); Ade et al. (2014). While these analyses generally find results consistent with the standard neutrino cosmology, they are difficult to interpret in terms of possible new interactions among neutrinos, as emphasized in Refs. Cyr-Racine and Sigurdson (2014); Oldengott et al. (2015). Other works Chacko et al. (2004); Beacom et al. (2004); Bell et al. (2006); Cirelli and Strumia (2006); Basboll et al. (2009); Cyr-Racine and Sigurdson (2014); Archidiacono and Hannestad (2014); Oldengott et al. (2015); Forastieri et al. (2015, 2017); Lancaster et al. (2017); Oldengott et al. (2017); Barenboim et al. (2019) have used more physical parameterizations that make the connection to the underlying particle nature of the neutrino interaction more transparent.

In particular, Ref. Oldengott et al. (2015) has developed a rigorous treatment of the evolution of cosmological neutrino fluctuations in the presence of neutrino self-interactions mediated by either a massive or massless new scalar particle. Using this framework, Ref. Oldengott et al. (2017) used CMB data to put constraints on the strength of neutrino self-interactions in the early Universe for the case of a massive mediator. These results largely confirmed earlier constraints from Refs. Cyr-Racine and Sigurdson (2014); Archidiacono and Hannestad (2014); Lancaster et al. (2017) obtained using an approximate (but nonetheless accurate) form of the neutrino Boltzmann hierarchy. Interestingly, these studies, which focused on four-neutrino interactions parametrized by a Fermi-like coupling constant $G_{\rm eff}$, found a bimodal posterior distribution for this latter parameter. While the first (and statistically dominant) posterior mode is consistent with the onset of neutrino free-streaming being in the very early Universe, the second posterior mode corresponds to a much delayed onset of free-streaming to $z_{\nu,\rm dec}\sim 8300$. In Ref. Lancaster et al. (2017), a previously unknown multi-parameter degeneracy involving the amplitude of scalar fluctuations, the scalar spectral index, the Hubble constant, and the neutrino self-interacting strength was identified as being responsible for the existence of this second posterior mode. While intriguing, the neutrino interaction strength favored by this mode is nearly ten orders of magnitude above the standard weak interaction. Taken at face value, this likely constitutes a very serious challenge from a model-building perspective.

Nevertheless, given the current tensions among terrestrial and atmospheric neutrino experiments described above, is the “interacting” neutrino mode hinting at the presence of new physics beyond the SM? The simplified interaction models used in Refs. Cyr-Racine and Sigurdson (2014); Lancaster et al. (2017); Oldengott et al. (2017); Barenboim et al. (2019) are likely capturing parts of a more realistic neutrino interaction scenario, hence leading to a somewhat suboptimal fit to the cosmological data. One aspect that has been neglected in studies of self-interacting neutrinos so far is the presence of neutrino mass. The impact of massive neutrinos on the CMB and matter clustering has been studied extensively in the literature (see e.g. Refs. Kaplinghat et al. (2003); Lesgourgues and Pastor (2006, 2014); Vagnozzi et al. (2017); Lattanzi and Gerbino (2018)). One of the aims of this paper is to understand how the effects of massive neutrinos on cosmological observables are modified when self-interactions are present in the early Universe.

Tensions are also growing between different late-time measurements of the Hubble constant $H_{0}$ Riess et al. (2016); Bernal et al. (2016); Birrer et al. (2019); Riess et al. (2019) and those based on CMB data Aghanim et al. (2018). Measurements of the amplitude of matter fluctuations at low redshifts (often parametrized using $\sigma_{8}$) from weak gravitational lensing and cluster counts are all consistently lower than that inferred from the CMB Heymans et al. (2012); Joudaki et al. (2018); Hikage et al. (2018). While the statistical significance of the deviation of each individual measurement is less than $3\sigma$, all recent measurements of the amplitude fluctuations in the local universe are below the Planck value. Physics beyond $\Lambda\mathrm{CDM}$ has been proposed to reconcile these tensions, such as early dark energy Mörtsell and Dhawan (2018); Poulin et al. (2018a, b, 2018), dark matter interactions Di Valentino et al. (2017, 2018), decaying dark matter Pandey et al. (2019); Vattis et al. (2019); Enqvist et al. (2015); Bringmann et al. (2018), modified gravity Renk et al. (2017); Khosravi et al. (2017), and new relativistic species D’Eramo et al. (2018), among others. However, these propositions often struggle to remedy both tensions simultaneously.

In this paper, we study how the presence of self-interacting massive neutrinos in the early Universe affect cosmological observables such as the CMB, with an eye on how these new effects could help relieve the current tensions among different datasets. In Sec. II, we describe the simplified neutrino interaction model used in this work. In Sec. III, we present the cosmological perturbation equations for massive self-interacting neutrinos. In Sec. IV, we describe the physical impacts that massive self-interacting neutrinos have on the CMB and the matter power spectrum. In Sec. V, we outline the data and method used in our cosmological analyses of self-interacting. The results from these analyses are presented in Sec. VI and discussed in Sec. VIII. We conclude and highlight future directions in Sec. IX.

II Neutrino Interaction Model

In this work, we focus on a simple framework that captures the most important cosmological aspects of realistic neutrino interaction models. We note however that building a successful model of neutrino self-interaction that respects the gauge and flavor structure of the SM likely requires the introduction of a light sterile species which mass-mixes with the active neutrinos and is itself coupled to a massive scalar or vector mediator (see e.g. Refs. Bringmann et al. (2014); Dasgupta and Kopp (2014); Hannestad et al. (2014); Cherry et al. (2014); Tang (2015); Kouvaris et al. (2015); Chu et al. (2015); Ghalsasi et al. (2017); Hasenkamp (2016)). The presence of these new interactions in the sterile sector suppresses the effective mixing angle between the active and sterile species at early times, ensuring that Big Bang nucleosynthesis (BBN) constraints are respected . At later times once the active-sterile oscillation rate becomes comparable to the finite temperature effective potential resulting from the new interaction, the mixing angle is no longer suppressed hence allowing the active and sterile sectors to partially thermalize with each other Archidiacono et al. (2015, 2016a, 2016b).

Diagonalizing the mass matrix of such a model leads to an effective interaction Lagrangian between the different neutrino mass eigenstates of the generic form

[TABLE]

where $g_{ij}$ is a (generally complex) coupling matrix, $\nu_{i}$ is a left-handed neutrino Majorana spinor, and the indices $i,j$ labeled the neutrino mass eigenstates. Here we have assumed a Yukawa-type interaction with a massive scalar $\varphi$, but note that the results presented in this work also apply if a massive vector is assumed instead. The Lagrangian given in Eq. (1) could also arise in models where neutrinos couple to a Majoron Gelmini and Roncadelli (1981); Simpson et al. (2018); Berlin and Blinov (2018).

In models where the new interaction arises through active-sterile mixing, the structure of the coupling matrix $g_{ij}$ would generally depend on the flavor content of each mass eigenstate. For instance, a mass eigenstate made of mostly active flavors will couple very weakly to the massive scalar $\varphi$, while an eigenstate being largely composed of the sterile species would couple more strongly to the mediator. In other models of neutrino interaction, the structure of the coupling matrix could be more arbitrary. In all cases though, $g_{ij}$ is subject to important flavor-dependent bounds Lessa and Peres (2007); Bakhti and Farzan (2017); Arcadi et al. (2019) arising from meson, tritium, and gauge boson decay kinematics.

In this work, we consider the simple case of a universal coupling $g_{\nu}$ between every neutrino mass eigenstate

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. While the universal coupling case is likely unrealistic for the reason outlined above, it does provide a simple benchmark to test the sensitivity of cosmological data to new neutrino physics.

We work in the contact-interaction limit in which the mass of the $\varphi$ mediator is much larger than the typical energy of the scattering event. In this case, one can integrate out this massive mediator and write the interaction as a four-fermion contact interaction. This is an excellent approximation at the energy scale probed by the CMB for $m_{\varphi}\gtrsim 1$ keV. In this limit, the squared scattering amplitude for a neutrino $\nu_{i}$ to interact with any other neutrino in the thermal bath is

[TABLE]

where we have defined the dimensionfull coupling constant $G_{\rm eff}\equiv|g_{\nu}|^{2}/m_{\varphi}^{2}$. Here, $s$, $t$, and $u$ are the standard Mandelstam variables. While our phenomenological model described by $G_{\rm eff}$ is unlikely to accurately capture all the complexity of novel neutrino interactions, it is nonetheless a useful framework to identify the interesting parameter space, as described in Ref. Ng and Beacom (2014).

Introducing new neutrino interactions has an impact beyond cosmology. For a low mass mediator ($<10$ MeV), SN 1987A Kolb and Turner (1987), Big Bang nucleosynthesis (BBN) Ahlgren et al. (2013); Huang et al. (2018), and the detection of ultra-high energy neutrinos at IceCube Ng and Beacom (2014); Ioka and Murase (2014); Cherry et al. (2016) provide some of the strongest constraints, with the latter bound having the potential of being the most stringent in the near future. Other limits Bilenky et al. (1993); Bardin et al. (1970); Bilenky and Santamaria (1999) coming from Z-boson decay do not directly apply at the energy scale probed by the CMB. Also, elastic collisions caused by the new interaction do not affect the time it takes for neutrinos to escape supernovae Manohar (1987); Dicus et al. (1989), although they could lead to interesting phenomena (see e.g. Refs. Davoudiasl and Huber (2005); Sher and Triola (2011); Fayet et al. (2006); Choi and Santamaria (1990); Blennow et al. (2008); Galais et al. (2012)). Finally, supernova cooling puts bounds on the coupling of majorons to SM neutrinos Kachelriess et al. (2000); Farzan (2003); Zhou (2011); Jeong et al. (2018), but the applicability of these likely depends on the details of the exact coupling matrix used.

III Cosmological Perturbations

In this section we summarize the key ingredients and simplifications entering our derivation of the Boltzmann equation governing the evolution of massive and self-interacting neutrino fluctuations, at first order in perturbation theory. Our computation uses two main approximations:

• •

Based on previous studies Cyr-Racine and Sigurdson (2014); Lancaster et al. (2017), we assume that neutrinos decouples while still in the relativistic regime. We thus neglect the presence of the small neutrino mass in the computation of collision integrals. As we shall see, our final results are consistent with this approximation.

• •

We assume that the neutrino distribution function remains exactly thermal throughout the epoch at which neutrinos decouple and start free-streaming. This thermal approximation (also called, relaxation time approximation) implies that the only possible neutrino perturbations are local temperature fluctuations. This approximation was shown to be very accurate in Ref. Oldengott et al. (2017) for the type of interaction we consider here.

Conformal Newtonian gauge is used throughout this section.

III.1 Neutrino distribution function and perturbation variables

We present a detailed derivation of the left-hand side of the Boltzmann equation for massive neutrino in Appendix B (see also Ref. Ma and Bertschinger (1995)). Our starting point is to expand the neutrino distribution function as

[TABLE]

where $\mathbf{x}$ denotes the spatial coordinates, $\tau$ is conformal time, and $\mathbf{p}$ is the proper momentum. The background (spatially uniform) neutrino distribution function is taken to be of a Fermi-Dirac shape

[TABLE]

where $p=|\mathbf{p}|$. In the ultra-relativistic regime, for which the thermal approximation implies that the only possible neutrino perturbations are local temperature fluctuations, the perturbation variable $\Theta_{\nu}$ admits the form

[TABLE]

where $\bar{T}_{\nu}$ is the background neutrino temperature, and $\delta T_{\nu}$ is its perturbation. It is therefore convenient to introduce the temperature fluctuation variables $\Xi_{\nu}$

[TABLE]

which is independent of $\mathbf{p}$ in the thermal approximation for massless neutrinos. However, the presence of a nonvanishing neutrino mass and the non-negligible momentum transfered in a typical neutrino-neutrino collision would in general introduce some extra $\mathbf{p}$-dependence to $\Xi_{\nu}$ Ma and Bertschinger (1995). This turns the Boltzmann equation of self-interacting neutrinos into a differentio-integral equation that is particularly difficult to solve exactly Oldengott et al. (2015). In practice though, the absence of energy sources or sinks coupled to the neutrino sector implies that the momentum dependence of the right-hand side of Eq. (7) should be vanishingly small at early times when neutrinos form a highly-relativistic tightly-coupled fluid. This allows us to neglect the momentum-dependence of $\Xi_{\nu}$ in the computation of the collision integrals, an approximation that was found to be accurate in Ref. Oldengott et al. (2017). We do retain, however, the momentum dependence of $\Xi_{\nu}$ in the left-hand side of the Boltzmann equation.

In this work, we only consider scalar perturbations and thus expand the angular dependence of the $\tilde{\Xi}_{\nu}$ variable (the Fourier transform on $\Xi_{\nu}$) in Legendre polynomials $P_{l}(\mu)$

[TABLE]

where $\mu$ is the cosine of the angle between $\mathbf{k}$ and $\mathbf{p}$. Before presenting the equation of motion for the neutrino multipole moments $\nu_{l}$, we discuss the structure of the collision integrals.

III.2 Collision term

The details of the collision term calculation for the $\nu\nu\rightarrow\nu\nu$ process is given in Appendix C. As explained above, the main simplification entering this calculation is the use of the thermal approximation in which we neglect the momentum dependence of the $\nu_{l}$ variables. Under this assumption, the collision term at first order in perturbation theory $C_{\nu}^{(1)}$ can be written as

[TABLE]

where the functions $A(x)$, $B_{l}(x)$, and $D_{l}(x)$ are given in Eqs. (84), (85), and (86), respectively. Here, we have adopted the notation $T_{\nu}\equiv\bar{T}_{\nu}$ to avoid clutter.

III.3 Boltzmann equation for self-interacting neutrinos

Substituting the collision term from Eq. (9) into Eq. (32) and performing the $\mu$ integral yields the equation of motion for the different neutrino multipoles $\nu_{l}$. They can be summarized in the following compact form

[TABLE]

where we have introduced the comoving momentum $\mathbf{q}\equiv a\mathbf{p}$, $q=|\mathbf{q}|$, $\epsilon=\sqrt{q^{2}+a^{2}m_{\nu}^{2}}$, a is the scale factor normalized to $a=1$ today, $\delta_{mn}$ is the Kronecker delta function, and $T_{\nu,0}$ is the current $(a=1)$ temperature of the neutrinos. The fact that the collision term is directly proportional to $\nu_{l}$ is a consequence of our use of the thermal approximation. We note that energy and momentum conservation ensures that $A+B_{0}-2D_{0}=0$ and $A+B_{1}-2D_{1}=0$, respectively.

As is standard in analyses of massive neutrino cosmologies, we shall consider our neutrino sector to be composed of a mix of massive and massless neutrinos. In the massless case ($q=\epsilon$), one can integrate Eq. (10) over the comoving momentum to yield a simpler neutrino multipole hierarchy Cyr-Racine and Sigurdson (2014); Lancaster et al. (2017)

[TABLE]

where

[TABLE]

and where we denoted the massless perturbations as $F_{l}$ to distinguish them from the massive neutrino variables $\nu_{l}$.

We implement these modified Boltzmann equations in the cosmological code CAMB Lewis et al. (2000). For computational speed, we precompute the functions $A$, $B_{l}$ and $D_{l}$ on a grid of $q/T_{\nu,0}$ values and use an interpolation routine to access them when solving the cosmological perturbation equations. As in standard CAMB, we use a sparse 3-point grid of $q/T_{\nu,0}$ values to evaluate the integrals required to compute the energy density and momentum flux of massive neutrinos. We have checked convergence of our scheme against a 5-point momentum grid and found negligible difference in the CMB and matter power spectrum in the parameter space of interest. We also precompute the coefficient $\alpha_{l}$ and tabulate them. We emphasize that energy and momentum conservation ensures that $\alpha_{0}=\alpha_{1}=0$, which we have checked with high accuracy.

For simplicity, we assume throughout this paper that the neutrino sector contains one massive neutrino, with the remaining neutrino species being massless. All neutrinos are assumed to interact with the same coupling strength $G_{\rm eff}$. We find that varying the number of massive neutrinos and number of mass eigenstates, while holding $N_{\mathrm{eff}}$ and $\sum m_{\nu}$ constant, has a very small impact on the CMB and matter power spectra for all values of $G_{\rm eff}$ consistent with the data used here. It is however possible that future data might be sensitive to the way $\sum m_{\nu}$ is spread among different mass eigenstates.

At early times, the large self-interaction rate of neutrinos renders the equations of motion for multipoles $l\geq 2$ extremely stiff. To handle this, we employ a tight-coupling scheme Cyr-Racine and Sigurdson (2011) in which multipole moments with $l\geq 2$ are set to zero at early times. Once the neutrino self-interaction rate falls to about a 1000 times the Hubble expansion rate, we turn off this tight-coupling approximation and allows power to flow to the higher multipoles. We have checked that this switch happens early enough as to not affect the accuracy of our results. After neutrino decoupling, once they become non-relativistic, we revert to the standard velocity-integrated truncated Boltzmann hierarchy as described in Ref. Lewis and Challinor (2002). We also modified the adiabatic initial conditions for the cosmological perturbations to take into account the absence of free-streaming neutrinos at early times. Finally, throughout this work, we use the standard BBN predictions to compute the helium abundance given the abundance of relativistic species and the baryon-to-photon ratio.

IV Effect on Cosmological Observables

IV.1 Cosmic microwave background

In the standard cosmological paradigm, free-streaming neutrinos travel supersonically through the photon-baryon plasma at early times, hence gravitationally pulling photon-baryon wavefronts slightly ahead of where they would be in the absence of neutrinos Bashinsky and Seljak (2004); Baumann et al. (2016); Choi et al. (2018). As a result, the free-streaming neutrinos imprint a net phase shift in the CMB power spectra towards larger scales (smaller $\ell$), as well as a slight suppression of its amplitude. Free-streaming neutrinos thus lead to a physical size of the photon sound horizon at last scattering $r_{*}$ that is slightly larger than it would otherwise be. This phase shift is thought to be a robust signature of the presence of free-streaming radiation in the early Universe Follin et al. (2015); Baumann et al. (2017); Choi et al. (2018).

The neutrino self-interactions mediated by the coupling constant $G_{\mathrm{eff}}$ delay the time at which neutrinos begin to free-stream. Fourier modes entering the causal horizon while neutrinos are still tightly-coupled will not experience the gravitational tug of supersonic neutrinos and will therefore not receive the associated phase shift and amplitude reduction. Compared to the standard $\Lambda$CDM model, neutrino self-interactions thus shift the CMB power spectra peaks towards smaller scales (larger $\ell$) and boost their fluctuation amplitude. This leads to a net reduction of the physical size of the photon sound horizon at last scattering $r_{*}$. As we shall see, this is the key feature of our model that helps reconcile CMB and late-time measurements of the Hubble constant $H_{0}$.

The left panels of Fig. 1 show the temperature CMB power spectra and their relative difference to a $\Lambda\mathrm{CDM}$ model for different values of $G_{\mathrm{eff}}$, $\sum m_{\nu}$, and $N_{\rm eff}$ to illustrate the effects of neutrino self-scattering in the presence of a non-vanishing mass term. Here, we keep $\Omega_{\rm m}$ fixed as $\sum m_{\nu}$ changes, and use the best-fit Planck TT+lowP+lensing $\Lambda\mathrm{CDM}$ values as our fiducial cosmology (Planck Collaboration et al., 2016b). The middle left panel of Fig. 1 displays the combined effect of changing both $G_{\rm eff}$ and $\sum m_{\nu}$. For the minimal sum of neutrinos masses $\sum m_{\nu}=0.06\,\mathrm{eV}$, an interaction strength of $G_{\mathrm{eff}}=10^{-4}\,\mathrm{MeV}^{-2}$ (solid blue line) has for only effect a slight increase of power at large multipoles. On the other hand, increasing the neutrino coupling strength to $G_{\mathrm{eff}}=10^{-2}\,\mathrm{MeV}^{-2}$ (solid red line) significantly boosts the amplitude of the TT spectrum and introduces a clear phase shift (identifiable from the oscillatory pattern of the residuals), which are the two telltale signatures of self-scattering neutrinos as described above.

Increasing the sum of neutrinos masses to $\sum m_{\nu}=0.23\,\mathrm{eV}$ (at fixed $\Omega_{\rm m}$) delays the time of matter-radiation equality. The delay slightly increases the amplitude of the TT spectrum near the first few acoustic peaks and dampens the spectrum at smaller scales (see dashed black line in Fig. 1). The resulting changes to the photon-baryon sound horizon at recombination and to the angular diameter distance to the surface of last scattering create a net phase shift towards low $\ell$ (Lesgourgues and Pastor, 2006), that is, in the opposite direction to that caused by increasing $G_{\mathrm{eff}}$. This opens the door for possible cancellations between the relative phase shift (as compared to $\Lambda$CDM) caused by neutrino self-scattering and that resulting from a large sum of neutrino masses. Such cancellation partially occurs in the middle left panel of Fig. 1 for $G_{\mathrm{eff}}=10^{-2}\,\mathrm{MeV}^{-2}$ as $\sum m_{\nu}$ is increased from $0.06$ to $0.23$ eV (dashed red line). Similarly, the boost in amplitude from $G_{\mathrm{eff}}$ can also compensate for the damping effects of increasing $\sum m_{\nu}$ at small scales (see e.g. the dashed blue line). Overall, we see that the effect of massive neutrinos and increased interaction strength are nearly additive111Indeed, combining the spectrum for {$G_{\mathrm{eff}}=10^{-2}\,\mathrm{MeV}^{-2}$, $\Sigma m_{\nu}=0.06\,\mathrm{eV}$} (solid red line) with that of the $\Sigma m_{\nu}=0.23\,\mathrm{eV}$ $\Lambda\mathrm{CDM}$ model (dashed black line) yields a spectrum similar to the model with {$G_{\mathrm{eff}}=10^{-2}\,\mathrm{MeV}^{-2}$, $\Sigma m_{\nu}=0.23\,\mathrm{eV}$} (dashed red line)., reflecting the fact that the physical processes associated with each of these properties take place at different times in the cosmological evolution.

The lower left panel of Fig. 1 displays the impact of increasing the energy density of the neutrino fluid, which we parametrized here through the standard parameter $N_{\rm eff}$, defined via the relation

[TABLE]

where $\rho_{\rm R}$ and $\rho_{\gamma}$ are the total energy density in radiation and in photons, respectively. The effects on the CMB of increasing $N_{\rm eff}$ have been well-studied in the literature (see e.g. Ref. Hou et al. (2013)) for the case of free-streaming neutrinos. For fixed values of the angular scale of the sound horizon, the epoch of matter-radiation equality, and the physical baryon abundance, it was found that the most important net impact of increasing $N_{\rm eff}$ was to damp the high-$\ell$ tail of the TT spectrum and to induce a phase shift towards larger scales (low-$\ell$). Interestingly, self-interacting neutrinos can partially compensate for these effects, hence pointing to a possible degeneracy between $G_{\rm eff}$ and $N_{\rm eff}$. An example of this can be seen in the dotted red line in the lower left panel of Fig. 1, where the excess of damping caused by $N_{\rm eff}=4.046$ (dotted black line) is compensated by suppressing neutrino free-streaming with $G_{\rm eff}=10^{-2}$ MeV*-2*.

$G_{\mathrm{eff}}$ affects the EE polarization power spectrum in a similar manner as the temperature spectrum. The right panel of Fig. 1 shows that the phase shift between the standard $\Lambda$CDM model and that with self-interacting neutrinos is more visible in this case due to the sharp, well defined peaks of the polarization spectrum (Baumann et al., 2016). This allows to directly see in which direction the spectrum is shifted compared to $\Lambda$CDM since the oscillations in the residuals lean in the direction of the phase shift, that is, there is a sharper drop off in the residuals in the direction that the spectrum is shifted. Once again, we clearly see that the absence of phase shift caused by a large value of $G_{\mathrm{eff}}$ can be partially canceled by increasing $\sum m_{\nu}$, in a nearly additive fashion. For the EE polarization spectrum, suppressing neutrino free-streaming can somewhat compensate the extra damping caused by a large $N_{\rm eff}$ (at fixed $\theta_{*}$, $z_{\rm eq}$, and $\Omega_{\rm b}h^{2}$; see lower right panel of Fig. 1).

IV.2 Matter power spectrum

The growth of matter fluctuations is sensitive to the presence of self-interacting neutrinos through the neutrinos’ impact on the two gravitational potentials $\phi$ and $\psi$. Indeed, neutrino self-interactions suppress the anisotropic stress of the universe, leading to $\phi-\psi=0$ before the onset of neutrino free-streaming. This contrasts with the $\Lambda$CDM case for which $\phi=(1+2R_{\nu}/5)\psi$ on large scales at early times for the adiabatic mode Ma and Bertschinger (1995), where $R_{\nu}$ is the radiation free-streaming fraction. This difference in the evolution of the potentials modifies the gravitational source term driving the growth of matter fluctuations. The equation describing the evolution of dark matter fluctuations can be written in Fourier space as Bashinsky and Seljak (2004)

[TABLE]

where

[TABLE]

and where $\delta_{\rm c}=\delta\rho_{\rm c}/\rho_{\rm c}$ is the standard dark matter energy density contrast in Newtonian gauge. Here, an overhead dot denotes a derivative with respect to conformal time $\tau$. The gauge-invariant variable $d_{\rm c}$ represents the fractional dark matter number density perturbation by unit coordinate volume. At late times, $d_{\rm c}$ is nearly equal to $\delta_{\rm c}$ and it is thus a useful quantity to understand the structure of the matter power spectrum at $z=0$. In the radiation-dominated epoch where $\dot{a}/a=\tau^{-1}$, the solution to Eq. (14) can be written Dodelson (2003)

[TABLE]

where $\phi_{\rm p}$ is the primordial value of $\phi$ on large scales. The integral appearing in Eq. (16) obtains most of its contribution when $k\tau\sim 1$. The changes to the growth of dark matter fluctuations can thus be understood by examining the behavior of the $\psi$ potential at horizon entry.

We compare the evolution of $\psi$ in the presence of self-interacting neutrinos with $G_{\rm eff}=10^{-2}$ MeV*-2* to that of standard $\Lambda$CDM in the left panel of Fig. 2. There, we track the evolution of three different Fourier modes: $k=10\,h/\mathrm{Mpc}$ which enters the horizon during the radiation dominated era while neutrinos are still tightly-coupled to each other, $k=0.3\,h/\mathrm{Mpc}$ which roughly corresponds to the scale entering the horizon when neutrinos begin to free-stream, and $k=10^{-3}\,h/\mathrm{Mpc}$ which does not enter the horizon until far after neutrino decoupling. We use here the same cosmological parameters as in Fig. 1. The resulting evolution of dark matter fluctuations for these three modes is shown in the right panel of Fig. 2.

When modes enter the horizon during the radiation-dominated era, the gravitational potential $\psi$ decays in an oscillatory fashion Dodelson (2003). The absence of anisotropic stress implies that $\psi$ starts its oscillatory decaying behavior from a larger amplitude. This boosts the amplitude of the envelope of the decaying oscillations as compared to $\Lambda$CDM, leading to an overall slower decay. While this at first increases the amplitude of dark matter fluctuations at horizon entry as compared to $\Lambda$CDM (see bottom right panel of Fig. 2), the subsequent oscillations of the integrand appearing in Eq. (16) lead to a net damping of the dark matter perturbation amplitude. Another way to think about this is that the slower decay of the potential $\psi$ in the presence of self-interacting neutrinos reduces the horizon-entry boost that dark matter fluctuations experience as compared to $\Lambda$CDM.

For modes entering the horizon at the time of neutrino decoupling, the potential $\psi$ begins decaying from its larger value with $R_{\nu}=0$ but rapidly locks into its standard $\Lambda$CDM evolution due to the onset of neutrino free-streaming. This case thus displays the quickest damping of the $\psi$ potential after horizon entry, which leads to a net boost of dark matter fluctuations as compared to $\Lambda$CDM. Indeed, these modes receive a positive contribution near horizon entry from the integral in Eq. (16), but without the subsequent extra damping due to the $\psi$ potential quickly converging to its $\Lambda\mathrm{CDM}$ behavior. The evolution of the $k=0.3\,h/\mathrm{Mpc}$ mode in Fig. 2 displays this behavior.

Finally, modes entering the horizon well-after the onset of neutrino free-streaming behave exactly like their $\Lambda\mathrm{CDM}$ counterparts, as illustrated by the $k=10^{-3}\,h^{-1}\mathrm{Mpc}$ mode in Fig. 2. Taking together the evolution of the different Fourier modes entering before, during, and after neutrino decoupling, we expect the matter power spectrum to have the following properties (at fixed neutrino mass). For large wavenumbers entering the horizon while neutrinos are tightly coupled, we expect the matter power spectrum to be suppressed compared to $\Lambda$CDM. As we go to larger scales and approach modes entering the horizon at the onset of free-streaming, we expect a “bump”-like feature displaying an excess of power as compared to $\Lambda$CDM. As we go to even larger scales, the matter power spectrum is expected to asymptote to its standard $\Lambda$CDM value.

These expectations are indeed realized as shown in Fig. 3. The middle panel shows the power spectrum ratios between the interacting neutrino models and $\Lambda$CDM. Focusing for the moment on the cases with $\sum m_{\nu}=0.06$ eV, we see that the matter power spectrum is damped at large wavenumbers and then displays a broad peak-like feature with an excess of power as compared to $\Lambda$CDM. The shape of this power excess is determined by the neutrino visibility function Cyr-Racine and Sigurdson (2014) encoding the details of neutrino decoupling. Increasing the sum of neutrino masses (at fixed $\Omega_{\rm m}$) leads to a damping of the matter power spectrum on small scales (Lesgourgues and Pastor, 2006; Lattanzi and Gerbino, 2018). This standard reduction of power is shown for $\Lambda$CDM as the thick black dashed line in Fig. 3. Interestingly, this small-scale suppression is also present for self-interacting neutrinos and occurs in addition to that caused by the slower decay of the gravitational potential $\psi$ discussed above. Thus, the matter power spectrum for massive self-interacting neutrinos is even more suppressed at large $k$ than in the standard $\Lambda$CDM case with massive neutrinos.

This fact might seem counterintuitive at first since the reduction of small-scale power from massive neutrinos is often refereed to as “free-streaming” damping. We see this moniker is somewhat of a misnomer since the damping is present whether or not neutrinos are actually free-streaming. Instead, the small-scale reduction of power is simply caused by the large pressure term that prohibits neutrino clustering on these scales. This pressure term is always there as long as neutrinos are relativistic, even when neutrinos are self-scattering. As was the case for the CMB, the effects of a non-vanishing sum of neutrino masses and large $G_{\rm eff}$ are largely additive. Comparing the $\sum m_{\nu}=0.23$ eV cases to that with $\sum m_{\nu}=0.06$ eV in Fig. 3 for both interacting neutrino models shown illustrates this well. Again, this near additivity reflects the fact that part of the effect comes from the behavior of dark matter fluctuations at horizon entry, while the rest is caused by the large pressure term of relativistic neutrinos on small scales.

The lowest panel of Fig. 3 shows the effect of increasing $N_{\rm eff}$ (at fixed $\theta_{*}$, $z_{\rm eq}$, and $\Omega_{\rm b}h^{2}$) on the matter power spectrum. For $\Lambda$CDM, the main impact is to increase the amplitude of Fourier modes that enter the causal horizon during radiation domination. This results from the larger radiation density and free-streaming fraction $R_{\nu}$ Bashinsky and Seljak (2004) at early times. Suppressing neutrino free-streaming for $N_{\rm eff}=4.046$ (dotted red line) nullifies this increase of power on small-scales, even leading to a net damping compared to $\Lambda$CDM for $k>10h$/Mpc. However, as the neutrinos start to decouple from one another, the larger radiation density leads to a higher amplitude feature on scales entering the horizon at that time.

We thus see that taken together, the joint effect of $G_{\mathrm{eff}}$, $\sum m_{\nu}$, and $N_{\rm eff}$ can lead to matter power spectra having a significantly different structure and shape than the standard $\Lambda$CDM paradigm.

V Data & Methodology

We use our modified versions of CAMB Lewis et al. (2000) and CosmoMC + Multinest Lewis and Bridle (2002); Feroz and Hobson (2008) to place constraints on $G_{\mathrm{eff}}$, $N_{\mathrm{eff}}$, and $\sum m_{\nu}$, as well as the standard cosmological parameters. We use nested sampling Skilling (2006) to ensure that we properly sample our posterior, which we expect to be multi-modal as in previous cosmological studies of self-interacting neutrinos Cyr-Racine and Sigurdson (2014); Lancaster et al. (2017); Oldengott et al. (2017).

We use a combination of CMB and low-redshift data sets in our analysis:

• •

TT: low-$\ell$ and high-$\ell$ CMB temperature power spectrum from the Planck 2015 release222Explicitly, we use the likelihood plik_lite_v18_TT for high-$\ell$ and commander_rc2_v1.1_l2_29_B at low-$\ell$. (Planck Collaboration et al., 2016a).

• •

EE, TE: low-$\ell$ and high-$\ell$ CMB E-mode polarization and their temperature cross-correlation from the Planck 2015 data release333Explicitly, we use the likelihood plik_lite_v18_TTTEEE for high-$\ell$ and lowl_SMW_70_dx11d_2014_10_03_v5c_Ap at low-$\ell$. (Planck Collaboration et al., 2016a). The 2015 polarization data is known to have residual systematics and results drawn using this dataset should be interpreted with caution. While our main conclusions will not make use of this dataset, we nonetheless present results including this dataset for completeness.

• •

lens: CMB lensing data from the Planck 2015 data release (Planck Collaboration et al., 2016c).

• •

BAO: Baryon Acoustic Oscillation (BAO) measurements from the 6dF Galaxy Survey constraining $D_{V}$ at $z=0.106$ (Beutler et al., 2011), Sloan Digital Sky Survey (SDSS-III) Baryon Oscillation Spectroscopic Survey (BOSS) data release 11 low-$z$ data measuring $D_{V}$ at $z=0.32$ and CMASS data measuring $D_{V}$ at $z=0.57$ (Anderson et al., 2014), and data from the SDSS Main Galaxy Sample measuring $D_{V}$ at $z=0.15$ (Ross et al., 2015)

• •

${\bf H_{0}}$: Local measurement444We note that the mean value of $H_{0}$ used in our analysis is slightly lower ($\sim 0.14\sigma$) than the value quoted in the published version of Ref. Riess et al. (2016) (ours corresponds to the value found in an earlier version of their manuscript). We do not expect this very small difference to impact our results in any way. of the Hubble parameter $H_{0}=73.0\pm 1.75\,\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$ at $z=0.04$ from Ref. (Riess et al., 2016).

We use the lite high-$\ell$ likelihood, which marginalizes over nuisance parameters, to reduce the number of free parameters in our analysis. We use the following data set combinations for our nested sampling analysis: ‘TT+lens+BAO’, ‘TT+lens+BAO+$H_{0}$’, ‘TT,TE,EE’, and ‘TT,TE,EE+lens+$H_{0}$’.

In Table 1 we list our adopted prior ranges. We place uniform priors on all these parameters, except for the Planck calibration parameter $y_{\rm cal}$, for which we use a Gaussian prior $y_{\rm cal}=1.0000\pm 0.0025$. For the analyses using the Planck polarization data, we also include a Gaussian prior on the optical depth to reionization given by $\tau=0.058\pm 0.012$ from Ref. Adam et al. (2016).

We use 2000 live points in our nested sampling runs, setting the target sampling efficiency to 0.3. We impose an accuracy threshold on the log Bayesian evidence of $20\%$, which ensures that our confidence intervals are highly accurate. We use the mode-separation feature of Multinest to isolate each posterior mode and compute their respective summary statistics.

VI Results

In this section, we first present the main highlights of our analysis, before discussing the physical properties of the two categories of interacting neutrino models that are favored by the data. We end this section with a brief discussion about which properties of interacting neutrino models help alleviate current tensions in cosmological data. Throughout this section, we quote and analyze results for the $\mathrm{TT+lens+BAO+}H_{0}$ data set combination unless otherwise specified. We discuss the impact of other data (including CMB polarization) in Sec. VIII and list parameter constraints for these other data set combinations in Table 6 and Table 7, in Appendix A.

VI.1 Highlights

Similarly to previous works Cyr-Racine and Sigurdson (2014); Lancaster et al. (2017); Oldengott et al. (2017), we find two unique neutrino cosmologies preferred by the data: a strongly interacting neutrino cosmology (hereafter denoted $\mathrm{SI}\nu$ model) characterized by $\mathrm{log}\left(G_{\mathrm{eff}}\,\mathrm{MeV}^{2}\right)=-1.35^{+0.12}_{-0.07}$ for the TT+lens+BAO+$H_{0}$ combination, and a moderately interacting neutrino cosmology (hereafter, $\mathrm{MI}\nu$ model) characterized by $\mathrm{log}\left(G_{\mathrm{eff}}\,\mathrm{MeV}^{2}\right)=-3.90^{+1.00}_{-0.93}$ for the same data set. Values of $G_{\rm eff}$ between these two modes are strongly disfavored by the data since they either prefer to have a phase shift that is largely consistent with free-streaming neutrinos, or no phase shift at all. We present constraints on cosmological parameters for the $\mathrm{SI}\nu$ and $\mathrm{MI}\nu$ modes in Table 2. While the MI$\nu$ cosmology was nearly indistinguishable from the $\Lambda$CDM scenario with massless neutrinos in previous work Lancaster et al. (2017), the addition of neutrino mass and $N_{\rm eff}$ in combination with the $H_{0}$ measurement from Ref. Riess et al. (2016) leads to a slight preference for a delayed onset of neutrino free-streaming. We expand more on this new development in subsection VI.3 below.

Cosmological parameters in the SI$\nu$ cosmology admit values that are significantly different from $\Lambda\mathrm{CDM}$:

The angular scale of the baryon-photon sound horizon at last scattering $100\theta_{*}=1.04604\pm 0.00056$ ($68\%$ C.L.) takes a value that is radically different ($>5\sigma$ away) than in the $\Lambda$CDM scenario, reflecting the absence of the free-streaming neutrino phase shift. 2. 2.

The large $N_{\mathrm{eff}}$ value $4.02\pm 0.29$ ($68\%$ C.L.) suggests the presence of an additional neutrino species, which might help reduce tensions between different neutrino oscillation experiments. 3. 3.

A smaller value of the baryon drag scale $r_{\mathrm{drag}}=138.8\pm 2.5\,\mathrm{Mpc}$ ($68\%$ C.L.) helps reconcile BAO with local Hubble constant measurements, leading to $H_{0}=72.3\pm 1.4\,\mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1}$ ($68\%$ C.L.). 4. 4.

The impact of self-interacting neutrinos on the growth of dark matter perturbations and a preferred suppressed spectrum of primordial scalar fluctuations lead to $\sigma_{8}=0.786\pm 0.020$ ($68\%$ C.L.).

To illustrate the ability of neutrino self-interactions to help resolve current cosmological tensions, we compare the $S_{8}\equiv\sigma_{8}\Omega_{\rm m}^{0.5}$ and $H_{0}$ 2D posteriors for the $\mathrm{SI}\nu$ model and $\mathrm{MI}\nu$ model with the base $\Lambda\mathrm{CDM}$ model in Fig. 4. We overlay bands for HSC constraints on $S_{8}$ Hikage et al. (2018) and local measurements of $H_{0}$ Riess et al. (2016). In order for our analysis to be independent from these measurements, we show posteriors for the TT+lens+BAO constraints for both the neutrino self-interaction models and $\Lambda\mathrm{CDM}$. Intriguingly, the strong neutrino self-interactions in the $\mathrm{SI}\nu$ model are able to independently produce the preferred values for $S_{8}$ and $H_{0}$, even without using these measurements in our analysis. The base $\Lambda\mathrm{CDM}$ model is unable to achieve these values, and the weak neutrino interactions of the $\mathrm{MI}\nu$ model can only achieve such values with weak significance.

In Table 3 we compute the $\Delta\chi^{2}$ values between the two neutrino self-interaction models and $\Lambda\mathrm{CDM}$. The data favor the strongly interacting neutrino model over $\Lambda\mathrm{CDM}$ with $\Delta\chi^{2}_{\mathrm{Total}}=-7.91$. This is a significant difference, even after accounting for the three extra parameters in the SI$\nu$ model (see subsection VII.2 for further discussion about this point). The preference for the self-interacting neutrinos comes from the local measurements of $H_{0}$, the high-$\ell$ TT data, and the BAO data.

In Fig. 5, we separate the posterior modes and plot their separate statistical distribution for the most salient parameters. For comparison, we also show the marginalized posteriors for the standard $\Lambda$CDM paradigm, as well as for its $N_{\rm eff}+\sum m_{\nu}$ two-parameter extension. In Fig. 6, we show the different covariances between the most relevant model parameters for three of the dataset combinations used in this work555Due to the presence of two posterior modes with different width, it is difficult to choose a smoothing scale that faithfully captures the intrinsic shape of the whole posterior while removing sampling noise. This particularly affects the $\mathrm{SI}\nu$ mode and results in a significantly reduced height which appears to visually suppress its statistical significance. See Fig. 11 in Appendix A for a figure with a smoothing scale more appropriate for the SI$\nu$ mode..

VI.2 Strongly interacting neutrino mode

The existence of the SI$\nu$ mode was first pointed out in Ref. Cyr-Racine and Sigurdson (2014), and further studied in Refs. Lancaster et al. (2017); Oldengott et al. (2017). As discussed there, the SI$\nu$ cosmology arises due to a multi-parameter degeneracy that opens up in CMB data when the onset of neutrino free-streaming is delayed until redshift $z\sim 8000$. This approximately coincides with the epoch when Fourier modes corresponding to multipole $\ell\approx 400$ enters the causal horizon (Lancaster et al., 2017), which lies somewhere between the first and second peak of the CMB temperature spectrum. We review below the properties of this alternate cosmology, emphasizing its differences with the standard $\Lambda$CDM model.

Sound horizon One of the most striking features of the SI$\nu$ model is the significantly larger value of the angular size of the sound horizon $\theta_{*}$. This is probably the most confusing aspect of our results since the angular size of the CMB sound horizon at last scattering is thought to be the best measured quantity in all of cosmology. To understand this apparent discrepancy, it is important to realize that the angular sound horizon is defined as $\theta_{*}\equiv r_{*}/D_{\rm A}$, where

[TABLE]

where $c_{\rm s}$ is the baryon-photon sound speed, $H$ is the Hubble rate, and $a_{*}$ is the scale factor at last scattering. We thus see that $\theta_{*}$ is purely defined in terms of background quantities, independent of the behavior of cosmological perturbations. In particular, it is independent of the gravitational tug that neutrinos exert on the photons.

Of course, when fitting CMB data we use the full temperature and polarization spectra computed from the Boltzmann equation which includes the effect of neutrinos. For the SI$\nu$ model, the absence of free-streaming neutrinos means that the CMB spectra do not receive the standard phase shift, and thus appear slightly displaced toward larger $\ell$ as compared to the corresponding $\Lambda$CDM spectra. In order to fit the data, we must compensate for this shift by increasing the value of $\theta_{*}$. Thus, the difference between the values of $\theta_{*}$ in the SI$\nu$ and $\Lambda$CDM models directly reflects the absence of the free-streaming neutrino phase shift in the former.

We note that it was a priori far from obvious that such a dramatic change in the angular size of the sound horizon was possible without introducing other artifacts that would significantly worsen the fit to CMB and BAO data. Our analysis shows that the larger value of $\theta_{*}$ is achieved by increasing $H_{0}$ and $\Omega_{\rm c}h^{2}$ above their $\Lambda$CDM values.

Primordial spectrum In addition to removing the CMB phase shift, suppressing neutrino free-streaming also increases the amplitude of the temperature and polarization spectra, as discussed in subsection IV.1. In the SI$\nu$ model, these changes are reabsorbed by modifying the primordial spectrum of scalar fluctuations parametrized by the amplitude $A_{\rm s}$ and spectral index $n_{\rm s}$. As was found in Refs. Cyr-Racine and Sigurdson (2014); Lancaster et al. (2017), lower values of both $A_{\rm s}$ and $n_{\rm s}$ are required to fit the temperature data in the SI$\nu$ mode. The difference between this alternative cosmology and $\Lambda$CDM is even more apparent if we compare the values of the parameter $A_{\rm s}e^{-2\tau}$ which directly determines the amplitude of the CMB temperature spectrum. As shown in Fig. 5, this amplitude parameter admits values that are radically ($>5\sigma$) different than in $\Lambda$CDM, again reflecting the large impact that suppressing neutrino free-streaming has on the CMB.

Neutrino properties The $\mathrm{SI}\nu$ model is consistent with having an entire additional neutrino species ($N_{\mathrm{eff}}=4.02\pm 0.29$, see Fig. 5), which has interesting implications for neutrino oscillation experiments. By comparing the SI$\nu$ cosmology with a more standard $\Lambda$CDM + $N_{\rm eff}$ + $\sum m_{\nu}$ model, we can understand how much of this preference is driven by the neutrino self-interaction. As shown in Fig. 5, the two-parameter extension of the $\Lambda$CDM cosmology already favors a larger $N_{\rm eff}$, but the introduction of strong neutrino self-interactions shifts the posterior to even larger values. To a certain extent, this shift is driven by the need to fit the large value of the local Hubble rate from Ref. Riess et al. (2016) by reducing the size of the sound horizon at the baryon drag epoch (see e.g. Ref. Aylor et al. (2019)). However, in the presence of free-streaming neutrinos, increasing $N_{\rm eff}$ also leads to a larger phase shift toward low $\ell$ which puts a limit on how much extra free-streaming radiation can be added before severely degrading the fit to CMB data. For the SI$\nu$ model, the absence of this phase shift allows for larger $N_{\rm eff}$, which leads to a smaller value of $r_{\rm drag}$ and, in turn, a larger Hubble constant. This is the key feature of the SI$\nu$ model that allows it to severely reduce the Hubble rate tension between CMB and late-time measurements, as we shall discuss in subsection VI.4.

The SI$\nu$ model also statistically prefers a nonzero value for the sum of neutrino masses. This preference was however already present at a less significant level ($<2\sigma$) in the $N_{\rm eff}$ + $\sum m_{\nu}$ extension of the $\Lambda$CDM scenario. In this latter case, nonzero neutrino masses arise from the need to suppress the amplitude of matter of fluctuations at late times (as measured here through CMB lensing) in the presence of a larger $N_{\rm eff}$ and $\Omega_{\rm m}$. For the SI$\nu$ model, the preference for a nonzero sum of neutrino masses is increased ($>2\sigma$) due to the even larger $N_{\rm eff}$ and $\Omega_{\rm m}$ values favored by this scenario.

We note that, in our analysis, the primordial helium abundance $Y_{P}$ is highly correlated with $N_{\mathrm{eff}}$ due to our use of the BBN consistency condition. Allowing $Y_{P}$ to take a different set of values in the SI$\nu$ scenario could lead to an even better fit to cosmological data.

Matter clustering Several competing effects act to set the amplitude of late-time matter fluctuations (as captured by the parameter $\sigma_{8}$) in the SI$\nu$ model. First, the large values of $N_{\rm eff}$ and $\Omega_{\rm m}$ (the latter necessary to keep the epoch of matter-radiation equality fixed) tend to boost the amplitude of matter fluctuation as discussed in subsection IV.2. On scales entering the horizon before the onset of neutrino free-streaming, this increase is counteracted by both a nonzero sum of neutrino masses and the reduction of the horizon entry boost for dark matter fluctuations in the presence of self-interacting neutrinos. Dark matter fluctuations entering the horizon during neutrino decoupling, which for the SI$\nu$ model are coincidentally those primarily contributing to $\sigma_{8}$, are however enhanced by the rapid decay of the gravitational potential on these scales. Finally, the lower amplitude and spectral index of the primordial scalar spectrum in the SI$\nu$ model tend to suppress power on scales probed by $\sigma_{8}$. Putting all of these effects together leads to a net lower value of $\sigma_{8}$, which, as discussed in the previous section, might be favored by some probes of late-time matter clustering. The overall shape of the matter power spectrum in the SI$\nu$ model will be further discussed in Sec. VIII.

VI.3 Moderately interacting neutrino mode

Within the $\mathrm{MI}\nu$ mode, the onset of neutrino free-streaming occurs before most Fourier modes probed by the Planck high-$\ell$ data enter the causal horizon. As such, the cosmological parameter values preferred by this mode are very similar to those from the $N_{\rm eff}$ + $\sum m_{\nu}$ extension of the $\Lambda$CDM scenario (see Fig. 5). The main difference here is that high-$\ell$ CMB modes do not receive the full amplitude suppression associated with free-streaming neutrinos due to the finite width of the neutrino visibility function. In other words, even though these high-$\ell$ modes enter the horizon after most neutrinos have started to free-stream, residual scattering in the neutrino sector still influences the amplitude of the CMB damping tail (see, e.g., the model with $G_{\rm eff}=10^{-4}$ MeV*-2* and $\sum m_{\nu}=0.06$ eV in Fig. 1). This increased small-scale power allows for a larger $N_{\rm eff}$, which, by reducing the baryon drag scale, leads to slightly larger Hubble constant. This shift is however quite small.

A surprising fact about the MI$\nu$ mode (also pointed out in Ref. Oldengott et al. (2017)) is that it shows a slight statistical preference for a nonzero value of $G_{\rm eff}$. As we can see in Fig. 6, this preference is nearly entirely driven by the local Hubble constant measurement of Ref. Riess et al. (2016). Indeed, removing this dataset from our analysis (blue contours) eliminates most of the preference for a nonzero value of $G_{\rm eff}$.

VI.4 Mediating Controversy: Effects on $H_{0}$ and $\sigma_{8}$

We show in 7(a) the impact of $G_{\rm eff}$, $N_{\rm eff}$ and $\sum m_{\nu}$ on the inferred value of the Hubble parameter. As described above, the large values of $N_{\rm eff}$ allowed in the presence of neutrino self-interactions reduce the size of the baryon drag scale, which allows a larger value of $H_{0}$ without damaging the fit to the BAO scale Aylor et al. (2019) and without introducing extra damping at large multipoles (see Fig. 1 for an illustration of this latter effect). In the SI$\nu$ model, this effect is compounded by the larger value of $\theta_{*}$ necessary to compensate for the absence of the free-streaming neutrino phase shift. This further slightly increases the value of $H_{0}$ necessary to fit the data, as can be seen by comparing the two modes in the left-most panel of 7(a).

It is worth noting that when $N_{\mathrm{eff}}$ is fixed at $3.046$, $H_{0}$ and $\sum m_{\nu}$ are usually negatively correlated (see e.g. Ref. Planck Collaboration et al. (2016b)). If both $N_{\mathrm{eff}}$ and $\sum m_{\nu}$ are allowed to vary, there is not a strong correlation between $N_{\mathrm{eff}}$ and $\sum m_{\nu}$ for CMB data alone. However, when $H_{0}$ or BAO data are added, $N_{\mathrm{eff}}$ and $\sum m_{\nu}$ become positively correlated (Hou et al., 2014) as shown in Fig. 6. The tight correlation between $N_{\mathrm{eff}}$ and $H_{0}$ then permits a positive correlation between $H_{0}$ and $\sum m_{\nu}$, seen in the third panel from the left of 7(a). Thus, instead of larger $\sum m_{\nu}$ being correlated with a smaller Hubble constant, here a larger sum of neutrino masses corresponds to a slightly larger $H_{0}$. Allowing the neutrinos to self-interact does not dramatically change the direction of this degeneracy, but does allow it to stretch to larger $H_{0}$ values.

$G_{\mathrm{eff}}$’s direct effects on matter clustering are scale dependent. As discussed in subsection IV.2, dark matter fluctuations that enter the horizon during neutrino decoupling receive a boost, while fluctuations that enter the horizon before neutrino decoupling are damped. For the $\mathrm{SI}\nu$ mode, neutrino decoupling is coincident with the modes entering the horizon that contribute most to $\sigma_{8}$, giving them a gravitational boost. However, low values of $A_{\rm s}$ and $n_{\rm s}$ must accompany a large $G_{\mathrm{eff}}$ for the $\mathrm{SI}\nu$, as discussed in subsection VI.2, which consequently damp these same scales. The combination of $G_{\mathrm{eff}}$’s effects thus leads to an overall decrease in matter clustering at scales probed by $\sigma_{8}$, seen in the $\mathrm{SI}\nu$ island in the left panel of 7(b).

The sum of neutrino masses is negatively correlated with $\sigma_{8}$ since massive neutrinos do not contribute to matter clustering for small scales where their pressure term is large (see third panel in 7(b)). Typically large $N_{\mathrm{eff}}$ boosts the dark matter fluctuations upon horizon entry, leading to a positive correlation between $N_{\mathrm{eff}}$ and $\sigma_{8}$. However, the positive correlation between $N_{\mathrm{eff}}$ and $\sum m_{\nu}$ when including BAO data causes $N_{\mathrm{eff}}$ and $\sigma_{8}$ to be negatively correlated (see second panel in 7(b)). This allows the interacting neutrino model to both be compatible with a large value of the Hubble constant and not overpredict the amplitude of matter fluctuations at late times.

VI.5 Impact of CMB polarization data

As can be seen in Fig. 6 and Fig. 7, the addition of EE polarization data tends to significantly reduce the statistical significance of the SI$\nu$ cosmology. This is in contrast with Ref. Lancaster et al. (2017) which found that polarization data slightly increased the significance of the SI$\nu$ mode (see also Ref. Oldengott et al. (2017)). The degradation of the fit for the SI$\nu$ model in our case is the result of (i) our use of the reionization optical depth prior from Ref. Adam et al. (2016) whenever we use polarization data, and (ii) our use of BBN calculations to predict the helium abundance for a given value of $N_{\rm eff}$. As we discuss in Sec. IX, it is likely that the fit could improve significantly by replacing this strong prior with the actual low-$\ell$ polarization data used to obtain it, and by letting the helium fraction float freely in the fit. As might be expected, the addition of the local Hubble constant measurement increases the statistical significance of the SI$\nu$ mode, as can be seen from TT,TE,EE+lens+$H_{0}$ data set in Fig. 11 in Appendix A.

VII Statistical Significance

In this section, we quantify the relative statistical significance of the two modes of the posterior, and compare the maximum likelihood values between our interacting neutrino models and standard extensions of the $\Lambda$CDM paradigm.

VII.1 Mode Comparison

To determine the statistical significance of the $\mathrm{SI}\nu$ mode relative to the $\mathrm{MI}\nu$ mode, we can compare their relative Bayesian evidence. It is defined as the parameter-averaged likelihood of the data

[TABLE]

where d is the data, $\mathcal{M}$ is the cosmological model, $\mathbf{\theta}$ are the parameters in model $\mathcal{M}$, and $\Omega_{\theta}$ is the domain of the model parameters. We use Multinest’s (Feroz et al., 2013) mode separation algorithm to compute the Bayesian evidence for each mode. In practice, this mode separation occurs near a neutrino coupling value of $\log_{10}\left(G_{\mathrm{eff}}\,\mathrm{MeV}^{2}\right)\approx-2.2$. This separation in parameter space defines $\Omega_{\theta}$ for each mode.

To compare the SI$\nu$ to the MI$\nu$ mode, we compute the following Bayes factor:

[TABLE]

We place a uniform prior on $\mathrm{log}_{10}\left(G_{\mathrm{eff}}\,\mathrm{MeV}^{2}\right)$ rather than a uniform prior on $G_{\mathrm{eff}}$ to avoid introducing a preferred energy scale. With our choice of prior, small values of $G_{\mathrm{eff}}$ can be thoroughly explored, which is particularly important since the actual Fermi constant governing neutrino interaction in the Standard Model takes the value $G_{\mathrm{F}}\sim\mathcal{O}\left(10^{-11}\,\mathrm{MeV}^{-2}\right)$. Taking a uniform prior on $G_{\mathrm{eff}}$ would greatly increase the statistical significance of the interacting mode (see Ref. Cyr-Racine and Sigurdson (2014)). We thus consider it conservative to adopt a uniform prior on $\mathrm{log}_{10}\left(G_{\mathrm{eff}}\,\mathrm{MeV}^{2}\right)$, but note that the statistical significance of the SI$\nu$ mode could be greatly enhanced by a different choice of prior.

The probability of the prior is equivalent for each mode (or cosmological model), so $\mathrm{Pr}\left(\mathcal{M}_{\mathrm{SI}\nu}\right)/\mathrm{Pr}\left(\mathcal{M}_{\mathrm{MI}\nu}\right)=1$. In Table 4 we show the Bayes factor for each data set combination we consider in this work. A Bayes factor less than unity indicates the data prefer the $\mathrm{MI}\nu$ mode for the specified parameter space. All values are below unity, indicating the data, on average, do not prefer the $\mathrm{SI}\nu$ mode. As expected though, incorporating the local Hubble rate measurement does increase the significance of the $\mathrm{SI}\nu$ mode.

A useful method to understand if the $\mathrm{SI}\nu$ mode is ever preferred and to further investigate the significance’s dependence on LSS data is to compare the maximum-likelihood value of each model:

[TABLE]

In Table 4 we show the maximum-likelihood ratios for the data set combinations in our analysis. Again, adding $H_{0}$ and CMB lensing data increases the likelihood of the strongly interacting mode. Intriguingly, the $\mathrm{SI}\nu$ mode has a larger maximum-likelihood (by a factor larger than 2) than the $\mathrm{MI}\nu$ mode for TT+lens+BAO+$H_{0}$ (see the unsmoothed posteriors in Fig. 11). It is reasonable that the Bayes factor for TT+lens+BAO+$H_{0}$ is below unity while the maximum-likelihood ratio is above unity since the former is a global, parameter-averaged statistic while the latter is based on a single set of best-case scenario parameters. This indicates that the parameter space for which strong neutrino interactions are preferred has a small volume.

It is also informative to look at the individual $\chi^{2}$ values for the different data sets. To compare the two modes, we list the $\Delta\chi^{2}=\chi^{2}_{\mathrm{SI}\nu}-\chi^{2}_{\mathrm{MI}\nu}$ values in Table 4. A positive $\Delta\chi^{2}$ value thus means that the MI$\nu$ mode is preferred, and vice versa. The $H_{0}$ and high-$\ell$ TT,TE,EE data show preference for the $\mathrm{SI}\nu$ mode for TT,TE,EE+lens+$H_{0}$, but this is compensated by a poorer fit to low-$\ell$ and CMB lensing data. For the TT+lens+BAO data combinations, the BAO and high-$\ell$ TT data display a slight preference for the $\mathrm{SI}\nu$ mode, which is again overshadowed by the low-$\ell$ data. We see that the slight preference for the SI$\nu$ mode with the TT+lens+BAO+$H_{0}$ data combination is largely due to improvement of the BAO and high-$\ell$ likelihoods.

VII.2 Comparison to $\Lambda\mathrm{CDM}$ and its extensions

Comparing how well each mode fits the data relative to $\Lambda\mathrm{CDM}$ and its common extensions tells us if these neutrino self-interaction models offer a viable improvement to current cosmological theory. For the TT+lens+BAO+$H_{0}$ data set, we list the $\Delta\chi^{2}=\chi^{2}_{\mathrm{SI}\nu}-\chi^{2}_{\Lambda\mathrm{CDM+ext}}$ values and the $\Delta\chi^{2}=\chi^{2}_{\mathrm{MI}\nu}-\chi^{2}_{\Lambda\mathrm{CDM+ext}}$ values for each observable in Table 5. Here, $\Lambda\mathrm{CDM+ext}$ refers to the $N_{\rm eff}+\sum m_{\nu}$ two-parameter extension of the $\Lambda\mathrm{CDM}$ cosmology. Comparison to plain $\Lambda$CDM was given in Table 3 above. For all data sets except the low-$\ell$ TT data, both modes offer a better fit to the data than $\Lambda\mathrm{CDM}+\mathrm{ext}$. In fact, the $\mathrm{SI}\nu$ mode has a total $\Delta\chi^{2}$ of $-3.33$, a significant difference. The improvement of the high-$\ell$ CMB data is notable since jointly fitting CMB and local $H_{0}$ data usually results in a worse fit to the CMB damping tail. For the SI$\nu$ model, this is somewhat compensated by a degradation of the low-$\ell$ fit.

What if the strong improvement in fit over $\Lambda\mathrm{CDM}$ is due to overfitting from the extra parameter we have added? To take this into account we compute the Akaike information criterion (AIC) (Akaike, 1974). The AIC takes into account how well the model fits the data and penalizes extra parameters, thereby discouraging overfitting. The AIC is defined as

[TABLE]

where $\chi^{2}_{\mathrm{Total}}=\chi^{2}_{\mathrm{low\,}\ell}+\chi^{2}_{\mathrm{high\,}\ell}+\chi^{2}_{\mathrm{lens}}+\chi^{2}_{H_{0}}+\chi^{2}_{\mathrm{BAO}}$, $\mathcal{L}$ is the maximum-likelihood, and $k$ is the number of fit parameters. Then we can write

[TABLE]

where $\Delta k$ is the difference in the number of parameters between the two models. The lower AIC between two models corresponds to the preferred model. Thus, for us, a negative $\Delta\mathrm{AIC}$ value indicates the data prefer the specified $\mathrm{I}\nu$ model over $\Lambda\mathrm{CDM}$, while a positive $\Delta\mathrm{AIC}$ value indicates the data prefer $\Lambda\mathrm{CDM}$ over the $\mathrm{I}\nu$ model.

We list the $\Delta\mathrm{AIC}$ values relative to $\Lambda\mathrm{CDM}+N_{\rm eff}+\sum m_{\nu}$ in Table 5. Here $\Delta k=1$, and the SI$\nu$ mode has a negative $\Delta\mathrm{AIC}=-1.33$, indicating a genuine statistical preference for the suppression of neutrino free-streaming in the early Universe for the TT+lens+BAO+$H_{0}$ data set. On the other hand, $\Delta\mathrm{AIC}=0.19$ for the $\mathrm{MI}\nu$ mode, indicating that the neutrino self-interactions do not add value to the fit beyond what is already provided by the $N_{\rm eff}+\sum m_{\nu}$ extension. Values of $\Delta\mathrm{AIC}$ between the I$\nu$ models and standard $\Lambda$CDM ($\Delta k=3$) are also given in Table 3. The fact that $\Delta\mathrm{AIC}$ values for the SI$\nu$ cosmology are similar ($-1.91$ versus $-1.33$) when comparing it to plain $\Lambda$CDM and $\Lambda$CDM $+N_{\rm eff}+\sum m_{\nu}$ means that suppressing neutrino free-streaming is the true driving factor behind the improvement of the fit. Thus, even after penalizing the self-interacting neutrino models for incorporating additional parameters, the TT+lens+BAO+$H_{0}$ data still significantly prefer the strongly interacting neutrino cosmology over $\Lambda\mathrm{CDM}$.

VIII Discussion

VIII.1 Cosmic microwave background

In Fig. 8, we plot the high-$\ell$ TT and EE power spectra residuals between the maximum likelihood parameters for each data set combination used and the best-fit Planck $\Lambda$CDM model. For the $\mathrm{SI}\nu$ mode (upper panels), the most striking feature of the residuals is the deficit of power at high multipoles ($\ell>1500$) as compared to $\Lambda$CDM for the TT+lens+BAO and TT+lens+BAO+$H_{0}$ data combinations. This is caused by the large value of $N_{\rm eff}$ and the resulting high helium abundance666We remind the reader that we use the standard BBN predictions to compute the helium abundance for given $N_{\rm eff}$ and $\Omega_{\rm b}h^{2}$ values. $Y_{P}$ for this category of models. This implies that the multi-parameter degeneracy that allows the SI$\nu$ cosmology to provide a decent fit to CMB temperature data at $\ell<1500$ could be broken by the addition of high-resolution CMB data (see e.g. Refs. Louis et al. (2017); Henning et al. (2018)). However, it is reasonable to assume that the BBN helium abundance is modified in the presence of the new neutrino physics we explore here, and that the deficit of power at large multipoles could be compensated by a smaller value of $Y_{P}$ Hou et al. (2013). We leave the study of the impact of a free helium fraction on interacting neutrino cosmologies to future works.

The EE polarization residuals shown in the right panel of Fig. 8 for the TT+lens+BAO and TT+lens+BAO+$H_{0}$ data combinations also display strong oscillations for the SI$\nu$ mode (upper panel). This implies that the shift in $\theta_{*}$ (and other parameters, see subsection VI.2) that was required to compensate for the absence of the free-streaming neutrino phase shift in the temperature spectrum does not fully realign the peaks of the polarization spectrum with the data. This is a consequence of the polarization data being more sensitive to the phase of the acoustic peaks Baumann et al. (2016). With the current size of the Planck error bars, this does not constitute an overwhelmingly strong constraint on the absence of a neutrino-induced phase shift, but it is possible that future CMB polarization data could entirely rule out this possibility.

The TT+TE+EE CMB-only data combination in the upper panels of Fig. 8 display an excess of power as compared to $\Lambda$CDM at nearly all scales, resulting in an overall poorer fit to the CMB data. At large multipoles, this is of course in contrast with the deficit of power that the TT+lens+BAO and TT+lens+BAO+$H_{0}$ fits display. Our use of the polarization-driven prior on the reionization optical depth from Ref. Adam et al. (2016) is largely responsible for this excess of power as compared to $\Lambda$CDM for the SI$\nu$ mode with the TT+TE+EE data set. Again, this shows that polarization data could in principle break the multi-parameter degeneracy that allows the SI$\nu$ cosmology to exist.

All dataset combinations we consider display an excess of power at $\ell<50$ for the SI$\nu$ mode. This is largely caused by the lower value of the scalar spectral index $n_{\rm s}$ which adds power on large scales. While error bars are large in this regime due to cosmic variance, the dip in power around $\ell\sim 20$-$30$ in the CMB temperature data tends to penalize any model displaying more low-$\ell$ power than $\Lambda$CDM. If this dip were to be explained by some other physics (from the inflationary epoch, say), then it is possible that the fit to the data using the SI$\nu$ cosmology could significantly improve.

It is important to emphasize how the suppression of neutrino free-streaming plays a very important role in the existence of the SI$\nu$ cosmology. To illustrate this, we plot in Fig. 9 the spectra corresponding to the best-fit $\mathrm{SI}\nu$ parameters for the TT+lens+BAO+$H_{0}$ data set but allow neutrino to free-stream at all times by setting $G_{\mathrm{eff}}=0$ (red dashed-dot line), along with the original TT+lens+BAO+$H_{0}$ best-fit $\mathrm{SI}\nu$ model (solid red) and a $\Lambda\mathrm{CDM}$ model with $\sum m_{\nu}=0.23\,\mathrm{eV}$ (black dashed), for reference. Here, the difference between the dashed-dot and solid red lines is entirely driven by the streaming property of neutrinos. Figure 9 reinforces our discussion from subsection VII.2 that the $G_{\mathrm{eff}}$ parameter plays a statistically significant role in improving the fit to the CMB data, beyond what is already provided by the $N_{\rm eff}+\sum m_{\nu}$ two-parameter extension of $\Lambda$CDM.

VIII.2 Matter clustering

We show in 10(a) the linear matter power spectrum residuals between the best-fit SI$\nu$ (and MI$\nu$) models and the corresponding $\Lambda$CDM models. The most striking feature for the SI$\nu$ mode is the overall red tilt of the matter power spectrum residuals for all data combinations shown. This tilt is due to the low preferred value of $n_{\mathrm{s}}$ for this mode. Despite this global shape difference with $\Lambda$CDM, the enhancement of matter fluctuations on scales entering the horizon at the onset of neutrino free-streaming discussed in subsection IV.2 causes the matter power spectrum to only slightly deviate from the CDM prediction on scales contributing the most to $\sigma_{8}$ ($0.02h\,\mathrm{Mpc}^{-1}\lesssim k\lesssim 0.2h\,\mathrm{Mpc}^{-1}$). For the TT+lens+BAO and TT+lens+BAO+$H_{0}$ data combinations, this difference is less than $5\%$ on these scales and results in a $\sigma_{8}$ value that is slightly lower than in $\Lambda$CDM, potentially bringing low-redshift measurements of the amplitude of matter fluctuations in agreement with CMB data as discussed in subsection VI.4.

Nevertheless, it is important to note that the localized feature in the matter power spectrum caused by the late onset of neutrino free-streaming nearly coincides with the BAO scale, that is, it is on scales where we have a large amount of data from, for example, spectroscopic galaxy surveys (see e.g. Ref. Anderson et al. (2014)). While an analysis that takes into account the full shape of the measured galaxy power spectrum at these scales is beyond the scope of this work, we note that both the SI$\nu$ and MI$\nu$ cosmologies only mildly deviate from the $\Lambda$CDM model near the BAO scale. On smaller scales, the $\mathrm{SI}\nu$ mode displays a net suppression of power which has implications for probes of small-scale structure such as the Lyman-$\alpha$ forest Iršič et al. (2017) and the satellite galaxy count surrounding the Milky Way Kim et al. (2018). It is an interesting possibility that the SI$\nu$ cosmology could help alleviate the small-scale structure problems Bullock and Boylan-Kolchin (2017) without introducing a nongravitational coupling between neutrinos and dark matter.

The MI$\nu$ residuals (lower panel) in 10(a) display an even richer structure than those shown in Fig. 3. Indeed, even in the case of relatively weak neutrino interactions, their impact on the matter power spectrum is significant, and potentially provide a different channel to constrain new physics in the neutrino sector. Since the dominant constraining power of the data used here comes from $k\sim 0.1\,h\,\mathrm{Mpc}^{-1}$, we observe that the $\mathrm{MI}\nu$ power spectra have values similar to $\Lambda\mathrm{CDM}$ near this scale. Outside the scales probed by $\sigma_{8}$, the linear matter power spectra deviate more significantly (up to $\sim 20\%$) from $\Lambda$CDM.

The lensing potential power spectrum in 10(b) shows a similar pattern to the matter power spectrum for the different best-fit models, as expected. The current large error bars of the Planck lensing measurements allow substantial freedom to the SI$\nu$ and MI$\nu$ cosmologies. As shown in Table 4, the lensing data prefer the MI$\nu$ mode for all data-set combinations, but we note that the $\mathrm{SI}\nu$ modes are typically within the error bars of the lensing data.

IX Conclusions

The presence of yet-unknown neutrino interactions taking place in the early Universe could delay the onset of neutrino free-streaming, imprinting the CMB and probes of matter clustering with distinct features. We have performed a detailed study of the impact of neutrino self interactions with a rate scaling as $\Gamma_{\nu}\sim G_{\rm eff}^{2}T_{\nu}^{5}$ on the CMB and the matter power spectrum, taking into account the presence of nonvanishing neutrino masses and of a nonstandard neutrino thermal history. Using recent measurements of the BAO scale, the local Hubble rate, and of the CMB, we find that a cosmological scenario (originally pointed out in Ref. Cyr-Racine and Sigurdson (2014)) in which the onset of neutrino free-streaming is delayed until close to the epoch of matter-radiation equality can provide a good fit to CMB temperature data while also being consistent with the Hubble constant inferred from the local distance ladder Riess et al. (2016).

This strongly interacting neutrino cosmology has the following properties:

• •

Using the data combination TT+lens+BAO+$H_{0}$, it displays a strong preference ($>3\sigma$) for an addition neutrino species ($N_{\mathrm{eff}}=4.02\pm 0.29$, $68\%$ C.L.). This can have important implications given the current anomalies in neutrino oscillation experiments. It also prefers a nonvanishing value of the sum of neutrino masses $\sum m_{\nu}=0.42^{+0.17}_{-0.20}$ eV ($68\%$ C.L.).

• •

It can easily accommodate a larger value of $H_{0}$ and smaller $\sigma_{8}$, hence possibly alleviating tensions between current measurements. Quantitatively, the data combination TT+lens+BAO+$H_{0}$ favors $H_{0}=72.3\pm 1.4\,\mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1}$ and $\sigma_{8}=0.786\pm 0.020$ at $68\%$ C.L.

It is remarkable that a cosmological model admitting parameter values that are so different (see Fig. 5) than in the standard $\Lambda$CDM paradigm can provide a better fit to the data at a statistically-significant level ($\Delta$AIC$=-1.91$). We believe that this is the most important lesson to be drawn from our work: While most analyses have focused on mild deformation from the standard $\Lambda$CDM scenario in trying to reconcile the current cosmological datasets, it is important to entertain the possibility that a radically different scenario (i.e. statistically disjoint in cosmological parameter space) could provide a better global fit to the data.

Despite the success of the strongly interacting neutrino cosmology in addressing tensions between certain cosmological data sets, there are several important obstacles that still tilt the balance towards the standard $\Lambda$CDM cosmology. First, the addition of polarization data seems to degrade the quality of the fit for the strongly interacting neutrino cosmology. We have traced this deterioration of the fit to our use of a Gaussian prior on the reionization optical depth from Ref. Adam et al. (2016). This prior was utilized as a way to capture the constraint on the optical depth from low-$\ell$ HFI Planck polarization data before the full likelihood is made available. It it likely that the Gaussian form of the prior leads to constraints that are too strong as compared to what the full likelihood will provide. Only a complete analysis with the legacy Planck data, once available, will allow us to determine whether this is the case. An important fact to keep in mind is that Ref. Lancaster et al. (2017) found that the addition of CMB polarization data (without an additional $\tau$ prior) tends to increase the statistical significance of the strongly interacting neutrino cosmology.

Second, the low values of the Bayes factor (see Table 4) consistently favor either very weakly interacting neutrinos or no interaction at all. This reflects the fact that strongly interacting neutrinos can only fit the data better for a narrow window of interaction strengths, while $\Lambda$CDM provides a decent (but overall less good) fit over a broader part of the parameter space. This is a fundamental feature of Bayesian statistics and it is unlikely to change in future analyses. This highlights the need to consider a portfolio of statistical measures to assess the quality of a given cosmological model.

Third, it might be difficult from a particle model-building perspective to generate neutrino self interactions with the strength required by the strongly interacting neutrino cosmology while not running afoul of other constraints on neutrino physics. A viable model might look similar to that presented in Ref. Cherry et al. (2014), but it remains to be seen whether the necessary large interaction strength can be generated while evading current constraints Lessa and Peres (2007) on new scalar particles coupling to Standard Model neutrinos. It is also possible that a successful self-interacting neutrino model could have a different temperature dependence than that considered in this work ($\Gamma_{\nu}\propto T_{\nu}^{5}$). This would change the shape of the neutrino visibility function (see Refs. Cyr-Racine and Sigurdson (2014); Lancaster et al. (2017)) and potentially improve the global fit to the data. We leave the study of different temperature scalings of the neutrino interacting rate to future works.

Our analysis could be improved in a few different ways. Given the computational resources we had at our disposal and the need to obtain accurate values of the Bayesian evidence, we used the “lite” version of the Planck high-$\ell$ likelihoods in our analysis. Since some of the assumptions that went into generating these likelihoods Planck Collaboration et al. (2015) might not apply to the interacting neutrino cosmologies, it would be interesting, given sufficient computing power, to reanalyze these models with the complete version of the likelihoods that include all the nuisance parameters. In particular, it is possible that some of the foreground nuisance parameters might be degenerate with the effect of self-interacting neutrinos. For simplicity, we have also assumed that the helium fraction is determined by the standard big-bang nucleosynthesis calculation throughout our analysis. Given the new physics and the resulting modified thermal history of the neutrino sector for the type of models we explore here, it reasonable to assume that the helium fraction would in general be different than in $\Lambda$CDM. While the details of the helium production within any interacting neutrino model are likely model-dependent, a sensible way to take these effects into account would be to let the helium fraction float freely in the fit to CMB data. We leave such analysis to future works.

Given the structure of the residuals between the best-fit interacting neutrino cosmologies and the $\Lambda$CDM model presented in Sec. VIII, it is clear that future high-$\ell$ CMB polarization and matter clustering measurements will play an important role in constraining or ruling out these models (see e.g. Park et al., prep). In particular, the overall red tilt of the matter power spectrum in the strongly interacting neutrino cosmology could have important consequences on both large and small scales. Since current anomalies in terrestrial neutrino experiments Aguilar-Arevalo et al. (2018, 2001) may indicate the presence of new physics in the neutrino sector, it is especially timely to use the complementary nature of cosmological probes to look for possible clues about physics beyond the Standard Model.

Acknowledgements.

We thank Kris Sigurdson and Roland de Putter for collaboration at early stages of this work. We also thank Jo Dunkley, David Spergel, and Lyman Page for comments on an early version of this manuscript, and Prateek Agrawal and David Pinner for useful conversations. C. D. K. acknowledges the support of the National Science Foundation award number DGE1656466 at Princeton University and of the Minority Undergraduate Research Fellowship at the Jet Propulsion Laboratory. F.-Y. C.-R. acknowledges the support of the National Aeronautical and Space Administration (NASA) ATP grant NNX16AI12G at Harvard University. This work was performed in part at the California Institute of Technology for the Keck Institute for Space Studies, which is funded by the W. M. Keck Foundation. Part of the research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.

Appendix A Results for all data sets

We display in Table 6 and Table 7 the $68\%$ confidence limits for the strongly interacting and moderately interacting neutrino modes, respectively. In Fig. 11, we show the marginalized posteriors for key cosmological parameters for a choice of smoothing kernel that represents more accurately the shape of the SI$\nu$ mode. In Fig. 12, we compare the marginalized posterior distribution of the SI$\nu$ mode for the four data set combinations considered in this work.

Appendix B Perturbation equations for interacting massive neutrinos

In this appendix, we derive the Boltzmann equation governing the evolution of the distribution function of massive self-interacting neutrinos which we denote by $f_{\nu}(\mathbf{x},{\bf P},\tau)$, where ${\bf P}$ is the canonical conjugate variable to the position $\mathbf{x}$, and $\tau$ is the conformal time. In the scenario considered here, neutrinos can exchange energy and momentum via 2-to-2 scattering of the type $\nu_{i}+\nu_{j}\rightarrow\nu_{k}+\nu_{l}$. The Boltzmann equation of neutrino species $i$ can be written as

[TABLE]

where $\lambda$ is an affine parameter that described the trajectory of the observer (see below) and $C_{\nu_{i}+\nu_{j}\rightarrow\nu_{k}+\nu_{l}}$ is the collision term for the process $\nu_{i}+\nu_{j}\rightarrow\nu_{k}+\nu_{l}$. In the conformal Newtonian gauge, the space-time metric takes the form

[TABLE]

where $a$ is the cosmological scale factor and $\phi$ and $\psi$ are the two gravitational potentials. We can define the affine parameter in terms of the four-momentum $P$ of an observer

[TABLE]

where $x=(\tau,\vec{x})$ is a four-vector parametrizing the trajectory of the observer. Using Eq. (24), we can then write

[TABLE]

where we have used the dispersion relation $g_{\mu\nu}P^{\mu}P^{\nu}=-m_{\nu}^{2}$. Here, we have defined $E=\sqrt{p^{2}+m_{\nu}^{2}}$, where $p=|\mathbf{p}|$ is the proper momentum, which is related to the conjugate momentum ${\bf P}$ via the relation $p^{2}=g_{ij}P^{i}P^{j}$. We note that Eq. (26) is valid to first-order in perturbation theory. As in other work in the literature, we choose to write the distribution function in terms of the proper momentum $\mathbf{p}$. This choice is valid as long as we also modify the phase-space volume element as $d^{3}P\rightarrow a^{3}(1-3\phi)d^{3}p$. The left-hand side of the Boltzmann equation takes the form Ma and Bertschinger (1995); Dodelson (2003)

[TABLE]

where $\mathcal{H}\equiv\dot{a}/a$ is the conformal Hubble parameter, a overhead dot denoting a derivative with respect to conformal time. We expand the neutrino distribution function as

[TABLE]

At early times, neutrinos form a relativistic tightly-coupled fluid with an equilibrium background distribution function $f_{\nu}^{(0)}(\mathbf{p},\tau)$ given by the Fermi-Dirac distribution. If the interactions mediated by the Lagrangian in Eq. (1) go out of equilibrium while neutrinos are relativistic, the background distribution function would maintain this shape, with a temperature red shifting as $T_{\nu}\propto a^{-1}$. As mentioned in Sec. III, we work under this approximation here and assume that the background distribution function maintains its equilibrium shape throughout the epoch of neutrino decoupling. In the absence of energy source or sink, and for the type of interaction we consider in this work, this is an excellent approximation Oldengott et al. (2017).

Substituting Eq. (28) in Eq. (27) and keeping terms that are first order in the perturbation variables, we obtain

[TABLE]

where the superscript $C^{(1)}_{\nu_{i}}$ denotes the part of the collision term that is first order in the perturbation variables $\Theta_{\nu_{i}}$. It is useful at this point to introduce the comoving momentum $\mathbf{q}\equiv a\mathbf{p}$ and comoving energy $\epsilon\equiv aE$. Going to Fourier space, Eq. (29) becomes

[TABLE]

where $\tilde{\Theta}_{\nu_{i}}$ is the Fourier transform of the perturbation variable $\Theta_{\nu_{i}}$, $\mathbf{k}$ is the Fourier conjugate of $\mathbf{x}$, $k=|\mathbf{k}|$, $\mu\equiv\hat{q}\cdot\hat{k}$, and $\hat{k}=\mathbf{k}/k$. In this work, we focus on (helicity) scalar perturbations and expand the angular dependence of the $\tilde{\Theta}_{\nu_{i}}$ variable in Legendre polynomials $P_{l}(\mu)$

[TABLE]

We note that this decomposition is always valid for scalar perturbations since they must be azimuthally symmetric with respect to $\mathbf{k}$, independently of the structure of the collision term. Substituting the above expansion in the first-order Boltzmann equation and and integrating both sides with $\frac{1}{2(-i)^{l}}\int_{-1}^{1}d\mu P_{l}(\mu)$ yields the hierarchy of equations

[TABLE]

where $\delta_{ij}$ is the Kroenecker delta and where we have suppressed the arguments of $f^{(0)}_{\nu_{i}}$ and $\theta_{l}$ for succinctness. We now turn our attention to the collision integral.

Appendix C Collision Integrals

We now compute the first-order collision term for neutrino scattering, $\nu_{i}(\mathbf{p}_{1})+\nu_{j}(\mathbf{p}_{2})\leftrightarrow\nu_{k}(\mathbf{p}_{3})+\nu_{l}(\mathbf{p}_{4})$. We start from the general expression Kolb and Turner (1990)

[TABLE]

where $|\mathcal{M}|_{\nu}^{2}$ here is the spin-summed (not averaged) matrix element for the scattering as defined in Eq. (II), $\mathbf{p}_{i}$ denotes the $i^{th}$ three-momentum, $p_{i}=|\mathbf{p}_{i}|$, and where

[TABLE]

and

[TABLE]

Using Eq. (28) and keeping only the first order term in the perturbation variable $\Theta_{\nu}$, we can rewrite the collision term as:

[TABLE]

where we have suppressed the ${\bf x}$ and $\tau$ dependence of the $\Theta_{\nu}$ variables to avoid clutter, and where we have used the symmetry $\mathbf{p}_{3}\leftrightarrow\mathbf{p}_{4}$ to simplify the integrand. Here, we have taken the background neutrino distribution function to have a relativistic Fermi-Dirac shape. As mentioned in Sec. III, we assume neutrinos decouple in the relativistic regime and thus neglect the small neutrino mass in the computation of the collision integrals, $g_{\mu\nu}P^{\mu}P^{\nu}\approx 0$ and $E\approx p$. We use the technique developed in Refs. Yueh and Buchler (1976); Hannestad and Madsen (1995); Oldengott et al. (2015) to perform the majority of the integrals. We first perform the $p_{4}$ integration using the identity

[TABLE]

where $H(x)$ is the Heaviside step function. The collision term then reduces to:

[TABLE]

where we have gathered all the terms dependent on the distribution functions inside $\tilde{F}$. To make progress, we have to choose a coordinate system. We take $\mathbf{p}_{1}$ to point in the z-direction, and $\mathbf{p}_{3}$ to lie in the x-z plane. and define the following angles:

[TABLE]

where $\mathbf{k}$ is the Fourier wavenumber of the perturbations and $\beta$ is the azimuthal angle for $\mathbf{p}_{2}$ to wrap around $\mathbf{p}_{1}$. The integration measure than takes the form

[TABLE]

The $\phi$ angle is the azimuthal angle for $\mathbf{p}_{3}$ to wrap around $\mathbf{p}_{1}$. Since we are only dealing with scalar perturbations here, we are free to redefine this angle at will since no physical quantity depends on it. Within this coordinate system, we can write $P_{1}\cdot P_{2}=-p_{1}p_{2}+p_{1}p_{2}\cos{\alpha}$, $P_{1}\cdot P_{3}=-p_{1}p_{3}+p_{1}p_{3}\cos{\theta}$, and $P_{2}\cdot P_{3}=-p_{2}p_{3}+p_{2}p_{3}(\cos{\alpha}\cos{\theta}+\sin{\alpha}\sin{\theta}\cos{\beta})$. We can now use the delta function to do the $\beta$ integration. Setting the argument of the delta function to zero and solving for $\cos{\beta}$ yields

[TABLE]

Performing the integration introduces a Jacobian

[TABLE]

In the above, $a_{\alpha}$, $b_{\alpha}$, and $c_{\alpha}$ are

[TABLE]

[TABLE]

[TABLE]

The above form of the collision term is useful when the $\theta$ integration needs to be performed first. In some instances, it will be easier to first perform the $\alpha$ integral first. In this latter case, the collision term can equivalently be written as:

[TABLE]

In the above, $a_{\theta}$, $b_{\theta}$, and $c_{\theta}$ are

[TABLE]

[TABLE]

[TABLE]

We remark that using Eq. (41) simplifies $P_{2}\cdot P_{3}$

[TABLE]

We note that since the matrix element for the type of interaction of interest (see Eq. (II)) only depends on a sum of Mandelstam variables or squares of Mandelstam variables, we can write

[TABLE]

depending on which of the $\theta$ or $\alpha$ integral we want to perform first. The coefficients are as follow:

[TABLE]

[TABLE]

We note that we can perform the $\cos{\theta}$ or the $\cos{\alpha}$ integration using the following results:

[TABLE]

[TABLE]

[TABLE]

To make further progress in performing the angular integration, we need to specify the angular dependence of the $\Theta_{\nu}(\mathbf{p}_{i})$ variables. As in Eq. (31), we expand their angular dependence in Legendre polynomials with respect to the angle between the vector $\mathbf{k}$ and the $\mathbf{p}_{i}$ vectors. Within our coordinate system, these angles are

[TABLE]

The following identity will be useful later in order to perform the remaining azymuthal integral (the “$\phi$” integral)

[TABLE]

We now consider separately the different terms in the perturbative expansion in $\Theta_{\nu}(\mathbf{p}_{i})$.

C.1 Terms involving $\Theta_{\nu}(\mathbf{p}_{1})$

This is the simplest case since $\Theta_{\nu}(\mathbf{p}_{1})$ can be carried outside the integrals. The azymuthal $\phi$ integral is trivial and yield an extra factor of $2\pi$.

[TABLE]

Performing the $\alpha$ integration first using Eqs. (54)-(56), we obtain

[TABLE]

Writing $\eta\equiv\cos{\theta}$, we have

[TABLE]

where we have use the definition

[TABLE]

The Heaviside step function $H(b_{\theta}^{2}-4a_{\theta}c_{\theta})$ determines the range of integration of both $\eta$ and $p_{3}$. It yields

[TABLE]

where

[TABLE]

We can then write

[TABLE]

Defining $x_{i}\equiv p_{i}/T_{\nu}$, we obtain

[TABLE]

C.2 Terms involving $\Theta_{\nu}(\mathbf{p}_{2})$

For the term involving $\Theta_{\nu}(\mathbf{p}_{2})$, we start from Eq. (42) and substitute the expansion from Eq. (31).

[TABLE]

where we have used the available freedom to redefine the azymuthal angle $\phi$. We can now perform the $\phi$ integral using the identity given in Eq. (58)

[TABLE]

Performing the $\cos{\theta}$ integral yields

[TABLE]

Similarly to the previous section, we define

[TABLE]

We use the Heaviside step function $H(b_{\alpha}^{2}-4a_{\alpha}c_{\alpha})$ to determine the range of integration for $\rho\equiv\cos{\alpha}$ and $p_{3}$

[TABLE]

where

[TABLE]

We thus obtain

[TABLE]

The remaining difficulty is the $p_{2}$ dependence of $\theta_{l}(p_{2})$. Since we are working in the thermal approximation in which the only possible neutrino perturbations are local temperature fluctuations, we note that the perturbation variable $\Theta(\mathbf{p}_{2})$ admits the form

[TABLE]

It is therefore convenient to introduce the temperature fluctuation variables (see Eq. (7) in main text) $\nu_{l}$

[TABLE]

which are independent of $\mathbf{p}_{2}$ in the ultra-relativistic limit. We note that as the neutrinos transition to the non-relativistic regime, the $\nu_{l}$ variables will develop a momentum dependence due to the presence of the mass term in the left-hand side of the Boltzmann equations. However, since we expect the neutrinos to self-decoupled in the relativistic regime, we can safely assume that $\nu_{l}$ is independent of the neutrino momentum. Substituting

[TABLE]

in Eq. (C.2) and writing down the answer in terms of $x_{i}$, we obtain

[TABLE]

C.3 Terms involving $\Theta_{\nu}(\mathbf{p}_{3})$

For the term involving $\Theta_{\nu}(\mathbf{p}_{3})$, we begin from Eq. (46) and substitute the expansion from Eq. (31)

[TABLE]

We can now perform the $\phi$ integral using the identity given in Eq. (58)

[TABLE]

Performing the $\cos{\alpha}$ integral yields

[TABLE]

Writing $\eta\equiv\cos{\theta}$ and using the same integration limits as in Eq. (63), we obtain

[TABLE]

Again, writing down the $p_{3}$ dependence of $\theta_{l}(p_{3})$ in terms of the $\nu_{l}$ variables and writing the integrals in terms of the dimensionless variables $x_{i}$, we get

[TABLE]

C.4 Total Collision Term

The complete collision term can then be written as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

where $x_{i}=p_{i}/T_{\nu}$.