A Precision Calculation of the Effective Number of Cosmological Neutrinos

The neutrino energy density of the Universe can be conveniently parametrized in terms of the so-called effective number of neutrinos, N_nu^eff. This parameter enters in several cosmological observables. In particular it is an important input in those numerical codes, like CMBFAST, which are used to study the Cosmic Microwave Background anisotropy spectrum. By studying the neutrino decoupling with Boltzmann equations, one can show that this quantity differs from the number of massless neutrino species for an additional contribution due to a partial heating of neutrinos during the electron-positron annihilations, leading to non thermal features in their final distributions. In this paper we review the different results obtained in the literature and perform a new analysis which takes into account, in a fully consistent way, the QED corrections at finite temperature to the photon and e^+- plasma equation of state. The value found for three massless active neutrinos is N_nu^eff=3.0395, in perfect agreement with the recommended value used in CMBFAST, N_nu^eff=3.04. We also discuss the case of additional relativistic relics and massive active neutrinos.


Introduction
At temperatures below the muon mass and above ∼ 10 MeV, the Universe is filled by a plasma of photons, electrons, positrons, and neutrinos, kept in thermodynamical equilibrium by the electroweak interactions. As the temperature drops below this value, the rate of weak interactions starts to be comparable with the Universe expansion rate, and neutrinos decouple from the electromagnetic γ, e ± plasma. For most practical purposes, it is accurate enough to consider the freeze-out of neutrinos as fully achieved at temperatures of about 2 − 3 MeV. In this limit neutrinos do not share any entropy release from e ± annihilations, once the temperature drops further, below the electron mass. Assuming that all the entropy produced by the annihilations is transferred to photons, their temperature T is increased with respect to the neutrino temperature by the well known factor T /T ν = (11/4) 1/3 (see, for example, [1]). Actually more accurate calculations show that neutrinos are still slightly interacting with the electromagnetic plasma and thus receive a small portion of the entropy from e ± annihilations [2]. Neutrinos with higher momenta are more heated, since, in the relevant range of energies, weak interactions get stronger with rising energy. This produces a momentum dependent distortion in the neutrino spectra from the equilibrium Fermi-Dirac behaviour.
A further, though smaller, effect on T /T ν is induced by finite temperature QED corrections to the electromagnetic plasma. In fact electromagnetic interactions modify the e ± and γ dispersion relations, and thus the energy density and pressure of the plasma. More precisely, the energy density is lowered so that the e ± annihilation phase releases less entropy with respect to the non-interacting particle limit calculation. Since most of this energy ends up into photons, this decrease results in a smaller T /T ν ratio [3].
While any direct observation of the actual distortions in the neutrino distributions is out of question, nevertheless we can hope that the increase of the total energy of the relic neutrinos may have a sizeable effect on the expansion rate of the Universe. The energy density stored in relativistic species, ρ R , is customarily given in terms of the so-called effective number of neutrino species [4,5], N eff ν , through the relation where ρ γ is the energy density of photons, whose value today is known from the measurement of the Cosmic Microwave Background (CMB) temperature. Eq. (1) can be also written as where ρ 0 ν denotes the energy density of a single specie of massless neutrino with an equilibrium Fermi-Dirac distribution, and ρ 0 γ is the photon energy density in the instantaneous neutrino decoupling approximation. The normalization of N eff ν is such that it gives N eff ν = 3 in the standard case of three families of massless neutrinos, in the limit of an instantaneous decoupling. In principle, N eff ν can receive contribution from other relativistic relics with energy density ρ X as well. In the following we will mostly restrict our analysis to the standard case, but we will further consider this more general framework in Section 4.
As we will discuss below, when considered separately, the non instantaneous neutrino decoupling gives ∆N eff ν ≡ N eff ν − 3 ≃ 0.034 [6][7][8], while QED effects contribute for about ∆N eff ν ≃ 0.011 [3]. Could the two effects be added linearly, they therefore would produce a final value N eff ν ≃ 3.045. Recently, these effects have been reconsidered in [9]. This work combines the results for the non instantaneous neutrino decoupling obtained by a numerical calculation in [10] and then replacing the ratio T /T ν = (11/4) 1/3 with the value obtained by considering QED corrections. This procedure may provide a good estimate, but it is worth pointing out that, since QED effects modify several points of the Boltzmann equations describing neutrino decoupling, it is advisable to study any possible interplay of the two effects. Actually a precision calculation of this kind, where both these effects are included at the same time, is still lacking, and it is the aim of the present paper. We will show that this interplay ends up in a 10% smaller total correction to N eff ν than what would be obtained by simply adding the two contributions. We also notice that the calculation of [10] seems to underestimate 1 the corrections from non instantaneous neutrino decoupling with respect to the result found in [6][7][8]. Our analysis is based on a numerical code described in [8], which has been modified to take into account QED finite temperature effects.
From the observational point of view the effects considered in this paper are too small to influence Big Bang Nucleosynthesis (BBN), since they produce a change in the 4 He mass fraction, ∆Y ( 4 He) ∼ 10 −4 [8,11,12], which is smaller than the actual theoretical, 5 10 −4 [12,13], and experimental, > ∼ 2 10 −3 [14], uncertainties on this quantity. This change is actually slightly larger when one takes into account the effect of flavour neutrino oscillations, as recently shown in [15]. More promisingly, they might be detected via future precision CMB anisotropy measurements at high multipoles, since cosmic variance prevents their resolution on scales probed by present balloon experiments. According to a recent analysis [16] (see also [13,[17][18][19][20][21]), CMB temperature measurements by the Planck satellite experiment will be able to provide a measure of the relativistic energy density with a precision of about ∆N eff ν ≃ 0.2 , without any strong prior on the other cosmic parameters. Furthermore, as discussed in [22], the situation is more promising if polarization measurements will be available, and some stronger priors are imposed on other cosmological parameters, enforced for example by independent measurements. In this case, ∆N eff ν would be determined with an accuracy comparable or, according to [22], even higher than the order of magnitude of the effects we are here considering. We notice that in the CMB data analysis the presence of a non vanishing ∆N eff ν is already considered with, for example, a recommended default value for three massless active neutrinos N eff ν = 3.04 in the CMBFAST code [23] . A precise check of this important parameter is therefore mandatory.
The paper is organized as follows. In Section 2 we summarize the set of equations and our numerical approach to compute the neutrino distortion due to their incomplete decoupling during the e ± annihilation phase. In Section 3 we also consider the corrections to the equation of state of e ± and γ plasma due to QED effects at finite temperature, showing how these effects modify the numerical computation of Section 2, and we present our results. We also discuss the general case of extra relativistic degrees of freedom and of massive active neutrinos in Section 4. Finally our conclusions are reported in Section 5.

Corrections from to the non instantaneous decoupling
The effects of non-instantaneous neutrino decoupling have been addressed in several studies, either based on analytical methods [2], on numerical analysis with some simplifying approximations [24][25][26][27], or finally on full numerical computations, which solve the Boltzmann equations until neutrinos are completely decoupled [6][7][8]10,28].
The procedure is easily summarized. Following the notation of [6,8], a convenient choice for the time variable is given by x ≡ m R, where R(t) is the scale factor of the Universe (chosen to have dimension of length) and m some reference mass, taken to be the electron mass as in [8]. We also introduce the dimensionless comoving momentum y ≡ p R , and the rescaled photon temperature z ≡ T R . In this notation, the Boltzmann equations for the neutrino distributions can be written in the form where H is the Hubble parameter, while I να represents the collisional integral, in momentum space y, for the single neutrino specie ν α , and is a functional of all neutrinos and electron/positron distributions 2 . For the processes of interest 2 Since e ± are kept in thermodynamical equilibrium with γ by the electromagnetic (see [6] for the detailed calculation), some of the integrations appearing in the I να can be analytically performed, and the collisional integrals can be reduced to two-dimensional integrals. In the range of temperatures we are interested in, electron neutrinos experience charged current interactions due to the presence of e ± in the thermal bath, while all active neutrino flavours interact via neutral current interactions. For this reason, the nonequilibrium corrections to the distribution function f νe are different from those of the other two neutrino species f νx (x denoting both µ and τ ).
The two Boltzmann equations for f νe and f νx must be supplemented by the continuity equation 4 , and ρ and P are the energy density and pressure of the γ, e ± , ν plasma. Eq. (4) can be rewritten as an evolution equation for the quantity z which gives the ratio T /T ν . In case of instantaneous neutrino decoupling, the asymptotic value of z, denoted by z 0 , results to be the well-known value (11/4) 1/3 .
As in Ref. [8] the unknown neutrino distributions are parametrized as where P i (y) are orthonormal polynomials with respect to the Fermi function weight By substituting Eq. (5) into Eqs. (3), one gets Since at high temperature the neutrinos are in thermal equilibrium, the initial condition for the coefficients is a α i = 0.
interactions, they have Fermi-Dirac distributions and share the same temperature of photons, T . We neglect the completely irrelevant e ± asymmetry. To perform the numerical computation it is necessary to truncate the infinite series in Eq. (5) at some finite value n, where the choice of n depends on the accuracy we do require on the results; it is found in [8] that n = 3 gives an accuracy of about 1% in the neutrino distortions. In this case, the asymptotic value found for z, denoted by z fin , slightly differs from the instantaneous decoupling result z 0 The final neutrino distribution functions show a non thermal behaviour due to the presence of non vanishing ripple terms. Having fixed n = 3, we rewrite the expression of Eq. (5) in the form and report the final values for the coefficients c α i in Table 1. By using Eq. (9) it is now immediate to compute the additional contribution to the neutrino energy density due to incomplete decoupling. Defining one obtains, for the values of Table 1, δρ νe /ρ 0 ν = 0.953% and δρ νx /ρ 0 ν = 0.399%. Finally, from the definition of N eff ν of Eq. (1), it is straightforward to get the following expression which gives the effective number for three massless active neutrinos for temperatures below the neutrino decoupling (∼ 1 MeV). The values found from the numerical analysis then give N eff ν = 3.0345, which fixes the additional contribution to be ∆N eff ν = 0.0345. Notice that from Eq. (11), ∆N eff ν takes two contributions. The term proportional to δz/z 0 accounts for the smaller profit that photon temperature gets from e ± annihilations, whereas the contributions due to δρ νe,x /ρ 0 ν are a measure of the non thermal behaviour of neutrino distribution functions. These two terms are indeed of the same size.
On the other hand we disagree with the results of [10], where by solving the same Boltzmann equations it is found δρ νe /ρ 0 ν = 0.607%, δρ νx /ρ 0 ν = 0.256% and δz/z 0 = −0.888·10 −3 , which finally gives N eff ν = 3.022. In both calculations [6,7] and [10], the integral-differential Boltzmann equations are solved through a grid in momentum space. The grid adopted in [10] extends to a larger momentum range than the one of [6,7]. The grid used in [6,7] is however denser in the interval 0.1 < y (m/MeV) < 20 . For neutrino distributions in thermal equilibrium, the 97.5% of the energy density comes from particles with momentum in the range 1 < y (m/MeV) < 10 . According to the analysis of [7], the different results of [10] may be due to a lack of precision in this most relevant interval. To conclude, we can only point out that the beautiful agreement of our findings with what obtained in [6,7], despite of the rather different numerical method, make us rather confident on our result for N eff ν .

