Analytic expressions for the background evolution of massive neutrinos and dark matter particles

TL;DR

This paper derives exact analytic formulas for the background evolution of massive neutrinos and dark matter particles, and validates them against numerical models, enhancing precision in cosmological computations.

Contribution

It introduces precise analytic expressions for neutrino and dark matter background properties, improving upon existing numerical methods in cosmological modeling.

Findings

Analytic expressions agree within 0.0001% with CLASS code.

Modifications do not affect computational performance.

Provides limits at relativistic and late times for neutrino equations of state.

Abstract

We provide exact analytic expressions for the density, pressure, average number density and pseudo-pressure for massive neutrinos and generic dark matter particles, both fermions and bosons. We then focus on massive neutrinos and we compare our analytic expressions with the numerical implementation in the CLASS Boltzmann code. We find that our modifications including the exact analytic expressions are in agreement to better than $10^{-4}\%$ with the default CLASS implementation in the estimation of the CMB power spectrum; our modifications do not have an impact on the performance of the code. We also provide several specific limits of our expressions at the relativistic regime, but also at late times for the neutrino equation of state.

Equations

$$n(a)$$

$$\rho (a)$$

$$P(a)$$

$$psP(a)$$

$$f_0(p)=\frac{g_s}{e^{\frac{E(p)}{T}}\pm 1},$$

$$f_0(p)$$

$$f_0(p)=\frac{g_s}{e^{\frac{p}{T_{\nu }(a)}}+1}.$$

$$n(a)$$

$$\rho (a)$$

$$P(a)$$

$$psP(a)$$

$$I_{n,k}\equiv {\int }_0^{\infty }dQ{Q}^n\frac{{\left(Q^2+a^2{M}^2\right)}^{k/2}}{e^{\frac{Q}{\tilde{T}}}+1},$$

$$I_{n,k}$$

$$\int_{0}^{\infty}z^{2\nu-1}\left(u^{2}+z^{2}\right)^{\alpha-1}e^{-\mu z}dz=\frac{u^{2\nu+2\alpha-2}}{2\sqrt{\pi}\Gamma(1-\alpha)}G_{1,3}^{3,1}\left(\begin{matrix}1-\nu\\ 1-\alpha-\nu,0,\frac{1}{2}\end{matrix}\bigg{|}\frac{\mu^{2}u^{2}}{4}\right),$$

$$\displaystyle I_{n,k}=\sum_{i=1}^{\infty}(-1)^{i+1}\frac{(aM)^{1+k+n}}{2\sqrt{\pi}\Gamma(-k/2)}G_{1,3}^{3,1}\left(\begin{matrix}\frac{1-n}{2}\\ -\frac{1}{2}(1+k+n),0,\frac{1}{2}\end{matrix}\bigg{|}\frac{x_{i}^{2}}{4}\right),$$

$$n(a)=\frac{6\pi g_s\zeta (3)}{a^3}{\tilde{T}}^3.$$

$$\rho (a)$$

$$\rho (a)=\frac{7{\pi }^5{g}_s}{30a^4}{\tilde{T}}^4.$$

$$\frac{d\rho (a)}{dM}$$

$$P(a)$$

$$P(a)=\frac{7{\pi }^5{g}_s}{90a^4}{\tilde{T}}^4=\frac{\rho (a)}{3},$$

$$psP(a)$$

$$psP(a)=\frac{7{\pi }^5{g}_s}{90a^4}{\tilde{T}}^4\equiv P(a)=\frac{\rho (a)}{3}.$$

$$\rho (a)$$

$$P(a)$$

$$n(a)=\frac{8\pi g_s}{a^3}{\tilde{T}}^3.$$

$$P(a)=\frac{4{\pi }^5}{45a^4}{g}_s{\tilde{T}}^4=\frac{\rho (a)}{3}.$$

$$k_{FT}(t)$$

$${\lambda }_{FT}(t)$$

$$⟨p⟩$$

$$=$$

$$\simeq$$

$${\lambda }_{FT}(t)$$

