Active-sterile neutrino oscillations in the early Universe with full collision terms

Sterile neutrinos are thermalised in the early Universe via oscillations with the active neutrinos for certain mixing parameters. The most detailed calculation of this thermalisation process involves the solution of the momentum-dependent quantum kinetic equations, which track the evolution of the neutrino phase space distributions. Until now the collision terms in the quantum kinetic equations have always been approximated using equilibrium distributions, but this approximation has never been checked numerically. In this work we revisit the sterile neutrino thermalisation calculation using the full collision term, and compare the results with various existing approximations in the literature. We find a better agreement than would naively be expected, but also identify some issues with these approximations that have not been appreciated previously. These include an unphysical production of neutrinos via scattering and the importance of redistributing momentum through scattering, as well as details of Pauli blocking. Finally, we devise a new approximation scheme, which improves upon some of the shortcomings of previous schemes.


Introduction
Light sterile neutrinos with masses around an eV have long been postulated as an explanation for the anomalous results seen in a number of short-baseline neutrino oscillation experiments (see, e.g., [1] for a review). Although intriguingly each experiment appears to point to a slightly different mass value, the preferred ranges of masses and active-sterile mixing angles are such that, within standard cosmology, the sterile neutrino state inevitably becomes thermalised in the early Universe. The presence of an additional thermalised neutrino species in the eV-mass range is however in conflict with observations of the cosmic microwave background (CMB) anisotropies and the large-scale structure (LSS) distribution: the best limit from the ESA Planck mission and other astrophysical observations on the number of thermalised neutrino species currently stands at N eff = 3.04±0.18 (68% C.I.) [2]. Thus, if the short-baseline anomalies are to be interpreted in terms of active-sterile neutrino oscillations despite their mutual tension, some mechanism to reconcile the sterile state with cosmology would be required.
Several such mechanisms to circumvent the CMB/LSS limits on sterile neutrinos have been investigated in the past, some of which employ a large lepton asymmetry [3][4][5][6][7] or interactions among the sterile neutrinos [8][9][10][11][12] to suppress the production of sterile neutrinos in the early Universe. What these two mechanisms have in common is that they both delay the production of sterile neutrinos until at the earliest the neutrino decoupling epoch (T ∼ MeV), after which the active neutrino states cannot be fully repopulated even if the activesterile oscillation probability should become large. This limits the total number of neutrinos populating the Universe, and hence avoids an unacceptably large energy density in relativistic particles.
To accurately track the sterile neutrino thermalisation process-in the presence of large lepton asymmetry, self-interactions, or otherwise-requires that we solve the quantum kinetic equations (QKEs), which describe the evolution of the density matrix of the active and sterile neutrino phase space distributions in the presence of flavour oscillations and scattering (forward and non-forward). The formal expressions for all components of the QKEs up to order G 2 F are well known [13,14]. However, as the collision integrals describing the non-forward scattering are generally quite complicated and nonlinear, it is customary to approximate them in numerical solutions of the QKEs. Common approximations include neglecting Pauli blocking and feedback from the real-time neutrino phase space distributions, as well as ignoring the electron mass in the evaluation of the scattering matrix elements [7,15].
In many ways these were reasonable approximations in their time. However, as observational precision improves, it also becomes necessary to examine more closely their precise impact on the observables. A conservative and naive estimate of the error due to neglecting the Pauli blocking factors, for example, is ∼ 10% [15] for active-sterile conversion before the neutrino decoupling epoch, which in view of Planck's sensitivity to N eff already appears to be pushing the limits of validity. As the conversion temperature approaches the decoupling epoch, the details of how the active neutrino states are repopulated and their oscillations to sterile states damped through collisions can only have an even larger impact; apart from thermalisation of the sterile neutrino and hence the N eff parameter, the detailed repopulation mechanism affects in principle also the active neutrino energy spectra, which could subsequently alter the weak reaction rates during big bang nucleosynthesis (BBN).
In this paper we include for the first time the full collision integral in the solution of the QKEs for active-sterile neutrino oscillations. For simplicity we restrict our attention to a two-neutrino model without lepton asymmetries or sterile neutrino self-interaction. We focus on oscillations between the electron and the sterile neutrinos described by a mass squared difference δm 2 and vacuum mixing angle θ, although we will keep the equations general and bear in mind also the cases of sterile neutrino oscillations with muon and tau neutrinos. We assume the sterile neutrino to be heavier than the active neutrinos as this is the most interesting case for mass squared differences above δm 2 ∼ 0.1 eV 2 : current cosmological limits on the sum of the neutrino masses are in the sub-eV range [2,16,17], and would be violated if the active neutrinos were heavier than the sterile state for the δm 2 values of interest.
The rest of the paper is organised as follows. We begin with an introduction to the QKEs in section 2, where we describe the repopulation of the active neutrinos and the damping of flavour oscillations from collisions. We discuss various approximations found in the literature, and devise new, improved approximations that incorporate more of the relevant physics. In section 3 we first test our implementation of the full collision term for convergence, before proceeding to systematically compare the various approximate solutions with the full result. We find that although the deviations are smaller than expected for most neutrino mixing parameters, there are some systematic biases that can be significantly reduced using our new approximation scheme. We give our conclusions in section 4. Throughout the paper we employ a unit system in which c = = k B = 1, and express all dimensionful quantities in powers of eV.

