lopscience

¡opscience.iop.org

Home Search Collections Journals About Contact us My IOPscience

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

This content has been downloaded from IOPscience. Please scroll down to see the full text. View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 59.98.174.246

This content was downloaded on 24/06/2016 at 18:50 Please note that terms and conditions apply.

/ournal of Cosmology and Astroparticle Physics

An IOP and SISSA journal

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

Steen Hannestad,ab Rasmus Sloth Hansen,a,c Thomas Tramd and Yvonne Y.Y. Wongc

"Department of Physics and Astronomy, Aarhus University, 8000 Aarhus C, Denmark b Aarhus Institute of Advanced Studies, Aarhus University, 8000 Aarhus C, Denmark cSchool of Physics, The University of New South Wales,

Sydney NSW 2052, Australia ^Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, United Kingdom

E-mail: sth@phys.au.dk, rshansen@phys.au.dk, thomas.tram@port.ac.uk, yvonne.y.wong@unsw.edu.au

Received June 22, 2015 Accepted July 6, 2015 Published August 11, 2015

Abstract. Sterile neutrinos are thermalised in the early Universe via oscillations with the active neutrinos for certain mixing parameters. The most detailed calculation of this thermal-isation 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.

Keywords: cosmological neutrinos, neutrino theory, physics of the early universe ArXiv ePrint: 1506.05266

Article funded by SCOAP3. Content from this work may be used

A f th rd^K termhof tte1Creative CTmoits A"1;"^®.0 LTT) doi:10.1088/1475-7516/2015/08/019

Any further distribution of this work must maintain attribution to the author(s) / ill and the title of the work, journal citation and DOI.

Contents

1 Introduction

2 Quantum kinetic equations

2.1 Repopulation and damping

2.2 Approximation schemes

2.2.1 Repopulation

2.2.2 Damping

2.3 Numerical implementation 2.3.1 Approximation schemes

2.4 Full collision term

3 Numerical results

3.1 Numerical convergence

3.2 Comparison of approximation schemes

3.3 Electron mass and Pauli blocking

3.4 Impact on distribution functions

3.5 Dependence on mixing parameters

4 Conclusions

14 14 17

A Derivation of the full collision terms 21

A.1 Neutrino-electron scattering in the s-channel 21

A.2 The massive case 23

A.3 The full collision terms 25

B Repopulation and damping coefficients in the A/S approximation 27

1 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 Neff = 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-7] or

interactions among the sterile neutrinos [8-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 active-sterile 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 G2F 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 Neff 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 Neff 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 5m2 and vacuum mixing angle d, 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 5m2 ~ 0.1 eV2: 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 5m2 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 = h = kB = 1, and express all dimensionful quantities in powers of eV.

2 Quantum kinetic equations

We consider oscillations between an active neutrino flavour va, where a = e, ^ or t, and a sterile flavour vs in an ensemble of neutrinos in the early universe. The density matrices p(k) encode the flavour content and coherence of the ensemble, and are conveniently expressed in terms of polarisation vectors (P0(k), P(k)), i.e.,

p(k)=(Pit) £<*))=1 /o (k)[Po (k)1+P(k) ■ a]- >

where f0(k) denotes a reference phase space density, which we take here to be the relativistic Fermi-Dirac distribution with vanishing chemical potential,1 and a = [ax,ay,az} are the Pauli matrices. In this convention the active and sterile neutrino phase space distributions are given, respectively, by

fva (k) = 1[P0(k)+ Pz (k)],

21 (2.1) fvs (k) = 2[Po(k) - Pz(k)],

and the QKEs that govern their evolution can be written as [13, 15]

P(k) = V(k) x P(k) + ^-^Z - D(k)PT(k),

R (k) /0(k) (2^) ^

po(k) = fi ^

Here, V(k) = VX(k)X + Vz(k)z contains the vacuum oscillation term as well as the matter potential from forward scattering, and the components are given respectively by

V»(k) = -2F sin20,

