Neutrino decoupling beyond the Standard Model: CMB constraints on the Dark Matter mass with a fast and precise $N_{\rm eff}$ evaluation

TL;DR

This paper introduces a fast, accurate method for modeling neutrino decoupling that includes BSM physics and applies it to constrain MeV-scale dark matter using CMB data, improving previous approaches.

Contribution

The authors develop a simplified differential equation approach for neutrino decoupling that accurately incorporates BSM species and finite temperature effects, enabling efficient constraints on dark matter.

Findings

Planck data rules out dark matter particles below 3.0 MeV with certain interactions.

The new method achieves 0.01 accuracy in $N_{eff}$ within the Standard Model.

CMB Stage-IV will test dark matter particles up to 15 MeV.

Abstract

The number of effective relativistic neutrino species represents a fundamental probe of the thermal history of the early Universe, and as such of the Standard Model of Particle Physics. Traditional approaches to the neutrino decoupling are either very technical and computationally expensive, or assume that neutrinos decouple instantaneously. In this work, we aim to fill the gap between these two approaches by modeling the neutrino decoupling in terms of two simple coupled differential equations for the electromagnetic and neutrino sector temperatures, in which all the relevant interactions are taken into account and which allows for a straightforward implementation of BSM species. Upon including finite temperature QED corrections we reach an accuracy on $N_{\rm eff}$ in the SM of $0.01$. We illustrate the usefulness of this approach to the neutrino decoupling by considering, in a model independent manner, the impact of MeV thermal dark matter on $N_{\rm eff}$. We show that Planck rules out electrophilic and neutrinophilic thermal dark matter particles of $m< 3.0\,\text{MeV}$ at 95\% CL regardless of their spin, and of their annihilation being $s$-wave or $p$-wave. We point out that thermal dark matter particles with non-negligible interactions with both electrons and neutrinos are more elusive to CMB observations than purely electrophilic or neutrinophilic ones. In addition, assisted by the accuracy of our approach, we show that CMB Stage-IV experiments will generically test thermal dark matter particles with $m \lesssim 15\,\text{MeV}$. We make publicly available the codes developed for this study at https://github.com/MiguelEA/nudec_BSM .

Figures

Tables

Table 1. Table 1 : Summary of the SM results on the number of effective relativistic neutrino species in the early Universe when different effects are taken into account. These numbers should be compared with the most precise evaluation in the SM: N eff SM = 3.045 superscript subscript 𝑁 eff SM 3.045 N_{\text{eff}}^{\rm SM}=3.045   [ 5 ] . The difference of ∼ 0.01 similar-to absent 0.01 \sim 0.01 results from not accounting for spin-statistics and the electron mass in the interaction rates.

Neutrino Decoupling in the SM	$T_{\gamma }/T_{\nu }$	$N_{\text{eff}}$
Instantaneous decoupling	1.4010	3
Instantaneous decoupling + QED	1.3998	3.011
MB collision term	1.3961	3.042
MB collision term + QED	1.3949	3.053

Table 2. Table 2 : Lower bounds on the masses of various thermal dark matter particles (in MeV) at 95% CL. ν 𝜈 \nu refers to neutrinophilic dark matter particles, and e 𝑒 e refers to electrophilic dark matter particles. Given the current disagreement between local and CMB measurements of the Hubble constant, a conservative bound would be to consider the smallest number from the Planck+BAO and Planck+BAO+ H 0 subscript 𝐻 0 H_{0} columns. Note that the bounds apply to both s 𝑠 s -wave and p 𝑝 p -wave annihilating relics.

Dark Matter		$N_{\mathrm{eff}}$ Constraints
Particle	g-Spin	Planck+BAO		Planck+BAO+$H_0$		Stage-IV Sensitivity
		$>2.66$		$>2.97$		$>2.99$
Majorana DM-$\nu$	2-F	-	8.5	-	6.0	-	14.3
Dirac DM-$\nu$	4-F	-	11.2	-	8.8	-	16.9
$\varphi$ DM-$\nu$	1-B	-	5.7	-	3.0	-	11.7
${\varphi }^{\star }$ DM-$\nu$	2-B	-	8.5	-	6.1	-	14.3
$Z^{\prime }$ DM-$\nu$	3-B	-	10.1	-	7.7	-	15.8
Majorana DM-$e$	2-F	6.4	-	12.8	-	14.1	-
Dirac DM-$e$	4-F	9.2	-	15.5	-	16.7	-
$\varphi$ DM-$e$	1-B	3.5	-	10.2	-	11.5	-
${\varphi }^{\star }$ DM-$e$	2-B	6.5	-	12.9	-	14.1	-
$Z^{\prime }$ DM-$e$	3-B	8.1	-	14.4	-	15.6	-

Table 3. Table 3 : Lower bound on the mass of a generic Majorana WIMP (in MeV) at 95% CL as a function of the annihilation final state ratio to electrons and neutrinos. The upper/lower rows correspond to electrophilic/neutrinophilic WIMPs. The only four slots with a range correspond to an exclusion range. Given the current disagreement between local and CMB measurements of the Hubble constant, a conservative bound will correspond to N eff > 2.66 subscript 𝑁 eff 2.66 N_{\rm eff}>2.66 and N eff < 3.57 subscript 𝑁 eff 3.57 N_{\rm eff}<3.57 . Note that the bounds apply to both s 𝑠 s -wave and p 𝑝 p -wave annihilating relics. Constraints for other branching fractions can be inferred by interpolation or by use of the  online code .

Annihilation Ratio	$N_{\mathrm{eff}}$ Constraints
	Planck		Planck+BAO+$H_0$		Stage-IV Sensitivity
$e:\nu$	$>2.66$		$>2.97$		$>2.99$
$1:{10}^6$	-	8.4	-	5.8	-	14.3
$1:{10}^5$	-	7.9	-	4.8	-	14.2
$1:{10}^4$	-	5.3	-	2.6	-	13.6
$1:{10}^3$	-	3.7	-	2.4	-	11.1
$1:{10}^2$	-	3.8	-	2.7	-	8.7
$1:9$	-	4.3	-	3.1	-	8.7
$9:1$	-	4.0	-	3.0	-	7.4
${10}^2:1$	-	3.1	-	2.4	-	5.4
${10}^3:1$	-	2.1	-	1.6	-	3.2
${10}^4:1$	-	1.0	1.7-11.6	0.8	1.6-13.3	1.3
${10}^5:1$	0.6-5.0	-	0.3-12.7	-	0.3-14.0	0.2
${10}^6:1$	6.2	-	12.8	-	14.1	-

Table 4. Table 4 : Lower bound on the mass of a generic WIMP (in MeV) from the measurements of the primordial Helium abundance Y p subscript 𝑌 𝑝 Y_{p} . The columns correspond to a neutral scalar, charged scalar, Majorana, neutral vector boson, and Dirac thermal dark matter particles. The rows correspond to different annihilation final states.

Annihilation Ratio	$Y_p=0.245\pm 0.006$ (95% CL) Constraints
$e:\nu$	$\text{B},g=1$	$\text{B},g=2$	$\text{F},g=2$	$\text{B},g=3$	$\text{F},g=4$
$1:{10}^6$	1.1	2.4	2.3	3.4	4.0
$1:{10}^5$	1.1	2.0	2.0	2.4	2.7
$1:{10}^4$	1.2	1.8	1.8	2.1	2.3
$1:{10}^3$	1.4	1.8	1.8	2.0	2.2
$1:{10}^2$	1.6	2.0	2.0	2.2	2.3
$1:9$	1.9	2.2	2.2	2.4	2.5
$9:1$	1.8	2.2	2.2	2.3	2.4
${10}^2:1$	1.5	1.9	1.9	2.1	2.2
${10}^3:1$	1.2	1.5	1.5	1.8	1.9
${10}^4:1$	0.7	1.1	1.1	1.4	1.6
${10}^5:1$	0.4	0.7	0.7	1.0	1.2
${10}^6:1$	0.4	0.5	0.5	0.7	0.9

Table 5. Table 5 : Same as Table  3 but for a neutral scalar boson dark matter particle.

Annihilation Ratio	$N_{\mathrm{eff}}$ Constraints
	Planck+BAO		Planck+BAO+$H_0$		Stage-IV Sensitivity
$e:\nu$	$>2.66$		$>2.97$		$>2.975$
$1:{10}^6$	-	5.7	-	3.0	-	11.7
$1:{10}^5$	-	5.4	-	2.5	-	11.6
$1:{10}^4$	-	3.8	-	1.8	-	11.2
$1:{10}^3$	-	3.0	-	1.9	-	9.3
$1:{10}^2$	-	3.3	-	2.3	-	7.5
$1:9$	-	3.8	-	2.8	-	7.8
$9:1$	-	3.6	-	2.7	-	6.7
${10}^2:1$	-	2.8	-	2.1	-	4.8
${10}^3:1$	-	1.7	-	1.3	-	2.7
${10}^4:1$	-	0.7	1.4-9.3	0.6	1.3-10.9	1.0
${10}^5:1$	0.5-2.6	-	10.1	-	11.4	-
${10}^6:1$	3.4	-	10.2	-	11.5	-