Quantum kinetic equations
We consider oscillations between an active neutrino flavour ν α , where α = e, µ or τ , and a sterile flavour ν s in an ensemble of neutrinos in the early universe. The density matrices ρ(k) encode the flavour content and coherence of the ensemble, and are conveniently expressed in terms of polarisation vectors (P 0 (k), P(k)), i.e., where f 0 (k) denotes a reference phase space density, which we take here to be the relativistic Fermi-Dirac distribution with vanishing chemical potential, 1 and σ = {σ x , σ y , σ z } are the Pauli matrices. In this convention the active and sterile neutrino phase space distributions are given, respectively, by 1) and the QKEs that govern their evolution can be written as [13,15] Here, V(k) ≡ V x (k)x + V z (k)ẑ contains the vacuum oscillation term as well as the matter potential from forward scattering, and the components are given respectively by where G F is the Fermi constant, M Z the mass of the Z boson, n i the number density normalised to the equilibrium value, g µ/τ = 1, and g e = 1 + 4 sec 2 θ W /(n νe + nν e ) with the Weinberg angle θ W . Finally, R α (k) and D(k) are respectively the repopulation and damping terms, which we describe in more detail in section 2.1. As we assume a vanishing lepton asymmetry, the same QKEs apply to both neutrinos and antineutrinos, and n να = nν α .

Repopulation and damping
The repopulation and damping terms are integrals over the matrix elements for annihilation and elastic scattering processes. Beginning with equation (24) of [13], which includes Pauli blocking and appears also in [18], we find these integrals to be

4) and
is the 1D Dirac delta function, the summation index i runs over all spectator neutrino flavours (i.e., ν β where β = α) and the electron, while j runs in addition over all their antiparticles as well as the oscillating neutrino and antineutrino. The terms V 2 are denotes the energy of particle i, and S|M | 2 (a(p), b(k)|c(p ′ ), d(k ′ )) is the squared matrix element for the forward process a(p)b(k) → c(p ′ )d(k ′ ), summed (but not averaged) over initial and final spins, and symmetrised over identical particles in the initial and the final state. If two ν α s are present in the initial state, then S|M | 2 must additionally be multiplied by 2 to account for the fact that ν α (k)ν α (p) → . . . and ν α (p)ν α (k) → . . . constitute two identical processes.
The first part of both the repopulation and the damping integrals (2.4) and (2.5) pertains to annihilation processes, while the rest describes scattering processes. The repopulation integral (2.4) incorporates Pauli blocking in the form of additional multiplicative factors of the form (1 − f j ) for every particle j in the final state of both the forward and reverse processes, and conforms with expectations. For the damping integral (2.5), however, Pauli blocking enters in a way that may not be entirely intuitive. Compared to the expression used 2 The prefactor Ni = 1/2Ei used here differs from the definition in [13] due to the normalisation of the Dirac spinor. Our choice of normalisation gives the completeness relations u(pi)ū(pi) = / pi + mi and v(pi)v(pi) = / pi − mi.
by McKellar and Thomson [13], (2.7) we find two modifications for each interaction process: one additional term (the "first term") and new multiplicative factors in the second term. In terms of the evolution of the density matrix ρ, Pauli blocking enters the equation of motion as a multiplicative matrix factor of the form δ ij − ρ ij (see equation (24) of [13]). The appearance of the first term can be traced to the off-diagonal part of this matrix factor, while the second term includes a factor from the matrix diagonal. Because the two corrections differ in sign, they cancel one another to some extent.
A naive introduction of Pauli blocking into the damping integral (2.5) might lead one to miss the first term, which would result in an underestimation of the damping. However, as it turns out, the negative correction contained in the second term dominates anyhow when using equilibrium distributions, so that the effect of Pauli blocking is similar to what would naively be expected, albeit smaller.
In appendix A we evaluate the repopulation and damping integrals (2.4) and (2.5) using the technique described by Hahn-Woernle, Plümacher and Wong [19]. With this approach we can reduce the nine-dimensional integral to three dimensions. Of these it is possible to perform one integral analytically, but the remaining two must be calculated numerically.

Approximation schemes
The customary way to treat the collision terms is to assume most of the particles to be in equilibrium with the background photons. This simplifies the integrals so that solving the final expression is less numerically demanding. Assuming that Pauli blocking can be neglected and that all species follow equilibrium distributions except for the one being repopulated, the repopulation and damping terms evaluate in what we shall call the equilibrium approximation [4,15,20] to where Γ = C α G 2 F xT 5 , with x = k/T , C e ≈ 1.27, and C µ,τ ≈ 0.92 [20][21][22]. While this expression makes intuitive sense in terms of bringing the distribution towards equilibrium and coincides with the relaxation time approximation commonly found in the Boltzmann transport literature, it can also be derived from a firmer basis [13,15].
An alternative approximation scheme is that introduced by Chu and Cirelli [7] (CC approximation), which is based on a second quantised approach [14]. Here, the combined collision (damping and repopulation) term is [6,7] where G s,a = diag γ α s,a , 0 , with (γ e s ) 2 = 3.06 and (γ µ,τ s ) 2 = 2.22 for scattering processes, and (γ e a ) 2 = 0.50 and (γ µ,τ a ) 2 = 0.28 for annihilations [7], k ≈ 3.15 T is the average momentum, and ρ 0 ≡ f 0 . The CC approximation was originally applied in [7] to the momentumaveraged quantum rate equations, and to adapt it for the momentum-dependent quantum kinetic equations we have had to scale equation (2.10) by a factor k/ k [6] to approximate the momentum dependence. Evaluating the matrix products, we finḋ which can be further recast as expressions for repopulation and damping, compatible with equation (2.2).