Corrections from QED at finite temperature
Finite temperature QED corrections modify in several points the calculations of neutrino decoupling described till now. First, through the change on the electromagnetic plasma equation of state, they affect Eq. (4). Moreover, since e ± masses are renormalized by a finite temperature term, the r.h.s. of Eq.s (3) and (7) should be modified correspondingly. Finally, the change of energy density modifies the expansion rate H in the Boltzmann equations (3). The effects of these corrections on neutrino decoupling temperature have been considered in [29]. Their analysis is however in the framework of the instantaneous neutrino decoupling limit, therefore they do not consider the non thermal effects described in the previous Section.
The change in the electromagnetic plasma equation of state can be evaluated by first considering the corrections induced on the e ± and photon masses.
They can be obtained perturbatively by computing the loop corrections to the self-energy of these particles. For the electron/positron mass, up to order α ≡ e 2 / (4 π) we find the additional finite temperature contribution [3] 3 where E k ≡ k 2 + m 2 e . While the first two terms of this expression depend on the plasma temperature T only, the last one depends on the e ± momentum p as well. However, by averaging this term over the equilibrium e ± distribution, one easily finds that it contributes for less than 10 % to δm 2 e . For this reason we will neglect it in the following (see also [12]).
The renormalized photon mass in the electromagnetic plasma is instead given, up to order α, by [29] The corrections (12) and (13) modify the corresponding dispersion relations as E 2 i = k 2 + m 2 i + δm 2 i (T ) (i = e, γ). This in turn affects the total pressure and the energy density of the electromagnetic plasma Expanding P with respect to δm 2 e and δm 2 γ , one obtains the first order correction Note that an important factor 1/2 must be introduced for not double counting the renormalization effect on the total pressure, which now reads P = P 0 +P int  Table 2 Values of the coefficients of Eq. (9) when QED effects are included. [3], P 0 being the value of the pressure for the non-interacting particle gas. The energy density is then obtained by using P in Eq. (15).
The presence of the additional contributions P int and ρ int modify the evolution equation for z contained in Eq. (4), which now reads where with The functions K ′ (ω) and J ′ (ω) stand for the first derivative of K(ω) and J(ω) with respect to their argument. Note that neutrinos affect the final value of z [6,7] (no-QED) [10] (no-QED) [8] (no-QED) our work (QED) 3.0340±0.0003 3.022 3.0345 3.0395 Table 3 The results of the different analyses.
through the terms df νx /dx, which are not vanishing only if neutrinos have a non thermal behaviour.
In case one neglects the finite temperature QED corrections, the functions G 1 (x/z) and G 2 (x/z) in (17) vanish, and one recovers the expression reported in [8]. Notice that the presence of G 2 (x/z) in the denominator of the r.h.s. of (17) makes, at least in principle, not correct to simply sum the neutrino contribution with the QED one.
Once both effects, neutrino incomplete decoupling and QED corrections to electromagnetic plasma, are included into the code of Ref. [8], we find the new results δρ νe /ρ 0 ν = 0.935%, δρ νx /ρ 0 ν = 0.390% and δz/z 0 = −1.841·10 −3 which give N eff ν = 3.0395. The values of the c α i coefficients can be found in Table 2. In Table 3 we report a comprehensive summary of our results. Comparing the last two columns of this table, we see that introducing the finite temperature QED corrections in the non instantaneous decoupling scenario, leads to a change on the effective number of neutrinos of 0.005 , which is a factor 2 less than what has been obtained in Refs. [3,12]. This is actually due to the interplay of incomplete decoupling of neutrinos and plasma effect (see the r.h.s of Eq. (17)), which therefore, at the level of accuracy we are considering, cannot be considered separately and then added linearly. Notice that the increase of N eff ν is mainly due to a variation of δz/z 0 , the changes on δρ να /ρ 0 ν being much smaller. Remarkably, the overall result N eff ν = 3.0395 turns out to be in excellent agreement with the recommended default value N eff ν = 3.04 used in the CMBFAST code [23].