Table 6. Table 6 : Same as Table  3 but for a charged scalar boson dark matter particle.

Annihilation Ratio	$N_{\mathrm{eff}}$ Constraints
	Planck+BAO		Planck+BAO+$H_0$		Stage-IV Sensitivity
$e:\nu$	$>2.66$		$>2.97$		$>2.975$
$1:{10}^6$	-	8.4	-	6.0	-	14.3
$1:{10}^5$	-	8.0	-	4.9	-	14.2
$1:{10}^4$	-	5.3	-	2.6	-	13.6
$1:{10}^3$	-	3.7	-	2.4	-	11.1
$1:{10}^2$	-	3.8	-	2.7	-	8.7
$1:9$	-	4.3	-	3.1	-	8.7
$9:1$	-	4.0	-	3.0	-	7.4
${10}^2:1$	-	3.1	-	2.4	-	5.4
${10}^3:1$	-	2.1	-	1.4	-	3.2
${10}^4:1$	-	1.0	1.7-11.7	0.6	1.6-13.3	1.3
${10}^5:1$	0.5-5.3	-	18.0	-	14.0	-
${10}^6:1$	6.3	-	12.9	-	14.1	-

Table 7. Table 7 : Same as Table  3 but for a neutral vector boson dark matter particle.

Annihilation Ratio	$N_{\mathrm{eff}}$ Constraints
	Planck+BAO		Planck+BAO+$H_0$		Stage-IV Sensitivity
$e:\nu$	$>2.66$		$>2.97$		$>2.975$
$1:{10}^6$	-	10.0	-	7.5	-	15.9
$1:{10}^5$	-	9.4	-	6.2	-	15.8
$1:{10}^4$	-	6.2	-	3.1	-	15.0
$1:{10}^3$	-	4.1	-	2.7	-	12.1
$1:{10}^2$	-	4.1	-	2.8	-	9.3
$1:9$	-	4.5	-	3.3	-	9.2
$9:1$	-	4.2	-	3.2	-	7.8
${10}^2:1$	-	3.3	-	2.5	-	5.8
${10}^3:1$	-	2.3	-	1.8	3.7-16.6	3.4
${10}^4:1$	-	1.1	1.9-13.0	0.9	1.8-14.7	1.4
${10}^5:1$	0.5-6.6	-	0.3-14.2	-	0.3-15.5	-
${10}^6:1$	7.9	-	14.4	-	15.6	-

Table 8. Table 8 : Same as Table  3 but for a Dirac dark matter particle.

Annihilation Ratio	$N_{\mathrm{eff}}$ Constraints
	Planck+BAO		Planck+BAO+$H_0$		Stage-IV Sensitivity
$e:\nu$	$>2.66$		$>2.97$		$>2.975$
$1:{10}^6$	-	11.0	-	8.6	-	16.9
$1:{10}^5$	-	10.3	-	7.1	-	16.8
$1:{10}^4$	-	6.9	-	3.5	-	15.9
$1:{10}^3$	-	4.3	-	2.8	-	12.9
$1:{10}^2$	-	4.3	-	3.1	-	9.7
$1:9$	-	4.7	-	3.4	-	9.5
$9:1$	-	4.4	-	3.3	-	8.1
${10}^2:1$	-	3.5	-	2.7	-	6.0
${10}^3:1$	-	2.4	-	1.9	-	3.6
${10}^4:1$	-	0.6	2.0-13.9	1.0	1.9-15.6	1.6
${10}^5:1$	0.6-7.5	-	0.4-15.2	-	0.4-16.6	0.3
${10}^6:1$	9.0	-	15.5	-	16.7	-

Equations

$$N_{\mathrm{eff}}\equiv \frac{8}{7}{\left(\frac{11}{4}\right)}^{4/3}\left(\frac{{\rho }_{\mathrm{rad}}-{\rho }_{\gamma }}{{\rho }_{\gamma }}\right).$$

$$N_{\mathrm{eff}}=3{\left(\frac{11}{4}\right)}^{4/3}{\left(\frac{T_{\nu }}{T_{\gamma }}\right)}^4.$$

$$\frac{\partial f}{\partial t}-Hp\frac{\partial f}{\partial p}=\mathcal{C}[f],$$

$$H=\sqrt{\frac{8\pi }{3}\frac{\rho }{m_{Pl}^2}},$$

$$\mathcal{C}[f_{\psi }]\equiv -\frac{1}{2E_{\psi }}\sum _{X,Y}\int \prod _{i}d{\Pi }_{X_i}\prod _{j}d{\Pi }_{Y_j}(2\pi {)}^4{\delta }^4(p_{\psi }+p_X-p_Y)\times$$

$$\left[|\mathcal{M}{|}_{\psi +X\to Y}^2{f}_{\psi }\prod _i{f}_{X_i}\prod _j\left[1\pm f_{Y_j}\right]-|\mathcal{M}{|}_{Y\to \psi +X}^2\prod _j{f}_{Y_j}\left[1\pm f_{\psi }\right]\prod _i\left[1\pm f_{X_i}\right]\right],$$

$$\frac{dn}{dt}+3Hn$$

$$\frac{d\rho }{dt}+3H(\rho +p)$$

$$\frac{d{T}_{\nu }}{dt}$$

$$\frac{d{T}_{{\nu }_e}}{dt}$$

$$\frac{d{T}_{{\nu }_{\mu }}}{dt}$$

$${\dot{\rho }}_{\mathrm{Tot}}=-3H({\rho }_{\mathrm{Tot}}+p_{\mathrm{Tot}}).$$

$$\frac{d{T}_{\gamma }}{dt}=-\frac{4H{\rho }_{\gamma }+3H\left({\rho }_e+p_e\right)+\frac{\delta {\rho }_{{\nu }_e}}{\delta t}+2\frac{\delta {\rho }_{{\nu }_{\mu }}}{\delta t}}{\frac{\partial {\rho }_{\gamma }}{\partial T_{\gamma }}+\frac{\partial {\rho }_e}{\partial T_{\gamma }}}.$$

$${\frac{\delta {\rho }_{{\nu }_e}}{\delta t}|}_{\mathrm{SM}}^{\mathrm{MB}}$$

$${\frac{\delta {\rho }_{{\nu }_{\mu }}}{\delta t}|}_{\mathrm{SM}}^{\mathrm{MB}}$$

$$F_{\mathrm{MB}}(T_1,T_2)=32(T_1^9-T_2^9)+56T_1^4{T}_2^4(T_1-T_2).$$

$$\frac{d{T}_{\gamma }}{dt}=-\frac{4H{\rho }_{\gamma }+3H\left({\rho }_e+p_e\right)+3H{T}_{\gamma }\frac{d{P}_{\text{int}}}{d{T}_{\gamma }}+\frac{\delta {\rho }_{{\nu }_e}}{\delta t}+2\frac{\delta {\rho }_{{\nu }_{\mu }}}{\delta t}}{\frac{\partial {\rho }_{\gamma }}{\partial T_{\gamma }}+\frac{\partial {\rho }_e}{\partial T_{\gamma }}+T_{\gamma }\frac{d^2{P}_{\text{int}}}{d{T}_{\gamma }^2}},$$

$$P_{\text{int}}=-{\int }_0^{\infty }\frac{dp}{2{\pi }^2}\left[\frac{p^2}{\sqrt{p^2+m_e^2}}\frac{\delta m_e^2(T)}{e^{\sqrt{p^2+m_e^2}/T}+1}+\frac{p}{2}\frac{\delta m_{\gamma }^2(T)}{e^{p/T}-1}\right],$$

$$\delta m_{\gamma }^2(T)$$

$$\frac{d{T}_{\nu }}{dt}$$

$${\frac{d}{d{T}_{\gamma }}\left[\frac{\delta {\rho }_{\nu }}{\delta t}\right]|}_{T_{\gamma }=T_{\nu }}$$

$$T_{{\nu }_e}^{\text{dec}}$$

$$N_{\mathrm{eff}}=2.99_{-0.33}^{+0.34}(95\%\text{CL, TT+TE+EE+lowE+lensing+BAO}).$$

$$N_{\mathrm{eff}}$$

$$\frac{d{T}_{\nu }}{dt}$$

$$\frac{d{T}_{\gamma }}{dt}=-\frac{4H{\rho }_{\gamma }+3H\left({\rho }_e+p_e\right)+3H\left({\rho }_{\chi }+p_{\chi }\right)+3H{T}_{\gamma }\frac{d{P}_{\text{int}}}{d{T}_{\gamma }}+\frac{\delta {\rho }_{{\nu }_e}}{\delta t}+2\frac{\delta {\rho }_{{\nu }_{\mu }}}{\delta t}}{\frac{\partial {\rho }_{\gamma }}{\partial T_{\gamma }}+\frac{\partial {\rho }_e}{\partial T_{\gamma }}+\frac{\partial {\rho }_{\chi }}{\partial T_{\gamma }}+T_{\gamma }\frac{d^2{P}_{\text{int}}}{d{T}_{\gamma }^2}},$$