Repopulation
Observe that in the CC approximation only the annihilation processes contribute to the repopulation term (2.12). This approximation is reasonable for the momentum-averaged quantum rate equations where the expression was first used, but is not strictly correct in the momentum-dependent quantum kinetic equations where elastic scattering processes is also expected to contribute to the equilibration of individual momentum bins. The equilibrium approximation (2.8) on the other hand does take into account scattering processes. However, in doing so, it also sacrifices the principle of detailed balance in that scattering processes can only redistribute momentum but not repopulate a momentum bin without reference to the population of other bins. A better solution would be to separate the two contributions, but this would require a more advanced treatment than a simple relaxation towards one equilibrium distribution.
Here, we develop a new repopulation approximation scheme that keeps the annihilation and scattering terms separate. We call this the A/S approximation. Looking first at annihilations, the full expression for ν ανα annihilation into a andā is an integral over (f a fā − f να fν α ), where we have ignored Pauli blocking factors. In the equilibrium approximation, it is assumed that f i ≈ f 0 for all particles i = a,ā,ν α but ν α . For a andā this is a good assumption. Forν α , however, we know that it overestimates fν α if nν α < 1. Thus, to compensate for the depletion ofν α , we adopt fν α ≈ n να f 0 (remembering that nν α = n να ) in the annihilation part of the A/S approximation, while keeping f i = f 0 for i = a,ā.
For the scattering processes on the other hand, we must take care to preserve the neutrino number density at all times. We accomplish this in two ways, depending on whether the scattering process involves energy and momentum exchange between the ν α population and an external bath. ν β , all of which are assumed to be in equilibrium with the photons, serve to drive the ν α population to a thermal distribution with temperature equal to the photon temperature T . We therefore model the process using a repopulation term of the relaxation form (2.8), but with the equilibrium distribution to which f να relaxes replaced by f eq, scat = 1 e x−ξ + 1 , (2.14) where µ = ξT is a pseudo-chemical potential. 3 The ξ parameter can be determined in practice from the condition that the scattering contribution toṅ να must be zero, or equivalently, that the third (not the second) kinetic moments of f eq, scat and f να must be equal (because the collision rate is proportional to momentum).

2.
Vanishing energy and momentum exchange. Scattering amongst ν α andν α conserves both the energy and number densities of the active oscillating neutrinos. Energy conservation implies that the ν α population could relax to a thermal distribution with a temperature T να different from the photon temperature T , i.e., Here, µ να and T να are determined from the combined condition that the scattering contributions toṅ να andṄ να must both be zero, where N να is the energy density normalised to the equilibrium value. The former constraint is equivalent to requiring equality between the third kinetic moments of f eq, να and f να , while the latter constraint calls for equality between the fourth moments.
Thus, the full repopulation term in the A/S approximation can now be expressed as (2.16) where C e,a = 0.180, C e,s = 0.718, C e,ν = 0.407, C µ/τ,a = 0.102, C µ/τ,s = 0.407, and C µ/τ,ν = 0.407. Full expressions for the coefficients can be found in appendix B. As it turns out, the separation of the scattering processes into vanishing and non-vanishing energy exchange with an external bath is not crucial for the accuracy of the approximation; it has been included here mainly for consistency. If it were to be abandoned, one could simply truncate equation (2.16) at the second term, and obtain a new coefficient for the scattering term by adding together the numerical values of C α,s and C α,ν given above.
Finally, note that equation (2.16) does not accommodate a large lepton asymmetry; if one were present, it would be necessary to reexamine the assumptions behind the equilibrium approximation (2.8) and all subsequent approximations that lead to (2.16).

Damping
The damping term in equation (2.5) is affected by Pauli blocking in two ways as already discussed in section 2.1. As the negative correction dominates when considering equilibrium distributions, D eq (k) in the equilibrium approximation (2.9) tends to overestimate the amount of damping. Conversely, Chu and Cirelli [7] included only the diagonal Pauli blocking terms when calculating the numerical coefficients for D CC (k) in equation (2.13), thereby underestimating the damping. 4 Here, we propose to evaluate the damping integral (2.5) again with the approximations f i ≈ f 0 for i = ν α ,ν α and f να = fν α ≈ n να f 0 . Due to Pauli blocking, these approximations lead to a damping term that is quadratic in n να : where C e,2 = −0.020, C e,1 = 0.569, C e,0 = 0.692, C µ/τ,2 = −0.020, C µ/τ,1 = 0.499, and C µ/τ,0 = 0.392, the negative C α,2 values reflecting the origin of the n 2 να term in the negative part of the Pauli blocking factors (expressions for the coefficients are given in appendix B). For convenience we shall continue to call this the A/S approximation, although, unlike in the repopulation treatment, there is no strict separation between annihilation and scattering; the coefficient C α,1 is predominantly from annihilation, while C α,0 comes mainly from scattering.
Lastly, we note that it is also possible to capture some of the k-dependence of the coefficients by omitting the k-integral in the computation of the damping coefficients (see appendix B), and instead fit the result with a function of the form where a, b, c, d, and g are constants determined by the fit. This kind of expression improves the agreement between the approximation and the full treatment insofar as it reproduces the sterile neutrino spectrum with greater success than the constant coefficients. It does not however lead to significant changes in ∆N eff and in the active spectrum which are often the most interesting quantities. For this reason we do not incorporate C α,fit in the A/S approximation.