Extra relativistic degrees of freedom and massive neutrinos
We now consider the possibility that, at the stage of neutrino decoupling, there are extra relativistic degrees of freedom, provided by some X field excitations, which are assumed to have a thermal distribution with some temperature T X . Their contribution ρ X to the total energy density of the Universe can be parametrized in terms of an additional contribution ∆N eff ν in the effective number of neutrinos, as defined in eq. (2) where and we have defined z X = T X R. The parameter g X depends on the spin (g X = 1 for a real scalar, g X = 7/4 for a Weyl spinor, etc.) as well as on the additional internal degrees of freedom of the X particles. Notice that if the X excitations have decoupled between µ ± and e ± annihilation phases we simply have N X = 4/7g X . For an earlier decoupling we have instead N X < 4/7g X .
The presence of ρ X , apart from introducing a new contribution to N eff ν , slightly affects the results obtained in the previous Sections, namely the relative change of neutrino energy density δρ νe,x /ρ 0 ν induced by incomplete neutrino decoupling, as well as the asymptotic photon temperature z fin . In fact since their energy density increases the expansion rate of the Universe, we expect δρ νe,x /ρ 0 ν to decrease with growing N X , and the ratio z 0 /z fin to become closer to unity, since neutrino decoupling process starts at earlier times. If we denote with δρ νe,x (N X )/ρ 0 ν and z fin (N X ) the new values for these parameters as functions of N X we therefore have Since we are interested to those contributions to N eff ν corresponding to species which are relativistic for temperatures in the MeV range, we can severely bound ∆N eff ν using the results from BBN, which leads to the conservative bound ∆N eff ν ≤ 1 (see for example [30] for a recent analysis). Using our numerical code we have evaluated how z fin (N X )/z 0 and δρ νe (N X )/ρ 0 ν change when N X varies in the range 0 < N X < 1. In this interval, with an accuracy of 10 −4 on N eff ν , we find with β = 0.0014 and where now both z 0 /z fin and δρ νe,x /ρ 0 ν in this expression are those reported in the last column of Table 3, i.e. with N X = 0. The changes of these parameters with N X is encoded in an additional term, which is weighted by the small parameter β.
It should be clear from what we said before that only species which are relativistic at the neutrino decoupling, down to the e ± annihilation stage, contributes for this extra term. It is known that the effective number of neutrinos at later stages, as for example at recombination, is much less constrained from data [16][17][18][19][20][21]30,31] and it is still possible the case that energy density in the form of relativistic plasma may be injected only well after the Big Bang Nucleosynthesis epoch. This implies that the value of N eff ν may well be rather different at the CMB and BBN epochs. In this case, Eq. (26) at recombination reads where N BBN X and ∆N CM B X are the contribution of species which are relativistic at the BBN and neutrino decoupling, and recombination epochs, respectively.
We finally consider the case of massive active neutrinos. This nowadays plausible scenario affects, in general, both the expectations for CMB anisotropy and large scale structure formation. As it is seems more and more clear from neutrino experiments on solar and atmospheric neutrinos, it is unlikely that neutrino mass differences may be greater than ∼ 0.1 eV [32], unless we enlarge the standard scenario introducing sterile neutrino states. At the same time, it is also quite well established from Tritium decay data that ν e mass is bound to be smaller or, at most, of the order of 1 eV. It is therefore clear that in this scenario all neutrino masses are completely negligible as far as their decoupling is concerned. However if their values is as large as ∼ 1 eV, they start to be relevant as the temperature approaches the range relevant for CMB, and the presence of a finite mass modify of course the neutrino contribution to N eff ν . It is interesting to consider how the effects of incomplete decoupling and QED thermal effects studied in the previous Section would now affect the neutrino energy density. This is conveniently parametrized by the (time dependent) quantity which is defined as in Section 2, but for a neutrino with a finite mass m α . Since for m α ∼ 1 eV active neutrino are fully relativistic for temperatures of the order of MeV , so we can still take at the e ± annihilation phase the results of Section 3, it is easy to see that δρ να (m)/ρ 0 να (m) is given by In Figure 1, we plot δρ να (m α )/ρ 0 να (m α ) for the reference choice m α = 1 eV, using the coefficients c α i of Table 2. The x range corresponds to a variation of the temperature from T ∼ 1 keV till T ∼ 10 −3 eV. We see that the effect of incomplete decoupling and QED plasma masses for e ± and γ on neutrino energy density decreases as the temperature becomes comparable with the value of the neutrino mass. Notice that, from (29), at very low values of T , when the neutrino energy is dominated by the mass term, δρ να (m α )/ρ 0 να (m α ) reaches an asymptotic value, which is given by the change in the neutrino number density due to the above effects, normalized with the number density for a pure thermal Fermi-Dirac distribution.