$$\frac{d{T}_{\nu }}{dt}$$

$$\frac{d{T}_{\gamma }}{dt}$$

$$\frac{d{T}_{\nu }}{dt}$$

$$\frac{d{n}_{\chi }}{dt}+3H{n}_{\chi }$$

$$\frac{\delta {\rho }_{\chi }}{\delta t}$$

Full text

1 Introduction

The number of effective relativistic neutrino species in the early Universe ($N_{\rm eff}$) constitutes a fundamental probe of the thermal history of the early Universe, and as such of the Standard Model of Particle Physics (SM). The agreement between the observations of the Cosmic Microwave Background (CMB) [1, 2], and the primordial element abundances [3, 4] with the SM prediction [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15] of three relativistic neutrino species represents a success of the Standard Model, and constitues a non-trivial test for some of its extensions in which light or massless species are present [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. In addition, CMB Stage-IV experiments [28] are expected to make this test considerably more stringent by achieving sensitivities to $N_{\rm eff}$ of $0.03$, and therefore improving by a factor of $5-6$ current measurements by the Planck satellite [1, 2].

Given the usefulness of the number of effective neutrino species as a tool to test the SM and its extensions, the neutrino decoupling process in the early Universe has been a subject of intense study [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. The state of the art modeling of neutrino decoupling [5, 6] includes all the relevant physics: exact collision terms, neutrino oscillations, and finite temperature QED corrections, which results in the SM prediction of $N_{\rm eff}^{\rm SM}=3.045$ [5]. However, although the state of the art treatment of neutrino decoupling in the SM is very accurate, it is considerably involved and computationally expensive since it makes use of the density matrix formalism and requires one to solve a very complex system of kinetic equations111The system consists on 200 stiff coupled integro-differential equations as a result of the binning (in momentum space) of the neutrino distribution functions [5, 6, 11].. On the other hand, studies that have focused on possible BSM implications for $N_{\rm eff}$ with species at the MeV scale have assumed that neutrinos decouple instantaneously and modeled the temperature evolution by assuming entropy conservation (see e.g. [29, 30, 31, 20, 21, 22, 23, 32, 33]). To the author’s knowledge, the sole exception to this approach is [34], where a BSM entropy transfer rate between neutrinos and electrons was considered but where the SM electron-neutrino interactions were not taken into account.

In this study, we carry out a first-principles calculation of $N_{\rm eff}$, in and beyond the Standard Model. We consider all relevant $\nu$-$e$ and $\nu$-$\nu$ interactions, and by assuming thermal equilibrium distribution functions for all the species, negligible chemical potentials, and Maxwell-Boltzmann statistics in the collision terms, we are able to model the neutrino decoupling process in terms of two simple coupled differential equations for $T_{\gamma}$ and $T_{\nu}$. We show that this approach has an accuracy of $0.01$ for $N_{\rm eff}$ in the SM, and what is more important: that it allows for a straightforward implementation of BSM states without the need of any a priori knowledge of the time at which neutrinos decouple (which is generically altered by BSM interactions). Therefore, this approach represents a simple – and yet accurate – treatment of neutrino decoupling in BSM theories. These two characteristics, accuracy and the ease of BSM species implementation, are crucial in order to make precise theoretical predictions for $N_{\rm eff}$ and the primordial element abundances as synthesized during Big Bang Nucleosynthesis (BBN). In particular, given the upcoming sensitivity of CMB Stage-IV experiments [28].

In this work, we illustrate the usefulness of our simplified approach to neutrino decoupling by revisiting the effect of MeV thermal dark matter (i.e. WIMPs [35, 36, 37]) on $N_{\rm eff}$ in a model independent manner. By comparing the theoretical predictions for $N_{\rm eff}$ from our approach with the measured values by the Planck satellite [1, 2] we set bounds on the masses of thermal dark matter particles with $s$-wave and $p$-wave annihilations to electrons and neutrinos. We find that current CMB observations rule out purely electrophilic and neutrinophilic thermal relics with masses $m<3.0\,\text{MeV}$ at 95% CL. In addition, we point out that thermal dark matter particles that have non-negligible interactions with both electrons and neutrinos are generically more elusive to CMB observations.

This paper is organized as follows. First, in Section 2 we provide a detailed derivation of the temperature evolution equations that model the relic neutrino decoupling process in the early Universe. This section includes justifications for the various approximations that allow us to capture the physics of the neutrino decoupling process in terms of just two simple differential equations. In Section 3, we review the current and upcoming constraints from CMB observations on $N_{\rm eff}$. In Section 4, we work out the time evolution of neutrino decoupling in the presence of MeV scale thermal dark matter particles and set constraints on their masses from CMB observations on $N_{\rm eff}$. In Section 5, we comment on the complementary bounds that could be inferred from BBN on the dark matter mass. We conclude in Section 6.

2 Neutrino Decoupling in the Early Universe

2.1 Simplifying approximations for Neutrino Decoupling

Our aim in this study is not to reproduce the very precise value of $N_{\rm eff}^{\rm SM}=3.045$ as obtained in Reference [5], but to formulate neutrino decoupling in a simplified form that allows one to capture the time dependence of the process while accounting for all possible interactions that can alter it. Which, in addition, thanks to several well justified approximations, allows for a fast, yet reliable and accurate, computation of $N_{\text{eff}}$ – as relevant for CMB observations – in the SM and beyond. Our formulation will make use of the integrated Boltzmann equation, and makes the following approximations:

• Approximation: the relevant particles follow thermal equilibrium distribution functions i.e.: Fermi-Dirac or Bose-Einstein for fermions and bosons respectively.
• Advantage: the distribution function of any species is simply determined by two parameters: its temperature $T$ and its chemical potential $\mu$.
• Justification: we will be considering particles that at high temperatures are in thermal equilibrium with the SM, and therefore corrections to the equilibrium distributions could only occur when the particles decouple from the plasma. The decoupling happens when the interaction rates are small, and therefore non-thermal spectral distortions are expected to be small. For the case of neutrinos in the SM, the non-thermal corrections amount to $\delta\rho_{\nu}/\rho_{\nu}<1\,\%$ [5, 6]. However, we note that if the neutrinos possess any other stronger BSM interactions, the non-thermal corrections may be even smaller. In addition, the electrons, positrons and photons present negligible non-thermal corrections due to the strong electromagnetic interactions.

• Approximation: Maxwell-Boltzmann statistics in the collision terms.
• Advantage: allows to fully integrate the 6 dimensional collision term for $1\leftrightarrow 2$ processes, and reduces the 9 dimensional phase space of $2\leftrightarrow 2$ annihilation processes to just 1.
• Justification: the statistical factors in the collision terms modify the rates by less than $15\%$ for thermal distribution functions (see e.g. [38]).

• Approximation: neglect the electron mass in the $e^{+}e^{-}\leftrightarrow\bar{\nu}_{i}\nu_{i}$ and $e^{\pm}\nu_{i}\leftrightarrow e^{\pm}\nu_{i}$ processes.
• Advantage: together with Approximations 1 & 2, permits the full integration of the collision terms, namely, one can write an analytic expression for the energy transfer rates.
• Justification: since the electrons are still relativistic when the electron-neutrino interactions freeze-out ($T\sim 2-3\,\text{MeV}$), the error of neglecting the electron mass in the rates is less than a few percent for the relevant temperatures [39].

• Approximation: neglect chemical potentials.
• Advantage: reduces the number of evolution equations.
• Justification: $\mu_{\gamma}=\mu_{e}=0$ is an excellent assumption due to the very strong scattering and annihilations in the electromagnetic sector. $\mu_{\nu_{e}}=\mu_{\nu_{\mu}}=\mu_{\nu_{\tau}}=0$ is also a fairly good approximation since the $e^{+}e^{-}\leftrightarrow\bar{\nu}_{i}\nu_{i}$ and $\bar{\nu}_{j}\nu_{j}\leftrightarrow\bar{\nu}_{i}\nu_{i}$ process damp any possible chemical potentials in the neutrino sector for the temperature range of interest for the process of neutrino decoupling. In addition, we will only consider scenarios in which no primordial leptonic asymmetry exists in the early Universe, i.e. $\mu_{\nu}=\mu_{\bar{\nu}}=0$.

• Approximation: the time evolution of neutrino decoupling will be characterized by just two variables $T_{\nu}=T_{\nu_{e}}=T_{\nu_{\mu}}=T_{\nu_{\tau}}$ and $T_{\gamma}=T_{e}$.
• Advantage: makes the time evolution of neutrino decoupling extremely simple.
• Justification: $T_{\gamma}=T_{e}$ is a very good approximation due to the strong scattering and annihilations between electrons, positrons and photons. $T_{\nu_{\mu}}=T_{\nu_{\tau}}$ is also a good assumption because these neutrinos have the same interactions in the SM. Furthermore, given the observed pattern of masses and mixings [40, 41, 42] the neutrino oscillations become active in the early Universe for temperatures $T=3-5\,\text{MeV}$ [43, 44] (see Appendix A.1 for more details). The neutrino oscillations act as to equilibrate the flavor distribution functions [44]. Thus, given the fact that the $\nu_{\mu}$ and $\nu_{\tau}$ are – or are very close to be – maximally mixed [40, 41, 42], $T_{\nu_{\mu}}=T_{\nu_{\tau}}$ is also a very good approximation, even in BSM scenarios. In addition, since $\theta_{12},\,\theta_{13},\,\theta_{23}\neq 0$, $T_{\nu_{e}}\simeq T_{\nu_{\mu}}\simeq T_{\nu_{\tau}}$ will be an approximate but almost realistic scenario as a consequence of the neutrino oscillations in the early Universe. Finally, we have explicitly checked that taking $T_{\nu_{e}}\neq T_{\nu_{\mu,\,\tau}}$ results in very similar values for $N_{\rm eff}$, and although in this work we will assume a common neutrino temperature, the online code allows for a different temperature evolution for each neutrino flavor.

Collectively, these approximations are well justified and, as we shall see, hugely simplify the time evolution of neutrino decoupling, reducing the system to just two ($dT_{\nu}/dt,\,dT_{\gamma}/dt$) ordinary differential equations. In what follows, we provide the derivation of the time evolution equations for the neutrino and electromagnetic sector temperatures in the early Universe.

2.2 Modeling of Neutrino Decoupling

Assisted by the well justified approximations adopted in Section 2.1 we can write down two simple evolution equations for $T_{\gamma}$ and $T_{\nu}$ that capture the time dependence of the process of neutrino decoupling and that take into account all possible interactions. This will allow for a fast, yet precise, computation of the effective number of relativistic neutrino species as relevant for CMB observations:

[TABLE]

Which, when only neutrinos and photons are present, and their distribution functions can be characterized by their temperatures, reads:

[TABLE]

2.2.1 Derivation of the temperature evolution equations

The starting point of our derivation is the Liouville equation in a flat, isotropic and homogeneous expanding Universe:

[TABLE]

where $f$ is the distribution function of a given species, $p$ is the momentum, $H$ is the Hubble rate:

[TABLE]

where $m_{Pl}=1.22\times 10^{19}\,\text{GeV}$, $\rho$ is the total energy density of the Universe, and $\mathcal{C}[f]$ is the collision term that takes into account any potential departure of $f$ from the equilibrium distribution $f^{\rm eq}$ as a result of any scattering/annihilation/decay process. Explicitly, the collision term for a species $\psi$ reads [7]:

[TABLE]

where $f_{\psi},f_{X_{i}},f_{Y_{j}}$ are the distribution functions of the particles $\psi,\,X_{i}$ and $Y_{j}$ respectively. The $-$ sign applies to fermions and the $+$ sign applies to bosons. $d\Pi_{X_{i}}\equiv\frac{g_{X_{i}}}{(2\pi)^{3}}\frac{d^{3}p}{2E}$ represents the phase space of a particle $X_{i}$ with $g_{X_{i}}$ internal degrees of freedom. The delta function ensures energy-momentum conservation, and $\mathcal{M}_{\psi+X\to Y}$ is the amplitude for the $\psi+X\to Y$ process.

We can integrate Equation (2.3) over momenta with $g\,d^{3}p/(2\pi)^{3}$ and $g\,E\,d^{3}p/(2\pi)^{3}$ as a measure to obtain:

[TABLE]

where $n$, $\rho$ and $p$ are the number, energy, and pressure densities of the given species respectively. $g$ represents the number internal of degrees of freedom. $\delta n/\delta t$ represents the number density transfer rate, and $\delta\rho/\delta t$ is the energy density transfer rate. Since we neglect chemical potentials, we can write down the temperature evolution of any species just in terms of the energy density evolution (2.7). Assuming $T_{\nu}=T_{\nu_{e}}=T_{\nu_{\mu}}=T_{\nu_{\tau}}$, and by using the chain rule, we can write down the temperature evolution equation for the neutrino fluid:

[TABLE]

where the first term in the r.h.s. represents the effect of the Universe’s expansion, and the second term takes into account non-adiabatic contributions due to the energy exchange between neutrinos and other particles (in the SM, at $T<10\,\text{MeV}$, electrons and positrons). If $T_{\nu_{e}}\neq T_{\nu_{\mu}}$, the temperature evolution equations read:

[TABLE]

We can obtain the evolution equation for the photon temperature from the continuity equation that relates the net change in the total energy density to the Universe’s expansion, the total energy density, and the total pressure [45]:

[TABLE]

We note that Equation (2.10) can be derived both from the fluid perspective ($\nabla_{\mu}T^{\mu}{}_{\nu}=0$) [45] and from the microphysics perspective by summing Equation (2.7) over all the species present in the Universe, since $\sum_{i}\int g_{i}\,E_{i}\,d^{3}p_{i}\,\mathcal{C}[f_{i}]\equiv 0$ (but note that in general $\sum_{i}\int g_{i}\,d^{3}p_{i}\,\mathcal{C}[f_{i}]\neq 0$).

Finally, from Equations (2.10) and (2.7) for the neutrinos we can, by using the chain rule, write down the temperature evolution for the electromagnetic plasma:

[TABLE]

Where again, the first terms in the r.h.s. ($4H\rho_{\gamma}+3H\left(\rho_{e}+p_{e}\right)$) take into account the adiabatic expansion together with the entropy release from electron-positron annihilation while the rest of the terms take into account the non-adiabatic contribution from the energy exchange with other particles (in the SM with neutrinos). This is simply a result of energy conservation: in the SM, the energy exchange rate of the electron-positron fluid is: $\delta\rho_{e}/\delta t=-\delta\rho_{\nu_{e}}/\delta t-2\,\delta\rho_{\nu_{\mu}}/\delta t$.

Therefore, by assuming negligible chemical potentials and thermal distributions for the neutrinos, electrons and photons, we are able to model the time evolution of neutrino decoupling in terms of just two coupled differential Equations (2.8) and (2.11) subject to the constraint (2.4).

To finally pose the system we are just left with calculating $\delta\rho_{\nu_{e}}/\delta t$ and $\delta\rho_{\nu_{\mu}}/\delta t$ from Equation (2.7). The processes that contribute to them, in the SM, are $e^{+}e^{-}\leftrightarrow\bar{\nu}_{i}\nu_{i}$, $e^{\pm}\nu_{i}\leftrightarrow e^{\pm}\nu_{i}$, $e^{\pm}\bar{\nu}_{i}\leftrightarrow e^{\pm}\bar{\nu}_{i}$, $\nu_{i}\nu_{j}\leftrightarrow\nu_{i}\nu_{j}$, $\nu_{i}\bar{\nu}_{j}\leftrightarrow\nu_{i}\bar{\nu}_{j}$, and $\bar{\nu}_{i}\nu_{i}\leftrightarrow\bar{\nu}_{j}{\nu}_{j}$. The exact collision terms have been calculated in the literature [9]. Here we use the matrix elements outlined in Table 2 of [7] and follow the integration method developed in [9] to obtain the relevant rates222The reader is referred to the appendices of [46] and [47] which very recently appeared in the literature and which explain in more detail and with explicit formulae how to perform the different integration steps.. We neglect statistical factors, the electron mass, and assume Maxwell-Boltzmann statistics for all the particles involved in each process to obtain analytic formulae for the rates333In previous versions of this paper we took the rates from Ref. [39]. However, we found that the formulae in Ref. [39] have missing factors of 2 and 8 for annihilation and scattering rates respectively. Such errors have also been reported in footnote 4 of Ref. [48].. Taking into account all the processes (see Appendix A.2), the energy exchange rates read:

[TABLE]

where we have defined:

[TABLE]

and where $G_{F}=1.1664\times 10^{-5}\,\text{GeV}^{-2}$ is the Fermi constant, and $s_{W}$ is sinus of the Weinberg angle, which given the observed values of $M_{Z}$ and $M_{W}$ reads $s_{W}^{2}=0.223$ [3].

2.2.2 Finite temperature QED corrections

We include the leading order QED finite temperature corrections for the electron and photon mass that result in a smaller pressure of the electromagnetic plasma and hence in an increase of $\sim 0.01$ on $N_{\rm eff}$ [49]. We use the finite temperature corrections at leading order in $\alpha$ from [49] and [50], and proceed as [13] to obtain the temperature evolution of the electromagnetic plasma including finite temperature corrections:

[TABLE]

where $P_{\text{int}}$ reads [13]:

[TABLE]

and the leading order QED finite temperature correction to the masses are [49, 50]:

[TABLE]

where $\alpha$ is the fine structure constant, and for $\delta m_{e}^{2}(T)$ we have neglected a momentum dependent contribution that amounts to less than 10% for an average electron [50, 13]. In addition, $\rho_{\rm int}=-P_{\rm int}+T_{\gamma}\frac{dP_{\rm int}}{dT_{\gamma}}$ must be included in the Hubble rate (2.4).

Note that in practice, since $P_{\rm int}$ is defined as a double integral which evaluation could be slow we simply tabulate it and its derivatives.

2.2.3 Estimate of the neutrino decoupling temperature in the Standard Model

The neutrino-electron interactions have a very characteristic energy dependence $\sigma\sim G_{F}^{2}E_{\nu}^{2}$, which means that the most energetic neutrinos will decouple later from the electromagnetic plasma. Namely, the neutrino decoupling is not an instantaneous process. However, from Equation (2.8) we can easily find an estimate for the temperature at which the average neutrino decouples from the electrons and positrons. By expanding $\delta\rho_{\nu}/\delta t$ in a Taylor series around $\Delta T=T_{\gamma}-T_{\nu}$ we find that:

[TABLE]

And therefore, we notice that the neutrinos will decouple from the plasma, i.e. will track the adiabatic expansion of the Universe, when the term accompanying $\Delta T$ becomes $\mathcal{O}(1)$:

[TABLE]

By considering the expressions in (2.12), $\partial\rho_{\nu}/\partial T=4\rho_{\nu}/T$, and that only electrons, positrons, photons and neutrinos are present in the early Universe we obtain:

[TABLE]

By numerically integrating the relevant collision terms we find that the neutrino decoupling temperature, including the Fermi-Dirac suppression of the rate, reads:

[TABLE]

And hence the decoupling temperatures are enhanced by less than $5\%$ with respect to the MB approximation when the statistical factors are considered.

2.2.4 $N_{\rm eff}$ in the Standard Model

By including the scattering and annihilation interactions between neutrinos and electrons – together with several well justified approximations – we have been able to model the process of neutrino decoupling in terms of just two Equations (2.14) and (2.8) subject to the constraint (2.4). Table 1 summarizes the obtained values for $N_{\text{eff}}$ and the final neutrino temperature in the SM under various assumptions. We note that in the case that we take into account all relevant interactions and finite temperature corrections, we obtain $N_{\text{eff}}^{\rm SM}=3.053$ which is remarkably close to the most precise determination in the SM $N_{\text{eff}}^{\rm SM}=3.045$ [5]. We attribute the difference between our result and that of Ref. [5], to be mainly due to the fact that we do not account for spin-statistics and the electron mass in the interaction rates. Although our modeling does not reproduce the very precise value of $N_{\text{eff}}^{\rm SM}=3.045$ obtained in [5], a pertinent comment is that our evaluation takes $\lesssim$ 10 s on an average computer, while the evaluation performed in [5] is computationally very expensive. This fact represents a clear advantage of our approach to explore the impact of BSM scenarios with light particles that can alter the process of neutrino decoupling.

2.3 Numerics

In order to solve the system of differential equations for $dT_{\gamma}/dt$ (2.14) and $dT_{\nu}/dt$ (2.8), in and beyond the SM, we use the Odeint package of Python that is itself based on the Fortran library Odepack. We require the relative and absolute accuracies of the integrator to be $10^{-8}$ and we calculate the pressure and the energy density of massive particles with relative and absolute accuracies of $10^{-8}$. We have explicitly checked that the continuity equation (2.10) is fulfilled at each point of the integration with a relative accuracy $<10^{-5}$ for the SM and all the BSM cases presented in this work. This ensures that i) no error has been committed while coding the Hubble rate and the expressions for the temperature evolution, and ii) the quoted values of $N_{\rm eff}$ have numerical relative accuracies $<10^{-5}$.