Numerical implementation
We employ a modified version of the public code lasagna 5 [23] to solve the QKEs assuming first the full collision term (2.4) and (2.5), and then in the various approximation schemes discussed above. The QKEs are solved on a fixed momentum grid, with the explicit requirement thatṖ −Ṗ = 0 so as to enforce a zero lepton asymmetry.

Approximation schemes
Implementation of the approximate collision terms is straightforward in the equilibrium and the CC approximation. In the A/S approximation, however, additional root-finding routines are required to evaluate the chemical potential µ of the normal scattering term (2.14), as well as the temperature T να and chemical potential µ να of the ν α ν α -scattering term (2.15). The chemical potential of the normal scattering term satisfies the equation which can be solved numerically. The parameter space of the ν α ν α -scattering term, on the other hand, is two-dimensional and subject to the conditions In order to solve for µ να and T να , we first reduce the problem to one dimension by constructing the ratio using equation (2.19) and (2.20). Here Li s (z) denotes the polylogarithm. Since the RHS depends on µ να /T να but not directly on T να , equation (2.21) can be solved immediately for µ να /T να using a one-dimensional root-finding algorithm.
Once we have established µ να /T να , equation (2.20) can be evaluated explicitly for the temperature T να . We find: where T is the photon temperature, and again we have x = k/T . Finally, we take the µ να and T να values obtained from equations (2.21) and (2.22), and further tune them by solving the discretised moments equations (2.19) and (2.20) using a 2D Newton's method initialised with these numbers. This last step ensures that the chosen µ να and T να values do indeed satisfy conservation of number and energy densities; the untuned values might not have this desired property because of the discretisation of momentum space. In practice, however, the amount of tuning is quite small.

Full collision term
For the full repopulation and damping terms (2.4) and (2.5), we use the following tricks to simplify the calculations. Consider first the repopulation integral, which we split into three parts: where the second subscript on the RHS labels the contributing processes. We use equilibrium distributions f eq for the electrons and positrons in the processes which should be a good assumption as these particles are tightly coupled via electromagnetic interactions. This assumption enables the pre-evaluation (as opposed to real-time evaluation while solving the QKEs) of one of the two energy integrals in each of equations (A.17), (A. 18) and (A. 20), so that the contribution of processes (2.23) and (2.24) to R α can be expressed as where the integration kernelsR α,x (k, k ′ , p) encode the kinematics of the interaction processes, and can be read off equations (A.17), (A.18) and (A.20). We note also that p ′ is fixed by energy and momentum conservation once a combination of k, k ′ , p has been specified.
In the limit of a massless electron, R α,e,s,i and R α,e,a,i of equation (2.26) scale trivially with temperature as ∝ T 4 . However, as the temperature drops below ∼ 1 MeV, the massless electron approximation also becomes increasingly tenuous. Reinstating a nonzero electron mass significantly complicates the temperature dependence of R α,e,s,i and R α,e,a,i ; we handle this by pre-evaluating equation (2.26) on a grid in (k, p, T ) (or (k, k ′ , T )), which we interpolate in real time using a 3D spline when solving the QKEs.
We use the same equilibrium approximation for the spectator neutrinos ν β and antineutrinosν β in the processes where α = β, and (−) ν β are assumed to be non-oscillating. Thus R α,β and the associated kernels R α,β,s,i and R α,β,a,i are, save for the e → β relabelling, formally given by equations (2.25) and (2.26) respectively, and the corresponding expressions forR α,x (k, k ′ , p) can be read off equations (A.13), (A.14) and (A.19). Since the massless neutrino limit always holds in our considerations, the temperature dependence of R α,β,s,i and R α,β,a,i is trivial, and the integrals need only be pre-evaluated on one 2D (k, p) (or (k, k ′ )) grid.
Finally, for those processes involving only the active oscillating neutrinos and/or antineutrinos, i.e., we evaluate the full double integrals (A.15) and (A.16) without approximation, since the momentum distributions of ν α andν α are a priori unknown quantities.
For the damping term, we see from equations (2.4) and (2.5) that it differs from the repopulation term only in the missing factors of f (k) and 1 − f (k). This means that the same simplifications of the repopulation integral discussed above and consequently all of the pre-evaluated integrals apply also to the damping term. For example, contribution of processes (2.23) and (2.24) to D α can be expressed as where R α,e,s,i and R α,e,a,i are exactly the pre-evaluated quantities of equation (2.26).

Numerical Results
The main quantity we wish to study is the change in the energy density of neutrinos due to a sterile neutrino population, This one quantity affects directly the Hubble expansion rate, making it possible to constrain ∆N eff using various cosmological observations.