$$f_X(p)=f_{eq}\left(p\frac{a(\eta )}{a({\eta }_D)},T_D\right)=\frac{g_s}{e^{\frac{\sqrt{p^2{a}^2/a_D^2+M^2}}{T_D}+1}},$$

$$n(a)$$

Full text

1 Introduction

Over the past decades Dark Matter (DM) has become a fundamental ingredient in the standard model of cosmology [1]. Although we know relatively little about its nature, it is clear that taking into consideration DM when modelling the Universe makes it possible to explain a wide variety of astrophysical observations [2, 3, 4, 5]. Nowadays, it is commonly believed that DM comprises beyond Standard Model particles which move slowly with respect to the speed of light and whose interaction with other particles does not go beyond gravity: the so-called Cold Dark Matter (CDM). However, in the Standard Model of particle physics there exist candidates with similar properties which can account for a fraction of the DM in the Universe. Neutrinos weakly interact with other particles and their speed of propagation is different at late- and early-times in the cosmic evolution: in the beginning their speed of propagation is very close to the speed of light and recently they became non-relativistic. This sort of DM is usually dubbed non-Cold Dark Matter (nCDM).

Even though there is compelling evidence for flavour neutrino oscillations which implies that neutrinos are massive particles [6, 7, 8, 9, 10, 11], current constraints do not fully determine their absolute mass scale [12, 13, 14]. Nevertheless, this situation is expected to change with upcoming galaxy surveys which will be able to measure the galaxy distribution on scales comparable to the horizon [15]. Since massive neutrinos suppress power on small scales [16], accurate measurements of the matter power spectrum will lead to a detection of their absolute mass scale thus reducing our ignorance of the abundance of DM in the Universe [17, 18, 15]. Furthermore, measurement of neutrino masses could give hints about new fundamental theories having the Standard model of particle physics as a low-energy limit.

Due to their weakly interacting nature, neutrinos obey a collisionless Boltzmann equation. However, since neutrinos are massive particles the evolution of their phase-space distribution function is not trivial [19]. In order to find the unperturbed density and pressure for neutrinos current implementations in Boltzmann solvers, such as CAMB222https://camb.info/ [20] and CLASS333http://class-code.net/ [21, 22], employ numerical methods. Shortcomings of the numerical approach include non-trivial weighting scheme to carry out the numerical integration, possible limited precision, increase of computing time, but more importantly hindering the understanding of the underlying physics. In this paper we show that a careful analytical treatment of the integrals makes it possible to overcome these difficulties. We provide explicit analytical solutions for the neutrino’s unperturbed density, pressure, number density, and pseudo-pressure. Our expressions agree with previous phenomenological attempts of analytical approximations444See, for instance, Ref. [23]. and also with the fully numerical implementation of the code CLASS. We have implemented our solutions in CLASS and verified that the fully numerical approach (current implementation in CLASS) and the fully analytical approach are in very good agreement. These changes in the code leave precision and computing time unchanged.

This paper is organized as follows. Firstly, in Section 2 we derive our main results, namely, analytical expressions for the background evolution of massive fermions and bosons that are either relativistic or non-relativistic at decoupling. Secondly, in order to compare with previous phenomenological attempts of analytical approximations we provide asymptotic expansions at late times for the quantities governing the neutrino background evolution in Section 3. Thirdly, in Section 4 we implement our analytical expressions for massive neutrinos in the code CLASS and compare with the current numerical implementation in the code. Finally, we conclude in Section 5.

2 Theoretical framework

In this section we will derive simple analytic expressions for several key quantities that are relevant for the background evolution of massive particles, such as the average number density $n(a)$, the density $\rho(a)$ and pressure $P(a)$ of a particle given its phase-space distribution. For the implementation in Boltzmann codes, it is also useful to calculate the derivative of the so-called pseudo-pressure, which we denote by $psP(a)$. All of these quantities are given by the following expressions:555Note that here and in what follows, we will use natural units in which $c=\hbar=k_{B}=1$.

[TABLE]