In all the results presented in this work, we solve the time evolution equations for $T_{\gamma}$ and $T_{\nu}$ starting from $T_{\gamma}=T_{\nu}=10\,\text{MeV}$ where the SM weak interactions are highly efficient ($\Gamma/H\sim 50$). The starting time for the evolution is simply $t=1/(2H)$, which in the SM corresponds to $t_{0}\sim 8\times 10^{-3}\,\text{s}$, and we evolve the system until $t_{\rm final}=5\times 10^{4}\,\text{s}$. This corresponds to temperatures $T_{\gamma}\sim 5\times 10^{-3}\,\text{MeV}$ at which the electrons and positrons have already annihilated away.

3 Current CMB measurements on $N_{\rm eff}$ and upcoming sensitivity

Both the polarization and temperature CMB spectra are sensitive to the number of effective neutrino species in the early Universe [3, 51, 52]. The most obvious effect of a change in $N_{\rm eff}$ is to modify the time of matter-radiation equality which is of pivotal importance for structure formation [53]. However, such an effect is not unique as can result form a variation in other cosmological parameters ($\Omega_{b}h^{2},\,\Omega_{\rm cdm}h^{2}$). The most genuine effect of an increase in $N_{\rm eff}$ is to augment the Silk damping scale [54], i.e. it amounts to a reduction of the power in the high-$\ell$ spectra [55, 56, 3]. This is how, Stage-IV generation of CMB experiments [28] – that aim to have high resolution at very small angular scales – expect to reach sensitivities of $\sim 0.03$ on $N_{\rm eff}$.