Numerical convergence
An important concern in the solution of the QKEs is detailed balance. If detailed balance is not fulfilled, at least so to a very good approximation, it will lead to unphysical excitation of the sterile neutrinos. As discussed in section 2.2.1, detailed balance depends on the implementation of and approximations applied to the repopulation term; it is violated, for example, by the equilibrium approximation already at the level of the equation. The full repopulation term, even in the presence of simplifications introduced in section 2.3, preserves detailed balance in principle. In practice, however, numerical solution of the QKEs using the full collision term requires that we sum over a set of discretised momentum bins at every time step. This discretisation can potentially violate detailed balance, unless the number of momentum bins is sufficiently large. In figure 1, we solve the QKEs for the normalised neutrino number density n ν and effective energy density N eff using the benchmark mixing parameter values δm 2 = 0.1 eV 2 and sin 2 2θ = 0.025, employing variously 40, 80, 100 and 150 bins; the outcomes are expressed relative to the 200 bin results. Beyond 80 bins, convergence towards the 200 bin results to better than 0.002 across the whole temperature range of interest is immediately manifest.
Comparing the different bin choices, the offset at high temperatures originates simply in the improved representation of the distribution function as we increase the number of bins. The additional discrepancy at T 10 eV is likely to have arisen from a very small deviation from detailed balance, since this is temperature regime at which the thermalisation process is most efficient. For the remainder of the analysis we use 100 bins which gives an absolute error of ∼ 0.001 for the benchmark mixing parameter values. This choice is a compromise between accuracy and speed, as the evaluation of the collision integrals scales with the number of momentum bins cubed: higher resolutions rapidly become prohibitively expensive in terms of CPU time.  The full collision term sources from a variety of scattering and annihilation processes of the oscillating active neutrinos with electrons, positrons, spectator active neutrinos, and the oscillating active neutrinos themselves. Since we have computed their individual contributions explicitly, we can now also study their individual effects on the sterile neutrino thermalisation process. This is a sanity check, but also serves to gauge the level to which our implementation of repopulation conserves energy and number densities and hence fulfils detailed balance. To this end, we consider in figure 2 two scenarios without annihilation, and compute the corresponding changes in the neutrino energy and number densities relative to the no-oscillation case for the benchmark mixing parameter values.
In the first scenario, we include all scattering processes (red solid and dot-dash lines) and find that the neutrino number density is conserved to an excellent degree, while the energy density drops below that of the no-oscillation case. This decrease can be understood as follows. Energy is removed from the oscillating active neutrino ν α population through conversion to sterile neutrinos beginning at the low end of the momentum distribution. Some of these low-momentum active states are refilled from higher-momentum states through ν α scattering with electrons, positrons, and the spectator neutrinos. This effectively reduces the total energy contained in the combined active and sterile neutrino population, where the surplus energy has been absorbed into the background of e + , e − , and The second scenario (blue dashed and dotted lines) consists of only scattering processes amongst the active oscillating neutrinos themselves. Number density conservation is again satisfied to a good degree. Energy density is likewise approximately conserved; this is expected, as the isolated ν α population has no contact with the background plasma with which to exchange energy. Observe that the degree of non-conservation is a larger here than in the first scenario. This is because the evaluation of the R α,α collision integrals (A.15) and (A.16)  for the ν α ν α processes are more sensitive to numerical inaccuracies owing to their nonlinear dependence on the distribution function fν α (see section 2.3). Increasing the number of momentum bins from the canonical 100 to 150 bins (thick black dashed and dotted lines) improves the situation, although it is clear that we still do not have perfect detailed balance. Nonetheless, the induced error is small enough that we can conclude it is under control.

Comparison of approximation schemes
The different approximation schemes have different strengths and weaknesses as already discussed in section 2.2. On the one hand, the equilibrium approximation has the advantage that the assumptions involved have been analysed quite extensively [15]. It also allows the scattering processes to push the distributions towards the equilibrium form as they should, albeit at the cost of sacrificing detailed balance. The CC approximation on the other hand enforces detailed balance as it prevents the same scattering processes from repopulating the active neutrino states. However, in doing so, the approximation also neglects possible refilling due to redistribution between different momentum states. With these considerations in mind, we expect the equilibrium approximation to overestimate the neutrino number and energy densities, and the CC approximation to underestimate them.
These trends are reinforced and possibly enhanced by the behaviour of the damping term in the two approximation schemes. By neglecting Pauli blocking, the equilibrium approximation overestimates the damping, while the CC approximation underestimates it through a missing positive Pauli blocking term discussed in section 2.2.2. The size of the damping term has direct consequences for ∆N eff as the production of the sterile neutrinos in our case is non-resonant, where the effective production rate, is directly proportional to the damping term Γ collision ∼ D and the in-matter mixing angle θ m [24]. Thus, the larger the damping, the higher the resulting ∆N eff , and vice versa. As shown in figure 3, the over-and underestimation of ∆n ν and ∆N eff in the equilibrium and CC approximations, respectively, are exactly what transpires in our numerical solutions of the QKEs for the benchmark mixing parameter values. In contrast, the A/S approximation introduced in this work takes the best of both worlds, and, as is evident in figure 3, comes much closer to the full collision treatment than both the equilibrium and the CC approximations. The A/S approximation is also expected to overestimate somewhat the damping relative to the full treatment, since oscillations generally push the distribution functions below the equilibrium values. This is however a very small effect, and the fact that the resulting ∆N eff and ∆n ν values tend towards different sides of the full solutions suggests that the A/S damping term is not the cause of any major systematic bias.