where $p$ is the physical momentum of the particles, $a$ is the scale factor, $E$ is the energy, while the distribution $f_{0}(p)$ is given by

[TABLE]

where $g_{s}$ is the degeneracy of the species, $T$ is the temperature of the particles and the $\pm$ corresponds to fermions/bosons respectively.666We ignore the chemical potential $\mu$ in our analysis as in all realistic particles, it is much smaller than the temperature. Moreover, current bounds on the common value of the neutrino degeneracy parameter indicate that a neutrino chemical potential can be safely neglected [24, 25, 26, 27, 28, 29, 30, 31].

As the Universe expands and cools down, the temperature will reach the decoupling temperature $T_{D}$ and all interactions will freeze out, so that the phase space distribution of Eq. (2.5) of a particle with mass $M$ will remain frozen [32, 11]

[TABLE]

Here $\eta$ is the conformal time, $f_{eq}$ is the distribution at thermal equilibrium, the subscript $D$ denotes decoupling, $a_{D}\equiv a(\eta_{D})$, and $a\equiv a(\eta)$. Thus, we will consider two separate cases for the distribution $f_{0}(p)$ at the decoupling temperature $T_{D}$:

The particles are relativistic, with energy $E(p)\sim p$; 2. 2.

The particles are non-relativistic, with energy $E(p)=\sqrt{p^{2}+M^{2}}$.

Note that this will only affect the distribution $f_{0}(p)$ and not the energy in the integrand, which can be allowed to be time-dependent. Following, we will present the analytical expressions for the thermodynamic quantities for both fermions and bosons that are relativistic and non-relativistic in Section 2.1 and Section 2.2 respectively.

2.1 Relativistic fermions and bosons at decoupling

In this section we are mainly interested in massive neutrinos and we will specifically focus on them, but our results are readily applicable to other massive relics that are relativistic at decoupling. Neutrino decoupling happened at $T_{D}\sim 1\,\textrm{MeV}$ or $z\sim 10^{10}$, so at that point neutrinos are still relativistic and their distribution can be written as

[TABLE]

Taking into account the expansion of the Universe, we see that the physical momentum $p$ will be redshifted and can be written in terms of the comoving momentum $Q$ as $p=Q/a$, where $a=\frac{1}{1+z}$ is the scale factor and $z$ is the redshift. After neutrino decoupling the temperature scales as $T_{\nu}(a)=T_{\nu,0}/a$ and $T_{\nu,0}\simeq\left(\frac{4}{11}\right)^{1/3}T_{\textrm{cmb}}$ is the neutrino temperature today with a value $T_{\nu,0}\sim 1.68\cdot 10^{-4}\,\textrm{eV}$. Therefore, the combination $Q/T_{\nu,0}$ will be constant and does not depend on the redshift, thus the distribution is frozen.

Defining $T_{\nu,0}$ as $T_{\nu,0}=\widetilde{T}$, the previous equations for the evolution variables can be written as

[TABLE]

In the previous equations all the integrals are of the form

[TABLE]

where $(n,k)$ are integers. In order to calculate $I_{n,k}$ analytically, we multiply the numerator and denominator with the term $e^{\frac{-Q}{\widetilde{T}}}$ and then we use the expansion $\frac{x}{1+x}=\sum_{\beta=1}^{\infty}(-1)^{\beta+1}x^{\beta}$ for $x\leq 1$, which in our case is possible as $e^{\frac{-Q}{\widetilde{T}}}\leq 1$ for all $Q\in[0,\infty)$, thus our series will always converge. Then, we have

[TABLE]

To solve the previous integral we use Eq. (3.389.2) from Ref. [33]

[TABLE]

where Re $\mu>0$, Re $\nu>0$, $|\text{arg}\,u<\frac{\pi}{2}|$ and G_{p,q}^{m,n}\left(\begin{matrix}a_{1},\ldots,a_{n}&a_{n+1},\ldots,a_{p}\\ b_{1},\ldots,b_{m}&b_{m+1},\ldots,b_{q}\end{matrix}\bigg{|}z\right) is the Meijer-G function. With this expression we find that