In this section we review the current $N_{\rm eff}$ measurements that we should use below to set constraints on MeV scale thermal dark matter particles, and comment on the expected sensitivity from Stage-IV CMB experiments on $N_{\rm eff}$ in order to forecast the reach of Stage-IV CMB experiments to them.

3.1 Planck 2018

The Planck satellite [1, 2] reports unprecedented precision measurements of $N_{\text{eff}}$. Various combinations of data sets have been considered in the Planck 2018 analysis within the $\Lambda$CDM framework, and for our study we shall use the Planck 2018 $N_{\text{eff}}$ constraints from the combination of TT+TE+EE+lowE+lensing+BAO:

[TABLE]

In addition, since the observed values of $H_{0}$ from CMB [2] and local observations [57, 58, 59] are in a $\sim 3\sigma$ tension – which could potentially be attributed to a BSM $N_{\rm eff}$ [60] – we also consider the 95% CL Planck 2018 constraint from the combination of Planck+BAO+$H_{0}$:

[TABLE]

Given the current disagreement between local and CMB measurements of $H_{0}$, we consider that an scenario will be conservatively ruled out by Planck 2018 at 95% CL if $N_{\rm eff}<2.66$ from TT+TE+EE+lowE or $N_{\rm eff}>3.57$ from TT+TE+EE+lowE+lensing+BAO+$H_{0}$.

Finally, we note that generically a change on $N_{\rm eff}$ is also accompanied by a change on the primordial Helium abundance [21, 22, 18] ($Y_{p}$), and the CMB is sensitive to both [1, 2]. However, the combination of Planck observations together with the measurements of the primordial abundances from [61] results in $N_{\rm eff}=2.99^{+0.43}_{-0.40}$ at 95% CL [2], which is very similar to that from Planck considering only variations of $N_{\rm eff}$ (3.1). Thus, we consider only the Planck constraints in (3.1) but note that they are robust upon modifications of $Y_{p}$ given the current measurements of the primordial element abundances [3, 4, 61].

3.2 CMB Stage-IV sensitivity

CMB Stage-IV experiments aim to measure $N_{\rm eff}$ with a precision of $\sim 0.03$ [28]. We consider the expected 95% CL sensitivity by multiplying $0.03$ by $2$. In addition, for consistency, we shall assume as the fiducial measurement of CMB Stage-IV experiments the one we obtain in the SM: $N_{\rm eff}=3.053$.

4 Neutrino Decoupling BSM: Thermal Dark Matter

A massive stable particle that, in the early Universe, annihilates to light species at a rate of $\left<\sigma v\right>_{\rm WIMP}\simeq 3\times 10^{-26}\,\text{cm}^{3}/\text{s}$ can account for the observed dark matter abundance in the Universe [62, 63], this is the so-called WIMP miracle [45]. Such a particle would decouple from the plasma while non-relativistic at temperatures of $T\sim m/20$. Therefore, WIMPs of $m\lesssim 20\,\text{MeV}$ will generically alter neutrino decoupling and hence impact $N_{\rm eff}$.

The aim of this section is to work out the temperature evolution equations for the electromagnetic plasma (2.14) and the neutrino sector (2.8) in the presence of a WIMP (which will serve as an illustration of how to do so in BSM scenarios in general). From the calculation in the SM (Section 2) the recipe is clear, if a BSM particle is in thermal equilibrium with either the electromagnetic or neutrino sector of the plasma, then the BSM particle will share its temperature. Generically, BSM particles have interactions with both electrons and neutrinos, so one should include the respective energy exchange rate between the two sectors in Equation (2.7) in order to properly track the temperature evolution.

This Section is structured as follows. First, in subsection 4.1 we work out the temperature evolution equations when WIMPs, that only annihilate to electrons or neutrinos, are in thermal equilibrium with either electrons or neutrinos respectively. We will compare our results with previous approaches that assumed an instantaneous neutrino decoupling, and we will use CMB constraints on $N_{\rm eff}$ to set a lower bound on the WIMP mass. In subsection 4.2 we will consider scenarios in which the WIMP annihilates to both electrons and neutrinos, and will compare the resulting $N_{\rm eff}$ with the observations by the Planck satellite and with the sensitivity of upcoming CMB observations.