Electron mass and Pauli blocking
Next we test the consequences of ignoring Pauli blocking, and of assuming m e = 0 in the full collision term. As shown in figure 3, the latter assumption makes no practical difference to ∆n ν and ∆N eff , since the conversion to sterile neutrinos takes place largely before electrons become nonrelativistic. Ignoring Pauli blocking, however, induces a ∼ 0.02 absolute error as T drops below ∼ 5 MeV, smaller than the naive expectation of ∼ 10% estimated from the Pauli blocking terms in the relevant momentum range [15].
Note that there is a small subtlety when ignoring Pauli blocking: detailed balance can be satisfied for Fermi-Dirac distributions only when the Pauli blocking is taken into account. Otherwise, ignoring Pauli blocking generally demands that we replace Fermi-Dirac with Maxwell-Boltzmann statistics in order to fulfil detailed balance. We would however like to continue using Fermi-Dirac distributions, and therefore choose to enforce detailed balance in a different way. For all processes of the form a(p) + ν α (k) ⇄ b(p ′ )c(k ′ ), we assume that This is also the assumption behind the equilibrium approximation, and is exact if f (E p ′ ), f (E k ′ ) and f (E p ) take on the Maxwell-Boltzmann equilibrium form. Not surprisingly we see in figure 3 that the equilibrium approximation solution follows the no Pauli blocking full solution quite well; the difference between them is due mainly to rounding errors in the numerical value of C α in the equilibrium approximation immediately below equation (2.9).

Impact on distribution functions
It is also interesting to compare the active and sterile neutrino distribution functions in the different approximations. The distribution function of the electron neutrino affects directly the weak reaction rates in BBN [6], while the division of ∆N eff between ν α and ν s have consequences for large-scale structure formation if δm 2 is significantly larger than the solar and atmospheric mass splittings [11,12]. Figure 4 shows this comparison for the benchmark mixing parameters.
The first observation is that the full collision treatment gives rise to f > f 0 for the high-momentum active neutrinos at T = 10 MeV. This is not a numerical artefact, but   originates in both the ν α ν α -scattering and the annihilation processes in the presence of a step-like feature at low momenta, such as that produced by oscillations here. Consider first  ν α ν α -scattering. Removing neutrinos from the low momentum modes causes the energy per neutrino to increase. This in turn triggers the number and energy conserving scatterings to push the distribution towards an equilibrium with a higher effective temperature, which automatically leads to f > f 0 . Annihilation processes on the other hand enhance f via a slightly different mechanism. For processes involving only states above the step, the equilibrium is maintained. For processes such as ν α (k high )ν α (p low ) ⇄ a(k ′ )ā(p ′ ), however, the reaction will be pushed to the left because there is a deficit ofν α (p low ) relative to a(k ′ ) and a(p ′ ). Thus, annihilation can also overproduce ν α at high momenta, leading to f > f 0 .
The A/S approximation mimics both the annihilation and the scattering effects to an extent, and is the only approximation scheme that manages to reproduce f > f 0 at T = 10 MeV, albeit only at very high momenta. The annihilation effect is captured by the n να factor in equation (2.16), while the separate treatment of ν α ν α accounts for the scattering effect. The latter constitutes the main role of the C α,ν term in the A/S repopulation approximation (2.16); in contrast, the C α,s term with the pseudo-chemical potential ξ as the sole degree of freedom acts in the wrong direction: it becomes negative when the active sector is depleted, giving rise to a lower value of f for any choice of momentum.  Notwithstanding its success at capturing some degree of the f > f 0 phenomenon, the approximate treatment of scattering remains the main source of discrepancy between the full solution and the A/S approximation at T = 10 MeV. This is evident in figure 5 where we compare the distribution functions obtained assuming (i) only scattering, and (ii) only annihilation. The scattering-only result also demonstrates that these processes do not replenish the active neutrino population, yielding f /f 0 ∼ 0.85 for most momenta, while the annihilations are very effective at bringing f /f 0 to 1.
Apart from these effects, we find that the equilibrium approximation gives too large a final distribution for both the active and the sterile neutrinos, as expected from the oversized R eq and D eq terms. The CC approximation, on the other hand, comes surprisingly close to the real active neutrino distribution, but is short by ∼ 5% for the sterile states. As repopulation does not directly affect oscillations into sterile states, we conclude that this offset must originate in the undersized damping term D CC , which as discussed in section 3.2 affects directly the effective sterile neutrino production rate (3.2).

