GeV-scale hot sterile neutrino oscillations: a derivation of evolution equations

Starting from operator equations of motion and making arguments based on a separation of time scales, a set of equations is derived which govern the non-equilibrium time evolution of a GeV-scale sterile neutrino density matrix and active lepton number densities at temperatures T>130 GeV. The density matrix possesses generation and helicity indices; we demonstrate how helicity permits for a classification of various sources for leptogenesis. The coefficients parametrizing the equations are determined to leading order in Standard Model couplings, accounting for the LPM resummation of 1+n<->2+n scatterings and for all 2<->2 scatterings. The regime in which sphaleron processes gradually decouple so that baryon plus lepton number becomes a separate non-equilibrium variable is also considered.

however for GeV-scale sterile neutrinos the rates that we are interested in do not change much in the temperature range between 130 GeV and 160 GeV [20]. Therefore, for the conceptual discussion, we can imagine to work in the "symmetric" phase of the electroweak theory. 2 In the temperature range 130 GeV < ∼ T < ∼ 10 5 GeV, all Standard Model interactions can be assumed to be in thermal equilibrium (this includes lepton chirality flipping processes through the electron Yukawa coupling). Then the state of the system is characterized by a temperature, T , by three lepton chemical potentials, µ a , and by the baryon chemical potential, µ B . Suppose now that we extend the Standard Model through right-handed sterile (gauge singlet) neutrinos, and use these to generate active neutrino masses through the see-saw mechanism. If the mass of the sterile neutrinos is ∼ GeV, then the Yukawa couplings are so small (∼ 10 −7 ) that sterile neutrinos do not equilibrate in this temperature range. Therefore they constitute a non-equilibrium ensemble, evolving "slowly" in a Standard Model background. The Yukawa interactions also imply that lepton number densities are not conserved, so these become slowly evolving variables as well.
The neutrino Yukawa interactions need not be aligned with Standard Model lepton generations. This leads to sterile neutrino oscillations. Because the neutrino Yukawa couplings are tiny, coherent oscillations may be maintained for a long period of time, and sterile neutrinos need to be described by a density matrix.
The momenta of the sterile neutrinos are changed by the same slow interactions as their number densities are. Therefore, kinetic equilibrium cannot be assumed and the density matrix displays a non-trivial momentum dependence.
There is one further slow variable to be tracked, namely the helicity of the sterile neutrinos. As massive particles, sterile neutrinos can carry both helicities. The two helicity states experience different interactions: basically, one state interacts with Standard Model leptons and the other with antileptons. Both states need to be included in the density matrix.
To summarize, we need a density matrix for each sterile neutrino momentum mode, labelled by k ≡ |k|. It turns out that to a good approximation the two helicity states have no direct overlap with each other (cf. discussion at the beginning of sec. 2.5). Therefore, for each k the set of non-equilibrium variables consists of two complex matrices, denoted by ρ (τ )IJ , where τ = ± labels helicity and I, J label generations. In addition the three active lepton asymmetries, denoted by n a , and the baryon asymmetry, denoted by n B , evolve slowly.
The goal of our study is to derive evolution equations for the non-equilibrium variables to O(h 2 ) in neutrino Yukawa couplings. In principal the general form of the equations is valid to all orders in Standard Model couplings (in practice certain small corrections are omitted along the way), however we subsequently evaluate the coefficients at leading order (cf. sec. 3). The physically interesting effects, which are of O(h 4 ) or O(h 6 ) [9], originate from the coupled dynamics of the "slow" oscillations, and will be addressed in a separate study.