4.1 Purely electrophilic or neutrinophilic WIMPs

A WIMP generically couples to all SM fermions [35, 36, 37], and therefore, if sufficiently light ($m<m_{\mu}$), the WIMP annihilation will proceed to neutrinos and electrons. However, for simplicity, in this subsection we will consider scenarios in which WIMPs have couplings to either electrons or neutrinos. The scenario in which WIMPs have couplings to both of these SM fermions will be the case of study of the next subsection.

4.1.1 Temperature evolution equations

If a WIMP ($\chi$) is solely coupled to the neutrinos, then by summing Equation (2.7) for the neutrinos and the WIMP, and by using the chain rule, the neutrino temperature evolution equation reads:

[TABLE]

while the $T_{\gamma}$ evolution equation is unaltered and given by Equation (2.14). Similarly, for a WIMP that solely interacts with electrons or photons, by summing Equation (2.7) for the WIMP, the photons and the electrons, the photon temperature evolution reads:

[TABLE]

while the neutrino temperature evolution is the same as in the SM and given by Equation (2.8). Where, of course, it is understood that the energy density of the WIMP contributes to the Hubble expansion (2.4).

4.1.2 Results and discussion

The energy density and pressure of a given WIMP are uniquely determined by its internal number of degrees of freedom $g$, its spin, and its temperature. In Figure 1 we show the resulting $N_{\rm eff}$ – for fermionic and bosonic dark matter particles – as a function of the dark matter mass for various choices of $g$ and spins. Electrophilic relics release entropy into the electromagnetic plasma and hence result in $N_{\rm eff}<N_{\rm eff}^{\rm SM}$ (see Equation (2.2)), while neutrinophilic relics release entropy into the neutrino sector and hence result in $N_{\rm eff}>N_{\rm eff}^{\rm SM}$. In Table 2 we provide a summary of the current constraints and future sensitivities for the masses of WIMPs with various spins, and that couple to either electrons or neutrinos. We note that, irrespectively of the spin of the WIMP, Planck observations rule out electrophilic and neutrinophilic thermal dark matter particles with $m<3.0\,\text{MeV}$ at 95% CL. In addition, we show that thermal relics purely annihilating to electrons or neutrinos with $m\lesssim 15\,\text{MeV}$ will be generically tested by upcoming CMB experiments.

We note that the obtained lower bound on the dark matter mass is independent on the low velocity behaviour of the annihilation process. This is because the quantity that controls the WIMP relic abundance is the thermally averaged annihilation cross section, and this should be $\left<\sigma v\right>\simeq 3\times 10^{-26}\,\text{cm}^{3}/\text{s}$ regardless of the behaviour at low velocities [62, 63, 45]. This makes our bounds particularly relevant for $p$-wave annihilating relics to electrons since CMB observations are not sensitive to them [2, 64]. This will be for instance the case of a Majorana dark matter particle annihilating to electrons via a dark photon exchange, or a Dirac dark matter particle annihilating to electrons via a dark Higgs exchange [65, 66, 67, 68]. Our bounds are of relevance for both $s$-wave and $p$-wave annihilating WIMPs to neutrinos [69, 70, 71, 72, 73] that remain elusive to searches at terrestrial neutrino experiments [74, 75, 76].

4.1.3 Comparison with previous literature

The impact on neutrino decoupling of MeV scale particles coupled to either the electromagnetic plasma or the neutrino sector has been studied in the past [29, 30, 31, 20, 21, 22, 23]. These studies considered the effect of such particles by assuming that neutrinos decouple instantaneously at a temperature of $T_{\nu}^{\rm dec}=2.3\,\text{MeV}$, and then they tracked the temperature of the neutrino and electromagnetic sectors by imposing entropy conservation [29].

Figure 2 displays a direct comparison of the neutrino temperature evolution from our approach (solid lines) and the previously used in the literature (dashed lines). We clearly notice that the traditional approach shows a somewhat artificial behaviour of the neutrino temperature evolution in the region $T=1-5\,\text{MeV}$ as a result of imposing entropy conservation. Our approach, on the other hand, describes the temperature evolution of the neutrino fluid with a smooth evolution that is attributed to the continuous entropy release of the BSM particle into the sector to which it interacts with.

We also appreciate a notable difference on $N_{\rm eff}$, which could be as large as $\sim 0.2$. At the moment, this difference is below current CMB experiments sensitivities, however, for CMB Stage-IV experiments that aim at sensitivity of $0.03$ for $N_{\rm eff}$, the comparison would result in biased results. The difference on $N_{\rm eff}$ from our approach and the traditional one that assumes an instantaneous decoupling is a result of two effects. First, in the Maxwell-Boltzmann approximation, the neutrinos decouple at $T\simeq 1.83\,\text{MeV}$ (2.19) and not at $T=2.3\,\text{MeV}$ (while in reality they do so at $T\simeq 1.91\,\text{MeV}$ (2.20)). The earlier the neutrinos decouple, the larger the effect on $N_{\rm eff}$ is from a massive relic coupled to either the electromagnetic or neutrino sector. And second, our approach takes into account finite temperature QED corrections that amount to $\Delta N_{\rm eff}=0.011$. We have explicitly checked that the bounds presented in Table 2 are $7$-$16\,\%$, $14$-$24\,\%$, \sim$$-5\% less stringent than those that would be obtained using the traditional instantaneous decoupling approximation for purely neutrinophilic WIMPs (where each number corresponds to Planck+BAO, Planck+BAO+$H_{0}$, and Stage-IV respectively). Similarly, for WIMPs coupled to electrons, our bounds are $24$-$50\,\%$, $18$-$36\,\%$, $21$-$38\,\%$ more restrictive than those obtained within the traditional instantaneous decoupling approximation. These modifications can be understood from the difference between the traditionally used neutrino decoupling temperature, $T_{\nu}^{\rm dec}=2.3\,\text{MeV}$, and the actual neutrino decoupling temperature within the MB approximation, $T\simeq 1.83\,\text{MeV}$ (2.19) – which is roughly $30\%$.

Finally, we note that although the differences between our approach and the commonly used in the literature are small (up to 25%) for $\Delta N_{\rm eff}$, the key possibility that our approach offers for BSM scenarios is that one can easily include energy exchange rates between sectors. Namely, we can consider what is the effect of WIMPs that simultaneously interact with electrons and neutrinos, which is generically what happens [35, 36, 37]. This is the case of study of the next subsection.

4.2 Thermal Dark Matter annihilating to electrons and neutrinos

The couplings of a WIMP with mass $m_{\chi}\lesssim 30\,\text{MeV}$ to electrons and neutrinos will vary in a model to model basis [35, 36, 37]. In this section, we consider the case of a WIMP with couplings to these two SM fermions in order to illustrate how can be affected by a generic WIMP. We note that this scenario has not previously been addressed in the literature and serves as a neat example of the possibilities that our simplified approach to the process of neutrino decoupling offers.

4.2.1 Temperature evolution equations

The WIMP ($\chi$) will share the temperature of the sector to which it annihilates the most. So, if the annihilation predominantly proceeds to neutrinos in the early Universe, the temperature evolution equations will read:

[TABLE]

Where these equations differ from (4.1) and (2.14) in that we have allowed for a non-vanishing energy transfer rate between the neutrino and electromagnetic components of the plasma, as a result of the dark matter particle interactions: $\delta\rho_{\chi}/\delta t$.

In a similar spirit, if the annihilation proceeds predominantly to electrons in the early Universe, the temperature evolution equations reads:

[TABLE]

To solve these systems of equations, we are just left with calculating the energy exchange rates as induced by the dark matter interactions. Unfortunately, these energy exchange rates are model dependent. This is because $\delta\rho_{\chi}/{\delta t}$ takes into account both scattering and annihilation interactions, but the condition of obtaining the correct relic abundance only fixes the annihilation rate to be $\left<\sigma v\right>\simeq 3\times 10^{-26}\,\text{cm}^{3}/\text{s}$ [62, 63]. The parametric dependence on the couplings and masses of the scattering and annihilation cross sections are similar. However, although the energy exchanged in a scattering process $\sim T$ can be substantially smaller to that in an annihilation $\sim m+3\,T$, the annihilation rate should be suppressed by a factor $\sim e^{-m/T}$ with respect to the scattering rate because the light species to which the WIMP scatters are relativistic. Here, we consider only the annihilation energy exchange rate but note that for any given model one should also include the scattering rates (yet, this will serve as a clear illustration of the phenomenology of this scenario). We calculate the energy exchange rate in a model independent manner by considering the usual WIMP number density evolution (this is simply Equation (2.6) for a stable particle) [45, 62]:

[TABLE]

where $n^{\rm eq}$ corresponds to the equilibrium number density. We can easily calculate the energy exchange rate by noting that the energy exchanged in the annihilation process will simply be $\sim 2m_{\chi}$ for a non-relativistic relic ($T<m_{\chi}/3$). Thus, the energy exchange rate is:

[TABLE]