[TABLE]

where for convenience we have set $x_{i}=\frac{iaM}{\widetilde{T}}$. Next we will provide the explicit expressions for each of the key background quantities mentioned earlier.

2.1.1 Average number density

The average number density $n(a)$ corresponds to the parameters $(n,k)=(2,0)$, so combining Eqs. (2.8) and (2.15) gives the well known result:

[TABLE]

2.1.2 The density

The density corresponds to the parameters $(n,k)=(2,1)$ and the final result can be found to be

[TABLE]

where $K_{\nu}(x)=H_{\nu}(x)-Y_{\nu}(x)$ is the Struve K function, $H_{\nu}(x)$ is the Struve H function and $Y_{\nu}(x)$ the usual Bessel Y function of the second kind [34]. In the relativistic limit, where $M=0$, we find

[TABLE]

The derivative $\frac{d\rho}{dM}(a)$, which is also useful in calculations in Boltzmann solvers, corresponds to $(n,k)=(2,-1)$ and is given by

[TABLE]

2.1.3 Pressure

The pressure corresponds to the set of parameters $(n,k)=(4,-1)$ and as a result we have

[TABLE]

In the relativistic limit, where $M=0$, we find

[TABLE]

as expected for relativistic particles.

2.1.4 Pseudo pressure

The pseudo-pressure corresponds to the set of parameters $(n,k)=(6,-3)$

[TABLE]

In the relativistic limit, where $M=0$, we find

[TABLE]

2.1.5 Results for massive bosons

In the case of bosons, the analytical expressions for the background are very similar to the ones found for fermions. The only difference is the factor $(-1)^{i}$ that appears in the sum which has to be replaced by $(-1)$. Hence, the density and pressure are

[TABLE]

and for the average number density we obtain the well known result

[TABLE]

In the relativistic limit where $M=0$, we find

[TABLE]

2.1.6 Free-streaming length

Similarly, we can also calculate the free-streaming length, i.e., the typical distance particles travel between interactions, which is defined via [35, 36]:

[TABLE]

where $v_{th}\equiv\frac{\langle p\rangle}{m}$ is the thermal velocity and $\langle p\rangle$ the average particle momentum. After the particles become non-relativistic we can calculate the average momentum, by using the results in the previous sections for non-relativistic massive particles, as follows

[TABLE]

Finally, we have that the free-streaming length is

[TABLE]

which is in good agreement with the result of Ref. [35, 36].

2.2 Non-relativistic fermions and bosons at decoupling

When we have massive fermions that are non-relativistic at decoupling ($M_{X}\gg T_{D}$) their distribution function after the freeze out or decoupling can be written as [32]

[TABLE]

where the subscript $D$ denotes decoupling and $a_{D}\equiv a(\eta_{D})$, $a\equiv a(\eta)$. Defining $\tilde{T}=T_{0}$, the comoving momentum $Q$ as $p=Q/a$ and the temperature parameter $T_{0}\equiv Ta\equiv T_{D}a_{D}$, following a similar approach as in Section 2.1 we can compute the average number density, the energy density and pressure as

[TABLE]

2.2.1 Average number density

To solve Eq.(2.35) we first perform a change of variables to hyperbolic functions. Then, using Eq. (3.547.2) from Ref. [33] we find

[TABLE]

where $K_{n}(z)$ is the modified Bessel function of the second kind and $y_{i}=\frac{ia_{D}M}{\tilde{T}}$.

2.2.2 The density

Using Eq. (7.6.1) from Ref. [37] we find after some algebraic manipulations that the density can be written as

[TABLE]

where $K_{n}(z)$ is the modified Bessel function of the second kind and $y_{i}=\frac{ia_{D}M}{\tilde{T}}$.

2.2.3 The pressure

Following the same procedure as with the density, we find that the pressure can be written as

[TABLE]

where again $K_{n}(z)$ is the modified Bessel function of the second kind and $y_{i}=\frac{ia_{D}M}{\tilde{T}}$.

2.2.4 Numerical results