Basic variables and equations of motion
In the temperature range considered (T > ∼ 130 GeV) the Higgs mechanism gives a contribution small compared with thermal masses, i.e. m W ∼ gv/2 ≪ gT , where g denotes the SU(2) gauge coupling. Then the masses of sterile neutrinos are directly given by the Majorana masses, assumed real and positive and denoted by M I , I ∈ {1, 2, 3}. The Lagrangian reads l a a Rφ h * Ia N I +N I h Iaφ † a L ℓ a , (2.1) whereφ = iσ 2 φ * is a Higgs doublet; a L , a R are chiral projectors; ℓ a = (ν e) T a is a left-handed lepton doublet of generation a; and h Ia are the components of the neutrino Yukawa matrix.
We consider the so-called ultrarelativistic regime, k ∼ πT ≫ M I , so that a free dispersion relation reads In this regime the vacuum mass is corrected by a thermal effect of O(h 2 T 2 ) [21]. Even though h is small, the thermal correction is relevant because it should be compared with the mass differences M 2 I − M 2 J and because the initial temperature may be high, T ∼ 10 5 GeV. In our formalism, thermal masses originate as a part of the O(h 2 ) corrections, cf. appendix A. Therefore, we treat the kinematics as in vacuum for the moment. The kinematic approximation of eq. (2.2) is frequently invoked in order to simplify the discussion.
The sterile neutrino field operator in the interaction picture can be written as N I (X ) = k 1 2ω k I τ u kτ I a kτ I e −iK I ·X + v kτ I a † kτ I e iK I ·X , (2.3) where k ≡ d 3 k (2π) 3 and K I · X ≡ ω k I t − k · x. In accordance with the Majorana nature of N I , the on-shell spinors are related by v = Cū T , where C is the charge conjugation matrix. The creation and annihilation operators, which are time-independent in eq. (2.3), satisfy the commutation relations {a kτ I , a † qσJ } = (2π) 3 δ (3) (k − q)δ τ σ δ IJ . Sterile neutrinos interact through the Yukawa terms in eq. (2.1). We rephrase the interactions through an interaction Hamiltonian, By j a andj a we denote Standard Model currents from eq. (2.1), j a ≡ a L j a ≡φ † a L ℓ a ,j a ≡l a a Rφ . (2.5) In order to understand the dynamics induced by H int , it is helpful to go over to the Heisenberg picture for a moment (cf. ref. [22] for an analogous discussion, and appendix A.1 for a detailed step-by-step argument). Then the canonical equation of motion for the annihilation operator, defined by expressing the field operator in the form of eq. (2.3) and accounting for any additional time dependences through a kτ I and a † qσJ , becomes An analogous equation is obtained for a † qσJ . The canonical anticommutator remains timeindependent. Similarly, the lepton asymmetries evolve as [23] iL a (t) = For the physical observables that we are interested in, the evolution rate is of O(h 2 ). We extract the rate from an expectation value of an operator like in eq. (2.7). In order to evaluate the expection value, we return to the interaction picture. Then, the time evolution of the density matrix is determined by H int . In particular, assuming that the full density matrix is known at some time t = 0, its time evolution is to first order in h given by (2.8) The physical rate can then be defined as (cf. e.g. ref. [24], and eq. (A.6) for an explanation of intermediate steps) The expectation value with respect to the density matrix ρ full (0) is denoted by ... ≡ Tr {... ρ full (0)}. At the end of the computation, this can be re-interpreted as having been evaluated with the density matrix at time t, since the difference between ρ full (t) and ρ full (0) is of O(h). This way, so-called secular terms can be avoided. Because of the different times scales related to the "slow" and "fast" variables, we can assume the full density matrix to have a block-diagonal form, ρ full = ρ N ⊗ ρ SM , where ρ N is the density matrix associated with the sterile neutrinos. The density matrix associated with the Standard Model degrees of freedom, ρ SM , is in equilibrium at a temperature T and is parametrized by (slowly evolving) chemical potentials µ a and µ B : (2.10) In the canonical formalism there are no other chemical potentials, however in the path integral formalism the hypercharge gauge field gets an expectation value in the presence of µ a , µ B = 0 which effectively generates an additional chemical potential for all fields coupling to the hypercharge field (cf. sec. 4). We note that the operator equation of motion in eq. (2.7) has the forṁ whereas the interaction Hamiltonian in eq. (2.4) can be written as where Y ≡ (t ′ , y). Inserting these structures into eq. (2.9) we get vanish, because within the Standard Model there are to O(h 0 ) direct correlations only between j a andj a . For the same reason, the sum over b is saturated by b = a. Furthermore, we can take the limit t → ∞, given that reaction rates proportional to ∼ h 2 are much slower (with a much larger time scale) than the fast Standard Model rates.