Which, for WIMPs that are in thermal equilibrium with neutrinos but that also annihilate to electrons ($\left<\sigma v\right>_{\chi\chi\to e^{+}e^{-}}\leq\left<\sigma v\right>_{\chi\chi\to\bar{\nu}\nu}$), in the MB approximation reads:

[TABLE]

where $K_{2}$ is a modified Bessel function of the second kind. For WIMPs that are in thermal equilibrium with the electromagnetic plasma but that also have interactions with neutrinos, in the MB approximation reads:

[TABLE]

Note that the annihilation cross section for WIMPs that annihilate to electrons and neutrinos is fixed to be $\left<\sigma v\right>_{\rm WIMP}=\left<\sigma v\right>_{e}+\left<\sigma v\right>_{\nu}\simeq 3\times 10^{26}\,\text{cm}^{3}/\text{s}$.

4.2.2 Results and discussion

In Figure 3 we show the temperature evolution of the neutrino sector as a function of the photon temperature in the presence of a WIMP that annihilates to neutrinos and electrons. The various lines correspond to different final state ratios of the annihilation process. We notice that the effect of light WIMPs that have sizeable interactions with both electrons and neutrinos ($\left<\sigma v\right>_{e,\,\nu}/\left<\sigma v\right>_{\rm WIMP}\gtrsim 10^{-5}$) is to maintain the neutrinos coupled to the electromagnetic plasma for smaller temperatures than the expected one in the SM, $T_{\nu}^{\rm dec}\simeq 1.8\,\text{MeV}$ (2.19). To have an understanding of why this is the case, we can actually estimate the neutrino decoupling temperature by following a similar approach to that in the SM (see Section 2.2.3), namely, we can compare the change in the energy exchange rate from the WIMP annihilation to the Hubble expansion:

[TABLE]

By considering Equation (4.7) and approximating $K_{2}\left[\frac{m}{T}\right]\simeq\sqrt{\frac{\pi T}{2m}}\,e^{-m/T}$ we can easily find the temperature at which the energy transfer between electrons and neutrinos as a result of the WIMP annihilations becomes inefficient:

[TABLE]

Where $g_{\star}$ corresponds to the number of relativistic effective degrees of freedom contributing to radiation. We notice that the decoupling temperature dependence with the annihilation cross section is only logarithmic. Plugging in representative numbers, we find:

[TABLE]

Thus, we realize that a light WIMP that mainly annihilates to electrons but that only annihilates to neutrinos once every $10^{3}$ times can still significantly alter the temperature at which neutrinos decouple. For a WIMP with $m_{\chi}=5\,\text{MeV}$, $\left<\sigma v\right>_{e}=(1-10^{-3})\left<\sigma v\right>_{\rm WIMP}$, and $\left<\sigma v\right>_{\nu}=10^{-3}\left<\sigma v\right>_{\rm WIMP}$, then $T_{\nu}^{\rm dec}\sim 0.8\,\text{MeV}$, which is substantially smaller than the temperature at which neutrinos decouple in the SM $T_{\nu}^{\rm dec}\sim 1.8\,\text{MeV}$ (2.19). Thus, we notice that the effect of WIMPs with masses $m\lesssim 20\,\text{MeV}$ is not only to release entropy into the system but they generically act as to delay the time of neutrino decoupling.

In Figure 4 we show the resulting $N_{\rm eff}$, as relevant for CMB observations, for a Majorana and a scalar thermal dark matter particle for various annihilation final state ratios to electrons and neutrinos. In Table 3 we provide a summary of the constraints – that can be inferred from CMB observations – for a generic Majorana dark matter particle that annihilates differently to electrons and neutrinos (which serves as a representative example of any dark matter particle). The reader is referred to Appendix A.5 for tables with the constraints for any thermal dark matter particle, and figures of $N_{\rm eff}$ as a function of the WIMP mass.

We notice that unless the WIMP annihilates to electrons or neutrinos less than once every $10^{5}$ times, the resulting $N_{\rm eff}$ is quite different to that of the case in which the WIMP exclusively couples to either electrons or neutrinos, as can be appreciated from Figure 1. In addition, the CMB bounds are milder for scenarios in which the WIMP couples non-negligibly to both electrons and neutrinos simultaneously. This is because, in such scenarios, neutrino decoupling is delayed as compared with the SM by the dark matter interactions. Which in turn means that the entropy released by the WIMP annihilation is shared between the neutrino and electromagnetic sectors, thus resulting in $N_{\rm eff}\sim N_{\rm eff}^{\rm SM}$. For the case of a Majorana dark matter particle, CMB observations rule out $m_{\chi}<0.6\,\text{MeV}$ at 95% CL regardless of its annihilation final state ratio to electrons or neutrinos. While for a neutral scalar, CMB observations cannot probe $m_{\phi}>m_{e}$ independently of the annihilation final state. Although we expect that Stage-IV CMB experiments will generically improve the bound on the dark matter mass for any generic WIMP, we have found that there is a window of rather light WIMPs $m\in m_{e}-1.1\,\text{MeV}$ with final state ratios to electrons and neutrinos of $\sim 10^{4}:1$ which will remain elusive even to high-sensitive CMB experiments.

The bounds derived for WIMPs that annihilate to both electrons and neutrinos are relevant for various physically well motivated scenarios. For instance, a dark matter particle charged under a gauged $U(1)_{B-L}$ symmetry will annihilate to electrons and neutrinos at a similar rate [77, 78]. Also, some neutrinophilic dark matter candidates do have small ($10^{-5}-10^{-4}$) but non-negligible annihilations to electrons, as in the case of dark matter within a $U(1)_{\mu-\tau}$ gauge symmetry [79] or in scenarios with mixing between active and sterile neutrino states [80]. Finally, on more general grounds, our bounds will be applicable to dark matter particles that annihilate through a $Z^{\prime}$ with axial couplings to electrons or neutrinos [81]. An axial coupling to electrons implies the existence of an axial coupling to neutrinos, and the other way around, as granted by the SM particle content and $SU(2)_{L}$ gauge invariance.

5 BBN constraints on thermal Dark Matter

Light particles in thermal equilibrium with the SM plasma can impact the synthesis of the primordial elements during BBN [29, 30, 20, 21, 22]. The presence of such particles will affect BBN mainly because [29, 16, 17, 18] i) the proton-to-neutron conversion rates are modified by a change in the neutrino-photon temperature evolution, and ii) the early Universe expansion is modified, in particular the time at which Deuterium starts to form will be modified (while the temperature is fixed to be $T_{\gamma}\simeq 0.85\times 10^{9}\,\text{K}$).

In order to draw precise statements about the effect of MeV scale thermal dark matter particles on BBN, one should implement the temperature evolution obtained in Section 4 in a modern BBN code like PArthENoPE [82, 83], AlterBBN [84, 85] or PRIMAT [86]. However, that is beyond the scope of this work and will be left for future studies. Nonetheless, since the Helium abundance ($Y_{p}$) depends mainly upon the proton-to-neutron conversion rates and the time at which Deuterium forms [16, 17, 18], we can actually provide somewhat reliable estimates for $Y_{p}$ when generic WIMPs are present in the early Universe (see Appendix A.4 for details about the calculation).

In Figure 5 we show the resulting Helium abundance as a function of the mass of a Majorana and a scalar WIMP depending on the annihilation final state ratio to electrons and neutrinos. In Table 4 we show the resulting bounds on various dark matter particles by comparing the obtained Helium abundance to the measured value of $Y_{p}=0.245\pm 0.06$ at 95%CL as recommended by the PDG [3]. We notice that current measurements of the Helium abundance can be sensitive to WIMPs with masses $m<0.4-4\,\text{MeV}$. Although we find that Helium BBN bounds are generically less constraining than those obtained from current CMB observations of $N_{\rm eff}$, the Helium BBN bounds can be slightly more constraining for WIMPs that annihilate to $e:\nu$ in ratios $\sim 10^{4}:1$, thus providing complementary bounds to those inferred from CMB observations.

6 Summary, conclusions and outlook

In this study we have developed a simplified treatment of neutrino decoupling in the early Universe. Thanks to several well justified approximations we have been able to model the neutrino decoupling process in terms of two simple coupled differential equations for $T_{\gamma}$ (2.14) and $T_{\nu}$ (2.8). The virtues of our modeling are that: it incorporates the SM energy exchange rates between the electromagnetic and neutrino sectors of the plasma, it allows for a fast and precise evaluation of $N_{\rm eff}$, and it makes straightforward the inclusion of BSM species. This differs from previous studies in the literature that either assume that neutrinos decouple instantaneously, or that include all the relevant physics (exact collision terms, QED corrections, and neutrino oscillations), but which are computationally very expensive due to the complexity of the resulting kinetic equations.

