{'fraction_non_alphanumeric': 0.03550626718233074, 'fraction_numerical': 0.010033704549337557, 'mean_word_length': 4.437837837837838, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Relic neutrino asymmetry evolution from first principles

The exact Quantum Kinetic Equations for a two-flavour active-sterile neutrino system are used to provide a systematic derivation of approximate evolution equations for the relic neutrino asymmetry. An extension of the adiabatic approximation for matter-affected neutrino oscillations is developed which incorporates decoherence due to collisions. Exact and approximate expressions for the decoherence and repopulation functions are discussed. A first pass is made over the exact treatment of multi-flavour partially incoherent oscillations.


I. INTRODUCTION
Taken together, the solar neutrino [1], atmospheric neutrino [2] and LSND [3] results require at least four light neutrino flavours if neutrino oscillations are to be their explanation. Consistency with the measured invisible width of the Z boson then demands that any additional light flavours be sterile with respect to ordinary weak interactions. With these important experimental results as partial motivation, the cosmological consequences of sterile neutrinos have recently been re-examined [4][5][6][7][8][9][10][11]. It has been shown [5] that oscillations between active and sterile neutrinos can create significant asymmetries between the number densities of relic neutrinos and antineutrinos. 1 These asymmetries then have extremely important consequences for Big Bang Nucleosynthesis (BBN), through their subsequent suppression of active to sterile neutrino oscillations [7,8], and through the effect on BBN nuclear reaction rates of a large ν e asymmetry [9,10].
Various approximations have been explored in the study of neutrino oscillations in the early universe, with neutrino asymmetry or lepton number evolution being an important aspect of these investigations [12][13][14]. The pioneering analyses neglected the thermal spread in relic neutrino momenta: the evolution of the neutrino ensemble was assumed to closely track the evolution of neutrinos having the mean momentum p ≃ 3.15T . Other approximations often utilised were the rate equation approximation also neglecting the thermal momentum distribution, or a Pauli-Boltzmann approach when the momentum spread was not neglected. In the case of active-sterile oscillations, the number densities of sterile neutrinos were often neglected if they were small. Pauli blocking, except for the BBN reactions, has usually been neglected. Various other more subtle approximations have also been made. For instance, in the presence of a nonzero neutrino asymmetry, the collision rate for neutrinos differs from that of antineutrinos. This difference has always been neglected. The repopulation of active neutrino distributions by weak interactions is another process often treated approximately.
The main purpose of this paper is to clarify the nature of these various approximations, and to discuss their regions of applicability more fully and systematically than has been done in the past. Our approach will take the known exact Quantum Kinetic Equations (QKEs) [14,15] for a two-flavour active-sterile system and then develop systematic approximations to them, focussing on the evolution of lepton number. Approximate treatments of lepton number evolution have in the past been vital in obtaining analytical insight. In fact, the discovery that partially incoherent active-sterile oscillations can create lepton number exponentially quickly was made through a treatment that neglected the momentum spread and treated the evolution in the centre of the MSW resonance very approximately [5]. The effect of removing these approximations was subsequently investigated, and a more exact treatment of lepton number creation obtained [7]. Nevertheless, the simple and approximate starting point proved essential in the discovery of this interesting effect, whose existence had hitherto been overlooked because of the complicated nature of more exact treatments. Another important motivation for studying approximation schemes is the computationally intensive nature of numerical analyses based on the exact QKEs. The neutrino momentum spread is computationally demanding, especially during the BBN epoch when neutrinos are out of thermal equilibrium.
Another goal of this paper is to address multi-flavour effects in the early universe from first principles for the first time. Previous analyses have either considered a two-flavour system in isolation, or have approximated the multi-flavour system via pairwise two-flavour subsystems [7,8,11]. (Note that some genuine three-flavour effects were considered for fully coherent oscillations in Ref. [10].) We will clarify how the "pairwise two-flavour approximation" arises from a more exact treatment. The three-flavour problem is particularly important, because one of the principal applications of lepton number creation lies in the suppression of the large angle ν µ − ν s oscillations that may solve the atmospheric neutrino anomaly [16]. The lepton number responsible for this suppression is quite possibly or plausibly created by small angle ν τ − ν s oscillations [7,8]. If so, then a full three-flavour analysis of the ν µ , ν τ and ν s system is warranted. This paper is structured as follows: Section II discusses exact and approximate treatments of the evolution of lepton number in the context of a simple two-flavour active-sterile system. Section III briefly discusses the important case where the strictly sterile neutrino is replaced by a mirror neutrino. Section IV is devoted to a study of multi-flavour effects, and Section V is a conclusion.

II. THE TWO-FLAVOUR ACTIVE-STERILE NEUTRINO SYSTEM: FROM FIRST PRINCIPLES TO USEFUL APPROXIMATIONS
The early universe contains, in part, an ensemble of neutrinos. 2 Their evolution in time is affected by three phenomena. The first is simply the expansion of the universe. The second is decohering collisions with the background medium. The third is, of course, the coherent process of oscillations governed by a matter-dependent effective Hamiltonian. The effect of collisions and oscillations is quantified using the Quantum Kinetic Equations. The QKEs generalise the Pauli-Boltzmann Equations to include quantum coherence between the particle species, in our case the active and sterile neutrinos (and active and sterile antineutrinos).
We will focus on the two-flavour system comprising one active and one sterile neutrino species in this section. We do so because two-flavour active-sterile oscillations lead to the rapid creation of lepton number when the oscillation parameters are in the correct region. The multi-flavour case is deferred to Section IV. Our task is to track the evolution of the neutrino and antineutrino ensemble, focusing in particular on the evolution of lepton number. We will take as our initial, high temperature, condition that the number densities of sterile neutrinos and antineutrinos are zero. This is expected simply because, by definition, sterile neutrinos do not feel electroweak (or strong!) interactions and so will decouple very early. 3 We will focus on the epoch between T = m µ ≃ 100 MeV and BBN (T ∼ 0.1 − 1 MeV), so that the background plasma consists of essentially only photons, electrons, positrons, neutrinos, antineutrinos and whatever amount of sterile neutrinos and antineutrinos get generated through oscillations.

A. The exact Quantum Kinetic Equations
We first write down and discuss the exact QKEs, partly by way of review and to define notation, and partly to discuss approximations to the decoherence and repopulation functions (to be defined shortly).
Consider a two-flavour active-sterile system composed of ν α and ν s , where α denotes either e, µ or τ . The 1-body reduced density matrices describing ν α − ν s and ν α − ν s oscillations are given, respectively, by where plus an analogous expression for antineutrinos. The functions P 0 and P are just the coefficients of a convenient expansion of ρ in terms of Pauli matrices σ and the identity matrix (P is often called the "polarisation"). The quantity p is the neutrino or antineutrino momentum. It will be understood that P 0 and P depend on time (or, equivalently, temperature via t ≃ m Planck /11T 2 ) as well as momentum. The diagonal entries of ρ (ρ) are relative number density distributions in momentum space of ν α (ν α ) and ν s (ν s ): where N eq (p, 0) is a reference number density distribution that we have set equal to the equilibrium Fermi-Dirac distribution, with zero chemical potential. This choice amounts to normalising ρ such that Tr(ρ) = 2 if thermal equilibrium holds and the chemical potential is zero. The antineutrino relative number densities are given similarly. Since the physical interpretation of ρ is related to number density distribution ratios, the expansion of the universe does not directly contribute to the time evolution of ρ. Note that N eq (p, 0) depends on time only through the expansion of the universe, and not on coherent or incoherent matter effects. The evolution of P 0 (p) and P(p) is governed by the QKEs and ∂ ∂t P 0 (p) = R(p).
These equations are obtained by evolving the full density matrix for all particles in the plasma forward in time by means of the S matrix, and then tracing over all degrees of freedom other than ν α and ν s . The essential difference between the Quantum Kinetic and Pauli-Boltzmann approaches is that the former time evolves amplitudes while the latter evolves probabilities. See, for example, Ref. [14] for a detailed derivation of the QKEs and further discussion. The antineutrino QKEs take an identical form, but with the substitutions We now discuss the form and significance of each of the terms above. 4 The function V(p) describes the quantally coherent part of the matter-affected evolution of the ν α − ν s subsystem. It is given by with The quantities ∆m 2 and θ 0 are the mass-squared difference and vacuum mixing angle for ν α − ν s oscillations, respectively. We define the mass eigenstate neutrinos ν a,b by with θ 0 defined so that cos 2θ 0 ≥ 0 and ∆m 2 ≡ m 2 b − m 2 a . The function V α is the effective matter potential [17]. To leading order, including the lowest order finite temperature term, it is given by [18] with the dimensionless functions a(p) and b(p) given by where T is temperature, G F is the Fermi constant, m W is the W -boson mass, ζ is the Riemann zeta function, A e ≃ 17, A µ,τ ≃ 4.9 and The quantities on the righthand side of this equation are the various asymmetries, given by where n f = N f (p)dp (19) are number densities (note that we will abbreviate L να by L α in this Section). The a(p) term in Eq. (14) is the leading order density-dependent (Wolfenstein) contribution to the effective potential, while the b(p) term is the leading order finite temperature term. (Note that the Wolfenstein term also depends on temperature, but only indirectly through the cosmological red-shifting of fermion number density distributions.) Observe that V depends on ρ through the dependence of a(p) on L α . As we discuss later, this makes the evolution of L α non-linear. It is important to notice that the dependence of V on L α is an O(G F ) effect (due to coherent forward scattering), rather than an O(G 2 F ) effect. The quantity η is a small term related to the cosmological baryon-antibaryon asymmetry. For antineutrinos, the corresponding function V is obtained from V by setting and replacing λ by λ. In other words, one simply replaces L (α) by −L (α) . The Mikheyev-Smirnov-Wolfenstein (MSW) [17,19] resonance conditions are for neutrinos and antineutrinos, respectively. Note that prior to the creation of lepton number, both resonance conditions are simply b(p) = cos 2θ 0 , which can only be satisfied if This condition plays an important role (see later). It is also important to realise that the resonance conditions of Eq.(21) at a given temperature T are met only for neutrinos and antineutrinos of a certain momentum p = p res . This is why the spread of neutrino momenta plays a significant role in the analysis [7]. The decoherence or damping function D(p) quantifies the loss of quantal coherence due to ν α collisions with the background medium. Its exact expression is derived in a general form in Ref. [14], and is given in the Appendix by Eq.(A1). When thermal equilibrium holds, the general expression for D reduces to a useful compact form, given by where Γ α (p) is the total collision rate for ν α 's with momentum p. By examining the momentum dependence of the specific weak collision processes operating between the neutrino decoupling temperature and 100 MeV, we see that where is the average momentum for a relativistic Fermi-Dirac distribution with zero chemical potential (indicated by the subscript "0"), and Γ α is the thermally averaged total collision rate for ν α . This last quantity can be expanded as a power series in L α , with y e ≃ 4, y µ,τ ≃ 2.9, z e ≃ 0.1 and z µ,τ ≃ 0.04. The y α term, is the thermally averaged collision rate at zero chemical potential, while the z α term is the first order correction due to a finite neutrino asymmetry. For antineutrinos, the corresponding expression is In previous studies [7,8], Eq.(24) was adopted as an approximate result for Γ α (p). Actually, it turns out to be exact, as we show in the Appendix.
It is important to observe that although D is in principle a dynamical quantity through its dependence on neutrino number densities and hence on ρ [see Eq.(A1)], it becomes approximately kinematic provided two conditions hold: (i) thermal equilibrium and (ii) lepton number being small. If these are so then D(p) ≃ y α G 2 F T 5 p/2 p 0 . The fact that D is most of the time just a given external function of temperature and momentum rather than a dynamically evolving quantity simplifies the numerical solution to the equations. Observe that the absence of thermal equilibrium is correlated with sufficiently weak interaction rates, so D is necessarily small in that case relative to V (for relevant choices of oscillation parameters). If, on the other hand, thermal equilibrium holds, but the neutrino chemical potentials are nonzero and evolving in time, then the dependence of D on ρ will only be important when the asymmetries are very large. For interesting choices of oscillation parameters, this only occurs at temperatures approaching the BBN epoch when collisions are again becoming insignificant. It turns out that in practice one may neglect all but the leading neutrino asymmetry independent term y α G 2 F T 5 p/2 p 0 in the applications we have in mind, a point we will discuss more fully later on.
Pauli blocking is neglected in the calculation of D, and also the repopulation function R discussed below. To incorporate this effect, appropriate factors of (1 − f ), where f (p) ≡ 2π 2 N(p)/p 2 , need to be inserted into the integrals. To estimate the size of the error introduced by neglecting Pauli blocking, we can evaluate f at the peak of the equilibrium momentum distribution, p ≃ 2.2T . Since f (2.2T ) ≃ 0.1 we estimate the error to be at most 10%.
We now discuss the repopulation or refilling function R(p) which controls the evolution of P 0 (p). Since the evolution of P 0 is governed by processes that deplete or enhance the abundance of ν α 's with momentum p. The rate of change of P 0 (p) clearly receives no contribution from quantally coherent ν α − ν s oscillations, because the two flavours simply swap. It therefore equals the rate at which ν α 's of momentum p are generated by scattering processes, minus the rate at which ν α 's of momentum p are scattered out of that momentum value. Its general form is where and V (i, j|i ′ , j ′ ) denotes an interaction matrix element (defined in the Appendix) for the process i + j → i ′ + j ′ .
As for the decoherence function D, the general expression for R in Eq.(30) simplifies when thermal equilibrium holds. For temperatures between neutrino decoupling and 100 MeV, it is given by where N eq (p, µ) is the Fermi-Dirac momentum distribution function with chemical potential µ, and µ α is the chemical potential for ν α . This expression follows from Eq.(30) when all number densities except for ν α are thermal, and ν α is approximately thermal (see the Appendix for the derivation). The chemical potential µ α is obtained through the lepton asymmetry using The physical interpretation of this expression for R is that all weak interaction processes involving ν α are tending to send its number density towards the equilibrium configuration. Equation (32) was adopted in previous works on a heuristic basis [10]. In fact, it is an exact result under the conditions stated above (see the Appendix). Observe that above the neutrino decoupling temperature, ν α is approximately in thermal equilibrium and so R is very small (this issue is discussed more carefully in the next Subsection). Below the neutrino decoupling temperature, the form given in Eq.(30) should in principle be used. For antineutrinos, R is obtained from R by replacing µ α by −µ α (when thermal equilibrium holds). Clearly, the dependence of the neutrino momentum distributions on the repopulation function will be most sensitive at low temperatures and when the asymmetry has evolved to large values. At high temperatures, the weak interaction rates are fast enough to thermalise the distributions rapidly, spreading out an asymmetry across the distribution (as described by the chemical potentials). At lower temperatures, this does not occur as efficiently, so the correct form for the repopulation function becomes more complicated [see Eq.(30)] if we want to track the actual form of the neutrino momentum distributions. We will discuss the low temperature (T ≃ a few MeV) regime later on in Subsection II E.

B. Evolution of neutrino asymmetry
We now focus on extracting the evolution of lepton number or neutrino asymmetry L α from the Quantum Kinetic Equations. Recall that this is a crucial quantity because the Wolfenstein a(p) term in the effective matter potential [see Eq. (14)] is proportional to a linear combination of fermion, and in particular neutrino, asymmetries. If this term is appropriately large, then active-sterile oscillations will be suppressed at low temperatures. Also recall that L e affects primordial light element abundances through weak interaction BBN reaction rates.
In principle, the QKEs (7) and (8) (plus the antineutrino equations) are solved to indirectly yield the evolution of the asymmetry defined by Eq. (18). It is, however, very important to obtain simpler approximate evolution equations for L α . There are two reasons for this. First, it is difficult to understand the qualitative behaviour of the evolution of the asymmetry from the in-principle QKE procedure per se. Analytic insight requires that the evolution of L α be directly considered and approximations introduced. Second, the direct numerical solution of the full QKEs is computationally demanding because of the neutrino momentum spread and the MSW resonance phenomenon.
We begin our study of approximate evolution equations by observing that Eqs. (5), (18) and (19) combine to produce as the time evolution equation for the α-like neutrino asymmetry L α . We have used the fact that N eq (p, 0)dp/n γ is time independent, and also conservation of lepton number through So far, no approximations have been made (except that possible tiny flavour changing neutral current effects for massive neutrinos have been neglected. 5 ) Using Eq.(29) and the corresponding expression for antineutrinos, Eq.(34) implies that To proceed, we need approximate solutions for P y and P y . These are obtained by solving the QKEs for P and the corresponding antineutrino equations. In expanded form, the QKEs (7) are which also require Eq. (8).
The qualitative character of the evolution goes through distinct phases. We first discuss the high temperature initial conditions. As stated previously, it is expected that the initial 5 These effects will always be small, with their exact size being model dependent.
abundance of sterile species is very close to zero. This implies that P z ≃ 1 at high T , so that N s ≃ 0. At temperatures above the neutrino decoupling temperature, the neutrinos are by definition in thermal equilibrium, so the repopulation function R equals zero according to Eq.(32). This means that P 0 is a constant which depends on the initial neutrino chemical potential and hence on the initial neutrino asymmetry. We will also take the initial neutrino asymmetry to be negligible (it could be of the order of the baryon asymmetry, for instance 6 ). If this is so, then P 0 ≃ 1.
The first phase occurs at high temperatures. With the asymmetry set to its initial value of zero, the function λ is dominated by the finite temperature term of the effective potential Replacing p by its thermal average, we see that the magnitude of this term increases with temperature as T 5 . The function β effectively goes to zero at sufficiently high temperature because it scales as T −1 . The decoherence function D, like λ, scales as T 5 . (Note that, numerically, |V α /D| ≃ 60 or 24 for α = e or µ/τ , respectively.) Setting P 0 to be constant, and neglecting β, the QKEs (37) become Therefore P z (t) is frozen at its initial value P z (t) ≃ 1, and no sterile neutrinos are produced. Writing the approximate QKEs (39) become which are trivially solved to yield Since, at high temperatures, D is large compared to the expansion rate of the universe we see that |P | and hence also P x and P y get exponentially damped to zero very rapidly from any initial values. (Observe that they rapidly oscillate about zero as they do so due to the large λ.) The system is therefore completely incoherent ("in a mixed state") at high temperatures, and frozen (Quantum Zeno or Turing Effect).
The third phase occurs at very low temperatures where the collisional terms D and b are negligible. The neutrino subsystem then evolves as a coherently oscillating subsystem which is also coupled through R to the background. The oscillations are matter-affected if a sufficiently large neutrino asymmetry L (α) has been created by this time. This is the regime relevant for the BBN epoch, for interesting choices of oscillation parameters. We will discuss this regime in a later Subsection.
Intermediate between these two limiting regimes lies the second phase. During this phase, the system emerges from its frozen initial state. The evolution is driven by an interplay between collisions, which are still rapid compared to the expansion rate, and oscillations, which begin to affect the neutrino number densities through a non-negligible β (recall that β is proportional to sin 2 2θ 0 where θ 0 is the vacuum mixing angle). This is the regime of the "static approximation" introduced in Ref. [7]. It is during this epoch that neutrino asymmetries will evolve to values orders of magnitude higher than the baryon asymmetry provided the oscillation parameters are in the appropriate region. We will now provide new analytical insight into this phenomenon by deriving in a new way the very useful approximate form for the asymmetry evolution equation considered in Ref. [7]. We will see that the static approximation is actually a combination of an adiabatic approximation for partially incoherent oscillations, and a small β expansion. This point was not crisply realised hitherto.
We now take the first crucial step, that continues to be a good approximation. With this approximation, the QKEs (37) simplify to the homogeneous equations or, in a self-evident matrix notation, The continued justification for Eq.(43) is as follows: Above the neutrino decoupling temperature, most of the ν α ensemble is in thermal equilibrium. From Eq.(32) we see that R is therefore zero or very small. One might be tempted to conclude that in fact R is exactly zero above the neutrino decoupling temperature. This, however, cannot strictly be the case. Suppose ν α 's of some given resonance momentum p res oscillate strongly to ν s 's. Instantaneously, N α (p res ) goes to zero or close to it, but P 0 (p res ) does not change due to these oscillations because of the generation of a nonzero N s (p res ). However, the absent momentum mode of ν α 's is quickly repopulated from the background medium, leading clearly to an overall nonzero value for ∂P 0 ∂t evaluated at the momentum p res . We would therefore expect that R is very small above the neutrino decoupling temperature, except when the oscillations are very strong. This means that the approximation in Eq.(43) will be a good one except possibly in the centre of an MSW resonance when oscillations might be rapid.
The word "might" is used because the oscillations will be non-adiabatic in some regions of parameter space. So, we proceed with the understanding that our approximation scheme might not be strictly valid in the centre of a resonance. We will have cause to revisit the centre of the resonance later on.
To solve Eq.(44), we first introduce the "instantaneous diagonal basis" through where U is a time-dependent matrix that diagonalises K, with k 1,2,3 being eigenvalues. In the instantaneous diagonal basis, Eq. (44) becomes This equation resembles, but is not the same as, the MSW evolution equation written in the instantaneous mass eigenstate basis with a varying background density. In that case, the adiabatic approximation sees the time derivative of the instantaneous mixing matrix set to zero. Our next important approximation after Eq. (43) is to analogously set so that the evolution equation becomes This is clearly within the family of adiabatic-like approximations. (We defer the discussion of the region of applicability of this approximation to a later subsection.) The difference with the usual adiabatic approximation for MSW evolution is the existence of quantum decoherence through a nonzero D. We will now show that Eq.(50) reproduces the approximate evolution equation derived in Ref. [7] using what was termed the static approximation. The static approximation of Ref. [7] is therefore revealed to be just the adiabatic approximation for a partially incoherent system of neutrinos undergoing matter affected oscillations (in the small β limit -see below). (Note that this type of idea was briefly discussed in Ref. [13].) Equation (50) can be formally solved to yield which in turn implies that The mixing matrix U −1 is obtained by placing the normalised eigenvectors κ i in the columns, where and The inverse matrix U is then composed of the row vectors One of the eigenvalues will be real and the other two will be a complex conjugate pair: In order to progress further, we need to calculate the eigenvalues. Since we are dealing with a 3 × 3 matrix, this is somewhat awkward algebraically. Fortunately, the small β limit is of great interest. The numerical details of the small β limit will be given in a later subsection. In this limit it is easy to show that the eigenvalues are approximately given by Since we are still in the regime where the decoherence parameter D is significant, it follows that Furthermore, Combining these results with Eqs.(52-56) we obtain that For the antiparticle system we obtain a similar expression, so the neutrino asymmetry evolution equation (36) becomes having used and the corresponding antineutrino expression. From Eq.(60) we see that lepton number is constant if the number density distributions of active and sterile species are equal. This is to be expected, because oscillations effectively do nothing in that case.
Since we are working within the R = R = 0 approximation, we now substitute N eq (p, µ) for N α and N eq (p, −µ) for N α . (The chemical potentials of ν α and ν α are equal and opposite because of the rapid ν α ν α ↔ e + e − process.) Since we are taking the initial lepton number to be small, we will also expand N eq (p, µ) to first order in to obtain This allows us to form the quantities N + α and N − α , which are given by The asymmetry evolution equation (60) becomes, where c ≡ cos 2θ 0 , s ≡ sin 2θ 0 , and is the collision rate with the chemical potential set to zero. The quantity ∆ is an O(L α ) correction term given by and δ is an additional O(L α ) correction term which arises from allowing D to be different from D, where x 0 is x as in Eq.(66) but with Γ 0 in place of Γ α . The evolution equation (65), including the correction term ∆ but excluding δ, was studied in detail in Ref. [7]. From Eq.(65) we see that the effect of having different collision rates for the neutrinos and antineutrinos is to introduce the additional correction, δ, as well as a small correction through x and x to the denominators. The ∆ and δ correction terms, as well as the O(L α ) corrections to x and x in the denominators may readily be incorporated into the numerical study of the evolution equation (65), however to a good approximation we may neglect these corrections. It turns out that δ only amounts to a small contribution such that |δ| |∆| ∼ 0.1.
The corrections to the denominators are very small, even at a resonance. At the initial MSW resonance where a ≃ 0, . The approximate evolution equation (65) has proven to be very useful for gaining analytical insight and indeed for performing numerical work. The δ = 0 version was studied in depth in Ref. [7]. The derivation presented above shows for the first time that the full glory of Eq.(65) follows from first principles (that is, the QKEs). In Ref. [7], a heuristic derivation based on a Pauli-Boltzmann-like approximation was used. A minor issue associated with the definition of x arises on closer examination. In Ref. [7], this quantity contains an additional term sin 2 2θ 0 that is absent in the present definition. It turns out that keeping higher order terms in β in the derivation of Eq.(60) does not in fact reproduce this missing term as one might hope. The mystery of the missing sin 2 2θ seems to require an explanation at a deeper level, as yet unprobed. Note also that when the sterile neutrino number density is negligible, Eq.(65) becomes a self-contained non-linear differential equation describing the evolution of the single momentum-independent variable L α . It is a considerable simplification over the full QKEs which are coupled differential equations for the eight momentum-dependent functions P, P 0 , P and P 0 .
In order to integrate Eq.(65) when the sterile neutrino number density is not negligible, an evolution equation for N s is needed. Using and the QKEs with the same approximations as above, we easily obtain For completeness, we now briefly recall the qualitative features of lepton number evolution which are revealed by Eq.(65). For further discussion see Ref. [7]. Initially, neutrino asymmetries are small and hence the term a is small. A critical factor in the subsequent evolution of L α is the overall sign of its derivative. (Note also the distinction between L α and L (α) .) If it is of the opposite sign to L (α) , then the asymmetry L α will evolve such that If, on the other hand, the derivative is of the same sign as L (α) , then the evolution of L α will be such as to increase L (α) . For ∆m 2 > 0, the former situation always obtains. No explosive creation of lepton number can occur. However, when ∆m 2 < 0 the situation is more complicated. (Observe that for small mixing, ∆m 2 < 0 means loosely that the sterile neutrino is less massive than the active neutrino). In this case the quantity a has the same sign as L (α) . The quantity b is always positive, but it is a decreasing function of time. The sign of the derivative is therefore controlled by (c − b). At a given temperature T , the function (c − b) equals zero for neutrinos and antineutrinos of momentum p c , where Neutrinos and antineutrinos with momenta less than p c make a positive contribution to the righthand side of Eq.(65), whereas those with momenta greater than p c make a negative contribution. The critical momentum p c is very small at high T but increases with decreasing T , which means that the derivative changes sign from negative to positive once T is small enough for p c to be near the peak of the Fermi-Dirac distribution. (For maximal mixing, c = 0 and the derivative is always negative. The creation of lepton number occurs due to small angle oscillations.) Since a is initially very small, the equation is also the MSW resonance condition for both neutrinos and antineutrinos. The denominator in the integrand of Eq.(65) is thus at a minimum for neutrinos and antineutrinos with p = p c ≃ p res . The resonance momentum p res is generally found by solving where the plus (minus) sign pertains to neutrinos (antineutrinos). Provided the vacuum mixing angle is small enough, the critical temperature T c at which the sign of the derivative changes is well approximated by taking p c ≃ 2.2T c which corresponds to the peak of the Fermi-Dirac distribution. From Eq. (74), it is given approximately by for ν e − ν s (ν µ,τ − ν s ) oscillations. (Numerically, the critical temperature is found to be slightly higher than this estimate [7]. Also, when the vacuum mixing angle is sufficiently large, a non-negligible sterile neutrino number density is produced, causing the actual T c to differ from this estimate.) Consider now the parameter space region chosen so that the mixing angle is small and so that T c occurs above the neutrino decoupling temperature. At T ≃ T c , explosive exponential growth of lepton number begins 7 because Eq. (65) is then of the form with the proportionality factor augmented by the MSW resonances in the denominator of the integrand. After a short time, lepton number has grown sufficiently for the quantity a to no longer be negligible. Non-linearity through the denominator then alters the qualitative character of the evolution from exponential growth into a slower phase governed approximately by a( p ) ≃ ± cos 2θ 0 , depending on whether the sign of the asymmetry is positive or negative. This equation is the resonance condition for neutrinos (antineutrinos) when |a| ≫ |b|.
As the temperature continues to decrease, the collision dominated epoch begins to give way to an epoch in which coherent oscillations take over as the dominant driver of neutrino evolution. By the time the neutrino decoupling temperature is reached, coherent oscillations certainly dominate over collisions. This is the regime relevant for the BBN epoch. Equation (65) no longer provides a good description of lepton number evolution. The evolution of lepton number is, of course, in principle obtained by solving Eqs. (7) and (8), using the general expression given in Eq.(30) for the function R. (Numerically, setting D, b ≃ 0 is a good approximation in this regime.) This in-principle procedure has, to our knowledge, not yet been attempted in practice. Instead, the low temperature calculations that have actually been performed have utilised the usual adiabatic approximation for fully coherent matter-affected 7 Note that while the mixing angle θ should be small, it must be large enough so that sin 2 2θ > ∼ 5 × 10 −10 (eV 2 /|∆m 2 |) 1/6 is satisfied, otherwise the partially incoherent oscillations are not strong enough to create lepton number. See Ref. [7] for more details.
oscillations, and the growth of the asymmetry has been calculated using an equation which determines how fast the resonance momentum moves through the distribution converting active into sterile neutrinos. The derivation of this equation from the QKEs will be presented, for the first time, in Subsection II E. Furthermore, the high-temperature expression for R as given in Eq.(B4) has been adopted as an approximation to the much more complicated and rigorous expression defined in Eq.(30). We may estimate the theoretical error that this approximate repopulation method introduces into the final result by comparing the value of the asymmetry obtained in the manner described above with that obtained by integrating over the neutrino momentum distributions to find the total number densities. If everything is consistent, these should of course agree. However, for the calculations presented in Ref. [10], we find discrepancies of up to 10 to 20%. The repopulation should not change the total asymmetry, only distribute it from the particular momentum state at which it is created, across the momentum distribution. The discrepancy, which must be due to inconsistencies in the repopulation process, is a good indication of the theoretical error in all of the calculations [9,10] performed thus far. This theoretical error is small enough for the conclusions so far reached to be essentially unchanged.

D. Region of applicability for the approximations
The validity of Eq.(60) rests on the assumptions that the quantity β is small relative to the size of D and λ, and that their rates of change with respect to time are negligible. We now examine the regions where these conditions are satisfied.
The small β approximation is sufficiently accurate provided that The "small β expansion" is more properly to be thought of as an expansion in the dimensionless quantity β/ √ D 2 + λ 2 . In order to explore the bound, we will replace p by its thermal average ≃ 3.15T .
Away from resonance, |λ| ≫ D, so the small β expansion holds provided that |β| ≪ |λ|. For T > T c , λ ≃ ∆m 2 b/2p, which means that the small β limit requires for α = e and µ/τ , respectively. We would like the small β limit to hold just prior to lepton number creation at T = T c . Putting T = T c , as given by Eq.(77), in Eq.(82) we see that |∆m 2 | cancels out of the inequality, leaving as the constraint. Note that the use of Eq.(77) presupposes that the mixing angle is sufficiently small for the sterile neutrino number density to be small, so this is not a serious constraint.
In the centre of the lepton-number creating resonance at T = T c , λ momentarily goes to zero. The small β expansion will be valid in the centre of the resonance, provided the more stringent constraint |β| ≪ D holds. This translates into the requirement that Using Eq. (77), we find that |∆m 2 | cancels out of the inequality, leaving the constraint on the vacuum mixing angle. We now discuss the applicability of the adiabatic-like approximation which assumes the complete negligibility of the term U ∂U ∂t −1 that appears in Eq.(48). Evaluated explicitly in the small β limit, this hitherto discarded quantity takes on the form to the lowest order in β. These non-vanishing diagonal and off-diagonal entries contribute to corrections to both the eigenvalues and eigenvectors of the matrix K d in Eq.(47). Thus our immediate task is to ensure that these corrections are small in the region of interest. The 33 entry of the matrix U∂U −1 /∂t represents a "first order" correction to the eigenvalue k 3 . Imposing the condition and rendering it into a more illuminating form, we find that where the momentum p has been replaced with the mean momentum p . Note that in deriving the above, we have set λ ≃ 0, i.e., the resonance condition, and we have used the relation t ≃ m P lanck /11T 2 . Thus Eq.(90) purports the failure of the static approximation for resonances occurring at temperatures below roughly 3 MeV. Together with the resonance condition c ≃ b ≃ 1, Eq.(90) effectively places a lower limit on |∆m 2 |, that is for ν e − ν s (ν µ, τ − ν s ) oscillations. Additional bounds arise from evaluating Eq.(89) alternatively at λ ∼ D, i.e., by requiring that The relative importance of the two terms above is not immediately obvious by inspection. Suppose the first term dominates in the region of interest. The resulting inequality is but a less severe version of Eq.(90). Supposing now that the second term is predominant and that there are no accidental cancellations such that the inequality roughly reduces to two separate conditions While the former is essentially another encrypted bound on the temperatures at which our approximations are valid, similar to that obtained earlier in Eq.(90), the latter amounts to demanding that at the mean momentum. This condition agrees remarkably well with that obtained in Ref. [7] and should be checked when integrating the static approximation equation (65) for selfconsistency. We shall not reproduce here the constraints resulting from the (U∂U −1 /∂t) 11 and (U∂U −1 /∂t) 22 corrections to the eigenvalues k 1 and k 2 respectively, as they appear to be much less stringent than those arising from earlier considerations. Indeed, this is to be expected given that the correction terms are only of the order β 2 where the small β limit is effective.
Turning our attention now to the effects of U∂U −1 /∂t on the instantaneous eigenvectors, we observe that, using a perturbative analysis, the "first order" corrections to Q where and we have neglected terms proportional to β 2 . Thus the single requirement that must be satisfied in this case in the small β limit is It turns out that the above condition does not in fact lead to more severe bounds than those following from other requirements.

E. Oscillation dominated asymmetry evolution
At low temperatures, such that D ≪ |β|, |λ|, the collision rate becomes much slower and the generation of lepton number is dominated by oscillations. Taking the D = 0 limit, we can demonstrate that the QKEs reduce to usual expression for MSW transitions.
We begin by setting D = 0 in Eq.(44). In this case, setting ∂U −1 /∂t = 0 exactly corresponds to the usual adiabatic approximation for MSW evolution. The eigenvalues now become and the matrix U −1 (= U † ) is given by With the initial conditions P x (0) ≃ 0, P y (0) ≃ 0, the solution for P z (t) is given by where c 2θm = cos 2θ m and The number density of sterile neutrinos is then given by Since the adiabatic transition probability for a neutrino propagating through a medium where the mixing angle changes from θ m (0) to θ m (t) is given by we have so, as expected, the adiabatic approximation for the QKEs is equivalent to the adiabatic limit of the MSW effect if we turn off the decohering collisions. The low temperature equation, used in previous work [9,10] to calculate the growth of asymmetry by determining how quickly the resonance momentum moves through the distribution, can be derived from the QKEs in the D = 0 limit.
Since we have b ≃ 0, the resonance condition c ± a ≃ 0 can only be attained for either the neutrinos or antineutrinos, depending on the sign of the asymmetry. Assuming for definiteness that L α > 0, which implies a > 0, only the antineutrino resonance will contribute to the growth of lepton number, so that where dP 0 dT ≃ 0 has been assumed. From Eq.(100) we find where we have used a/c = p/p res and da dT = − aT pres d dT pres T . This results in where δp res is the resonance width, since the last term in Eq.(107) ensures the integral is strongly peaked about the resonance momentum. Averaging the integrand over the width of the resonance leads to so that in the limit of zero resonance width we recover where the N α − N s term is interpreted as evaluated just before the resonance passes a particular momentum value. This term expresses the difference between the number of ν α and ν s which are converted at the resonance, while the factor d dT pres T accounts for how quickly the resonance momentum moves through the distribution converting active to sterile neutrinos. See Refs. [9,10] for applications of this type of equation. The derivation given above for Eq.(109) is the first justification of it from first principles.

III. MIRROR NEUTRINOS INSTEAD OF STERILE NEUTRINOS
We now briefly discuss how the foregoing must be altered to deal with mirror instead of strictly sterile neutrinos. Mirror neutrinos form the neutral lepton sector of the mirror fields postulated in the Exact Parity Model (EPM) [20]. The EPM restores exact parity invariance to the microscopic world by parity-doubling the Standard Model. Interestingly, if neutrinos and mirror neutrinos have mass, and if the two sector mix, then in the absence of intergenerational mixing the mass eigenstates are maximal superpositions of ordinary and mirror neutrinos. Maximal oscillations of ν µ with its mirror partner ν ′ µ can explain the atmospheric neutrino anomaly, while maximal oscillation of ν e with its mirror partner ν ′ e can similarly solve the solar neutrino problem. The maximal mixing of ν α with ν ′ α is a consequence of the unbroken parity symmetry. This form of parity symmetry is arguably one of the most credible explanations for the maximal mixing of muon-neutrinos observed in atmospheric neutrino experiments. The LSND anomaly can be accomodated by switching on small intergenerational mixing with the appropriate mass hierarchy.
The issue of whether or not the EPM explanation of the neutrino anomalies is consistent with Big Bang Nucleosynthesis is an important one. We assume that in the very early stages of the Big Bang ordinary matter predominates over mirror matter. (See Refs. [21] for speculations about how this could arise.) The question is then whether interactions between the ordinary and mirror sectors at later times overcreates mirror matter, spoiling BBN. Focussing on ordinary-mirror neutrino mixing, we note that a sufficiently large ( > ∼ 10 −5 ) preexisting neutrino asymmetry can always be invoked to suppress ordinary-mirror oscillations, thus leading to standard BBN despite the presence of the mirror sector [4]. However, it is more interesting to suppose that the neutrino sector of the EPM "saves itself" through the production of large neutrino asymmetries in the manner discussed above. In order to do this analysis, the QKEs have to be modified to take into account the interactions of mirror neutrinos amongst themselves and with other mirror species such as the mirror photon and the mirror electron. See Ref. [8] for the proof that the EPM is consistent with BBN in an interesting region of parameter space.
Mirror neutrinos feel mirror weak interactions. The interaction strength is given by G F , since exact parity symmetry imposes equality of coupling constants and gauge boson masses between the sectors. We may classify the additional terms in the QKEs through powers of G F : coherent forward scattering induces O(G F ) terms, while incoherent scattering induces terms of O(G 2 F ) and higher. Provided that the number densities of mirror neutrinos remains low, all terms other than the O(G F ) coherent forward scattering effects can be neglected. If the number densities of mirror neutrinos become appreciable, then many of the complexities inherent in active-active scattering become an issue.
So, provided one is working in a parameter space region where the mirror neutrino number densities remain small, it is consistent to alter the QKEs of the strictly sterile neutrino case by the simple substitution, for ν α − ν ′ β oscillations (α, β = e, µ, τ ). The quantity L ′(β) is simply the mirror version of L (β) . This substitution fully incorporates the coherent forward scattering of mirror neutrinos off mirror neutrinos induced by mirror weak interactions.
It is important to understand why these O(G F ) mirror neutrino self-interactions are as important as O(G F ) neutrino self-interactions. The term a(p) in the effective matter potential is equal to the sum of terms directly proportional to the difference in number densities between neutrinos and antineutrinos, and the difference in number densities between mirror neutrinos and mirror antineutrinos. These number density differences will be of the same order of magnitude simply because mirror neutrinos and antineutrinos are produced, via partially incoherent oscillations, from ordinary neutrinos and antineutrinos. Since we know that the a(p) term dominates the b(p) term once significant neutrino asymmetries have been created, and since coherent effects anyway become more and more important as the temperature approaches the BBN epoch, the O(G F ) mirror neutrino self-interactions are of critical importance. The fact that they are completely non-negligible distinguishes the cosmology of mirror neutrinos from that of strictly sterile neutrinos [22].

IV. MULTI-FLAVOUR OSCILLATIONS
The purpose of this section is to begin the study of multi-flavour effects. In particular, we wish to determine the conditions under which it is appropriate to break up a multi-flavour situation into effective two-flavour subsystems, where two-flavour expressions can be applied. An understanding of when multi-flavour effects may arise is important, because while twoflavour subsystems are a useful approximation, a realistic situation will always involve three or more flavours. A multi-flavour system may involve features that do not occur in simpler two-flavour systems. We would expect, for instance, genuine multi-flavour effects to arise in systems in which there are "overlapping" resonances. With a general multi-flavour system, there is also the possibility of explicit CP violation due to complex phases in the neutrino mixing matrix.
Partially incoherent oscillations in a multi-flavour system have not previously been studied from first principles. In all prior work on lepton asymmetries, multi-flavour systems were dealt with in terms of two-flavour subsystems.
The system we shall consider here, chosen both for definiteness and importance, consists of ν µ maximally mixed a sterile neutrino ν s , and a heavier ν τ which has small mixing with ν µ and ν s . This particular system is well motivated by the atmospheric neutrino anomaly, with the parameter space region sin 2θ µs ≃ 1 and ∆m 2 µs ∼ 10 −2 − 10 −3 eV 2 leading to a resolution of the anomaly via ν µ → ν s oscillations. The ν τ,µ,s system was studied in Ref. [7], where the role of the heavier tau neutrino, taken to have a mass in the range where ∆m 2 τ s = 1−1000eV 2 , is to generate an L τ asymmetry via oscillations with the sterile neutrino as per the dynamics discussed in the previous Section. The L τ asymmetry generates a large Wolfenstein term in the effective potential for the ν µ − ν s subsystem, which then acts to suppress oscillations between the muon and sterile neutrinos. For a range of parameters, this prevents the maximal ν µ −ν s oscillations from bringing the sterile neutrinos into equilibrium, thus circumventing the supposedly stringent BBN bounds on the mixing of sterile and active neutrinos. This phenomenon is vital for reconciling the ν µ → ν s solution to the atmospheric neutrino problem with cosmology. Note, however, that in previous analyses the ν τ − ν s and ν µ − ν s resonances were dealt with in a two-flavour manner; possible three-flavour effects were ignored. Given that the resonances are well separated, this is expected to be a good approximation. However, in order to verify this is the case, we need to study it within a full three-flavour framework.
We begin by writing down the evolution equations for a three-flavour partially incoherent system using the density matrix formalism. While it is in a sense straightforward to generalise the Quantum Kinetic Equations to a three-flavour system, qualitatively different complexities are encountered due to the presence of active-active neutrino oscillations in addition to the more readily understood active-sterile oscillations. These active-active complications result in equations from which it is difficult to obtain analytical insight. 8 However, under the momentum averaged approximation the equations simplify significantly, allowing progress to be made. Although going to the momentum averaged approximation will introduce some errors, it is probably a reasonable approach. In any case, we adopt it here because it makes possible a first attempt at a first-principles analysis of a three-flavour system. (The differences between including and excluding the momentum spread for twoflavour active-sterile oscillations were discussed at length in Ref. [7]. The qualitative features of the evolution were found to be preserved in the mean momentum approximation.) The three-flavour momentum averaged density matrix ρ is parameterised in terms of the SU(3) Gell-Mann matrices λ i , i = 1, . . . , 8 such that The neutrino number densities are related to the diagonal entries of ρ : The evolution of a momentum averaged density matrix is described by the Quantum Rate Equations (QREs) rather than the Quantum Kinetic Equations. Generalising the twoflavour QREs derived in Ref. [14], which involves the lengthy but straightforward procedure of substituting Eq.(112) into the general QKE expression given by Eq. (22) of Ref. [14], we find the three-flavour QREs describing the considered system to be, with the evolution of P 0 , which is related to the total number density of ν τ , ν µ and ν s , given by where R α is the repopulation function [14] R α = 1 (n eq ) 2 j=e,νe,νµ,ντ where h e = 1 4 , h ν = 1 and the values for the collision rates Γ(jj → ν α ν α ) can be found in Refs. [14,23]. The n's are normalised such that n eq να = n eq and n eq e = 2n eq . The G terms arise from the mixed active-active part of the density matrix, with G e = 0.26G 2 F T 5 and G νe = 0.51G 2 F T 5 . The quantity D = 1 2 Γ α is the damping parameter, with Γ α as given by Eq. (27), while D ′ is the equivalent parameter for the active-active ν τ − ν µ oscillations: The parameter C in Eq. (114) is a damping-like term which couples the neutrino and antineutrino density matrices through the off-diagonal ρ µτ and ρ µτ elements, and only arises in the case of mixed active-active neutrinos. For ν τ − ν µ oscillations [14], The quantity H is given by where dp ≡ 1 (2π) 3 d 3 p. This is an order G 2 F T 5 quantity, however it will be small compared to D and C, provided that ρ τ µ is small.
We define where f ijk are the SU(3) structure constants. Finally, V is given by with where ω αβ are the (vacuum) energy eigenvalues in the flavour basis, with U being the mixing matrix relating the neutrino mass and flavour eigenstates. V αβ are the effective potential terms which, as well as the usual diagonal potential terms (V αα = ∆m 2 2p (−a + b)), also include off-diagonal contributions given by The appearance of off-diagonal terms in the effective potential is a consequence of activeactive oscillations. They depend on the quantities P 4 and P 5 which parameterise the mixed ν τ and ν µ states in the density matrix. Note, however, that the number densities of ν τ and ν µ are equal (to order L). As we shall see, this results in P 4 and P 5 being suppressed, with respect to P 1 & P 2 and P 6 & P 7 which parameterise the mixed ν τ,s and ν µ,s subsystems, respectively.
For simplicity, the ν τ − ν µ mixing angle will be set to zero. Observe, however, that there will still be a small effective mixing between ν τ and ν µ indirectly through ν s . We parameterise the mixing matrix U as where φ is the ν τ − ν s mixing angle, and the ν µ − ν s mixing angle has been fixed at π 4 . Unlike the two-flavour case, the neutrino and anti-neutrino density matrices are coupled through the C term in Eq.(114), so that we must consider a system of 16 coupled differential equations for P 1 , . . . , P 8 , P 1 , . . . , P 8 . As in the two-flavour case, we adopt the dP 0 dt ≃ 0 approximation. This allows Eq.(114) to be expressed in the form: where P = (P 1 , ...P 8 , P 1 , ..., P 8 ) and K is the matrix given by where M and C are the submatrices and M is obtained from M via the replacement λ → λ, σ → σ, ǫ → ǫ, D → D, D ′ → D ′ , and H → H. We use the notation, δ I = Imδ, δ R = Reδ, H I = ImH, H R = ReH, and where ∆m 2 = ∆m 2 τ s and r = ∆m 2 µs /∆m 2 τ s ≪ 1. The parameters β and λ take roughly the same form as for the two-flavour active-sterile system, and are related to the ν τ − ν s oscillations. The quantities δ and γ are the equivalent of β for the ν τ − ν µ and ν µ − ν s oscillations, respectively, while σ and ǫ are the ν τ − ν µ and ν µ − ν s analogues of λ.
To solve these equations, P 2 , P 4 , P 5 and P 7 must be determined as functions of time.
A first inspection of Eqs. (127), (128) and (129) gives the impression that the various oscillation modes are coupled together in a non-trivial manner. However, under the assumptions made, the solution reveals that the three oscillation modes effectively decouple to first order in the small parameters β, δ, γ and H.
To solve these sixteen coupled equations, the matrix will be diagonalised in the small β, δ, γ, H limit, in analogy to the two-flavour case. This is done by treating the β, δ, γ, H terms as a perturbation to the zeroth order solution obtained with β = δ = γ = H = 0. The parameters β and γ may be of approximately the same order of magnitude, depending on the mass-squared differences and mixing angles. The size of δ depends on the size of the offdiagonal potential term V τ µ which, as with H, is dependent on the values of P 4 and P 5 . Our perturbative solution relies on the assumption that P 4 and P 5 are small. This is expected for the physical reason of decoherent damping through D ′ , except at low temperatures (say in the neutrino decoupling regime). The solution obtained with this assumption may then be checked for self-consistency.
For the purposes of the perturbation, we treat β, δ, γ, H as being of the same order. Again, as with the two-flavour case, we adopt the adiabatic-like approximation that U ∂U ∂t −1 ≃ 0 where U is the time dependent matrix which instantaneously diagonalises the matrix in Eq.(127). Solving for the appropriate eigenvectors to first order in β, δ, γ, H yields P 2 (t) = − βD D 2 + λ 2 n ντ − n νs P 0 n eq , P 4 (t) = Y 4 n ντ − n νµ P 0 n eq + Z 4 n ν τ − n ν µ P 0 n eq , P 5 (t) = Y 5 n ντ − n νµ P 0 n eq + Z 5 n ντ − n νµ P 0 n eq , where Y 4,5 , Z 4,5 are order δ terms and are functions of the parameters D ′ , σ, σ and C, arising from the diagonalisation of the P 4 , P 5 , P 4 , P 5 submatrix We see that P 5 (and similarly P 4 ) is suppressed not only because δ is small but also since n ντ − n νµ = O(L). This ensures that V τ µ and H remain small, consistent with the assumptions made. The proportionality of P 4 and P 5 to n νµ − n ντ is one manifestation of the unimportance of oscillations between species of roughly equal number densities.
These expressions in Eq.(132) resemble those that would be obtained by breaking the system down into two two-flavour subsystems ν τ + ν s and ν µ + ν s . The evolution equations for L τ and L µ are coupled only via the effective potentials and n νs . In the initial stages of lepton number creation, when the sterile neutrino number density can be neglected, the equations are in essence coupled only through the L-dependent effective potentials. The equations for lepton number become with the P 4 and P 5 terms now neglected. These types of equations, generalised in the obvious way to incorporate the thermal momentum distribution, have been used in practical calculations [7,8]. The above analysis is the first to fully justify their use on the basis of first principles (admittedly in the Quantum Rate Equation approximation). The various oscillation modes would seem only to be coupled via higher order terms in β, γ, δ, H. In fact the restrictions we imposed on the mixing angles can be relaxed, and the only difference between Eq.(134) and the two-flavour equations used in Refs. [7,8] are the small corrections (in terms of the mass-squared differences and mixing angles) to the parameters β, λ, γ and ǫ.
This two-flavour subsystem solution will not be valid when the resonances of different subsystems are close enough to overlap. For example, the ν τ − ν s and ν µ − ν s resonances begin to overlap if λ ≃ ǫ. In this case some of the eigenvalues of the matrix in Eq.(128) become degenerate and hence the perturbative method used to obtain the solution may break down.
The study of the ν τ,µ,s system, with ν µ → ν s parameters set by the atmospheric neutrino anomaly, has recently been re-examined by Foot [11]. In this work, the "pairwise two-flavour approximation" is used, but with each two-flavour active-sterile subsystem treated using the full two-flavour QKEs.

V. CONCLUSION
Active-sterile neutrino oscillations will play an important role in early universe cosmology if light sterile neutrinos exist. The combined solar neutrino, atmospheric neutrino and LSND anomalies require at least one light sterile flavour. The partially incoherent neutrino oscillations that occur in the early universe are described using the Quantum Kinetic Equations. These equations provide an in-principle means of calculating the cosmological consequences of light sterile neutrinos, with Big Bang Nucleosynthesis being an important concern. However, their qualitative physical consequences are not always easy to extract, and they are computationally demanding. Approximation schemes are therefore welcome.
The main insights of this paper are twofold: 1. An adiabatic approximation for the partially incoherent oscillations described by the QKEs has been developed. This clarifies the origin of the "static approximation" discussed in Ref. [7], revealing it to be conceptually related to the usual adiabatic approximation of matter-affected oscillations. In the absence of collisions, we have explicitly shown that it reduces to the usual adiabatic approximation. When collisions are important, we have used it to rederive some very useful approximate evolution equations for lepton number in the small β limit. The evolution of lepton number is of particular importance because of the central role played by the Wolfenstein term in the effective matter potential in suppressing sterile neutrino production prior to BBN. The systematic approach to approximating the QKEs introduced in this paper suggests further development which could improve, in a hopefully practical way, on the adiabatic approximation for partially incoherent oscillations.

2.
A first principles treatment of partially incoherent three-flavour oscillations has been attempted for the first time. Focussing on a ν τ,µ,s subsystem with maximally mixed ν µ and ν s , important for resolving the atmospheric neutrino anomaly, the three-flavour Quantum Rate Equations were explicitly written down. The QREs are obtained from the QKEs by approximating the evolution of the neutrino momentum distribution by the evolution of neutrinos having the mean momentum. By employing a similar adiabatic-like approximation in the small β, γ, δ, H limit (three-flavour static approximation) to this case, it was demonstrated that the three-flavour system separated into two two-flavour subsystems, ν τ + ν s and ν µ + ν s , coupled only through the dependence on the Wolfenstein term on family lepton numbers, and through the common ν s ensemble. The "coupled two-flavour subsystem" approach had previously been used in studies which concluded that the maximal ν µ → ν s solution to the atmospheric neutrino anomaly was consistent with BBN provided that the ν τ −ν s oscillation parameters lay in a particular region. The conclusion that there is no cosmological objection to the ν µ → ν s solution to the atmospheric anomaly can thus be made with even more confidence.
An exciting new era in fundamental physics has begun, ushered in by the beautiful atmospheric neutrino results of SuperKamiokande. We await with great interest news from the Sudbury Neutrino Observatory [24], and further data from SuperKamiokande, regarding the existence or otherwise of light sterile neutrinos. Should they exist, then a new synergy between the microscopic and macroscopic worlds will be revealed through the evolution of partially coherent active-sterile (or ordinary-mirror) neutrino oscillations.
If we assume that all weakly interacting species are in thermal equilibrium with zero chemical potential, then we obtain which, if we set the momentum k equal to its thermal average k 0 , reduces to the standard expression in terms of the thermally averaged collision rate D = 1 2 Γ α . To allow for a nonzero lepton asymmetry, we will assume that ν α and ν α may have non-zero chemical potential, while all other particle species will be assumed to have zero chemical potentials. This results in small correction terms to Eq.(A5) where y e ≃ 4.0, y µ,τ ≃ 2.9, u α ≃ 0.72, v e ≃ 1.0, v µ,τ ≃ 0.8, w α ≃ 0.33, x e ≃ 0.46, and x µ,τ ≃ 0.36.