We compared the analytical expressions for the average number density, density and pressure Eqs. (2.38)-(2.40), with the numerical integration of Eqs. (2.35)-(2.37) for some realistic values of WIMP particles like neutralinos [38, 39] $z_{D}\sim 10^{13},m_{X}\sim 25\,\text{GeV},T_{D}\sim 1\,\text{GeV}$ and found that the agreement is better than $10^{-18}\%$ for only 4 iterations.

2.2.5 Results for massive bosons

In the case of bosons, the analytical expressions for the background are very similar to the ones found for fermions, the only difference being the factor $(-1)^{i+1}$ that appears in the sum which has to be replaced by $1$. Hence

[TABLE]

In the non relativistic limit $M_{X}\gg T_{D}$, we can find semi-analytical expressions valid for fermions and bosons for the average number density (2.35), density (2.36) and pressure (2.37) in the following way

[TABLE]

where we have to consider the series expansion $\sqrt{Q^{2}+a^{2}_{D}M^{2}}\sim a_{D}M+\frac{Q^{2}}{2a_{D}M}+O(Q)^{4}$ in the exponential and also assume that $e^{\frac{\sqrt{Q^{2}+a^{2}_{D}M^{2}}}{\tilde{T}}}\gg 1$. Again $K_{n}(\beta)$ is the modified Bessel function of the second kind and $\beta=\frac{a^{2}M}{4a^{2}_{D}T_{D}}$.

3 Asymptotic expansions at late times

The Struve K function $K_{\nu}(z)$ is a particular solution of the inhomogeneous Bessel differential equation

[TABLE]

and it admits the following asymptotic expansion for large values of the argument $z$ with fixed $\nu$ [34]:

[TABLE]

which can be used to obtain asymptotic expansions for the quantities in the previous section. Specifically, we find

[TABLE]

Keeping the zero-order terms for the density and the pressure gives an approximation at late times for the equation of state $w\equiv\frac{P}{\rho}$ as

[TABLE]

which is accurate to a few percent at late times $z<10$. This expression is also in excellent agreement with the ansatz of Ref. [23] that at late times the equation of state scales as $w(a)\sim 1/a^{2}$. Moreover, Eq. (3.7) also provides us with the exact numerical coefficient $\frac{5\zeta(5)}{\zeta(3)}\frac{\widetilde{T}^{2}}{M^{2}}$.

4 Numerical results and implementation in CLASS

Here we present numerical comparisons between our analytic results for massive neutrinos, see Sec. 2.1, and numerical calculations of the quantities based on double precision calculations from CLASS, arbitrary precision calculations in Mathematica and the CEPHES library777https://www.netlib.org/cephes/index.html that we used to implement the Struve K functions in C.

A lower limit on the neutrino mass of approximately $m_{\nu}\sim 0.06\textrm{eV}$ is settled by the existence of three-flavour oscillations (Refs. [8, 9, 10]), independently of their nature (Dirac or Majorana) and this is the value that we will use in what follows.

First, we compare the implementation of the Struve K functions in CEPHES with Mathematica’s arbitrary precision calculations. The results of this comparison are shown in Fig. 1, where we present the percent difference of the implementation in the CEPHES library vs. the arbitrary precision code of Mathematica for $K_{0}(x)$ (solid black line) and $K_{1}(x)$ (dashed black line) for $x\in[10^{-3},10^{6}]$. We find that in both cases, on average the agreement between the two codes is on the order of $\sim 10^{-12}\%$ for both functions, thus we are confident in our numerical implementation in what follows.

Next, we compare our numerical implementation of the analytical expressions for the neutrino density and pressure given by Eqs. (2.17) and (2.20), with the numerical integration done in CLASS. For this comparison we assumed no relativistic species and only 1 massive neutrino of mass $m_{\nu}=0.06\,\textrm{eV}$, while keeping all other parameters in CLASS in their default values. The results of the comparison are shown in Fig. 2, where we present the percent difference between the default version of CLASS and our analytical expressions for the density (left) and the pressure (right) for 10, 50 and 100 terms (black, green and blue lines) of the analytical expressions given by Eqs. (2.17) and (2.20). We find that keeping 50 terms in the expansion yields an accuracy of $10^{-4}\%$ on average for the density and pressure, without affecting the computational performance.