Dependence on mixing parameters
So far we have tested the various approximation schemes using the benchmark mixing parameter values δm 2 = 0.1 eV 2 and sin 2 2θ = 0.025. In reality, however, the errors caused by these approximations are also sensitive to the parameter values. Figure 6 shows ∆N eff at T = 0.1 MeV as a function of δm 2 and sin 2 2θ computed using the full collision term. Our results are generally quite similar to previous calculations using the equilibrium approximation [4,21]. Note however that for large mixing angles and small mass squared differences, the contours deviate slightly from the straight lines found in [4,21]. To highlight this deviation, we show in figure 7 the differences in ∆N eff incurred respectively by the equilibrium, CC, and A/S approximations relative to the full solution in the (δm 2 , sin 2 2θ)-parameter space.
As expected, the deviations always occur, independently of the approximation scheme, in a diagonal band corresponding to the transition region from ∆N eff = 0 to ∆N eff = 1. Beyond this common feature, however, the different approximations incur the largest deviations at different parameter values.
The equilibrium approximation follows the result of the full calculation quite well for δm 2 values above 0.01eV 2 , but overestimates ∆N eff by more than 0.1 at δm 2 < 0.001 eV 2 . We can understand this deviation by looking at the conversion temperature. The temperature of maximal conversion is proportional to (δm 2 ) 1/6 [21], so that low values of δm 2 generally correspond to low conversion temperatures. If the conversion temperature is sufficiently high (T ≫ 1 MeV), repopulation is rapid and ∆N eff is limited only by how fast ν s can be produced through oscillations and collisions [24]. The effective production rate is given by equation (3.2), and production ceases as soon as the collision rate becomes too low. If most of the conversion occurs at low temperatures, however, collisions become too inefficient to sustain the population of the active sector and consequently for equation (3.2)-which assumes instantaneous repopulation-to hold, thereby causing the real ∆N eff contours to deviate from straight lines in relation to δm 2 and sin 2 2θ in figure 6. The equilibrium approximation errs in its overestimation of the repopulation rate, yielding almost straight ∆N eff contours even in the low δm 2 , high sin 2 2θ region.
In contrast, the top right panel of figure 7 shows that the CC approximation generally underestimates ∆N eff , but works somewhat better at δm 2 0.01 eV 2 . For high δm 2 values the underestimation is due mainly to the undersized D CC damping term which diminishes the sterile neutrino production rate as already discussed in section 3.2. For low δm 2 values, however, the agreement becomes better (a deviation of 0.02 at δm 2 = 10 −4 , sin 2 2θ = 10 −0.5 ) because the deficiency of sterile neutrinos is compensated by an overproduction of active  neutrinos (similar to the overpopulation of the active sector discussed above in relation to the equilibrium approximation). This overproduction comes about as follows. Although technically the CC repopulation integral incorporates only annihilation processes, numerically the constant 2(γ e a ) 2 /3.15 = 0.317 is rather larger than its annihilation counterpart in the A/S approximation, C e,a = 0.177. Thus, the CC approximation does inadvertently contain some degree of repopulation due to scattering, which drives up the active neutrino repopulation rate relative to the true rate.
Finally, the bottom panel of figure 7 shows the A/S approximation. Here, the agreement is generally much better than either the equilibrium or the CC approximation, although there remains a small discrepancy of ∼ 0.01 in the low δm 2 , high sin 2 2θ corner. This is not surprising, as the region is characterised by large spectral distortions from conversion at low temperatures, and is thus most sensitive to how exactly we handle repopulation in the active sector. It is nonetheless remarkable that for most of the parameter region, the A/S approximation is able to reproduce the full results at a precision of ∼ 0.001 through a fairly simplistic description of the collision term. This deviation is in fact comparable to the numerical error expected of the full treatment due to our choice of momentum resolution (see section 3.1). Therefore, figure 7 not only validates the A/S approximation, but through comparison with a physically intuitive modelling of collisions, also confirms that the full collision term has been implemented correctly

Conclusions
We have calculated in this work the full collision term for active-sterile neutrino oscillations in the early universe, and for the specific case of (ν e , ν s )-oscillations implemented it in the computation of ∆N eff from sterile neutrino thermalisation. In particular we have included for the first time a nonzero electron mass and Pauli blocking in the collision integrals. The former turns out to have a negligible impact on the final results, while ignoring the latter gives rise to noticeable discrepancies. Nonetheless, our full treatment confirms previous analyses based on approximations of the collision terms that a 1 eV sterile neutrino coupled with a mixing angle sin 2 2θ ∼ 0.1 produces ∆N eff ≈ 1; the discrepancies arising from approximations are at the level that affects only precision calculations.
We then proceeded to perform a systematic comparison of the full collision treatment with two approximation schemes found in the literature. We find that the commonly used "equilibrium approximation" [15] reproduces ∆N eff at better than 0.04 for δm 2 > 0.01 eV 2 , but incurs large errors (> 0.1) for very small mass squared differences due to the unphysical repopulation of the active sector by elastic scattering. The "CC approximation" of Chu and Cirelli [7] is discrepant up to 0.04 for δm 2 > 0.0001 eV 2 , but improves for low δm 2 values because of a fortuitous cancellation between an underestimated damping term and an overestimation of repopulation from annihilation processes.
Recognising that the equilibrium and CC approximations neglect different physical effects, we have devised a new approximation scheme, the A/S approximation, in which scattering and annihilation contributions to repopulation are treated separately, and the damping term includes Pauli blocking. As the scheme is better able to capture the physics of repopulation and damping, it also pushes the error in ∆N eff down to 0.001 for most of the parameter region, although minor deviations (∼ 0.01) remain in the δm 2 < 0.001 eV 2 region, where the sterile neutrino conversion temperature is lowest and hence spectral distortions from oscillations are expected to have the largest effect.
The connection between low δm 2 deviations and low temperature spectral distortions also has implications for an inverted mass spectrum, i.e., where the active state is heavier than the sterile state, as well as for the various mechanisms designed to reconcile eV-mass sterile neutrinos with cosmology. In the case of an inverted mass spectrum, sterile neutrino thermalisation proceeds via a resonance, which, depending on the adiabaticity of the resonance and hence the mixing angle, can cause more disturbance to the active neutrino spectrum than in the non-resonant case. Likewise, mechanisms that block the production of eV-mass sterile neutrinos typically work by delaying the thermalisation to low temperatures, which by construction also makes them more sensitive to the approximations employed for the collision terms.
At the current level of observational precision-σ(∆N eff ) ≈ 0.2 from Planck [2]-our analysis shows that the CC and the A/S approximations can be reliably applied to most active-sterile oscillation scenarios, whereas the equilibrium approximation appears to be approaching its boundary of validity if the sterile neutrino conversion temperature is too low. In the future, large-volume galaxy surveys are expected to improve the sensitivity to ∆N eff to ∼ 0.03 [25]. Thus, should hints for a 1 eV sterile neutrino persist in the laboratory, a collision treatment more precise than either the equilibrium or the CC approximation can offer would become necessary. Short of evaluating the full collision terms, which is numerically costly or possibly even infeasible in some cases, the A/S approximation developed in this work appears to be a most convenient alternative.