Time evolution of a sterile neutrino density matrix
Let us apply the formalism above to the time evolution of a sterile neutrino "density matrix". We define it asρ where V is the volume, taken to be infinite at the end of the computation. This operator now plays the role of O(t) in eq. (2.9). Making use of the equation of motion in eq. (2.6), eq. (2.9) leads to 2-point correlators of the Standard Model currents j a andj a , defined in eq. (2.5). These can be expressed in terms of Wightman functions, where α, β ∈ {1, ..., 4} represent spinor indices. The operators j a ,j a have a non-trivial commutator with the lepton number operator L a appearing in the density matrix ρ SM , cf. eq. (2.10); as a consequence, following a text-book derivation, the Wightman functions can be expressed in terms of the spectral function, with relations depending on the index a through the chemical potential that is carried by active leptons: Here P = (ω, p) and n F (ω) ≡ 1/(e ω/T + 1) is the Fermi distribution.
In the expectation value following from eq. (2.9), the spectral function is bracketed by the on-shell spinors u and v. The expression can be simplified by making use of the fact that, with the exception of processes involving Yukawa couplings, chirality is preserved by Standard Model interactions at high temperatures. Omitting higher-order contributions involving the Yukawas, 3 chiral invariance implies that the spectral function is proportional to the Dirac matrices γ µ [21]. Then we can use the properties of the charge conjugation matrix C to show thatū Furthermore, chiral invariance implies that amplitudes betweenū and u conserve helicity (helicity states are defined in eq. (3.2), and examples of non-zero matrix elements are shown in eq. (A.22)),ū kσJ a L ρ a a R u kτ I ∝ δ στ ,ū kσJ a R ρ a a L u kτ I ∝ δ στ . (2.23) We make use of this important simplification in the following. The integration over the time t ′ in eq. (2.9) can also be simplified. In an equation like eq. (2.13) we are faced with integrations of the types Here I 1 multiplies terms bracketed betweenū and u, and I 2 those bracketed betweenū and v. Given that ω k I ≈ k + M 2 I /(2k) ∼ 3T , I 2 contains a rapid oscillation, similar to the "fast" Standard Model variations; these rapid oscillations will be omitted. In I 1 the oscillation rate is suppressed by Majorana mass differences; this is a "slow" process and needs to be kept. The large-t value of I 1 can be defined as a limit: Assuming that φ 1 is slowly varying around ω ≈ ω k J , the principal value part is antisymmetric around this point and corresponds to a higher time derivative correction; it amounts to a modification of ω k J through a thermal mass (this is shown in appendix A). We postpone the inclusion of this "dispersive" or "virtual" correction for the moment, focussing first on "absorptive" or "real" effects. For those, we need Re . Inserting the time integral as well as eqs. (2.17)-(2.20) into eq. (2.9), we find that absorptive time evolution is parametrized by the slowly evolving coefficientŝ where K J ≡ (ω k J , k). Noting that ρ a is real (cf. appendix B of ref. [25] for a general discussion), the evolution equation reads where the "equilibrium" terms containing n F originate from {a kσJ , a † kτ L }/V = δ τ σ δ LJ . The terms containingΓ + and n F (ω − µ a ) represent scatterings involving leptons, whereas those withΓ − and n F (ω + µ a ) represent the contributions of antileptons. Physically,Γ + (aτ )IJ describes the rate at which the in-medium wave function of the state (kτ J) gets projected in the direction of (kτ I), andΓ − (aτ )IJ does the same for the charge-conjugated process. It may be noted that the right-hand side of eq. (2.29) vanishes in equilibrium, i.e. if the density matrix is diagonal and all lepton chemical potentials vanish. Its general form is, however, valid both near and far from equilibrium.
It can also be observed that there is no equilibrium term with τ = σ: helicity non-diagonal correlations decrease to zero with time. In particular, if we start from an initial density matrix which is helicity-diagonal, it stays so. However, the values of the coefficientsΓ ± do depend on the helicity index τ (cf. sec. 3).
For formal considerations, it is convenient to have coefficients which are independent of time. This also offers for a simple way to include the dispersive thermal mass corrections mentioned above. This can be achieved by redefining the coefficients in eqs. (2.27) and (2.28) and the density matrix aŝ The evolution equation for ρ then obtains an additional term of O(h 0 ), where the part of O(h 2 ) has the same structure as in eq. (2.29) but with time-independent coefficients (Γ ± ). We remark that in text books density matrices are usually defined in terms of "states" rather than "operators", and then the free evolution has the sign of the Liouville -von Neumann equation, i∂ t ρ = [H 0 , ρ] + ... . In the present paper we defined the sterile neutrino "density matrix" through operators, cf. eq. (2.14). Correspondingly the free time evolution in eq. (2.31) has the same sign as appears in operator equations of motion, cf. eq. (2.6). If desired the sign difference could be eliminated by reversing the ordering of indices [26], however in practice it is inconsequential so we do not bother.

Time evolution of lepton asymmetries
The derivation of sec. 2.3 can be repeated for lepton asymmetries. The starting point is the equation of motion in eq. (2.7), and we take n a (t) ≡ L a (t)/V to play the role of O(t) in eq. (2.9). Apart from neutrino Yukawa interactions, at high temperatures lepton asymmetries are also violated by sphaleron processes. The observables not affected by the latter are the linear combinations n a − n B /3. For these the final result can be expressed in close analogy with eq. (2.29), The coefficientsΓ ± are identical to those in eq. (2.29). There are again two terms, reflecting the fact that lepton asymmetry can increase through the production of leptons or the disappearance of antileptons. Eq. (2.32) is odd in charge conjugation, and only the helicitydiagonal components of the sterile neutrino density matrix contribute.