Then we also compare the results of the CMB power spectrum for our implementation and that of the default version of CLASS. The results of the comparison are shown in Fig. 3, where we present the percent difference in the CMB power spectrum for 10 terms in the expansion (black line), 50 terms (green line) and 100 iterations (blue line). We find that keeping 50 terms in the expansion yields an accuracy of $10^{-4}\%$ on average for the $C_{\ell}^{TT}$ of the CMB spectrum, without having an impact on the performance of the code.

We also test the approximation for the equation of state $w(z)$ of the neutrinos at late times, given by Eq. (3.7). The comparison for one massive neutrino of mass $m_{\nu}=0.06\,\textrm{eV}$ is shown in Fig. 4, where we present the percent difference in the equation of state $w(a)$ between the numerical results (solid black line) and the approximation of Eq. (3.7) (dashed line), for which $w(a)\sim a^{-2}$. As can be seen in the inset plot, at late times ($z<10)$ the agreement is better that $1\%$, thus validating the ansatz of Ref. [23]. Finally, we have also checked and confirmed that using a slightly larger neutrino mass, such as $m_{\nu}=0.15-0.30\,\textrm{eV}$, does not affect the precision of our comparison with CLASS.

5 Conclusions

In this paper we presented simple but exact analytical expressions for the background evolution of the density $\rho(a)$, the pressure $P(a)$ and the average number density $n(a)$ for massive particles, both fermions and bosons. In both cases we considered the case when the particles are either relativistic or non-relativistic at the time of decoupling. We find that for non-relativistic massive particles the expressions are somewhat more cumbersome due to the presence of a double sum but in principle these results could be useful in future studies of dark matter candidates, such as WIMPs or any of the hypothetical superpartners of the leptons (sneutrino, etc.).

We also specifically tested our expressions, given by Eqs. (2.17) and (2.20) for the density and pressure respectively, in the case of massive neutrinos that are still relativistic at decoupling ($z\sim 10^{10}$), assuming one neutrino with mass of $m_{\nu}=0.06\,\textrm{eV}$. We implemented our analytical expressions in the Boltzmann code CLASS and found that by keeping 50 terms in the sum, e.g., in Eqs. (2.17) and (2.20), it is possible to achieve better than $10^{-4}\%$ accuracy with respect to the default implementation in CLASS. Our modifications in the code do not have an impact in the computational performance and avoid the involved quadrature integration scheme at the background level. Our analytical expressions provide validation for the current numerical implementations in public Boltzmann codes. By comparing CMB angular power spectra, we find the agreement between our analytical approach and the current numerical implementation is better than $10^{-4}\%$.

The main advantage of our approach is that our expressions are both exact and analytic, thus they can also provide useful intuition about the behavior of the background quantities for massive particles and how they affect the CMB. Moreover, our analytical expressions allow us to compute quantities such as the entropy density $s=\left(\rho+P\right)/T$ or the conserved number $Y=n/s$. For instance, it is possible to derive the exact behavior of the neutrino equation of state $w(a)$ at late times $(z<10)$ and show it behaves as $w\sim a^{-2}$ to better than $1\%$, in agreement with the ansatz of Ref. [23], thus demonstrating how fast massive neutrinos can become non-relativistic.

Numerical Analysis Files: The numerical codes used by the authors in the analysis of the paper and our modifications to the CLASS code will be released upon publication of the paper on the websites of the EFCLASS code here and here.

Acknowledgments

The authors would like to thank Juan García-Bellido and Julien Lesgourgues for useful discussions. They also acknowledge support from the Research Projects FPA2015-68048-03-3P [MINECO-FEDER], PGC2018-094773-B-C32 and the Centro de Excelencia Severo Ochoa Program SEV-2016-0597. S.N. also acknowledges support from the Ramón y Cajal program through Grant No. RYC-2014-15843.