Relic neutrino decoupling with flavour oscillations revisited

We study the decoupling process of neutrinos in the early universe in the presence of three-flavour oscillations. The evolution of the neutrino spectra is found by solving the corresponding momentum-dependent kinetic equations for the neutrino density matrix, including for the first time the proper collision integrals for both diagonal and off-diagonal elements. This improved calculation modifies the evolution of the off-diagonal elements of the neutrino density matrix and changes the deviation from equilibrium of the frozen neutrino spectra. However, it does not vary the contribution of neutrinos to the cosmological energy density in the form of radiation, usually expressed in terms of the effective number of neutrinos, N_eff. We find a value of N_eff=3.045, in agreement with previous theoretical calculations and consistent with the latest analysis of Planck data. This result does not depend on the ordering of neutrino masses. We also consider the effect of non-standard neutrino-electron interactions (NSI), predicted in many theoretical models where neutrinos acquire mass. For two sets of NSI parameters allowed by present data, we find that N_eff can be reduced down to 3.040 or enhanced up to 3.059.


Introduction
A generic prediction of the hot big bang model is the existence of a relic background of neutrinos, in number almost as abundant as the relic photons that form the cosmic microwave background (CMB). These neutrinos were produced in the early universe, when its temperature T was large enough so that weak interactions were effective and neutrinos were in thermal contact with charged leptons and the rest of the primeval plasma. But below the so-called neutrino decoupling temperature, at T dec ∼ 2 − 3 MeV, these elusive particles do not interact any longer and propagate freely. In the simple but accurate approximation known as instantaneous decoupling, neutrinos keep the energy distribution of a relativistic fermion, unchanged but for the effect of redshift of physical momentum, and do not share the entropy release from the annihilation of electron-positron pairs into photons that occurs at temperatures below T dec . In such a case, entropy conservation leads to the well-known ratio of temperatures of relic photons and neutrinos, T γ /T ν = (11/4) 1/3 , and to a contribution of neutrinos to the cosmological energy density of radiation given, in terms of the effective number of neutrinos, by N eff = 3 [1].
The fact that neutrino decoupling and e ± annihilations are processes quite close in time leads to the existence of relic interactions between electrons, positrons and neutrinos at cosmological temperatures smaller than T dec . These processes are more efficient for neutrinos with larger momenta, leading to non-thermal distortions in the neutrino spectra at the percent level and a slightly smaller increase of the comoving photon temperature, as calculated in [2][3][4][5][6][7][8][9][10][11][12][13][14]. The decoupling process of neutrinos was studied in these previous analyses with increasing precision, in particular with respect to the numerical solution of the kinetic equations for the neutrino distribution functions, including finite temperature QED radiative corrections and, in some papers, flavour neutrino oscillations. For instance, in [12] it was shown that neutrino mixing is important for fixing the final flavour spectral distortions and their small effect on primordial nucleosynthesis, but oscillations do not modify the contribution of neutrinos to the radiation energy density, that was found to be N eff = 3.046. Other recent works have obtained, in absence of neutrino oscillations but including QED corrections, values between 3.044 [13] and 3.052 [14].
The value found for N eff is thus only 0.04 − 0.05 larger than the expected number in the instantaneous decoupling approximation, but it fixes the contribution of neutrinos in the standard case that must be included in the cosmological model. Only very recently the analysis of cosmological data has achieved a level of precision that is clearly below one unit in N eff . In particular, the data on CMB anisotropies from the Planck satellite, in combination with baryon acoustic oscillations (BAO) measurements, lead to an allowed range of N eff = 3.15 ± 0.23 (at 68% CL, data from Planck TT+lowP+BAO), that is restricted to N eff = 3.04 ± 0.18 if the complete polarization likelihood from Planck is used (see [15] for all details on the different analyses and combinations of data). Forthcoming cosmological data from future CMB experiments, large-volume galaxy surveys, etc, are expected to improve the sensitivity on N eff , possibly down to values such as 0.02 − 0.04.
Prompted by the increased accuracy on the measurement of N eff , in this paper we revisit the process of neutrino decoupling in the early universe with the aim of checking some of the approximations assumed in previous works. In particular, here we solve the Boltzmann equations for the neutrino density matrix with full collision terms, for both the diagonal and off-diagonal terms, avoiding for the latter the so-called damping factors used in other analyses. We also consider the present best-fit values of neutrino mixing parameters and investigate whether the results depend on the two possible ways of ordering the neutrino masses, still allowed by current oscillation experiments. In addition, we take the opportunity to update the analysis performed in [16] and consider how neutrino decoupling is affected by the presence of non-standard neutrino-electron interactions.
This paper is organized as follows. In Sec. 2 we present the Boltzmann equations that we need to solve in order to follow the evolution of the neutrino spectra, for both the standard scenario and in the presence of non-standard neutrino-electron interactions. We discuss the approximations with respect to previous works, and comment on the main computation and technical issues. Our main results on the decoupling of cosmological neutrinos are described in Sec. 3, where we also include the values obtained for the final neutrino spectra and the corresponding contribution to N eff for each case. Finally, we present our conclusions.
2 Neutrino decoupling including three-flavour oscillations