Simplified form of evolution equations
We noted in the context of eq. (2.29) that non-diagonal helicity components of the density matrix decouple from the diagonal ones, and in eq. (2.32) that only the diagonal ones contribute to lepton asymmetries. Therefore, we omit the non-diagonal components in the following. Moreover, for easier inclusion of thermal mass corrections, we implement the redefinition in eq. (2.30), and denote The evolution equations in eqs. (2.29) and (2.32) can be simplified if we expand the righthand sides to first order in the lepton chemical potentials, assuming µ a , µ B ≪ πT . Within perturbation theory the presence of µ a , µ B = 0 implies that the temporal component of the hypercharge gauge potential develops an expectation value, guaranteeing the hypercharge neutrality of the plasma; this expectation value is conventionally referred to as the hypercharge chemical potential, denoted by µ Y (cf. sec. 4.1). In this limit the coefficients in eqs. (2.27) and (2.28), redefined through eq. (2.30), have the forms (cf. sec. In principle there are also coefficients proportional toμ B , appearing like those proportional toμ Y , however these vanish at leading order (cf. sec. 3) and are omitted here already. The functions Q, R and S, estimated in sec. 3, are found to be real. To a reasonable approximation they are also symmetric in I ↔ J, however this symmetry is broken by the "soft" 1 + n ↔ 2 + n scatterings evaluated in sec. 3.2. Roughly speaking, the coefficient Γ + (aτ )IJ describes the amplitude T J|I 0 , where |... T implies that the state evolves within a medium. Even though 0 J|I 0 = 0 I|J * 0 , it is possible that T J|I 0 = T I|J * 0 . The physical meaning of the equations can be made more transparent by taking the helicitysymmetric and antisymmetric parts of ρ (τ ) as the basic variables. Correspondingly we define Furthermore, in order to streamline the equations, we make use of the kinematic simplification in eq. (2.2). This implies that momenta k < ∼ M I ≪ πT are not treated properly, however their contribution to lepton asymmetries is strongly phase-space suppressed (M I ∼ 10 −2 πT ). Let us first inspect the equation for lepton asymmetries, eq. (2.32). Inserting eqs. (2.2), (2.34) and (2.35) into eq. (2.32) we obtain denote a symmetrization or antisymmetrization over the helicity-flipping and conserving contributions (cf. sec. 3), and a symmetrization over the generation indices, respectively. The first term on the right-hand side of eq. (2.37) is a "washout term", decreasing any lepton asymmetry towards zero. It agrees with the corresponding term derived from different considerations in ref. [23]. The second and third terms are "source terms", generating a lepton asymmetry. They display a product of structures manifesting Sakharov-type conditions, namely deviation from thermal equilibrium and CP violation. One of the sources originates from a helicity-symmetric ρ and the other from a helicity-asymmetric one. The parts proportional to chemical potentials in eqs. (2.39) and (2.40) are of second order in the language of linear response theory, containing a product of two deviations from equilibrium (chemical potentials and a non-thermal density matrix). We include them in the equations, given that the density matrix may deviate significantly from equilibrium.
Analogous equations are obtained for the density matrices. Displaying them as a pair of complex Hermitean matrices with generation indices, we obtaiṅ The coefficient matrices read (2.47) The coefficients C ± generate non-thermal density matrices if lepton asymmetries are present.
The coefficients D ± are "washout terms" in the sense that they drive the system towards equilibrium, but they are often called "production rates", because normally ρ + ≪ ½n F (k).
Then D + produces ρ + and D − produces ρ − . In the limit of a single generation D + agrees with a term derived from different considerations in refs. [22,27]. The Hermitean matrices H 0 and ∆ 0 on the first lines of eqs. (2.42) and (2.43) contain the vacuum energies but also the thermal mass corrections (cf. appendix A.2), The first term in eq. involving the Fermi distribution, should also contain the eigenvalues of the system described by H 0 and ∆ 0 as arguments. This is, however, a higher-order effect in the ultrarelativistic regime k ∼ πT ≫ M I . We conclude by remarking that a set of equations similar to eqs. (2.37), (2.42) and (2.43) was obtained in ref. [15], however without the inclusion of the hypercharge chemical potentialμ Y (and hence of the IR sensitivities discussed in sec. 3) and under the assumption that terms of orderμ a ρ − could be dropped, as they represent a double deviation from equilibrium. Moreover helicity conserving contributions were neglected together with the generation (IJ) dependence of the helicity flipping ones. As discussed in sec. 3 and confirmed numerically in sec. 3.6, the latter assumptions hold well for M I ≪ T . Specifically, upon settingμ Y → 0, (2) , dropping terms proportional toμ a ρ − , and recalling the discussion below eq. (2.31) concerning the sign of the commutator terms, we can reproduce eqs. (2.16) and (2.19) of ref. [15].

Rate for fermion number non-conservation
It is well-known that if all Majorana masses are set to zero in eq. (2.1), then the theory has an additional conserved charged, which we may call the "fermion number". Indeed the Majorana spinor can be replaced by a chiral Dirac spinor in this case, and the conserved charge then counts the total asymmetry in right-handed and left-handed leptons. Keeping instead the Majorana character intact, the fermion number is defined as the sum of the helicity asymmetry of the sterile neutrinos and the lepton asymmetry of the Standard Model sector. It is easy to demonstrate that eqs. (2.37) and (2.43) respect this symmetry for where the coefficients read  For determining the dependence on τ , the form of the spinor u kτ I is needed. We can write where the spinors satisfy τ =± η τητ = 1 2 (½ + γ 0 ). The precise form of η τ depends on the representation chosen for the Dirac matrices. Examples for the standard and Weyl Figure 1: Examples of 1 + n ↔ 2 + n scatterings contributing to the sterile neutrino spectral function (the spectral function is a cut, i.e. "amplitude squared" of such processes, convoluted with the appropriate distribution functions). Sterile neutrinos are denoted by a double line, whereas arrowed, dashed, and wiggly lines correspond to Standard Model fermions, scalars, and gauge fields, respectively.
representations are Here |τ is a helicity eigenstate, i k i σ i |τ = τ |k||τ , where σ i are the Pauli matrices.

1 + n ↔ 2 + n scatterings
Physically, a non-vanishing spectral weight originates as a result of scatterings. The operator j a in eq. (2.5) couples directly to two fields, the Higgs doublet and active lepton doublets, and if no further particles are involved we call the process a 1 ↔ 2 scattering (cf. fig. 1). Such scatterings give no contribution in the massless limit, because there is no phase space for the on-shell process. In the presence of thermal masses and Majorana masses, one of the kinematic channels may open up. If we count the masses as being of order M I ∼ m φ ∼ m ℓ ∼ gT , then this contribution, suppressed by the small masses, is parametrically of the same order as that of 2 ↔ 2 scatterings (cf. fig. 2). Given that the masses are small and that thermal momenta are of order k ∼ πT , the computation of 1 ↔ 2 scatterings can be simplified by considering ultrarelativistic kinematics, cf. eq. (2.2). However, there is also a complication, namely that soft scatterings which do not modify the kinematics are not suppressed and need to be resummed to all orders. This so-called Landau-Pomeranchuk-Migdal (LPM) resummation was worked out in ref. [28]. We need to generalize these results, because in ref. [28] a sum was taken over the two helicity states, and because only the diagonal elements of the density matrix were considered, and because the lepton chemical potentials were set to zero. 5 It is not too difficult to generalize the results of ref. [28] to apply to the spectral function. Adopting the notation in sec. 3.1 of ref. [20] and noting that all thermal masses are even in µ a , because they represent elastic scatterings through gauge exchange which can be both off fermions and off antifermions, the resummed "helicity-conserving" and "helicity-flipping" wave functions are denoted by g J and f J . The latter is a p-wave (vector) object, and it flips the helicity from that carried by Standard Model leptons. The spectral function can be expressed in terms of these as Here the chemical potentials are µ La ≡ µ a −µ Y /2 and µ H ≡ µ Y /2, where µ Y is the hypercharge chemical potential (the expression of µ Y in terms of the µ a is recalled in sec. 4).
Taking subsequently projections such as in eq. (3.1) (cf. eq. (A.22)); expanding as in eq. (2.2); and employing k as a variable instead of ω k J , whereby terms suppressed by O(M 2 I /k 2 ) are omitted, we find For the latter chiral projection in eq.  6 The Hamiltonian is also written in terms of k, where Γ(y) is given in eq. (3.3) of ref. [20], m 2 ℓ ≡ lim µ a →0 m 2 ℓ,a is given in eq. (3.24), and (the spectral function is a cut, i.e. "amplitude squared", convoluted with the appropriate distribution functions). The notation is as in fig. 1.
where m H is the physical Higgs mass. The wave functions are obtained from For the numerical solution we adopt the procedure described in ref. [29]; at T < ∼ 160 GeV, when we are in the Higgs phase, it needs to be modified as explained in ref. [20].

Hard 2 ↔ 2 scatterings
Unlike the 1 ↔ 2 processes, the 2 ↔ 2 scatterings, illustrated in fig. 2, are not phase-space suppressed. Therefore they can be computed at leading order in an expansion in M 2 I /k 2 , i.e. with massless right-handed neutrinos. Then only one helicity state contributes, 7 and we can write (cf. eq. (A.22)) For the latter chiral structure in eq. (3.1), δ τ,+ gets replaced with δ τ,− . Furthermore we can replace K J through K ≡ (k, k), whereby the result is independent of the indices I and J. The 2 ↔ 2 scatterings contain two logarithmic IR divergences, related to soft lepton exchange and scattering off soft Higgs particles, respectively. Following ref. [30], we handle the former by first carrying out a massless computation (this subsection), and subsequently correct the soft exchange contribution through an appropriate resummation (sec. 3.4). We denote the unresummed contribution from hard momentum exchange by ρ 2↔2,hard a (K). The latter divergence concerns terms proportional to µ Y and was not present in ref. [30], however an analogous procedure of "subtraction" and "correction" can be adopted (sec. 3.5). The phase space regions from which the divergences originate are illustrated in fig. 3.  Here (q, q 0 ) parametrizes the fourmomentum of the exchanged particle. There are logarithmic infrared divergences associated with the fermionic contributions Φ sf and Φ tf , from soft lepton exchange (exchanged particle has (q, q 0 ) ≈ (0, 0)) and soft Higgs scattering (exchanged particle has (q, q 0 ) ≈ (k, k) whereas external scatterer is soft). The divergences can be resummed as explained in secs. 3.4 and 3.5, respectively.
Accounting for the processes shown in fig. 2, we obtain (3.13) Here the chemical potentials read 3 , and dΩ n→m denotes the phase space integration measure. Furthermore p i ≡ |p i | denote incoming and k i ≡ |k i | outgoing momenta. The matrix elements read s u , (3.14) where s, t and u are the Mandelstam variables. Generalizing the techniques of ref. [30] and parametrizing by q ± ≡ (q 0 ± q)/2 the 4momentum of an exchanged particle, all but two of the phase space integrals can be carried out, yielding (3.16) Introducing the functions processes with bosonic and fermionic s-channel exchange lead to The corresponding t-channel contributions read

Resummation of soft t-channel lepton exchange
As already mentioned the massless matrix elements and phase space integrals lead to logarithmic IR divergences. A well-known divergence originates from fermionic t-channel exchange around (q + , q − ) ≈ (0, 0) where the integrand can be approximated as The divergence is regulated by a thermal mass, denoted by m ℓ,a , that the lepton obtains through its interactions with the Standard Model plasma: Computing the contribution of soft momenta q ⊥ ∼ m ℓ requires a Hard Thermal Loop (HTL) resummed computation [30]. Fortunately, this computation remains practically identical in the presence of a chemical potential, so we just briefly state the results. Following the presentation in sec. 4.1 of ref. [20], the resummed computation yields two separate ingredients. One is a "subtraction term" which removes the IR divergence in eq. (3.23) from the naive computation: The second ingredient is the correctly computed IR contribution. For this we obtain (3.26) In total, the 2 ↔ 2 contribution can then be expressed as

Resummation of soft s and t-channel Higgs scattering
There is another IR divergence which at leading order only affects chemical potential dependence, specifically S (τ ) defined in accordance with eq. (2.34). It originates from the fact that when expanded in µ H = µ Y /2, the bosonic distribution functions multiplying Φ sf and Φ tf in eq. (3.16) diverge as ∼ ±µ H T /(k −q 0 ) 2 at one corner of the integration range, cf. fig. 3. When the integration is defined as a principal value, most terms cancel between s and t-channel contributions. However, a small remainder is left over.
To be specific, the distribution functions appearing in the 2 ↔ 2 contributions can be expanded as in eqs. (3.6)-(3.8), whereas eqs. (3.19) The problem originates from the fact that Φ sf (0) and Φ tf (0) are not equal when approaching (q, q 0 ) = (k, k) from the s and t-channel sides, respectively. In order to cure the problem, we may subtract the divergent terms from the integrand of eq. (3.16) (denoting by ∆ terms within the curly brackets), where the function χ reads χ(k) = T l 1f (k) + T 2 k [ π 2 4 + l 2b (k) − l 2f (k)]. Subsequently, the schannel subtraction is reflected into the t-channel domain by q 0 → 2k − q 0 , q → 2k − q, whereby the logarithms and χ drop out. The remainder is integrated after noting that the t-channel integration domain in fig. 3 originates from the energy-conservation constraint δ(q 0 − k + ǫ φ k−q ), where a soft (i.e. |k − q| < ∼ gT ) Higgs energy is ǫ φ k−q ≡ (k − q) 2 + m 2 φ . Working out the t-channel integration range in the presence of m φ > 0 yields the correct IR contribution from soft Higgs scattering to S 2↔2 (+) , (3.34)

General setup
The left-hand side of eq. (2.37) contains charge densities, whereas on the right-hand sides of eqs. (2.37), (2.42) and (2.43) chemical potentials appear. In order to close the set of equations, the chemical potentials need to expressed in terms of the charge densities. To leading order the results are given by eqs. (4.3) and (4.7), whose derivations we wish to briefly review. Charge neutrality of the plasma poses a non-trivial constraint on the relation between chemical potentials and number densities. In the temperature range of interest we can to a good approximation assume the electroweak symmetry to be restored. Then charge neutrality concerns the hypercharge field. Defining a corresponding chemical potential as µ Y ≡ ig 1 B 0 , where B 0 is the hypercharge field in the imaginary-time formalism, a simple way to proceed is to first express the pressure (minus the free energy density) in terms of µ Y , µ a and µ B [31]. To leading order in Standard Model couplings, treating all particles as massless (masses are included in eq. (A.6) of ref. [20]), we obtain where χ F (m) ≡ k −2n ′ F (ω k ) and χ B (m) ≡ k −2n ′ B (ω k ) are so-called susceptibilities, with the special values χ F (0) = T 2 /6, χ B (0) = T 2 /3. Hypercharge neutrality corresponds now to ∂p/∂µ Y = 0, and the conserved charge densities are obtained as ∂p/∂µ a , ∂p/∂µ B .

T > 130 GeV
At T > 130 GeV the baryon chemical potential is eliminated through the sphaleron constraint µ B = − 1 3 a µ a , so that µ a couples to the strictly conserved quantity L a − 1 3 B. In a path integral formalism, the presence of µ a = 0 implies that the perturbative minimum lies at a non-zero value of µ Y , as determined by eq. (4.1). Minimizing with respect to µ Y we obtain and     The numerical uncertainties of these expression are about 20%, owing mostly to large O(α s ) corrections from the QCD coupling [23] and to IR sensitive effects from the Higgs [32]. The baryon asymmetry is given by [33]

T ∼ 130 GeV
When T ∼ 130 GeV, the sphaleron processes become slow. Consequently, baryon plus lepton asymmetry needs to be added as a dynamical non-equilibrium variable. The quantities for which equations of motion can be written are n a − n B /3 and n B + a n a . Coupling chemical potentials to these slow variables we can read off the original chemical potentials µ a and µ B : aμ a n a − n B 3 +μ B+L n B + a n a = a μ a +μ B+L µ a These values of µ a and µ B are inserted into eq. (4.1). The pressure is minimized with respect to µ Y like before, leading to Furthermore, taking partial derivatives of eq. (4.1) with respect toμ a andμ B+L , we obtain n a − n B /3 and n B + a n a as functions of the chemical potentials. Inverting these relations, we get

Evolution of baryon plus lepton asymmetry
Suppose that we start the evolution of the system from a high temperature, T ≫ 130 GeV, and are given some initial values of the lepton symmetries n a −n B /3, for instance n a −n B /3 = 0 ∀a.
To solve the evolution equations (2.37), (2.42) and (2.43), we first need to determine the chemical potentials µ a . These can be obtained from eq. (4.3). The baryon asymmetry is known as a "side product" of the evolution, from eq. (4.4).
The situation changes when the sphaleron processes become slow. We can switch to this setting at some temperature T 0 > 130 GeV at which we know the initial values of n a − n B /3 and n B + a n a from the computation described above (note that n B + a n a is in general non-zero, and can be determined from eq. (4.4)). The corresponding chemical potentials can be determined from eq. (4.7). The other chemical potentials are obtained from eqs. (4.5) and (4.6), and can then be inserted into eqs. (2.37), (2.42) and (2.43).
Obtaining the evolution equation for n B + a n a is non-trivial, given that the fluctuations of n a − n B /3 and n B + a n a are correlated, as exemplified by eq. (4.7). However the starting point is again an operator equation of motion analogous to eqs. (2.6) and (2.7). This time it takes the form of the anomaly equation, where we have introduced n G = 3 as the number of Standard Model generations; the factor 2 accounts for baryons and leptons; and c is the topological charge density. In principle this operator could be inserted into eq. (2.9), but it is not easy to express the first order time evolution of the density matrix in a useful way [33]. However, we can assume that to leading order inμ B+L andμ a , the topological charge density is only correlated with itself. Moreover, a general argument concerning correlated fluctuations [23] shows that we can writė where X a are general slowly evolving charges. By Ξ we have denoted a susceptibility matrix; its inverse multiplied by the charges yields the corresponding chemical potential. Specifically, Ξ cb = ∂ 2 p/∂μ c ∂μ b and Ξ −1 (B+L)b (X b /V ) =μ B+L . It remains to compute the symmetric correlator in eq. (5.2) for the operator on the righthand side of eq. (5.1). We denote x c(t, x) ≡Ṅ CS (t) where in the classical limit N CS is the Chern-Simons number. It is conventional to shift the time interval to run between zero and positive times; making use of the time-reversal symmetry of the anticommutator, and arguing furthermore that the dynamics is dominated by classical configurations which show linearly growing diffusive behaviour at large times, we can write 9 Here the infinite-volume limit is implicitly understood. The quantity in eq. (5.3) is precisely the one estimated with classical lattice gauge theory simulations in ref. [17]. 9 Somewhat more precisely, limt→∞ To summarize, recalling the factor 2n G from eq. (5.1) and the factor 1/(2T ) from eq. (5.2), the evolution equation for baryon plus lepton asymmetry obtains the simple form Hereμ B+L is a linear combination of all slowly evolving charges, as given by eq. (4.7).

Summary and outlook
In this paper we have presented a "field-theoretic" derivation of evolution equations of a coupled system consisting of lepton asymmetries, the baryon asymmetry, and a sterile neutrino density matrix. The basic equations are (2.37), (2.42), (2.43), and (5.4). Numerical values of the coefficients parametrizing these equations can be found in sec. 3.6 and in ref. [17]. On the right-hand sides of the equations, various chemical potentials appear; to close the system, the chemical potentials need to be expressed in terms of the lepton and baryon asymmetries, which can be achieved as re-iterated in secs. 4.2 and 4.3.
Prior to our work, many studies have appeared in which similar evolution equations have been derived (cf. e.g. refs. [4][5][6][7][8][9][10][11][12][13][14][15] and references therein). The main novelties of our investigation are the full inclusion of both helicity-flipping and conserving contributions (or, in the language of sec. 2.6, fermion-number conserving and violating effects); the inclusion of all chemical potentials and gauge field expectation values induced by them; a consistent leadingorder computation of all coefficients parametrizing the equations, both in the "symmetric" and in the "Higgs" phase; as well as a formulation general enough to permit for the treatment of the regime in which the sphaleron processes gradually switch off. We have also gone beyond linear response theory in the treatment of the sterile neutrino density matrix, permitting for both its small and large deviations from equilibrium. Even though we do not expect any of these improvements to change the previous results by orders of magnitude, many of them may play a role if a numerical precision at or below the 20% level is desired.
A numerical solution of the evolution equations within the background of an expanding universe, with a sphaleron rate [17] and equation of state [18,19] inserted from lattice studies, poses a non-trivial technical challenge, to which we hope to return in the near future.
For some qualitative insight, consider the coefficients producing or equilibrating sterile neutrinos, eqs. (2.46) and (2.47). At leading order in chemical potentials, the coefficients are determined by Q + and Q − , respectively. Here Q + contains a sum over helicity-flipping and conserving contributions, and Q − their difference, cf. eq. (2.41). The part parametrized by Q + generates a helicity-symmetric density matrix, and Q − a helicity asymmetry. Both yield a parametrically similar contribution to lepton asymmetry, cf. eqs. (2.37), (2.39) and (2.40). These effects are present even in the massless limit when helicity-conserving (fermion-number violating) contributions are absent. In the massless limit the total lepton asymmetry equals minus the total helicity asymmetry integrated over momenta, cf. sec. 2.6.
As a final comment, we note that in a recent paper [25] a coupled set of evolution equations was derived, within linear response theory, for the spin-averaged phase space distribution of one sterile neutrino species and for the total lepton asymmetry. Conceptually, this situation can be obtained from our framework by making two of the sterile neutrinos heavy so that they represent "fast variables"; integrating them out; averaging over the helicity components of the light sterile neutrino; and restricting to leading order in deviations from equilibrium. In practice we cannot proceed to that limit because different approximations are needed for treating fast and slow variables. Nevertheless, it would be interesting to understand whether analogues of the (small) "non-factorizable" contributions of O(h 4 ) that were found in ref. [25] could originate in our system, if our derivation were extended up to the O(h 4 ) level.
where a k and a † k are annihilation and creation operators, and j, j † are currents with which they interact. In the interaction picture (denoted with the subscript I), a k and a † k evolve with time: The time dependences here correspond to those in eq. (2.3). Like in the discussion below eq. (2.5), we now go over to the Heisenberg picture (denoted with the subscript H). Then As explained just below eq. (A.4), the operators in the latter terms can be identified asρ kj andρ ik , respectively, up to corrections of O(h). The ensemble averages of j and j † can be identified as advanced, retarded, and Wightman correlators: where we made use of time-translation invariance. Thereby At this point we approximate ω i ≈ ω j ≡ ω ≫ |ω i −ω j | within the Fourier transforms, whereby lim t→∞ t 0 dt ′ e iω (t ′ −t) Π < (t ′ − t) + e iω (t−t ′ ) Π < (t − t ′ ) = Π < (ω) , (A.14) The function Π < = 2n B ρ ω is real, whereas iΠ R,A have both a real and an imaginary part: iΠ R = i Re Π R − ρ ω , iΠ A = i Re Π R + ρ ω . The real parts (proportional to the spectral function, denoted here by ρ ω in distinction to the density matrix ρ) yield the absorptive effects discussed in the main text. Focussing now on the dispersive imaginary parts and carrying out a substitution like in eq. This matrix represents the "standard" energy correction for the system of eq. (A.1). Its generalization to the case of a Majorana fermion emerges through the first term in the dispersion relation in eq. (A. 19) and ultimately leads to eq. (A.25).

A.2. Modified dispersion relation for ultrarelativistic sterile neutrinos
Returning to the full system, consider the structure leading to the last term on the first row of eq. (2.29) as an example. 10 Before restricting to the absorptive part, this term reads The integral over t ′ can be carried out by making use of eq. (2.26). Subsequently, we are faced with a spectral representation which can be identified as the real and imaginary parts of the retarded correlator: The latter term leads to the absorptive behaviour in eq. (2.29), and we now focus on the first term. The retarded correlator is an analytic continuation of the Euclidean correlator, which for the operators in eq. At finite temperature the sum-integral is proportional to two independent Lorentz-tensors, / K and γ 0 . After the analytic continuation k n → −i(k 0 +i0 + ), with K = (k n , k) and K = (k 0 , k), so that Π E (K) → Π R (K), we can write For the opposite chiral projections, the roles of the helicity states are exchanged. In any case, for k ≫ M L , only the contribution proportional to β is needed. 10 Thanks to its diagonal structure the first term on the first row, containing the Fermi distribution, cancels against a contribution from the corresponding term on the second row, once we work up to leading order in the ultrarelativistic approximation ω k I ≈ ω k J . This is the same phenomenon which rendered the first row of eq. (A.13) into the purely real Π < of eq. (A.14).