Acknowledgments
We acknowledge use of computational resources provided by the Danish e-Infrastructure Cooperation. Y 3 W thanks the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work.

A Derivation of the full collision terms
A.1 Neutrino-electron scattering in the s-channel Consider first the case of ν α + e − ⇄ ν α + e − . The repopulation term due to this process takes on the form (A.1) where the matrix element is [26,27] with A e = 2 sin 2 θ W + 1, A µ,τ = 2 sin 2 θ W − 1, B e = 2 sin 2 θ W , and B µ,τ = 2 sin 2 θ W . Note the change of notation here in the appendix: p and p now denote respectively the 4-and 3-momentum, while in the main text we use p to indicate the magnitude of the 3-momentum. The matrix element has three different dependences on the momentum, and therefore requires three different parameterisations to simplify the integral. We label the first dependence (k · p)(k ′ · p ′ ) the s-channel, the second (k · p ′ )(k ′ · p) the u-channel, and the last dependence (k · k ′ ) the t-channel, where for each channel it is convenient to define a 3-momentum variable, given respectively by which replaces one of the integration variables of equation (A.1). Then, to evaulate the integral for each channel we simply follow the technique described by Hahn-Woernle, Plümacher and Wong in [19].
With this choice, we find the quantities and consequently for the first term of the matrix element (A.2).
Since the matrix element (A.3) now depends only on energies and the magnitude of q, we can use the Dirac delta functions in the integral (A.1) to integrate out the directional dependencies. Evaluating first the p ′ -integral as [19] to further simplify equation (A.4) to [19] (A.5) In a similar way, we rewrite where in the second line we have also changed the integration variable from p to q. Then, applying equation (A.5) and (A.6) to the integral (A.1), we find, after performing the trivial angular integrations and averaging over the direction of the incoming neutrinos ( d cos η/2), the s-channel contribution R α,s,e,s = 1 where F denotes all the distribution functions, and the arguments of the two Dirac delta functions can be read off equations (A.5) and (A.6) respectively.
The integrals over cos θ and cos η involve only Dirac delta functions. Taking the cos ηintegral as an example and noting that integration limits can be equivalently expressed as step functions, we find The two step functions can be reinterpreted as limits on |q|, and together they confine |q| to the interval Integrating over d cos θ gives the same formal result save for the replacements k → k ′ and p → p ′ . Thus the combined integration limits on |q| can be written equivalently as max ||k| − |p||, ||k ′ | − |p ′ || ≤ |q| ≤ min |k| + |p|, |k ′ | + |p ′ | , (A. 8) and we note that |p ′ | is determined from |p|, |k|, |k ′ | by imposing energy conservation. Applying the limits (A.8) to the integral (A.7) then yields R α,s,e,s = 1 as our reduced repopulation integral from the s-channel. The u-and t-channel integral reduction proceeds in a similar manner, using v and w respectively as an integration variable.

A.2 The massive case
The reduced integral (A.9) is but one of three contributions to the repopulation of ν α arising from neutrino scattering with electrons. The full collision term, including scattering and annihilation processes with e ± , ν β , has in total 14 such terms to be evaluated (see table 1). Fortunately, however, these 14 terms can all be recast into one of the standard s-, t-, and u-forms, and thus can be handled in ways similar to that discussed above.
As the reduction procedure concerns only kinematics and does not involve the actual matrix element besides the initial classification of the momentum dependence into s-, u-or Table 1. Matrix elements for all relevant reactions involving ν α , with x W = sin 2 θ W = 0.23864 [28], S is a symmetrisation factor of 1/2 for each pair of indistinguishable particles in the initial and the final state, and |M | 2 has been summed but not averaged over initial and final spins. For the process with two ν α s in the initial state, we have further multiplied the matrix element by 2 to account for the fact that ν α (k)ν α (p) → . . . and ν α (p)ν α (k) → . . . constitute two identical processes. Where there is a choice of ±, the plus signs are for α = e, and the minus signs for α = µ, τ . The corresponding matrix elements forν α can be obtained by the exchange (k · p)(k ′ · p ′ ) ↔ (k · p ′ )(k ′ · p) for the elastic scattering processes.
t-forms, we shall keep the calculation as general as possible and allow for the possibility that all initial and final states are massive. Then, the number of independent reductions to be performed is only three, which yield: (A.12) After inserting the matrix element S|M x | 2 and momentum distributions F , these reduced integrals are valid for any 2 → 2 process.

A.3 The full collision terms
The matrix elements for all elastic and inelastic processes involving ν α at temperatures T m µ are summarised in table 1. These have been computed at various times by several different groups [26,27,29], but can be easily obtained from first principles in the four-fermion limit. Using these matrix elements and the expressions (A.10), (A.11) and (A.12), we can now determine the contribution of each process to the repopulation integral. The results are as follows.
The integrals over |q|, |v| and |w| can be evaluated analytically, and it turns out that they fall into two different functional forms: The remaining two integrals over E p and E k ′ must be performed numerically, although as discussed in section 2.3 judicious assumptions about certain distribution functions in the integrand make it possible to pre-evaluate some of the integrals only once, rather than evaluating them in real time simultaneously with the numerical solution of the QKEs.

B Repopulation and damping coefficients in the A/S approximation
We summarise here the full expressions for the dimensionless repopulation and damping coefficients in the A/S approximation discussed in section 2.2. The quantity A ≡ 2π/ dΠ k kf 0 (k) is a normalisation factor.