Boltzmann equations
The process of neutrino decoupling in the early universe takes place at a temperature of the order of MeV, when weak interactions are no longer effective to keep neutrinos in thermal contact with electrons, positrons, and, indirectly, with photons. It has been shown in previous works (see e.g. [17]) that neutrino oscillations become effective at similar temperatures, if the mixing parameters have the values required to explain the data from solar, atmospheric, reactor and accelerator neutrino experiments. In order to take into account the effects of both interactions and oscillations, we consider the evolution of the neutrino density matrix p , where the diagonal terms are the occupation numbers f να of flavour neutrinos with momentum p and the off-diagonal terms, nonzero in the presence of mixing, are described by the real parameters a i , b i and c i (i = 1, 2). Since we neglect a potential neutrino asymmetry, quite constrained due to the measured value of the mixing angle θ 13 [18,19], neutrinos and antineutrinos share the same density matrices. The equations of motion for the neutrino density matrices are the corresponding set of Boltzmann equations in an expanding universe [20] (see also [21]) where H is the Hubble parameter, G F is the Fermi constant and m W the W boson mass. The commutator term includes the vacuum oscillation term proportional to M F , the masssquared matrix in the flavour basis that is related to the diagonal one in the mass basis diag(m 2 1 , m 2 2 , m 2 3 ) via the neutrino mixing matrix U , as shown in [12]. Here we consider the following values for the neutrino mixing parameters (assuming CP conservation) where s 2 ij ≡ sin 2 θ ij , that correspond to the best-fit values from the global analysis of neutrino oscillation experiments in [22] (using the values found in more recent analyses [23,24] does not modify our results). The two choices correspond to the two ways of ordering neutrino masses: normal (NH) and inverted (IH) hierarchy.
The second term inside the commutator in eq. (2.2) corresponds to neutrino forward scattering in the primeval medium, proportional to G F and the diagonal matrix E, that represents the energy densities of charged leptons. For the range of temperatures we are interested in, only the contribution of electrons and positrons is relevant (i.e. the only nonzero component of E is E 11 = ρ ee ≡ ρ e − + ρ e + ). We neglect other refractive terms that appear in the general form of eq. (2.2): the usual term proportional to the charged-lepton asymmetries (always smaller than either the vacuum or the E term) and the contribution of neutrinoneutrino interactions, that vanishes for zero neutrino-antineutrino asymmetry.
The last term in eq. (2.2) includes the effect of non-forward neutrino interactions, proportional to G 2 F . The most general form of I[ p (t)], for each of the interaction process with neutrinos, is a matrix with a collision integral for each of the components of p [20]. For instance, the collision term from the annihilation process ν(p 1 ) +ν(p 2 ) ↔ e − (p 3 ) + e + (p 4 ) is ,ē (4) ) + 2(p 1 · p 2 )m 2 e F RL ann (ν (1) ,ν (2) , e (3) ,ē (4) ) + F LR ann (ν (1) ,ν (2) , e (3) ,ē (4) ) , (2.4) where m e is the electron mass and we have defined Here f i = f e − (p i ) andf i = f e + (p i ) are the distribution functions of charged leptons and G a is a matrix of couplings (a = L, R, left-handed L or right-handed R). In the absence of non-standard neutrino interactions, G L,R are diagonal with components with θ W the weak mixing angle. For the scattering process with electrons we have and for the scattering with positrons In both cases we use the definition where f i = f (p i ) is the distribution function of the electron or positron, depending on which particle we are considering for the scattering. The collision terms for the neutrino-neutrino processes are similar to eqs. (2.4), (2.8) and (2.9), but with more complicated expressions that are non-linear in the neutrino density matrices [20]. As in previous analyses of neutrino decoupling in the early universe [6,7,[10][11][12][13][14], in our calculations the collision terms for the diagonal components αα (the distribution functions of flavour neutrinos f να ) are solved numerically after they are analytically reduced to twodimensional integrals, following the process described in [7]. However, previous works have approximated the off-diagonal collision terms of I αβ [ p ] in eq. (2.2) as damping factors with an expression −D αβ (p, t) αβ , where α = β and the D functions can be obtained under some approximations in a similar way as in [25][26][27]. With this prescription, one assumes that each time neutrinos interact in a process that distinguishes among flavors they collapse into weak-interaction eigenstates. In our analysis we relax that approximation and deal with the off-diagonal collision terms in the same way as for the diagonal variables for the processes that involve neutrinos and e ± . We do not consider the full off-diagonal terms for weak reactions involving only neutrinos, such as νν ↔ νν or νν ↔ νν, that play a less important role in the process of neutrino heating (see e.g. the detailed discussion in [14]). This assumption also allows us to reduce the computing time.
The Boltzmann equations for the neutrino density matrices must be solved simultaneously with the continuity equation for the total energy density of radiation ρ, where ρ and P are the total energy density and pressure, respectively, of the relativistic plasma, that includes the electromagnetic components γ and e ± (in equilibrium with temperature T γ ) and the three neutrino states. This equation gives the time evolution of the photon temperature T γ . Finally, for the cosmological period of interest, finite temperature QED corrections induce effective electron and photon masses, that in turn modify the equation of state of the electromagnetic plasma and lead to small changes (via δm 2 e ) in the collision rates of the processes involving e ± . These modifications are included in our calculations as described in [11,28,29].

Non-standard neutrino-electron interactions
Neutrino decoupling in the early universe would be modified in the presence of interactions beyond the weak processes present in the Standard Model (SM) of particle physics. Since neutrinos are massless in the framework of the SM, many extended models have been proposed were neutrinos acquire mass. Most of these models naturally lead to new non-standard interactions (NSI) involving neutrinos, whose value strongly depends on the particular model.
In this paper we follow refs. [16,30,31] and assume that new physics induces NSI only through the four-fermion operators (νν)(f f ), where f is a charged lepton or a quark. Moreover, for the decoupling of cosmological neutrinos only the NSI that involve electrons are important. They are described, together with the standard weak interactions, by the effective Lagrangian which, in addition to the SM term [16], contains the contribution from NSI where α, β = e, µ, τ and P = L, R = (1 ∓ γ 5 )/2 are the chiral operators. The NSI parameters ε P αβ can lead to a flavour-changing contribution (α = β) or induce a breaking of lepton universality (α = β). Their values are constrained by data from various laboratory experiments, as recently reviewed in [32], but due to possible cancellations the limits depend on the number of NSI parameters that are simultaneously included.
The effect of NSI on relic neutrino decoupling was studied in [16], where it was analytically shown how NSI parameters lead, in general, to a larger interaction between neutrinos and e ± , although in some particular regions of ε they could reduce the collision rate (see Fig. The Boltzmann equations in eq. (2.2) are modified due to NSI interactions as follows where the matrix is no longer diagonal and contains the effect of non-standard interactions via the combinations of NSI parameters The presence of NSI also modifies the statistical factors of eqs. (2.5) and (2.10) in the collision integrals. The coupling matrices adopt the form (2.17) In our analysis we will restrict the calculation of neutrino decoupling with NSI and full collision terms, as described in sec. 2.1, to two combinations of ε parameters with ε P ee and ε P τ τ , that are still allowed by present data [32,33] . We do not consider the contribution of ε that affect muon-neutrino interactions because they are strongly suppressed [34].

Computation and technical issues
We solve the system of Boltzmann equations (2.2) and the continuity equation (2.11) with the same comoving variables as in [12], where we have chosen as arbitrary mass scale the electron mass, p is the neutrino momentum, T γ is the photon temperature and a is the cosmological scale factor, normalized so that a(t) → 1/T γ at large temperatures. The expressions of the Boltzmann equations in terms of these comoving variables are listed in appendix A.
The system of kinetic equations to be solved is integro-differential due to the presence of the neutrino density matrices in the collision terms. In previous analyses of neutrino decoupling, this system was solved either using a discretization in a grid for the neutrino momenta y i as in [6-9, 12, 14] or expanding the non-thermal distortions of the neutrino distribution functions in moments as in [10,11,13]. Here we choose the first method and use a grid of values for the neutrino momenta y i in the range [0.01, 20]. We found good convergence for 100 bins. The system of integro-differential equations is solved using the free FORTRAN77 library ODEPACK, which includes several solvers for differential equations.
We start the numerical computation at a value x in = m e /T 0 γ , where T 0 γ = 10 MeV, while neutrinos are still in good thermal contact with the electromagnetic plasma via weak interactions and flavour oscillations are suppressed by the medium. Therefore, the initial values of the off-diagonal components of (y i ) are zero, while the diagonal terms are [exp(y i /z in )+1] −1 , with z in = 1.00003 the initial value of the dimensionless photon temperature. This number is found solving the continuity equation (2.11) with neutrinos fully coupled to electrons and positrons [9]. The system is solved up to a final value x fin = 30 when the neutrino distribution functions and z have reached their asymptotic values.

Standard case with flavour oscillations
First, we consider the case with standard weak interactions, with and without neutrino oscillations. In the presence of neutrino mixing, we are mainly interested in the possible effects of including the full expressions for the off-diagonal collision terms in the kinetic equations. On the other hand, once the collision terms are fixed, we also want to check whether there are differences between the two options for the neutrino mass hierarchy or the present best-fit values of the mixing parameters lead to modifications with respect to the results in [12]. We show in Fig. 1 the evolution of the flavour neutrino spectra for a particular value of the neutrino momentum (y = 5). The corresponding distortions f να /f eq , where f eq = [exp(y) + 1] −1 , are shown as a function of the photon temperature or the cosmological expansion (x). The behaviour of this evolution has been discussed in previous works (see e.g. [7,12]). At large temperatures (T γ 2 MeV) neutrinos are still interacting with e ± and their energy spectra keeps an equilibrium form with T γ . Later the cosmological expansion renders less efficient the weak processes and neutrinos decouple from the electromagnetic plasma in a momentumdependent way. The residual interactions lead to spectral distortions for neutrinos, which are larger for the electronic flavour. Neutrino oscillations, effective after T γ 3 MeV when the medium potential is diluted by the expansion, reduce the difference in the spectral distortions of electron neutrinos with respect to the other flavours.
The results shown in Fig. 1 when oscillations are taken into account were found with the full collision terms (the matrix I[ y ]) for neutrino-electron processes, avoiding the damping functions for the off-diagonal elements. In order to have an idea of the difference, we show in Fig. 2 the evolution of the collision terms for the real parts of the off-diagonal components αβ in eq. (2.1), for a particular neutrino momentum (y = 5). For each case (a 1 , b 1 , c 1 ), the result of the full integral is compared with the corresponding damping term. One can see that, while the overall behaviour is similar for each case, the full calculation presents a smoother evolution starting from zero at large temperatures. It is also interesting that the largest (absolute) values of the collision terms with the full integrals appear later than in the simple damping prescription. Since collisional damping leads to a loss of flavour coherence, a difference in the off-diagonal terms can modify the evolution of neutrino oscillations and eventually change the final distortions in the neutrino spectra.
The final values of the momentum-dependent distortions in the spectra of flavour neutrinos are shown in Fig. 3. Their non-thermal character is evident and can be understood from the fact that more energetic neutrinos are interacting with e ± at smaller temperatures. Frozen f Figure 4. Same as Figure 3 for the mass neutrino states and NH. The spectral distortions of the mass states for the IH case overlap with those for NH at this scale.
each flavour state is enhanced with respect to the instantaneous decoupling approximation. The results are given in Table 1 in terms of δρ να ≡ (ρ να − ρ ν 0 )/ρ ν 0 , where ρ ν 0 is the energy density of one neutrino state that does not share the entropy release from the e ± annihilations. In all cases the enhancement is at the sub-percent level, which indicates the well-known result that the non-thermal distortions for cosmological neutrinos are very small. The main effect is a slightly larger contribution of neutrinos to the cosmological energy density of radiation ρ r , usually given, as in the case of any other contribution of relativistic particles other than photons, in terms of the effective number of neutrinos, that can be also written as Here z 0 = (11/4) 1/3 and z fin are the asymptotic values of the dimensionless photon temperature in the instantaneous decoupling approximation and for each case in our calculations, respectively. The final values of z fin for each case and the calculated ones for N eff are given in Table 1.
Our main results correspond to the cases NH and IH in Table 1, computed with the full collision terms. Both cases lead to almost equal values of the enhanced neutrino energy densities and z fin , and to the same final effective number of neutrinos, N eff = 3.045, which is almost the same as 3.046, found in [12] and usually taken as reference for the standard value of N eff . Such a difference (0.001) is smaller than the accuracy quoted in [12] and corresponds to what we estimate for our own numerical calculations (performed with a completely independent code). We conclude that the value of N eff does not depend on the ordering of neutrino masses. However, the NH and IH cases are not exactly equivalent; there is a resonance for IH when neutrino oscillations become effective (T γ ∼ 3 MeV) that leads to an interchange among the mass states, although they evolve in a similar way after the resonance for both mass hierarchies. The difference is therefore not appreciable for the range of temperatures after e ± have annihilated and neutrinos are completely decoupled, as can be seen in Table 2 and Fig. 4, where our main results are shown for the mass neutrino eigenstates ν 1,2,3 . The corresponding neutrino distribution functions are the relevant energy spectra for all cosmological calculations at temperatures much smaller than the neutrino decoupling value (MeV), both when neutrinos are still relativistic or when they start to behave as non-relativistic particles.
It is interesting to check if, fixing the mass ordering to NH, the differences shown in Fig. 2 between the off-diagonal collision terms calculated with the full integrals or the simple damping prescription, have any effect on the final values that characterize relic neutrino decoupling. We find that the value of N eff is, within the accuracy of our numerical calculations, the same in both cases (3.045), while there exist small differences in the values of the energy density distortions δρ να , that are slightly smaller (larger) for muon or tau (electron) neutrinos when the damping approximation is not used. However, one can conclude that for practical purposes, in particular concerning the final value of N eff , the inclusion of the full collision integrals in the evolution equation of the off-diagonal components of p is not necessary.
We have also included in Table 1 our results for the case when neutrino oscillations are not considered. As shown in previous analyses, the energy densities of the neutrino flavour states in absence of mixing are quite different with respect to those found with neutrino oscillations, but the final value of N eff is again 3.045. This can be compared with the values found in previous works for zero neutrino mixing but including QED corrections: 3.046 [12], 3.044 [13] and 3.052 [14]. While some of the small differences among these values could be due to different methods for the numerical evaluation, we believe that a new, specific study of the effects of finite temperature QED radiative corrections on relic neutrino decoupling could be useful in order to understand the approximations assumed.

Non-standard neutrino-electron interactions
Here our aim is not to explore the whole space of NSI parameters allowed by laboratory data, but to show the effects of the new interactions on relic neutrino decoupling with a couple of examples. We consider two sets of diagonal parameters ε P ee and ε P τ τ that are not ruled out by experimental data, as described in [33]. The two chosen sets are ε R τ τ = −ε L τ τ = 0.37; and ε R ee = −0.42 and ε L ee = −0.09, setting in both cases the rest of NSI parameters to zero, and including neutrino oscillations with a normal ordering of neutrino masses.
We chose these values because they correspond to allowed NSI parameters that are the furthest from the standard case, where all ε = 0. In addition, with this selection of values we want to illustrate the two possible effects on neutrino decoupling: the first case (ε τ τ = 0) leads to an increase in the effective number of neutrinos with respect to N eff = 3.045, while the case with ε ee = 0 gives a smaller value of N eff . An explanation of this effect in terms of Case the effective decoupling temperature of neutrinos can be found in [16] (see in particular their Figs. 3 and 4).
We show in Fig. 5 the evolution of the neutrino spectral distortions f να /f eq as a function of the photon temperature or our comoving variable x, for both sets of NSI parameters and for a particular neutrino momentum (y = 5). The asymptotic values of the flavour neutrino distortions as a function of y are depicted in Fig. 6. The final values of the dimensionless photon temperature, neutrino energy distortions and N eff are listed in the lower rows of Tables 1 (flavour neutrino states) and 2 (mass neutrino states). The presence of non-zero NSI parameters modifies the neutrino spectral distortions and leads to small changes in the final value of N eff : +0.014 for the case ε τ τ = 0 and −0.005 when ε ee = 0. NSI ee NSI ττ Figure 6. Same as Figure 3 for the two cases with non-standard neutrino interactions, including flavour oscillations and masses in the NH case.

Conclusions
In this paper we have revisited the decoupling process of neutrinos in the early universe in the presence of three-flavour oscillations, both with standard weak interactions and with non-standard neutrino interactions with electrons. We have calculated the evolution of the momentum spectra of neutrinos and found their final contribution to the cosmological energy density in the form of radiation, given by the effective number of neutrinos, N eff . This problem requires the numerical solution of the corresponding momentum-dependent kinetic equations for the neutrino density matrix, that we have performed including the full collision terms for neutrino-electron processes. In the case of the off-diagonal components, these collision integrals replace the damping terms usually assumed in previous analyses. We find that, while this improved calculation does change the evolution of these off-diagonal elements, it does not modify the final value of the effective number of neutrinos, given by N eff = 3.045. The impact of including the off-diagonal collision terms calculated with the full integrals or the simple damping prescription could be more important for other cases, such as active-sterile neutrino oscillations [35] or very low-reheating scenarios [36], where the spectral flavour distortions could be significantly larger than in the standard case.
We have also considered whether the value of N eff depends on the neutrino mass ordering (normal or inverted), which is still unknown from oscillation experiments. We find very small differences in the final distortions of the flavour and mass neutrino spectra for the two cases, but the same final value of N eff = 3.045, in agreement with previous theoretical calculations [12] and fully consistent with the latest analysis of Planck data [15].
In absence of neutrino mixing the process of neutrino decoupling is less complicated. We have also considered this case, currently non-standard, because it can be compared with other previous analyses. Concerning the effective number of neutrinos, we again find N eff = 3.045, a similar result to the ones obtained in refs. [12][13][14]. It is possible that some of the small differences could arise from the different methods for the numerical evaluation, but we believe that a future specific study of the effects of finite temperature QED radiative corrections on relic neutrino decoupling could be useful. In our analysis we have not calculated the small effects of neutrino heating on primordial nucleosynthesis, and in particular on the production of deuterium or 4 He, but it has been discussed in refs. [7,10,12,14].
Fixing the contribution of neutrinos to the radiation energy density in the standard case that exceeds N eff = 3 is an important input for the cosmological model. Other relativistic relics or non-standard neutrino properties could lead to a different value of N eff , usually larger than 3.045. For instance, in our analysis we have also considered the presence of non-standard neutrino-electron interactions. For two sets of NSI parameters allowed by present data, we find that the final value of N eff can be reduced down to 3.040 or enhanced up to 3.059.
The small departures of N eff from 3 that we have found in this paper, either in the standard case of neutrino decoupling with flavour oscillations or with NSI, are beyond the current level of observational precision fixed by the analysis of Planck data (alone or in combination with other cosmological observations), of the order of σ(N eff ) 0.2 [15]. However, forthcoming data from future CMB experiments or large-volume galaxy surveys could improve the sensitivity on N eff down to values such as 0.02 − 0.04 (see e.g. [37][38][39]).

A Boltzmann equations and collision terms in comoving variables
In terms of the comoving variables in eq. (2.18) one can write the Boltzmann equations (2.2) as where the bar over some quantities indicate that they are written in the comoving variables of eq. (2.18), such as the dimensionless energy densityρ = ρ(x/m e ) 4 . On the other hand, the continuity equation (2.11) can be translated into an equation for dz/dx as in [11], including the contributions of finite temperature QED corrections.