Vz<k> = -IFcos2fl - 45^1kT4("'" + (M>

where gf is the Fermi constant, mz the mass of the Z boson, m the number density normalised to the equilibrium value, gM/T = 1, and ge = 1 + 4 sec2 /(nVe + n^e) with the Weinberg angle . Finally, Ra(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 nVa = .

1 Henceforth we shall reserve f0 to denote exclusively the relativistic Fermi-Dirac distribution with zero

chemical potential, and use feq to indicate an equilibrium distribution of a nonspecific form.

2.1 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

The prefactor N = 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)u(pi) = pi + mi and v(pi)v(pi) = pi - mi.

Ra(k) = J dnk'dnp'dnp se(kp|k'p') (2.4)

x E V2[v«(fc),v/«(p)li(k'), i(p')] [/i(Efc')f(Ep)(1 - fVa (k))(1 - f9a (p))

- fva (/(P)(1 - fi(Ek' ))(1 - f(Ep))] + ^ V2[va(k), j(p)K(k'), j(p')] [fva(k)fj(Ep)(1 - fva(k))(1 - fj(Ep))

fVa (k)fj (Ep)(1 - fva (k' ))(1 - fj (Ep' ))] ,

D(k) = n J dnk'dnp'dnp SE(kp|k'p') (2.5)

x V2[Va(k),V/«(p)|i(k'), i(p')] [fi(Ek' )fi(Ep' )(1 - fVa (p))

+ fVa (P)(1 - fi(Ep' ))(1 - fi(Ek' ))] O

+ E V2[va(k),j (p)K(k' ),j(p')] [fva (k' )fj (Ep' )(1 - fj (Ep))

+ fj (Ep)(1 - fj (Ep' ))(1 - fva (k'))] , where we have used the shorthand dnp = d3p/(2n)3, SE(kp|k'p') = S(1)(Ek+Ep-Ek' -Ep') is

the 1D Dirac delta function, the summation index i runs over all spectator neutrino flavours (i.e., vg where P = a) and the electron, while j runs in addition over all their antiparticles as well as the oscillating neutrino and antineutrino. The terms V2 are

V 2[a(p),6(k)|c(p' ),d(k')] = (2n)3S(3)(k + p - k' - p' )N2N62N2N,2S |M |2(a(p), b(k)|c(p' ),d(k')),

where N = i/1/2Ei,2 Ej 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)6(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 vas are present in the initial state, then S|M|2 must additionally be multiplied by 2 to account for the fact that va(k)va(p) ^ ... and va(p)va(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 - fj) 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 by McKellar and Thomson [13],

D(k) = n J dnfc'dUp'dUp se(kp|k'p') (2.7)

x ^ V2K(k),z7«(p)|i(k'),¿(p'f (p) + ^ V2[v«(k), j(p)K(k'),j(p')]fj(Ep), i j

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 p, Pauli blocking enters the equation of motion as a multiplicative matrix factor of the form Sj — Pj (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, Pliimacher 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.

Req(k) = r fo — 2 (Po + Pz) J , (2.8)

Deq(k) = |r, (2.9)

where r = xT5, with x = k/T, Ce « 1.27, and C^>r « 0.92 [20-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]

g2 t5 k

pcoii(k) = —^-^y ({G2, p — p0} —2Gs(p—p0)Gs + {Gi p—p0} +2Ga(p—p0)Ga) ,

(2.10)

2.2 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 re-populated, the repopulation and damping terms evaluate in what we shall call the equilibrium approximation [4, 15, 20] to

where Gs,a = diag (y?«, 0), with (y!)2 = 3.06 and (Ys'r)2 = 2.22 for scattering processes, and (Y«)2 = 0.50 and (Ya'T)2 = 0.28 for annihilations [7], (k) œ 3.15 T is the average momentum, and p0 = /0 1. The CC approximation was originally applied in [7] to the momentum-averaged 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 find

P (k) 1G2 T5 k /4(Ya)2(p«« - pL) [(y?)2 + (Yaa)2] P«A (2 n)

pcoll(k) = -2GFT (k) ^[(Y?2)2 + (Ya)2] Psa 0 )> (2.11)

which can be further recast as expressions for repopulation and damping,

i2 ^5 X 2

Rcc(k) = 2GFT5^(Yaa)^/o - J(Po + Pz)) , (2.12)

Dcc(k) = 1 GF T5 -X- [(Y?2)2 + (Yaa)2 ] , (2.13)

2 F 3.15

compatible with equation (2.2).

2.2.1 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 v^a« annihilation into a and a is an integral over (fafa — fVa faa), where we have ignored Pauli blocking factors. In the equilibrium approximation, it is assumed that f « f0 for all particles i = a, a, aa but va. For a and a this is a good assumption. For aa, however, we know that it overestimates faa if naa < 1. Thus, to compensate for the depletion of aa, we adopt faa « f0 (remembering that naa = ) in the annihilation part of the A/S approximation, while keeping f = f0 for i = a, 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 va population and an external bath.

1. Nonzero energy and momentum exchange. Scattering processes involving electrons, positrons and spectator neutrinos , 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 /Va relaxes replaced by

/eq, scat = gX—g + 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 /eq, scat and /Va must be equal (because the collision rate is proportional to momentum).

2. Vanishing energy and momentum exchange. Scattering amongst and Va 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 Tv

different from the photon temperature T, i.e.,

/eq, = gfc/Tva -Mva /Tva +1 ' (2.15)

Here, and TVa are determined from the combined condition that the scattering contributions to and must both be zero, where is the energy density normalised to the equilibrium value. The former constraint is equivalent to requiring equality between the third kinetic moments of /eq, Va and /Va, 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 Ra/S (k) = GFxT5 | Ca>a (/o - (P0 + Pz)) + Ca,s (Zeq, scat - f (P0 + Pz))

( /eq,Va - y (P0 + Pz)

3A pseudo-chemical potential appears with the same sign for both particles and anti-particles, whereas the

normal chemical potential has opposite signs for particles and anti-particles.

(2.16)

where Ce,a = 0.180, Ce,s = 0.718, Ce,v = 0.407, C^/t>0 = 0.102, C^/T>s = 0.407, and C^/T>v = 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 and 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).

2.2.2 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, Deq(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 Dcc(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 ~ f0 for i = va,aa and fVa = faa « nVaf0. Due to Pauli blocking, these approximations lead to a damping term that is quadratic in nVa:

da/s = 2GFxT5 «Ca,2 + nvaCa,i + C«^) , (2.17)

where Ce,2 = —0.020, Ce,i = 0.569, Ce,0 = 0.692, C„/r>2 = —0.020, C„/r>1 = 0.499, and C„/r,0 = 0.392, the negative Ca,2 values reflecting the origin of the n^ 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 Ca>1 is predominantly from annihilation, while Ca,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

C«,fit = a + 6e-cx--—, (2.18)

2.3 Numerical implementation

We employ a modified version of the public code lasagna5 [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 P — P = 0 so as to enforce a zero lepton asymmetry.

2.3.1 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 p of the normal scattering term (2.14), as well as the temperature TVa and chemical potential of the vava-scattering term (2.15). The chemical potential of the normal scattering term satisfies the equation

ek/T-„/T + 1 = 2/ dnfc k f0(P0 + P),

4The collision terms of the CC approximation have been presented in [7] in their integrated form, without details of how exactly they have been computed. However, reverse engineering suggests that they arise from the

integral -K f dUk duk' dnp' dnv Se (kp\k'p')J2 i V 2[va(k),Va(p)\i(k'),i(p' )]/eq(p)/eq(k)(l - /eq(p'))(l - feq(k' )) for the annihilations, and n/dHkdHk'dnp'dnp Se(kp\k'p')J2j V2[va(k),j(p)\va(k'),j(p')]feq(p)feq(k)(1 -feq(p'))(1 — feq(k')) for the scatterings, both normalised by f dnk/o(k).

5Available at https://github.com/ThomasTram/lasagna_public.

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 ANej and in the active spectrum which are often the most interesting quantities. For this reason we do not incorporate Ca>flt in the A/S approximation.

which can be solved numerically. The parameter space of the vava-scattering term, on the other hand, is two-dimensional and subject to the conditions

dnk k ek/Tva-.1 /Tva +1 = 1 / dnk k + ), (2-19)

dnk ^ ek/Tva/Tva +1 = ^ ^ ^^ + PZ)- (2-20)

In order to solve for and TVa, we first reduce the problem to one dimension by constructing the ratio

(/dEfc k2/o(Po + Pz))4/5 = (/dnfc k2/(efc/Tvq-.yq/Tva + 1))4/5

/ dnfc k/o(Po + Pz) / dnfc k/(ek/TVa-.va /Tva +1) ' Q

r / Mvq M4/5 Ivj

o 1/5 -24Li5 —e Tva

= (2n2)--V J\ , (2.21)

—6Li4 e Tv« J

using equation (2.19) and (2.20). Here Lis(z) denotes the polylogarithm. Since the r.h.s. depends on /TVa but not directly on TVa, equation (2.21) can be solved immediately for /TVa using a one-dimensional root-finding algorithm. Once we have established /TVa, equation (2.20) can be evaluated explicitly for the

( \V5 / dx x4/o(Po + Pz )/2

y —24Li5 (-eTva j y

T, (2.22)

where T is the photon temperature, and again we have x — k/T. Finally, we take the and TVa 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 TVa 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.

2.4 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:

Ra — Ra,e + R

where the second subscript on the r.h.s. labels the contributing processes. We use equilibrium distributions /eq for the electrons and positrons in the processes

v«(k)e±(p) ^ v«(k')e±(p'), (2.23)

v«(kK(p) ^ e-(k')e+(p'), (2.24)

temperature TVa. We find: o

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 Ra can be expressed as

«„(*) =(i - / (*))(/ ^m/(*') + /^Wi -/ (P»)

- /(k) (y dEfc'Ra,e,s,2(1 - /(k')) +J dEp/(p)R«;e)a;^ ,

Ra ,e ,s , 1 (k, k') = J dEp(Ra , s>e- + R« , s, e+ )/eq(Ep)(1 - /eq(Ep)),

Ra, e,s,1 (kj k')

Ra ,e,a,1 (k,p)

Ra, e,s,2(k, k')

Ra ,e,a,2(k, p)

; , a , ej eqV-1-7 feVJ eqV-1-7p'

^pl-^-a^e- -1 <<a,s,e+/ ■ZeqV-'^pA-L JeqV-^p'

Ra,e,a,2(k,p) = J dEfc/Ra,a,e(1 - /eq(Efcz))(1 - /eq(Ep')),

(2.25)

(2.26)

where the integration kernels Ra,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, Ra,e,s,i and Ra,e,a,i of equation (2.26) scale trivially with temperature as a T4. 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 Ra,e,s,i and Ra,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 o in real time using a 3D spline when solving the QKEs.

We use the same equilibrium approximation for the spectator neutrinos vg and antineutrinos z/g in the processes

Va(k)^(p) ^ Va(k')(-g (p'), Va(k)z/a(p) ^ Vg(k')z/g(p'),

where a = P, and(-g are assumed to be non-oscillating. Thus Ra,g and the associated kernels Ra,g,s,i and Ra,g)a,i are, save for the e ^ P relabelling, formally given by equations (2.25) and (2.26) respectively, and the corresponding expressions for Ra,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 Ra,g,S)i and Ra,g)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.,

Va(k)(-,a (p) ^ Va(k')(-a (p'),

we evaluate the full double integrals (A.15) and (A.16) without approximation, since the momentum distributions of va and are a priori unknown quantities.

3 Numerical results

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 /(k) and 1 — /(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 Da can be expressed as

Da,e(k) — ( [ dEfc'Ra,e,s,l/(k') + / d£pR«;e,a,1 (1 — /(p)))

7 7 . (2.27)

J d£fc' Ra>e>s>2(1 — / (k')) + J dEp/(p)R«,e,aJ ,

where Ra,e,s,i and Ra,e,a,i are exactly the pre-evaluated quantities of equation (2.26).

The main quantity we wish to study is the change in the energy density of neutrinos due to a sterile neutrino population,

ANeff — Neff ,a + Neff, s — 1 — ^^f dk/ (k) + (k))k3 — 1. (3.1)

This one quantity affects directly the Hubble expansion rate, making it possible to constrain ANeff using various cosmological observations.

3.1 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 and effective energy density Neff using the benchmark mixing parameter values ¿m2 — 0.1 eV2 and sin2 20 — 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.

£ 0.006

dm2 = 0.1eV2, sin220 = 0.025

---- 150 momentum bins --100 momentum bins ■ ■ ■ 80 momentum bins _ --- 40 momentum bins

.........................

5 0.004

T [MeV]

Figure 1. Convergence test using the full collision term. We compare the difference in the neutrino number and energy densities between using 200 momentum bins and 40 (blue long dashed), 80 (red dot-dashed), 100 (green dashed) and 150 bins (cyan dotted).

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 va 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 va 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 va 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 Ra,a collision integrals (A.15) and (A.16)

Sm2 = 0.1eV2, sin220 = 0.025

ANeff, no annihilation

AWeff, only va,va processes

ANefl, only va,va processes, 150 bins

Anv, no annihilation

Anv, only va, va processes

Anv, only va, va processes, 150 bins

15 T [MeV]

Figure 2. Different contributions to the collision terms and their effects on the neutrino number and energy densities relative to the no-oscillations case. The red solid and dot-dash lines denote a scenario with only scattering and no annihilation. The blue dashed and dotted lines represent scattering only amongst the oscillating active neutrinos va; the thick black dashed and dotted lines denote the same scenario, but now computed using 150 momentum bins instead of the canonical 100.

for the processes are more sensitive to numerical inaccuracies owing to their nonlinear dependence on the distribution function (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.

3.2 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 ANeff as the production of the sterile neutrinos in our case

2 O M U1

is non-resonant, where the effective production rate,

r — 2sin2(20m)rcollision, (3.2)

is directly proportional to the damping term rcollision ~ D and the in-matter mixing angle [24]. Thus, the larger the damping, the higher the resulting ANeff, and vice versa.

As shown in figure 3, the over- and underestimation of Anv and ANeff 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. It is not obvious what to expect for the A/S approximation, and the fact that the resulting ANeff and Anv values tend towards different sides of the full solutions suggests that there is no major systematic bias.

3.3 Electron mass and Pauli blocking

/ (Ep )/(Ek ) — / (Ep)/(k) — /eq(Ep)(/eq(k) — / (k)). (3.3)

This is also the assumption behind the equilibrium approximation, and is exact if /(Ep'), /(Ek) and /(Ep) 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 Ca in the equilibrium approximation immediately below equation (2.9).

3.4 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 ANej between and have consequences for large-scale structure formation if ¿m2 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 / > /0 for the high-momentum active neutrinos at T — 10 MeV. This is not a numerical artefact, but originates in both the vava-scattering and the annihilation processes in the presence of a

Next we test the consequences of ignoring Pauli blocking, and of assuming me — 0 in the full collision term. As shown in figure 3, the latter assumption makes no practical difference to Anv and ANeff, 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) + va(k) ^ b(p')c(k'), we assume that

0.8 0.6

¡^ 0.4 <

0.0 0.04

% 0.02

^ 0.00 I

' -0.02

Full collision term Assuming me = 0 Ignoring Pauli blocking Equlibrium approx. Chu & Cirelli approx. A/S approx.

0.8 0.6

I °-4 0.2

0 5 10 15 20 25 30

N * * N -

0 5 10 15 20 25 30

T [MeV]

- Full collision term Assuming me = 0 _ Ignoring Pauli blocking Equlibrium approx. -Chu & Cirelli approx. A/S approx. -

, -----r

e -0.02

15 T [MeV]

Figure 3. Comparison of the full treatment (black solid lines) with the equilibrium approximation (red dot-short dash), the CC approximation (blue long dashed), and the A/S approximation (purple dot-long dash) introduced in this work. Setting me = 0 (cyan dotted) has almost no effect, while ignoring Pauli blocking in the full expression (green dashed) gives almost the same result as the equilibrium approximation.

2 O M U1

0.991 1.003

0.991 1.003

Relative active distribution

10.0MeV

T = 1.0MeV

0.8 0.6 0.4 0.2

0.0 1.0 0.8 0.6 0.4 0.2 0.0 0.82 0.78 0.74 0.70 0.66

Relative sterile distribution

Full collision term Equlibrium approx. Chu & Cirelli approx. A/S approx.

x = k/T

Figure 4. The active and sterile neutrino momentum distributions computed from the full treatment (black solid lines), the equilibrium approximation (red fit-short dash), the CC approximation (blue dashed), and the A/S approximation (purple dot-long dash) for the benchmark mixing parameters ¿m2 — 0.1eV2 and sin (20) — 0.025 at three different temperatures. The distributions have been normalised to the relativistic Fermi-Dirac distribution /0.

step-like feature at low momenta, such as that produced by oscillations here. Consider first vava-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 / > /0. Annihilation processes on the other hand enhance / via a slightly different mechanism. For processes involving only states above the step, the equilibrium is maintained. For processes such as va(khigh)^a(plow) ^ a(k')a(p'), however, the reaction will be pushed to the left because there is a deficit of Va(plow) relative to a(k ') and a(p'). Thus, annihilation can also overproduce va at high momenta, leading to / > /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 / > /0 at T — 10 MeV, albeit only at very high momenta. The annihilation effect is captured by the factor in equation (2.16), while the separate treatment of ^ accounts for the scattering effect. The latter constitutes the main role of the term in the A/S repopulation approximation (2.16); in contrast, the 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 / for any choice of momentum.

Notwithstanding its success at capturing some degree of the / > /0 phenomenon, the

2 O M U1

0.95 0.90 f 0.85 0.80

1.002 1.000 0.998 0.996 0.994

Scatterings

5 10 15

Annihilations

Full collision term A/S approx.

x = k/T

Figure 5. The active momentum distributions at T = 10MeV calculated using a collision term incorporating only scattering processes (top) and only annihilations (bottom). The difference between the full treatment (black solid) and the A/S approximation (purple long-dot dash) is much larger for the scatterings than for the annihilations. We have used the benchmark mixing parameters ¿m2 = 0.1eV2 and sin (20) = 0.025, and the distributions have been normalised to the relativistic Fermi-Dirac distribution /q.

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 ///0 ~ 0.85 for most momenta, while the annihilations are very effective at bringing ///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 Req and Deq 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 Dee, which as discussed in section 3.2 affects directly the effective sterile neutrino production rate (3.2).

3.5 Dependence on mixing parameters

So far we have tested the various approximation schemes using the benchmark mixing parameter values ¿m2 = 0.1 eV2 and sin2 20 = 0.025. In reality, however, the errors caused by these approximations are also sensitive to the parameter values.

Figure 6 shows ANeff at T = 0.1 MeV as a function of ¿m2 and sin2 20 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 ANff incurred

Figure 6. ANeff at T = 0.1 MeV using the full collision term.

respectively by the equilibrium, CC, and A/S approximations relative to the full solution in the (¿m2, sin2 20)-parameter space.

As expected, the deviations always occur, independently of the approximation scheme, in a diagonal band corresponding to the transition region from ANej = 0 to ANej = 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 ¿m2 values above 0.01eV2, but overestimates ANeff by more than 0.1 at ¿m2 < 0.001 eV2. We can understand this deviation by looking at the conversion temperature. The temperature of maximal conversion is proportional to (¿m2)1/6 [21], so that low values of ¿m2 generally correspond to low conversion temperatures. If the conversion temperature is sufficiently high (T » 1MeV), repopulation is rapid and ANeff is limited only by how fast 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 ANeff contours to deviate from straight lines in relation to ¿m2 and sin2 20 in figure 6. The equilibrium approximation errs in its overestimation of the repopulation rate, yielding almost straight ANeff contours even in the low ¿m2, high sin2 20 region.

In contrast, the top right panel of figure 7 shows that the CC approximation generally underestimates ANeff, but works somewhat better at ¿m2 < 0.01 eV2. For high ¿m2 values the underestimation is due mainly to the undersized Dec damping term which diminishes the sterile neutrino production rate as already discussed in section 3.2. For low ¿m2 values, however, the agreement becomes better (a deviation of 0.02 at ¿m2 = 10-4, sin2 20 = 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

2 O M U1

Neff, eq. approx.

(0.1MeV) - Neff, fnii(0.1MeV)

log(sin2(20))

0.160 0.120

0.040 0.020 0.010 0.001 -0.001 -0.010 -0.020 -0.040

Neff, A/S approx. (0.1MeV) - Neff, fuii(0.1MeV)

-4 -3 -2 log(sin2(20))

0.160 0.120 0.080 0.040 0.020 0.010 - 0.001 0.001 0.010 0.020 0.040

Figure 7. Deviations in the ANeff(T = 0.1 MeV) values computed using various approximations from the full solutions. Top left: the equilibrium approximation. Top right: the CC approximation. Bottom: the A/S approximation. Note that the colour scale is linear, but the contour levels are not.

2 O M U1

to the equilibrium approximation). This overproduction comes about as follows. Although technically the CC repopulation integral incorporates only annihilation processes, numerically the constant 2(Ya)2/3.15 = 0.317 is rather larger than its annihilation counterpart in the A/S approximation, Ce,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 5m2, high sin2 20 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

4 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 (ve, vs)-oscillations implemented it in the computation of ANej 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 appro2ximations of the collision terms that a 1 eV sterile neutrino coupled with a mixing angle sin2 20 ~ 0.1 produces ANeff « 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 ANeff at better than 0.04 for ¿m2 > 0.01 eV2, 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 ¿m2 > 0.0001 eV2, but improves for low ¿m2 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 ANff down to 0.001 for most of the o parameter region, although minor deviations 0.01) remain in the ¿m2 < 0.001 eV2 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 ¿m2 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 — ff(ANeg) « 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 ANeff 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. Y3W 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

Ra,s,e = iEk / d3k'd3p'd3p ¿(4)(k + p - k' - p') ff (A.1)

2yn° J Ek' EpEp'

x |M|2(v«(k),e-(p)K(k'), e-(p')) [fa(E^)/e(Ep)(1 - fa(Efc))(1 - f(Ep))

-fa (Efc )/e(Ep)(l - /va (Efc' ))(1 - /e(Ep' ))]

where the matrix element is [26, 27]

|M|2 = 32GF (4;(k ■ p)(k' ■ p') + (k ■ p')(k' ■ p) - ■ k')) , (A.2)

with Ae = 2 sin2 0W + 1, AM)T = 2 sin2 0W - 1, Be = 2 sin2 0W, and BM)T = 2 sin2 0W. 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

q = p + k = p' + k', v = p — k' = p' — k, w = k — k' = p' — p,

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, Pliimacher and Wong in [19].

Taking the s-channel as a worked example, we use the 3-momentum variable q as a reference and define around it a coordinate system

q = |q|(0,0,1),

k = (0, sin n, cos n), k' = (cos 0 sin 0, sin 0 sin 0, cos 0).

With this choice, we find the quantities

s = (p + k)2 = (p' + k')2 = (Ep + Efc)2 - |q|2, 2

, i/i s - m2 k ■ p = k ■ p = —2—,

|q - k'|2 = |q|2 + |k'|2 - 2q ■ k' = |q|2 + E2, - 2|q|Efc/ cos 0,

|q - k|2 = |q|2 + |k|2 - 2q ■ k = |q|2 + E2 - 2|q|Efc cos n,

and consequently

|Ms|2 = 8 GFA^((Ep + Efc)2 - |q|2 - m2)2 (A.3)

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,

f ¿<4)(k + p - k' - p' ) = / d3p' dE/(E2 - ^ ^"P. 0(Ep - m.)

J 2Ep J 2^J|p'|2 + m2

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]

x ¿(Efc + Ep - Ek - Ep)¿3(k + p - k' - p') (A.4)

¿(Efc + Ep - Ek - ^|q - k'|2 + m2) V|q - k'|2 + m2

0(Efc + Ep - Ek - m.),

we use the property

(x))=f £o (x.)|

¿(x - Xi)

, 2 , l2 2, ¿(Ei -*/|pi|2 + m2) ¿(Ei W|pi|2 + m2) ¿(E2 - |pi|2 - m2) =-/ + V

and hence

V |pi|2 + m2 2^ |pi|2 + m2

to further simplify equation (A.4) to [19]

J dp¿(4)(k + p-k'-p')

= ¿((Efc+Ep - Efc')2 - |q-k' |2-m2) 0(Efc+Ep - Efc'-m.) = ¿ ((Efc+Ep - Efc' )2 - |q|2-E2 +2|q|Efc' cos 0 - m2) 0(Efc+Ep - Efc'-m.)

i ¿ /cos0-E2'-(E+Ep-Ek')2 + N2+m2) 0(Ek+Ep-Efc'-m.)

2|q|Efc' V 2|q|Efc.

In a similar way, we rewrite

2Ep = / d3p dEp ¿(Ep2 - |q - k|2 - m2) 0(Ep - m.)

f 3 1 ( E2 - Ep + |q|2 + m2\ ,

= J d3qdEp¿1cosn- k 22|q|Ek-" 0(Ep-^

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 n/2), the s-channel contribution

Ra s e s —of d cos nd cos 9d|q| dEk' dEp ¿(cos 9 — ... )^(cos n — •••)

27J (A.7)

x |Ms|2F9(Efc + Ep — Ev — me)9(Ep — me),

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 9 and cos n involve only Dirac delta functions. Taking the cos n-integral as an example and noting that integration limits can be equivalently expressed as step functions, we find

||k| — |p|| = |Efc — ^Ep2 — m21 < |q| < Efc + ^Ep2 — m2 = |k| + |p|.

Integrating over d cos 9 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

a , s , e ,s = 27"3E2 dEk' / dEp/ d|q||Ms|2F

2'n°Efc J0 Jmax(me,k'-fc+me) J (A.9)

x 9(|q| — max(||k| — |p||, ||k'| — |p'||))9(min(|k| + |p|, |k'| + |p'|) — |q|)

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 va arising from neutrino scattering with electrons. The full collision term, including scattering and annihilation processes with e^-i,(-g, 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

(1 d( u t Ek — Ep2 + |q|2 + m2\

Ld(cos n)' rs n--—) 0

J, E2 — Ep2 + |q|2 + m^ JE2k — Ep2 + |q|2 + m2 + \

= 2|qE +l) • ^

The two step functions can be reinterpreted as limits on | q| , and together they confine | q| to

the interval O

Reaction (a = ,0) S |M |2

Va(k)V/(p) ^ Va(k')V/(p') 32GF(k • p)(k' • p')

Va (k)V/ (p) ^ Va(k')%(p') 32GF(k • p')(k' • p)

Va(k)Va(p) ^ Va(k')Va(p') 32GF2(k • p)(k' • p')

Va(k)^a(p) ^ Va(k')^a(p') 32GF4(k • p')(k' • p)

Va(k)e—(p) ^ Va(k')e— (p') 32GF ((2xw ± 1)2(k • p)(k' • p') + 4xW(k • p')(k' • p) —(2xw ± 1)2xwm2(k • k'))

Va (k)e+(p) ^ Va(k')e+ (p') 32GF ((2xw ± 1)2(k • p')(k' • p) + 4xW(k • p)(k' -p') —(2xW ± 1)2xWm2(k • k'))

Va(k)^a(p) ^ V/(k')^/3(p') 32GF(k • p')(k' • p)

Va(k)^a(p) ^ e—(k')e+ (p') 32GF ((2xw ± 1)2 (k • k')(p • p') + 4xW (k • p')(k' • p) —(2xW ± 1)2xWm2(k • p))

Table 1. Matrix elements for all relevant reactions involving va, with xW = sin2 9W = 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 vas in the initial state, we have further multiplied the matrix element by 2 to account for the fact that va(k)va(p) ^ ... and va(p)va(k) ^ ... constitute two identical processes. Where there is a choice of ±, the plus signs are for a = e, and the minus signs for a = t. The corresponding matrix elements for Va 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:

Ra.s—

a, s—channel —

dEp / d|q| S|Ms|2F

Ra.i—

a, t—channel

Ra, u—channel

27n3Efc|k| ,/mfc' "" Jmax(mp,fc' — fc+mp) ' J ' (A.10) x 0(|q|-max(||k|-|p||, ||k'|-|p'||))0(min(|k| + |p|, |k'| +|p'|)-|q|),

o7 dEk dEW d|w| S |Mt |2F

2 n°Efc|k| Jmy Jmax(mp,fc' — fc+mp) J (A.11)

x 0(|w|-max(||p|-|p'||, ||k|-|k'||))0(min(|p| + |p'|, |k| + |k'|)-|w|),

dEpj d|v| S|Mu|2F

-27 3 E k dEfc' ^p

2 n Efc|k| Jmy Jmax(mp,fc'—k+mp)

x 0(|v| —max(||k| — |p'||, ||k'| — |p||))0(min(|k| + |p'|, |k'| + |p|)-|v|).

(A.12)

After inserting the matrix element S|Mx|2 and momentum distributions F, these reduced integrals are valid for any 2 ^ 2 process.

2 O M U1

A.3 The full collision terms

The matrix elements for all elastic and inelastic processes involving va at temperatures T < 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.

1. (p) ^ (p'):

rEk +EP

2(2n)3Ek 7o

'max(Q,Ek/ —Ek)

' max( | Ek-EP|,|2Efc,— Ep—Ek |)

(Ep + Efc )2 — |q|2

(A.13)

f (Ek )f (Ep )(1 — f (Efc ))(1 — f (Ep)) — f (Efc )f (Ep)(1 — f (Ek' ))(1 — f (Ep ))

2. va(fc)z/ß(p) ^ Va(k')vß(p'):

a>s>ß 2(2n)3E,2 ./o

rmin(2Ek +Ep-Ek ,Ek' +EP)

'max(Q,Ek/ —Efc)

|Efc'-Ep|

(Ep — Efc' )2 — |v|2

(A.14)

f (Ek' )f (Ep )(1 — f (Efc ))(1 — f (Ep)) — f (Efc )f (Ep)(1 — f (Ek' ))(1 — f (Ep ))

INJ O M U1

3. Va(k)va(p) ^ (k')va(p'):

/^2 rœ

(2n)3E, Jo

rEfc +EP

'max(Q,Ek/ —Ek)

/max(|Efc-Ep|,|2Efc'-Ep-Efc |)

(Ep + Ek )2 — |q|2

(A.15)

f (Ek' )f (Ep' )(1 — f (Ek ))(1 —f (Ep )) —f (Ek )f (Ep)(1 — f (Ek' ))(1 — f (Ep ))

4. Va(fc)z/a(p) ^ (k')Va(p'):

_ 2GF r dE r

Ra,s,a — ~ ^tty) I dEk'

(2n)3E2 Jo

(Ep — Ek )2 — |v|2

max(0,Ek' —Ek) 2

/•min(2Ek +Ep —Ek ,Ek' +EP) |Efc' — Ep|

(A.16)

f (Ek )f (Ep )(1 — f (Ek ))(1 —f (Ep )) —f (Ek )f (Ep)(1 — f (Ek ))(1 — f (Ep ))

5. v«(k)e (p) ^ v«(k')e (p'):

2(2n)3E2 ,

[ dIqIв(IqI -max(|Ek-|p||, |Ek'-|p'||))e(min(Ek + |p|,Ek' + |p'|)-|q|)

(A.17)

x (2xw ± 1)2((Ep+Ek)2-|q|2-m2)2

+ J d|v|e(|v|-max(|Ek-|p'||, Ek '-|p||))e(min(Ek + |p'\,Ey + |p|)-|v|) x 4xW ((Ep - Ek)2 I v 12-m2)2

+ У d|w|e(|w|-max(||p|-|p' ||, |Ek - Ek' |))e(min(|p| + |p' |,Ek+Ek' )-|w|) x 4m2(2xw ± 1)xw ((Ek - Ek' )2-|w|2)

f (Ek )f (Ep )(1-f (Ek ))(1-f (Ep))-f (Ek )f (Ep )(1-f (Ek ))(1-f (Ep ))

6. v«(k)e+(p) ^ v«(k')e+(p'):

e+ —

o/o \3T?2 dEk' dEp

2(2П) Ek -/o Jmax(me,Ey-Efc +me)

(A.18)

d|q|в(Iq|-max(|Ek-|p||, |Ek'-|p'||))e(min(Ek + |p|,E^ + |p'|)-|q|)

x 4xW((Ep+Ek)2-|q|2-m2)2

+ / d|v|e(|v|-max(|Ek-|p'||,Ek-|p||))e(min(Ek + |p'\,Ey + |p|)-|v|)

x (2xw ± 1)2((Ep - Ek' )2 I v 12-m2)2

+ J d|w|e(|w|-max(||p|-|p'||, |Ek-Ek'|))e(min(|p| + |p'|,Ek+Ek')-|w|) x 4m2(2xw ± 1)xw ((Ek - Ek' )2-|w|2)

f(Ek)f(Ep)(1 -f(Ek))(1 -f(Ep))-f(Ek)f(Ep)(1 -f(Ek))(1 -f(Ep))" .

2 О M

7. v«(k)z/«(p) ^ ve(k')z/g(p'):

2(2n)3Ek2 Л

œ /-Efe +Ep r min(2Efc +Ep-Efc' ,Ek' +Ep)

dEp dEk d|v|

(Ep-Ek' )2-|v|2

|Efc'-Ep I

(A.19)

f (Ek' )f (Ep )(1-f (Ek ))(1-f (Ep)) - f (Ek )f (Ep )(1-f (Ek ))(1-f (Ep ))

8. v«(k)z7a(p) ^ e-(k/)e+(p/):

G r<x rEk +Ep-me

= 9(9n)F3E 2 dEp / (A-20)

2(2n)3Efc ./min(0,2me-Efc)

d|q|0(|q|-max(|Efc-Ep|, ||k'|-|p'||))0(min(Efc + Ep, |k'| + |p'|)-|q|) x (2xw ± 1)4xwm2(|q|2-(Ep+Efc)2)

+ J d|v|0(|v|-max(|Efc-|p' ||, ||k' |-Ep |))0(min(Efc + |p' |k' | + Ep)-|v|) x 4xW((Ep-Ek)2 |v|2-m2)2

+ f d|w|0(|w|-max(|Ep-|p' ||, |Efc-|k' ||))0(min(Ep + |p' |,Efc + |k' |)-|w|)

x (2xw ± 1)2((Efc-Ek)2 |w|2-m2) f (Ek )f (Ep )(1 - f (Efc ))(1 -f (Ep))-f (Efc )f (Ep)(1 - f (Ek ))(1 - f (Ep ))

/7/ 2\ 2 2 2 3 X

dx(a - x ) = ax--ax +---h constant,

dx(a - x2) = ax —3 + constant.

Ca>0 = A J dnfcdnfc'dnpdnp 5e(MfcV) ^ V2[v«(k),i/«(p)|i(k/), )]fo(p)fo(k),

= A ^ dnfcdnfc'dnpdnp 5e(kp|kV) ^ V2[v«(k),j(p)Mfc'),j(P)]fo(p)fo(k), Ca,v = A f dnfcdnfc'dnpdnp 5e(M*y) ^ V2[v«(k),j(p)K(fc'),j(P)]fo(p)fo(k),

We remind the reader again that |p' | is not a free parameter, but is constrained by energy conservation.

The integrals over |q|, |v| and |w| can be evaluated analytically, and it turns out that they fall into two different functional forms:

2\2 2 2 3 . x5 , ,

The remaining two integrals over Ep and Ek 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 = 2n/J dnkkfo(k) is a normalisation factor.

vdnp ¿e(fcp|fcV) Y] V

J = Va,Va

Ca,0 = Ay dnkdnk'dnpdnp ¿e(kp|k'p')/o(k)

x [j>2K(k),v7«(p)|i(k'),f(p')]/o(k')/o(p ')

+ £ V2K(k),j (p)K(k' ),j(p')]/o(p)(1 — /o(p')) C,1 = A/**dnv« ¿E(MkV)/o(*)

X [ £ V2[Va(k), i>„(p)|!(k'), i(P)]/o(p)(1 — /o(p) — /o(k))

+ X) V2 [Va (k), j (p)K(k' ),j(p')]/o(k')(/o (p') — /o(p)) + £ V2[v«(k),j(p)|v«(k'),j(p')]/o(p) ,

C„, = A/**dnv« ¿E<kp|kV)/o(*)

X ]T V 2[Va(k), j (p)|va(k '),j(p' )](/o (k' )/o(p') — /o(k' )/o(p) — /o(p)/o(p')).

See equations (2.16) and (2.17) for the implementation of these coefficients in the A/S repopulation and damping terms.

References

[1] C. Giunti, M. Laveder, Y.F. Li and H.W. Long, Pragmatic View of Short-Baseline Neutrino Oscillations, Phys. Rev. D 88 (2013) 073008 [arXiv:1308.5288] [inSPIRE].

[2] Planck collaboration, P.A.R. Ade et al., Planck 2015 results. XIII. Cosmological parameters, arXiv: 1502. 01589 [inSPIRE].

[3] R. Foot and R.R. Volkas, Reconciling sterile neutrinos with big bang nucleosynthesis, Phys. Rev. Lett. 75 (1995) 4350 [hep-ph/9508275] [inSPIRE].

[4] S. Hannestad, I. Tamborra and T. Tram, Thermalisation of light sterile neutrinos in the early universe, JCAP 07 (2012) 025 [arXiv:1204.5861] [inSPIRE].

[5] A. Mirizzi, N. Saviano, G. Miele and P.D. Serpico, Light sterile neutrino production in the early universe with dynamical neutrino asymmetries, Phys. Rev. D 86 (2012) 053009

[arXiv: 1206. 1046] [inSPIRE] .

[6] N. Saviano et al., Multi-momentum and multi-flavour active-sterile neutrino oscillations in the early universe: role of neutrino asymmetries and effects on nucleosynthesis,

Phys. Rev. D 87 (2013) 073006 [arXiv:1302.1200] [inSPIRE].

[7] Y.-Z. Chu and M. Cirelli, Sterile neutrinos, lepton asymmetries, primordial elements: How much of each?, Phys. Rev. D 74 (2006) 085015 [astro-ph/0608206] [inSPIRE].

[8] S. Hannestad, R.S. Hansen and T. Tram, How Self-Interactions can Reconcile Sterile Neutrinos with Cosmology, Phys. Rev. Lett. 112 (2014) 031802 [arXiv:1310.5926] [inSPIRE].

[9] B. Dasgupta and J. Kopp, Cosmologically Safe eV-Scale Sterile Neutrinos and Improved Dark Matter Structure, Phys. Rev. Lett. 112 (2014) 031803 [arXiv:1310.6337] [inSPIRE].

[10] N. Saviano, O. Pisanti, G. Mangano and A. Mirizzi, Unveiling secret interactions among sterile neutrinos with big-bang nucleosynthesis, Phys. Rev. D 90 (2014) 113009 [arXiv:1409.1680] [inSPIRE].

[11] A. Mirizzi, G. Mangano, O. Pisanti and N. Saviano, Collisional production of sterile neutrinos via secret interactions and cosmological implications, Phys. Rev. D 91 (2015) 025019 [arXiv:1410.1385] [inSPIRE].

[12] X. Chu, B. Dasgupta and J. Kopp, Sterile Neutrinos with Secret Interactions — Lasting Friendship with Cosmology, arXiv:1505.02795 [inSPIRE].

[13] B.H.J. McKellar and M.J. Thomson, Oscillating doublet neutrinos in the early universe, Phys. Rev. D 49 (1994) 2710 [inSPIRE].

[14] G. Sigl and G. Raffelt, General kinetic description of relativistic mixed neutrinos, Nucl. Phys. B 406 (1993) 423 [inSPIRE].

[15] N.F. Bell, R.R. Volkas and Y.Y.Y. Wong, Relic neutrino asymmetry evolution from first principles, Phys. Rev. D 59 (1999) 113001 [hep-ph/9809363] [inSPIRE].

[16] M. Archidiacono, E. Giusarma, S. Hannestad and O. Mena, Cosmic dark radiation and neutrinos, Adv. High Energy Phys. 2013 (2013) 191047 [arXiv:1307.0637] [inSPIRE].

[17] M. Costanzi, B. Sartoris, M. Viel and S. Borgani, Neutrino constraints: what large-scale structure and CMB data are telling us?, JCAP 10 (2014) 081 [arXiv:1407.8338] [inSPIRE].

[18] A.I. Akhiezer and S.V. Peletminski, Methods of Statistical Physics, Pergamon Press, (1981).

[21] K. Enqvist, K. Kainulainen and M.J. Thomson, Stringent cosmological bounds on inert neutrino mixing, Nucl. Phys. B 373 (1992) 498 [inSPIRE].

[22] J.M. Cline, Constraints on almost Dirac neutrinos from neutrino - anti-neutrino oscillations, Phys. Rev. Lett. 68 (1992) 3137 [inSPIRE].

[23] S. Hannestad, R.S. Hansen and T. Tram, Can active-sterile neutrino oscillations lead to chaotic behavior of the cosmological lepton asymmetry?, JCAP 04 (2013) 032

[arXiv: 1302. 7279] [inSPIRE] .

[24] K. Kainulainen, Light Singlet Neutrinos and the Primordial Nucleosynthesis, Phys. Lett. B 244 (1990) 191 [inSPIRE].

[25] T. Basse, O.E. Bjaelde, J. Hamann, S. Hannestad and Y.Y.Y. Wong, Dark energy properties from large future galaxy surveys, JCAP 05 (2014) 021 [arXiv:1304.2321] [inSPIRE].

[26] J. Flaig, Neutrino-Oszillationen in thermischer Umgebung, Diploma Thesis, Max-Planck-Institute for Physics and Astrophysics (1989).

[27] S. Hannestad and J. Madsen, Neutrino decoupling in the early universe, Phys. Rev. D 52 (1995) 1764 [astro-ph/9506015] [inSPIRE].

[28] particle Data Group collaboration, K. Olive et al., Review of Particle Physics, Chin. Phys. C 38 (2014) 090001.

[29] A.D. Dolgov, S.H. Hansen and D.V. Semikoz, Nonequilibrium corrections to the spectra of massless neutrinos in the early universe, Nucl. Phys. B 503 (1997) 426 [hep-ph/9703315] [inSPIRE].

[19] F. Hahn-Woernle, M. Pliimacher and Y.Y.Y. Wong, Full Boltzmann equations for leptogenesis including scattering, JCAP 08 (2009) 028 [arXiv:0907.0205] [inSPIRE].

[20] K. Kainulainen and A. Sorri, Oscillation induced neutrino asymmetry growth in the early universe, JHEP 02 (2002) 020 [hep-ph/0112158] [inSPIRE].