In this study, we have illustrated the usefulness of our approach to neutrino decoupling in the early Universe to make precise statements about the value of $N_{\rm eff}$ beyond the Standard Model. In particular, we have studied the effect of thermal dark matter (WIMPs) on the process of neutrino decoupling. We have revisited the case of electrophilic and neutrinophilic thermal relics, and for the first time explored scenarios in which the WIMP annihilates simultaneously to electrons and neutrinos. For the latter case, we have shown that WIMPs – in addition to releasing entropy into the system – act as to generically delay the time of neutrino decoupling. By confronting the predicted values of $N_{\rm eff}$ from our treatment, with the observed values as reported by the Planck Collaboration in their 2018 analysis, we have set lower bounds on the masses of various types of WIMPs. Our analysis shows that Planck rules out both $s$-wave and $p$-wave electrophilic and neutrinophilic thermal dark matter candidates with $m<3.0\,\text{MeV}$ at 95% CL. We have pointed out that WIMPs with non-negligible interactions to both electrons and neutrinos are considerably more elusive to CMB observations than purely electrophilic and neutrinophilic dark matter particles. In addition, assisted by the precision of our approach to neutrino decoupling, we have shown that upcoming CMB observations will be sensitive to WIMPs of $m\lesssim 15\,\text{MeV}$.

Although we have calculated the temperature evolution of the neutrino and electromagnetic sectors of the plasma for $T<20\,\text{MeV}$, we have not addressed in detail the impact of BSM sectors in the synthesis of the primordial element abundances during BBN. We have estimated the impact of light thermal dark matter on $Y_{p}$, and found that that particles with $m<0.4-4\,\text{MeV}$ are expected to be constrained by current Helium abundance measurements. However, the actual implementation of the temperature evolution in a modern BBN code is beyond the scope of this work, and is left for future studies.

Our approach to neutrino decoupling relies on very useful and well justified simplifying assumptions, which make the description of the neutrino decoupling process straightforward. Some improvements that would however lead to a more involved system of equations or non-analytical collision terms are: the inclusion of neutrino oscillations (which could be done with the methods developed in [43]), including non-negligible chemical potentials for the neutrinos (which could be relevant for certain BSM scenarios), and to include spin statistics in the neutrino collision terms (which would require 2-d integrations at each time step).

Finally, we note that the possibilities of our simplified approach to neutrino decoupling are countless given the variety of extensions of the Standard Model with massless or light particles (for instance, see [87]).

Acknowledgments

The author thanks Dan Hooper, Gordan Krnjaic and Mathias Pierre for collaboration on related topics. The author is grateful to Sam Witte for a careful read of the draft, and also thanks Chris McCabe, Nuria Rius, Sergio Pastor and Sam Witte for their very helpful comments on earlier versions of the manuscript. The author is supported by the European Research Council under the European Union’s Horizon 2020 program (ERC Grant Agreement No 648680 DARKHORIZONS).

This paper is dedicated to my mother (in memoriam), Antonia Abenza Ortiz, and to my partner, Andrea González Montoro, for their endless love, support and encouragement.

Appendix A Appendices

A.1 Neutrino oscillations in the primordial Universe

Here we provide an estimate for the temperature at which neutrinos begin to oscillate in the early Universe. The reader is referred to [43, 44, 7] for detailed treatments.

We can estimate the temperature at which neutrinos start to oscillate by comparing the characteristic time for oscillation, scattering, and expansion of the Universe. It is clear that neutrino oscillations will be active when $t_{\nu}^{\text{os}}<t_{\rm exp}<t_{\nu}^{\text{scat}}$, i.e. when the oscillations take place within the lifetime of the Universe and are not decohered by scattering interactions. The oscillation timescale can easily be estimated from the 2 neutrino oscillation formula [3]:

[TABLE]

where $E$ is the neutrino energy, $\Delta m^{2}$ is the mass squared difference, and $\theta_{\beta\alpha}$ corresponds to the mixing angle. Therefore, at the time at which the argument of the second sinus becomes $\mathcal{O}(1)$, the oscillations will be active unless the scattering interactions damp them. Taking $E\to\left<E\right>=3\,T$ we can find all the relevant timescales:

[TABLE]

Figure 6 shows the evolution of the different timescales as relevant for neutrino decoupling. Hence, we can clearly appreciate that oscillations become active for $T\sim 3-5\,\text{MeV}$ [43, 44].

We also note that, even if BSM interactions are present, in order no to modify significantly $N_{\rm eff}$, the BSM-neutrino rates will become smaller than the Hubble rate at $T\sim\text{MeV}$ and therefore neutrino oscillations will be active regardless of possible BSM interactions given the strong energy dependence of the oscillation rate.

A.2 Neutrino collision terms in the SM in the Maxwell-Boltzmann approximation

In this appendix we outline the results of the calculation of the collision terms for annihilations and scatterings assuming Maxwell-Boltzmann statistics for all the particles involved in the process, neglecting the statistical factors and assuming that all the species are massless. These results differ from those presented by the author in previous versions of the paper. The results outlined in previous versions of this work were based on the formulae presented in [39] which we have found to have missing factors of 2 and 8 for the respective annihilation and scatterings interactions. Such errors have also been reported in footnote 4 of [48].

We have been able to fully calculate the collision terms following the phase space reduction method developed in [9]. We refer the reader to the appendices of [46] and [47] which very recently appeared in the literature and which explain in more detail and with explicit formulae how to perform the different integration steps.

By assuming Maxwell-Boltzmann statistics for all the particles involved in the process, neglecting the electron mass, neglecting statistical factors, and by considering the matrix elements reported in Table 2 of [7] we arrive to the following annihilation collision term:

[TABLE]

where $G_{F}$ is the Fermi constant, $g_{L}^{\nu_{e}}=\frac{1}{2}+s_{W}^{2}$, $g_{R}^{\nu_{e}}=s_{W}^{2}$, $g_{L}^{\nu_{\mu,\,\tau}}=-\frac{1}{2}+s_{W}^{2}$, $g_{L}^{\nu_{\mu,\,\tau}}=s_{W}^{2}$, and $s_{W}$ is the sinus of the Weinberg angle. Expression (A.3) precisely matches the one obtained from Equation B.2 in [11]. The scattering rate between neutrinos and electrons can similarly be obtained to give:

[TABLE]

If we were to consider only the forward part of the annihilation and scattering processes ($\bar{\nu}\nu\to{e^{+}e^{-}},\,\nu e^{\pm}\to\nu e^{\pm}$) between electrons and neutrinos we will recover Equation 63 in [7]. The momentum averaged rates are:

[TABLE]

By following a similar procedure, we compute the neutrino-neutrino collision terms, and we perform the integration over momenta for all the relevant processes which yield:

[TABLE]

With these results, we are able to write the total number density and energy density transfer rates in the SM within the Maxwell-Boltzmann approximation:

[TABLE]

where we have defined:

[TABLE]

Similarly, the number density transfer rates read:

[TABLE]

where $\delta n_{i}/\delta t\equiv\sum_{j}\left<\sigma v\right>_{ij}(n_{i}^{2}-n_{j}^{2})$.

A.3 Thermodynamics formulae

Here we list the formulae used to compute the energy density, pressure density, and the variation of the energy density with respect to the temperature for the relevant species:

[TABLE]

Note that in practice, we truncate the integral for massive particles at $E=100\,T$, and set $\rho_{\chi}=p_{\chi}=\partial\rho_{\chi}/\partial T=0$ for $T<m_{\chi}/30$.

A.4 Helium abundance

Here we outline the equations used for our calculation of the primordial Helium abundance. The Boltzmann equation that governs the number of neutrons in the early Universe is [18]:

[TABLE]

where $X_{n}\equiv\frac{n_{n}}{n_{n}+n_{p}}$. We can easily estimate the Helium abundance by noting that nearly all the present neutrons at the time at which Deuterium is no longer dissociated by photons ($T_{D}=0.85\times 10^{9}\,\text{K}=0.073\,\text{MeV}$) will form Helium [18]. This means that:

[TABLE]

In the early Universe, there are three processes that can efficiently convert protons into neutrons and the other way around. They are $n+\nu_{e}\leftrightarrow p+e^{-}$, $n+e^{+}\leftrightarrow p+\bar{\nu}_{e}$ and $n\leftrightarrow p+e^{-}+\bar{\nu}_{e}$. The rate for neutron to proton conversion reads [8]:

[TABLE]

where $z_{\gamma}=m_{e}/T_{\gamma}$, $z_{\nu}=m_{e}/T_{\nu}$, $q=(m_{n}-m_{p})/m_{e}$, and $m_{n}-m_{p}=1.2933\,\text{MeV}$ [3] and $K\simeq(1.636\,\tau_{n})^{-1}$ and we take $\tau_{n}=880.2\,\text{s}$ [3]. The proton-to-neutron rates are obtained by $\Gamma_{pn}=\Gamma_{np}(-q)$.

A.5 Supplemental material for generic thermal Dark Matter candidates

Here we show the corresponding $N_{\rm eff}$ as a function of the mass of a generic Majorana, Dirac, neutral scalar, charged scalar, and neutral vector boson dark matter particle. We also provide tables with the corresponding bounds given the observed values of $N_{\rm eff}$ from current CMB observations and the sensitivity from upcoming CMB experiments. Constraints for other branching fractions can be inferred by interpolation or by use of the online code.