Conclusion
In this paper we have considered how the two effects due to incomplete neutrino decoupling and QED corrections to the e ± and photon plasma equation of state affect the effective neutrino degrees of freedom, a crucial parameter for many cosmological observables.
The main result of our analysis, which have been carried out by numerically solving the set of Boltzmann equations describing neutrino decoupling, is a value for N eff ν = 3.0395, for the standard case of three active neutrino flavours. The issue is certainly not new, but we hope our study will contribute to reach a better accuracy in the theoretical determination of relativistic energy density of the Universe. In particular we have stressed that there is quite an interesting interplay between the two effects we have considered, so that at the level of accuracy of 0.005 on N eff ν , their corrections to the neutrino energy density cannot be naively summed as they were fully independent.
We have also considered the less standard scenario of extra species contributing to N eff ν , and how they affect the calculation of the thermal distortion of neutrino distribution, as well as the final value for photon temperature after the e ± annihilation phase.
Despite of the smallness of the corrections we have been concerned with in this paper, nevertheless a careful analysis of data on the CMB anisotropy spectrum at large multipoles may be able in the near future to reveal their effects. Estimated sensitivity of Planck satellite experiment to relativistic energy density is at least of the order of ∆N eff ν ≃ 0.2, but it is strongly improved, ∆N eff ν ≃ 0.005 [22], when including polarization measurements and strong priors.