Precision study of GeV-scale resonant leptogenesis

Low-scale leptogenesis is most efficient in the limit of an extreme mass degeneracy of right-handed neutrino flavours. Two variants of this situation are of particular interest: large neutrino Yukawa couplings, which boost the prospects of experimental scrutiny, and small ones, which may lead to large lepton asymmetries surviving down to T<5 GeV. We study benchmarks of these cases within a"complete"framework which tracks both helicity states of right-handed neutrinos as well as their kinetic non-equilibrium, and includes a number of effects not accounted for previously. For two right-handed flavours with GeV-scale masses, Yukawa couplings up to $|h| \sim 0.7 \times 10^{-5}$ are found to be viable for baryogenesis, with $\Delta M/M \sim 10^{-8}$ as the optimal degeneracy. Late-time lepton asymmetries are most favourably produced with $\Delta M/M \sim 10^{-11}$. We show that the system reaches a stationary state at T<15 GeV, in which lepton asymmetries can be more than $10^3$ times larger than the baryon asymmetry, reach flavour equilibrium, and balance against helicity asymmetries.


Introduction
An extension of the Standard Model through two or three generations of right-handed neutrinos, which account for the observed active neutrino mass differences and mixings, offers for a simple explanation of the baryon asymmetry in the present universe [1]. The Euclidean Lagrangian is where h is a Yukawa matrix, M a Majorana mass matrix, ℓ L a left-handed lepton doublet, andφ = iσ 2 φ * a conjugated Higgs doublet. After a singular value decomposition and field rotation we may assume M = diag(M 1 , M 2 , M 3 ), where M I ≥ 0. In the following we focus on the minimal case that effectively only two generations (with masses M 1 , M 2 ) play a role; this is sufficient for explaining all known active neutrino properties.
In its classic implementation [1], leptogenesis assumes that M I ≫ 200 GeV, so that righthanded neutrinos become non-relativistic and fall out of equilibrium at a time when baryon number violating interactions through sphaleron processes are still in thermal equilibrium [2]. If the Majorana masses are furthermore assumed to be "hierarchical", only the lightest among them plays a substantial role in leptogenesis. This prototypical example has been studied to great detail by now, including the effect of radiative corrections (cf., e.g., refs. [3][4][5][6]). The drawback of this scenario is that it is not falsifiable: leptogenesis depends on high-energy parameters which cannot be uniquely fixed in low-energy experiments (cf., e.g., ref. [7]). Falsifiability can be boosted by making the right-handed neutrinos light, whereby they could become accessible e.g. to B-factory type experiments. The price to pay is that a certain degree of mass degeneracy is then needed. We refer to this framework [8][9][10][11][12] as "low-scale resonant leptogenesis". The near-degeneracy can be argued to be "natural" in the sense that it may originate from a slightly broken symmetry (cf., e.g., ref. [13]). The neutrino Yukawa couplings can be tuned relatively large, perhaps making the framework particularly well suited for experimental detection. The purpose of the current paper is to scrutinize the parameter space of this scenario, following many recent investigations [14][15][16][17][18][19][20][21][22][23][24][25][26][27].
Right-handed neutrino oscillations become efficient when the oscillation rate equals the Hubble rate, i.e. around the temperature where M ≡ (M 1 + M 2 )/2, ∆M ≡ M 2 − M 1 , and k osc is the comoving momentum at temperature T osc . Baryon asymmetry generation through sphaleron processes stops at T sph ∼ 130 GeV [28]. If we make ∆M very small, the dynamics relevant for baryogenesis takes place at temperatures just above T sph [23]. Electroweak crossover is at T ∼ 160 GeV [29,30], and therefore we may find ourselves on the side of the Higgs phase 1 in this situation. Our study is based on a quantum field theoretic formalism that we have developed in a series of previous papers [14,19,25], drawing upon earlier investigations [31][32][33][34][35]. The system is characterized by a number of slow equilibration rates, which are mediated by neutrino Yukawa couplings and are of magnitude ∼ |h| 2 g 2 T /π, where g 2 ≡ 4πα is a generic Standard Model coupling, as well as by a slow flavour oscillation rate, which is of magnitude ∼ |M 2 2 − M 2 1 |/k. The slow rates imply that right-handed neutrinos are neither in chemical, nor in kinetic, nor in helicity, nor in flavour equilibrium, and need to be tracked through density matrices. The equilibration rates contain both "direct" and "indirect" contributions, with the former referring to 1 ↔ 2 and 2 ↔ 2 decays or scatterings and the latter to rates experienced by off-shell left-handed neutrinos, which subsequently "oscillate" into right-handed neutrinos thanks to the presence of the Higgs mechanism at T < ∼ 160 GeV.
The plan of this paper is as follows. After reviewing the overall theoretical framework in sec. 2 and the parametrization of a charge-asymmetric ensemble in the presence of a Higgs mechanism in sec. 3, we discuss the structure of the indirect contribution in sec. 4, keeping consistently track of both helicity states. All ingredients appearing in the rate coefficients are computed in sec. 5, generalizing previous results in order to account for both helicity states and the presence of chemical potentials. The direct contributions are discussed in sec. 6, again resolving existing results to the chemical potential assignments relevant for the Higgs phase. The resulting system is solved numerically in an approximate form in sec. 7, in order to identify relevant corners of the parameter space. A more precise solution is presented in sec. 8, for a benchmark with large Yukawa couplings, and in sec. 9, for a benchmark with small ones. We conclude in sec. 10, and relegate some technical details to four appendices.

Overview of the framework
We start by summarizing the form of the master equations that were derived in ref. [19] from operator equations of motion and from arguments based on a separation of time scales. The variables considered are the yield parameters for lepton minus baryon asymmetries, Y a − Y B /3, and the helicity-symmetrized and antisymmetrized density matrices for right-handed neutrinos, ρ ± (k). The cosmological evolution is conveniently tracked through a variable x ≡ ln(T max /T ), and momentum through the co-moving variable k T ≡ k [s(T )/s(T min )] 1 3 , where s is the entropy density. The yield parameters evolve as where the first structure on the right-hand side may be called a washout term and the latter structures source terms. The trace goes over the flavour indices and n F denotes the Fermi distribution. To O(h 2 Ia , µ 2 ) the coefficients read where h Ia are Yukawas coupling a sterile neutrino of flavour I to a lepton of generation a; µ i ≡ µ i /T are rescaled chemical potentials; and rate coefficients Q, R, S (to be defined in sec. 4, cf. eqs. (4.4) and (4.5)) are normalized as Q ≡ Q/(3c 2 s H), where H is the Hubble rate and c 2 s the speed of sound squared. The superscripts ± indicate a symmetrization/antisymmetrization over helicity, and {IJ} indicates a symmetrization over flavour indices. Right-handed neutrino density matrices evolve as where v ≃ 246 GeV is the Higgs expectation value, κ ± (a)IJ and δ ± (a)IJ are given in eq. (4.8), and β ± (a) is in eq. (6.14). There is also an evolution equation for baryon plus lepton asymmetry, specified above eq. (A.12) and parametrized by the Chern-Simons diffusion rate Γ diff [28]. 2 2 Compared with refs. [19,25], we have displayed a subscript (a) in Q, R, S because these coefficients can depend non-linearly on lepton chemical potentials µ a in the broken phase; we have inserted a superscript (i) in S because a larger set of chemical potentials plays a role; and, most importantly, we have included all the mass corrections relevant for the broken phase, parametrized by the coefficients β ± , κ ± and δ ± .
To close the set of equations, the yields appearing on the left-hand side of eq. (2.1) and the chemical potentials appearing on the right-hand sides of eqs. (2.1) and (2.5) need to be related to each other. This "static" relation can be established as n i = ∂p/∂µ i , where the µ i -dependence of the pressure p is specified in sec. 3 and in more detail in appendix A.
As seen from eqs. (2.2)-(2.4) and (2.8)-(2.10), the microscopic information needed for solving the rate equations is contained in the mass corrections β, κ, δ and in the rate coefficients Q, R, S, which at high temperatures are functions of the temperature T , the momentum k, and the right-handed neutrino masses M I . At low temperatures T < ∼ 160 GeV, when we find ourselves in the Higgs phase, the coefficients become more complicated, depending also on v and on various particle masses. In the class of gauges in which the Goldstone modes and the gauge fields do not couple to each other, the coefficients can be expressed as [14] (2.11) Here the direct contributions refer to 1 + n ↔ 2 + n and to 2 ↔ 2 processes also present in the symmetric phase, whereas the indirect contributions are proportional to v 2 , and originate from the "oscillation" of left-handed (active) neutrinos into right-handed (sterile) ones. The direct contributions were derived in ref. [19], but require a modification with respect to their chemical potential dependence in the Higgs phase (cf. sec. 6). The indirect contributions require a lengthier analysis, as we need to generalize the results of ref. [14] to include dependences both on helicity and on various chemical potentials. After specifying the chemical potentials (sec. 3), we thus first turn to the indirect contributions (cf. secs. 4 and 5).

Parametrization of the asymmetric ensemble
As shown in eqs. (2.2)-(2.4) and (2.8)-(2.10), we aim to compute the coefficients entering the rate equations to leading non-trivial order in chemical potentials. Having non-zero chemical potentials at T < ∼ 160 GeV implies that the Higgs field and both neutral components of the gauge potentials develop expectation values. The Feynman rules pertinent to this situation are non-standard and somewhat subtle; moreover sign conventions can be a source of trouble. We summarize in this section the conventions and Feynman rules that are needed later on. With the density matrix where L a ≡ x [l La γ 0 ℓ La +ē Ra γ 0 e Ra ] is the lepton number for generation a, the part of the Euclidean action containing the kinetic terms for ℓ La is  Table 1. Effective chemical potentials carried by Standard Model particles in the chiral limit (µ q ≡ µ B /3). In the symmetric phase v ≪ T , we impose eq. (3.5), so only µ Y plays a role. Deep in the broken phase v ≫ T , when fermion masses and the chiral anomaly lead to rapid transitions between the two chiral states, we impose eq. (3.6), guaranteeing that both chiral states have the same chemical potential. The same applies to Goldstone modes and the corresponding gauge fields. The intermediate regime v ∼ T is more delicate and the assignments above are only suggestive (cf. the text). No chemical potential is indicated for right-handed neutrinos, which are not necessarily in chemical equilibrium.
The covariant derivative acting on ℓ La reads where B µ is the hypercharge field and σ a are the Pauli matrices. Gauge field backgrounds (we employ Euclidean conventions for B µ , A a µ ) are denoted by The resulting chemical potentials for ν La , e La and for other particles are collected in table 1. Now, in the "symmetric phase", where the Higgs mechanism is not operative, the SU L (2) gauge symmetry is intact, so within a perturbative treatment we should have In contrast, in the "broken phase", fermion masses induced by Yukawa couplings, as well as the chiral anomaly, violate chirality. If we assume that these reactions are in chemical equilibrium and that a quasiparticle description is viable, we should assign the same chemical potential to both chiral states. 3 According to table 1, this implies that In contrast a large chemical potential can be assigned to the electromagnetic field (≡ µ Q ), which means that we may write Here s ≡ sin(θ) denotes a temperature-dependent weak mixing angle (cf. eq. (B.4)).
In the following, we keep both µ A and µ Y non-zero, with the motivation of having expressions that can be extrapolated both to v ≪ T and v ≫ T . Furthermore this helps to illustrate the challenges that arise in the regime v ∼ T , |µ Z | ∼ |µ Q |. We are interested in determining rate coefficients and mass corrections up to linear order in chemical potentials. With the choice of eq. (3.6), terms linear in µ Z arise from 1-loop "tadpoles" mediated by Z 0 exchange; the corresponding value of µ Z is given in eq. (A.7).
Next, consider fluctuations around the minimum. The covariant derivative acting on the Higgs field is given by We write the (fluctuating parts of the) Higgs doublet and gauge potentials as (3.9) We also denote W − µ ≡ W + * µ and Z ′ µ ≡ (g 1 B µ − g 2 A 3 µ )/ g 2 1 + g 2 2 . Feynman gauge fixing is adopted because it simplifies the power counting relevant for the ultrarelativistic regime [14]. The gauge constraints are chosen to contain components of the background fields, (3.12) 3 Put another way, only by assigning the same chemical potential to both chiral states do we obtain simple propagators for massive particles (top, bottom, Higgs, W ± , Z 0 ). If we violate this condition, which happens in the regime v ∼ T , chemical potentials should probably be treated as "insertions" within perturbation theory, rather than being resummed into propagators. We have not undertaken this rather cumbersome treatment. At the same time the violation of eq. (3.6) induces a certain free energy cost in the landscape parametrized by v, µ Y and µ A , and this has been fully accounted for, as explained around eq. (3.17) and in appendix A.
With this gauge fixing, the quadratic part of the charged sector is where p n denotes a bosonic Matsubara frequency and P ≡ (p n , p). It is observed that with eq. (3.6) (or, more generally, to linear order in µ A + µ Y ), the gauge propagator obtains a simple form (here Σ Pδ (P ) ≡ 1): (3.14) Similarly, φ + can be assigned the chemical potential µ φ + = (µ Y − µ A )/2 as given in table 1. An analogous consideration can be carried out in the neutral sector. The mass splitting between the Higgs field h and the neutral Goldstone φ 3 complicates matters, so that the quadratic part now reads The coupling of the temporal gauge field component to the scalars disappears for µ A +µ Y = 0, and the Z propagator reads For the neutral scalar field φ 0 , a simple propagator parametrized by In order to fix the values of µ A and µ Y , we need to extremize the corresponding effective potential [36]. The effective potential equals minus the pressure. Since the chemical potentials are small compared with the temperature, only the leading non-trivial order is needed, and we can indeed treat chemical potentials as insertions. Restricting to leading order in Standard Model couplings, the result can be represented as a smooth interpolating function which has correct leading-order limits at πT ≪ m W and πT ≫ m W [14]: Here the susceptibilities are defined as where n F and n B are the Fermi and Bose distributions. The neutrino masses m νa serve as a symbolic indicator of the origin of the contribution. The relations between chemical potentials and lepton and baryon asymmetries following from eq. (3.17) are given in appendix A. Corrections, which are of O(g), have so far only been determined for the symmetric phase [37,38].

General structure of the indirect contribution
As can be inferred from eqs. (2.2)-(2.4) and (2.8)-(2.10), the rate coefficients Q are related to C-even and R, S to C-odd processes. In the symmetric phase, R, S could be determined from a Taylor expansion in chemical potentials. In contrast, the dependence on chemical potentials is non-linear in the broken phase, so we need to generalize the definitions. At O(h 2 Ia ), all rate coefficients can be related to the Euclidean 2-point correlator of the operators to which the right-handed neutrinos couple: where k n is a fermionic Matsubara frequency, X = (τ, x), andK · X = (k n − iµ a )τ + k · x. In the language of the canonical formalism, the expectation value is taken with respect to the density matrix in eq. (3.1). In perturbation theory, µ a , µ B = 0 induce expectation values for gauge field zero modes, which effectively act as additional chemical potentials (cf. sec. 3). The central object is the spectral function corresponding to eq. (4.1). It is the imaginary part of the retarded correlator Π R (K), 4 which in turn is an analytic continuation of Π E (K): The rate coefficients are obtained by taking matrix elements of ρ a (K), a R are chiral projectors, and u kτ I is an on-shell spinor for sterile flavour I in the helicity state τ = ±. As we work at O(h 2 Ia ) in neutrino Yukawa couplings and the mixing of active and sterile neutrinos was already accounted for within the reduction of the non-equilibrium problem into the correlators in eq. (4.3), the mixing can be omitted in the definition of the on-shell spinors u kτ I .
In the equations of ref. [19], another version of eq. (4. 3) also appears, in which the chiral projectors and helicity states are interchanged (a L ↔ a R , τ → −τ ) and the four-momentum is simultaneously put to −K J . In the chiral limit, flipping the helicity is compensated for by exchanging the chiral projectors, however changing the sign of K J does have an effect. Specifically, without chemical potentials the real (imaginary) part of the neutrino self-energy is odd (even) in K J → −K J , whereas a single insertion of a chemical potential reverses these properties. This implies that the substitution K J → −K J corresponds to µ → −µ, and we can write The dependence of Q, R, S on the flavour index a vanishes in the symmetric phase, where a Taylor expansion in chemical potentials is viable.
Let us now focus on the indirect contribution in the language of ref. [14], obtained by replacingφ by its vacuum expectation value where ∆ is an analytic continuation of the inverse neutrino propagator and u ≡ (1, 0) is the plasma four-velocity. 5 Suppressing chiral projectors, let us write ∆ as Assuming that the self-energy is proportional to either / K or / u [39], we find (4.10) The self-energy Σ can be parametrized as where the coefficients are defined as real. All the coefficients in eq. (4.11) are proportional to g 2 . We may expect that the coefficient a can be omitted, given that it is subleading compared with the tree-level term / K in eq. (4.9), but for completeness we keep it for the moment and verify that it indeed does not contribute in eqs. (4.8), (4.17), or (4.18).
Inserting eq. (4.11) into eq. (4.10), using the momentum K J as needed in eq. (4.3), and counting M 2 Here, in an expansion in M J /k, the coefficients read (4.14) Terms up to O(g 4 T 2 ) have been retained in β K and up to O(g 6 T 3 ) in β u ; this is because β K is weighted by a coefficient of O(g 2 T ) in eq. (4.16). 5 The real part of ∆ −1 also plays a role, leading to a "dispersive" correction as elaborated upon in appendix A of ref. [19]. Following an analysis similar to that leading to eqs. (4.17) and (4.18), we find that this amounts to the terms ∝ v 2 in eqs. (2.6) and (2.7), with where the coefficients are from eq. (4.11), and µ ± indicates a symmetrization/antisymmetrization with respect to chemical potentials. In the degenerate ( whereas in the temperature regime T > ∼ 30 GeV most relevant for us, the helicity-flipping factors " 1 2 " dominate. Equivalent mass corrections, apart from the chemical potential dependence, were obtained in ref. [18]. The matrix elements of eq. (4.12), needed in eq. (4.3), read Inserting eqs. (4.13) and (4.14) and working to leading order in g 2 , we find the helicityconserving and helicity-flipping coefficients Parametrically, the helicity-flipping rate Ω indirect (a+) is suppressed by O(g 2 ) with respect to the helicity-conserving rate Ω indirect (a−) , and does not contain the possibility of resonant enhancement (the latter observation conforms with refs. [18,40]). On the other hand, we find that in general 2kΓ K > Γ u at high temperatures, cf. sec. 5. This anticipates the situation in the symmetric phase, where Ω (a−)IJ is suppressed by ∼ M I M J /(gT ) 2 with respect to Ω (a+)IJ [19].

Determination of rate coefficients for the indirect contribution
We now turn to the determination of the coefficients b, Γ u and Γ K that are defined through eqs. (4.9) and (4.11) and that parametrize the indirect contribution to masses and rate coefficients through eqs. (4.8), (4.17) and (4.18), respectively. Three different regimes are considered (remaining always in the broken phase): "high temperatures", 500 GeV > ∼ πT ≫ m W ; "intermediate temperatures", πT ∼ m W ; and "low temperatures", 15 GeV < ∼ πT ≪ m W .
The starting point, eq. (4.9), involves a specific analytic continuation, and some care is needed for implementing it properly in the broken phase. We first note that in eq. (4.2) the combination k n − iµ a is analytically continued to −i[k 0 + i0 + ], whereas in eqs. (4.6) and (4.7) the combination −k n + iµ a is analytically continued to i[k 0 + i0 + ]. The latter can be re-interpreted as k n + iµ a analytically continued to the "advanced" frequency −i[k 0 − i0 + ], and subsequently taken with an inverted sign of k 0 . In other words, denoting by Φ(K, ...)| k n +iµ a →−i[k 0 −i0 + ] the advanced self-energy before the last sign inversion, and factoring out gauge couplings corresponding to Z 0 or W ± exchange, we can write Here ... stands for masses and chemical potentials pertinent to the channel in question. Now, there is a complication with this setup, arising because in the broken phase most particles feel a gauge field background, parametrized through µ Y and µ A via eq. (3.4), cf. table 1. Whenever possible it is very convenient to "resum" this gauge field background into the corresponding propagators. But then we must make sure that the relationship corresponding to chemical equilibrium, is respected in any 1 ↔ 2 reaction. Thus the Matsubara frequencies of the corresponding particles should readk n = k n + iµ 3 ,q n = q n + iµ 1 , andp n = p n + iµ 2 , withk n =q n +p n , and the analytic continuation needed for computing Φ(K, ...) with resummed propagators reads replacing the unresummed analytic continuation

Real part of the active neutrino self-energy
Let us first consider the real part of the advanced self-energy, parametrized by the function b in eq. (4.11). Like in eq. (5.1), there are two gauge channels, and in addition there is a term linear in µ Z , originating from a 1-loop Z 0 -boson tadpole contribution. 6 With the sign conventions of eq. (4.9), this implies that where the arguments show the masses and chemical potentials appearing in the loop, and µ Z is given in eq. (A.7). According to table 1 and the definitions in eq. (3.7), We omit the appearance of µ Z inside the function E, because this contribution is suppressed by ∼ α/π compared with the explicit appearance of µ Z in the last term of eq. (5.4).
In order to determine E, the inverse neutrino propagator of eq. (4.9) can be computed with the gauge propagators of eqs. (3.14) and (3.16). For Φ of eq.
Here ǫ 1 ≡ |k − p|, ǫ 2 ≡ p 2 + m 2 , and the chemical potentials are related by eq. (5.2). After carrying out the Matsubara sum, taking the real part of the advanced propagator, and recalling the conventions in eqs. (4.9), (4.11) and (5.4), we obtain (È ≡ principal value) Eq. (5.6) can be simplified at high and low temperatures. For πT ≫ m, we find In each structure only the leading term in an expansion in m/(πT ) is shown. The µindependent part corresponds to an "asymptotic" lepton thermal mass [39]. For πT ≪ m, The µ-independent part is equivalent to the classic result from ref. [49]. After inserting µ Z from eq. (A.7) and recalling that the Fermi constant reads agrees with the function −c as given in eq. (3.13) of ref. [34].

Widths at high temperatures: 2 ↔ 2 scatterings with soft gauge exchange
In the high-temperature regime, the determination of the active neutrino width requires a resummed computation [14], which profits from light-cone sum rules [41][42][43][44]. The leading contribution originates from scatterings mediated by Bose-enhanced soft gauge bosons. In order to determine this contribution, the gauge boson propagator needs to be Hard Thermal Loop (HTL) resummed [45][46][47][48]. Parametrically, HTL effects are important when m W < ∼ gT , i.e. v < ∼ T . As elaborated upon around eq. (3.6), the inclusion of chemical potentials is complicated in this regime beyond linear order. Nevertheless, we can show that chemical potentials are not expected to play a role at linear order, because of a general symmetry property of the soft contribution (see below).
In terms of Φ of eq. (5.1), the HTL-resummed result has the form 7 wherek n ≡ k n + iµ 3 ,p n ≡ p n + iµ 2 , and the chemical potentials are related through eq. (5.2). In Feynman gauge the gauge propagator can be expressed as where m depends on the gauge channel; the self-energies Π T,E can be found in appendix B of ref. [14] and their relevant limiting values in eqs. (B.2) and (B.3); and the projectors read . Inserting the projectors into eq. (5.9) we obtain In order to proceed, we write the resummed propagators in a spectral representation: Then the Matsubara sum can be carried out. Subsequently we take the cut and keep the channel leading to a soft contribution from momenta p 0 , p ≪ k, πT . Setting K → −K as in eq. (5.1), and taking the sign relevant for / Σ in eq. (4.9), we find where ρ free is the spectral function corresponding to 1/(P 2 + m 2 ). Next, let us denote the structures relevant for eqs. (4.17) and (4.18) by (5.14) The corresponding contributions to Φ are denoted by Φ HTL (τ ) , τ = ±. Carrying out the angular integral in eq. (5.13) and setting D → 4, we get where the energy conservation constraint |k − p| + p 0 = k in eq. (5.13) permitted us to write As a final step, we again focus on the contribution from the soft domain p, p 0 ≪ k, πT . Then we can drop terms suppressed by p/k or p 0 /k from eq. (5.15). According to eq. (5.16), we can subsequently write Furthermore, the leading-order contribution originates from the Bose-enhanced structure n B (p 0 − µ 2 ) ≈ T /(p 0 − µ 2 ) ≫ 1. The oddness of the p 0 -integrand implies that µ 2 only contributes at O(µ 2 2 ) and can be omitted. Thereby eq. (5.15) becomes where the spatial momentum is p ≈ p ⊥ + p 0 e k . The integral over p 0 can now be carried out. It is illustrative to first consider the term involving ρ free . Expressing the spectral function as a discontinuity of the "resolvent" R, we are faced with an integral of the type Noting that is actually regular near the real p 0 -axis, and definingR(p 0 ) ≡ 1/(p 2 ⊥ + m 2 ), I can be reexpressed as a complex integral, We proceed with the help of the residue theorem. There is a contribution from the pole at p 0 = 0, amounting to π/(p 2 ⊥ + m 2 ). In addition there is a contribution from arcs that can be sent to |p 0 | → ∞, yielding −π/(p 2 ⊥ + m 2 ). Summing together, ρ free does not contribute. Now, let us inspect the other terms in eq. (5.17). For those involving p 2 ⊥ /(p 2 ⊥ + p 2 0 ), the arcs at |p 0 | → ∞ do not contribute because of the additional suppression by ∼ 1/p 2 0 . On the other hand there is an additional pole at p 0 = ±ip ⊥ , but this does not contribute either, because Π T and Π E coincide for p = 0. Therefore only the pole at p 0 = 0 has an effect; this is the content of the sum rule obtained in refs. [41,42]. As recalled in appendix B, in this limit Π T vanishes and Π E is replaced by a mass parameter m 2 E . Finally, for the term −ρ E + 2ρ T in eq. (5.17), there is a contribution from both the pole at p 0 = 0 and the arcs at |p 0 | → ∞ [43,44]. As elaborated upon in appendix B, at the far-away arcs Π E vanishes and Π T is replaced by a mass parameter m 2 E /2. Combining the terms we obtain For the neutral sector, the temperature-dependent weak mixing angle needs to be evaluated in the proper momentum domain. Inserting the prefactors from eq. (5.1) and making use of the angles θ,θ,θ and the thermally modified masses m W , m Z , m Q , mW , mZ, mQ defined in appendix B, we get

Widths at intermediate temperatures: Born 1 → 2 decays
In the intermediate temperature range πT ∼ m W , no resummations are necessary and the inverse neutrino propagator is given by eq. (5.5). Its imaginary part, called the Born rate, can be expressed in terms of logarithms and dilogarithms denoted by . Parallelling the splitup in eqs. (5.1) and (5.4), we can write For πT ∼ m W , v ∼ πT /g ≫ T , so according to eq. (3.7) we can omit µ Z in comparison with µ Q , and set µ ν La → µ a , where it is understood that the results can be expanded to first order in chemical potentials. We note that in the high-temperature limit, when k ∼ πT and m 2 ≪ kT , the helicity-flipping interaction rate,Γ Born u + 2kΓ Born K , is larger than the helicity-conserving one,Γ Born u . In the low-temperature limit, both become exponentially suppressed.

Widths at low temperatures: Fermi 2 ↔ 2 scatterings and 1 → 3 decays
The low-temperature limits of Γ u and Γ K originate from 2 ↔ 2 scatterings and 1 → 3 decays among light fermions. To address these, we recall that at πT ≪ m W , weak gauge bosons can be integrated out and the physics described by the Fermi model. The four-fermion coupling is proportional to G F , and the rates of 2 ↔ 2 scatterings to G 2 F . In this situation, the advanced inverse neutrino propagator can be written as where the coefficients c i and the chemical potentials µ i are listed in table 2 (given that we are at T ≪ v, we can again set µ Z → 0). The different structures can be compactly expressed within the imaginary-time formalism: Here Σ {...} indicates a sum-integral over fermionic Matsubara momenta; Σ Pδ (P ) ≡ 1;P i ≡ (p ni + iµ i , p i ); and Euclidean conventions are used for Dirac matrices.
After carrying out the Matsubara sums, substitutingk n → −i[k 0 − i0 + ], setting K → −K, identifying the self-energy / Σ according to eq. (4.9), and taking the imaginary part, we get Here the structures are analytic continuations of eqs. (5.31)-(5.33). Restricting to those channels that are kinematically allowed in the massless limit, 8 we obtain If only the neutrino channels are open, the sum evaluates to +2. In this regime the equilibrium distribution n F (k T ) should also be replaced by n F ( k 2 T + M 2 J ).
where we have gone over to Minkowskian conventions, with p i ≡ |p i | and P i ≡ (p i , p i ). The Dirac traces D i appearing in eq. (5.35) can be easily taken: We refer to the structures containing P 1 · P 3 as "t-channel" and to those containing P 1 · P 2 as "s-channel" contributions. Through a renaming of integration variables, together with a permutation of chemical potentials, the three channels in eq. (5.35) can be transformed into the appearance of the first channel. As a final step, the phase space can be reduced to a convergent two-dimensional integral representation. For the widths in eqs. (4.17) and (4.18), we thereby obtain Here the t and s-channel integrals read where l 1f is defined in eq. (5.26) and p ± ≡ (p 0 ± p)/2. The integrands are supposed to be expanded to leading order in chemical potentials; the coefficients appearing after this expansion are collected in appendix C. At zeroth order in chemical potentials, eq. (C.2) reproduces eqs. (5.34-37) of ref. [14].

Determination of rate coefficients for the direct contribution
Let us turn to the direct contribution, which adds up to the indirect contribution according to eq. (2.11). At low temperatures, the largest indirect contribution is helicity-conserving, cf.
eq. (4.17), with the helicity-flipping channel in eq. (4.18) lacking the possibility of resonant enhancement. For the direct contribution the roles are interchanged.
In the ultrarelativistic regime m ≪ πT , 1 ↔ 2 reactions are phase-space suppressed. If m ∼ gT , this implies that 1 ↔ 2 rates are of the same order as unsuppressed 2 ↔ 2 rates. The 1 ↔ 2 processes are also substantially modified by soft higher-order scatterings, i.e. by 1 + n ↔ 2 + n processes with n ≥ 1, which therefore need to be summed to all orders, via a procedure known as Landau-Pomeranchuk-Migdal (LPM) resummation [31]. At low temperatures, when m ∼ πT , the phase-space suppression is not present, and it is sufficient to consider Born level 1 → 2 decays. In the following we consider 2 ↔ 2, resummed 1 + n ↔ 2 + n, and Born 1 ↔ 2 processes in turn.

High temperatures: 2 ↔ 2 scatterings with lepton or scalar exchange
The direct contribution from 2 ↔ 2 scatterings was originally determined in ref. [32], and subsequently resolved into helicity channels and generalized to include chemical potentials relevant for the symmetric phase in ref. [19]. Two separate resummations were needed in the presence of chemical potentials. In the broken phase, the chemical potentials and masses need to be re-adjusted, so that the results of ref. [19] change moderately.
The 2 ↔ 2 contribution originates from scatterings with hard momenta p i ∼ πT and is not phase-space suppressed. Therefore it can be evaluated in the massless limit, i.e. restricting to a term of the type β u (2k) in the language of eq. (4.16). The numerator of eq. (4.3) becomes Consequently we can write the contribution from hard momenta as where the coefficients c i and the associated chemical potentials are listed in table 3. The phase space integrals have forms analogous to eq. (5.35), and are collected in appendix D. The contribution from hard 2 ↔ 2 scatterings, eq. (6.2), needs to be resummed in two ways in order to render it IR finite. This can be implemented by subtracting the problematic terms, and subsequently adding them in a resummed form: Here "L" and "H" refer to scatterings mediated by soft lepton exchange and taking place off soft Higgs bosons, respectively. Table 3. The channels, coefficients and chemical potentials (cf. table 1) that appear in eq. (6.2). The field Z ′ , which is a linear combination of the physical Z and photon fields, is defined below eq. (3.9). When eq. (3.6) is satisfied, the chemical potentials of chiral partner states coincide.
Considering first the lepton exchange contribution, we define a thermal lepton mass as [39] The IR-sensitive contribution originates from the t and u-channel terms, c t1 Ξ t1 (+) + c u1 Ξ u1 (+) . The logarithmic divergence from small momenta can be subtracted with The resummed term, obtained by using a HTL propagator for the soft lepton, reads Turning to scatterings off soft Higgs bosons, the problem arises from expanding n B (k − p 0 − µ 2 ) in eq. (D.1), n B (k − p 0 − µ 1 ) in eq. (D.2), and n B (p 0 − k + µ 3 ) in eq. (D.3), to first order in chemical potentials, yielding ±µ i /(k − p 0 ) 2 . Then there is a logarithmic divergence from momenta p 0 ≈ k. Inserting the coefficients from table 3, the problematic terms can be subtracted with The resummed result is obtained by putting the integration domains together, whereby most terms cancel, and integrating the remainder over a domain regularized by a scalar mass. At this point, we recall the discussion around eqs. (3.15), (3.16), namely that neutral scalars cannot be treated as being in chemical equilibrium if µ Z = 0. Therefore we borrow an argument from the parametric regime v ≫ T , and impose eq. (3.6). Then The result for the symmetric phase is recovered by setting µ Q → µ Y and m W → m φ . In fact, apart from the values of running couplings, eq. (6.8) represents the only difference of the symmetric and broken phase values of the direct 2 ↔ 2 contribution. When πT ≪ m W , the 2 ↔ 2 contributions determined by using massless propagators need to switched off. We have done this by multiplying Ω 2↔2,direct (a+)IJ by a phenomenological factor κ(m W ), defined as 6.2. High temperatures: ultrarelativistic 1 + n ↔ 2 + n scatterings and decays The treatment of direct 1 + n ↔ 2 + n scatterings requires LPM resummation, a procedure that was first worked out for right-handed neutrinos in ref. [31]. Some chemical potentials were included in ref. [17]. These results were resolved into helicity channels and generalized to include all chemical potentials relevant for the symmetric phase in ref. [19]. In the broken phase, the assignment of chemical potentials and masses needs to be reconsidered. As discussed in ref. [14], the LPM contribution originates from four components of a wave function, describing different annihilation channels. We express this as where the superscript α ∈ {H, Z, W } enumerates the components. According to table 1 and the definitions in eq. (3.7), the chemical potentials read . Because of issues discussed above eq. (6.8), we impose eq. (3.6), omitting contributions from µ Z . The resummed terms read The s and p-wave functions g (α) and f (α) satisfy the matrix equations Here m ℓ is from eq. (6.4), whereas the matrix Γ 4×4 is given in eq. (3.20) of ref. [14].
Here, denoting ǫ 1 ≡ |k − p|, ǫ 2 ≡ p 2 + m 2 α and µ 3 ≡ µ 1 + µ 2 , the basic structure reads . (6.16) Carrying out the Matsubara sum, taking the imaginary part, and restricting to the kinematics K 2 < m 2 relevant for low temperatures, we find 17) where now ǫ 1 = |k + p|. After performing the angular integral, the result can be decomposed as β K / K + β u / u , so that matrix elements can be taken according to eqs. (4.15) and (4.16). Thereby we obtain for M ≪ k expressions similar to eqs. (5.28) and (5.29), except that the roles of the helicity channels have swapped places: Here the channels have the masses and chemical potentials given in eq. (6.15), and polylogarithms are defined in eqs. (5.25)-(5.26). Given that πT < ∼ m W implies v ≫ T , we make use of eq. (3.6) and set µ Z → 0 in the chemical potentials. If m 2 < K 2 , we omit this contribution.

Summary: putting everything together
In secs. 6.1-6.3 we have discussed the different "direct" contributions to the rate coefficients. Let us now specify how these are put together and then combined with the "indirect" ones. The full direct contribution to the rate coefficients in the broken phase can be expressed as Ω direct (aτ )IJ = Ω 2↔2,direct Here the 2 ↔ 2 part is from eqs. (6.3) and (6.9). The function I represents an interpolation between the two different 1 + n ↔ 2 + n computations, cf. eqs. (6.10), (6.18) and (6.19), in analogy with the procedure discussed in ref. [14]. 10 The chemical potential dependence is expanded to linear order, 10 More precisely, the LPM contribution is overtaken by the Born contribution at the smallest k/T and T /GeV when the Born contribution is smaller than the LPM one. The reason is that in these regimes the practical determination of the LPM contribution, making use of ultrarelativistic kinematics, becomes unreliable and overestimates the correct result.
This defines the coefficients Q direct (aτ ) , R direct (aτ ) and S (i)direct (aτ ) that are subsequently summed together with the indirect contributions.
As far as the indirect contributions go (cf. secs. 4 and 5), we invoke an expansion similar to eq. (6.21), but this time for the coefficients appearing in the active neutrino self-energy, cf. eq. (4.11): .

Approximate solution and overall parametric dependences
In order to gain insight on the behaviour of the equations specified in sec. 2, we first consider an approximate solution, similar to that followed in most of the literature. The idea is to assume that all components of the density matrix are in kinetic equilibrium, with ρ ± (k T ) ≡ ρ ± (x) n F (k T ), x = ln(T max /T ). If we subsequently integrate both sides of eq. (2.5) over k T , we end up with a coupled set of equations for the lepton asymmetries and the variablesρ ± , parametrized by integrals of the rate coefficients, weighted by n F (k T ) or n F (k T )[1 − n F (k T )]. 11 More precisely, the interpolation I makes use of the Fermi contribution at T ≤ m W /π and then freezes its value (in units of T ). The HTL contribution overtakes the Fermi one, once it exceeds the frozen value. In the rare case that the rapidly growing Fermi contribution is still smaller than the HTL one at T = m W /π, we continue to follow it until the two cross, and go over to the HTL one at higher temperatures.
The latter integrals can be carried out once and for all. As inputs for this we employ the values of Q, R, S, β, κ, δ Taylor-expanded to linear order in chemical potentials. It is appropriate to remark that, based on the analysis in ref. [25], it is not clear a priori whether a momentum-averaged solution can be accurate. First of all, the density matrices found in ref. [25] have the characteristic feature that they kinetically equilibrate very fast at small momenta, and remain close to their vanishing initial values at large momenta (a similar finding had been made in ref. [11]). Consequently, ρ ± 12 (k T )/n F (k T ) are peaked at around k T ∼ 0.5T in fig. 3 of ref. [25], rather than being constant. Second, the different momentum modes add up incoherently in the source terms for the lepton asymmetries, so that Y a − Y B /3 show much less oscillations than the momentum-averaged recipe suggests. In spite of these differences, we find that the momentum-averaged recipe performs reasonably well, with errors < ∼ 50% in many cases, even if differences of O(10) can also be found (cf. figs. 2 and 3).
For the numerical solution itself, we remark that the system contains a "charge" that is almost conserved at high temperatures, sometimes referred to as the "fermion number", and defined as the sum of the helicity asymmetries of right-handed neutrinos and the lepton asymmetries of the Standard Model particles. It is important to make sure that the integration algorithm respects this symmetry, as otherwise a non-zero fermion number generated inadvertently by numerical inaccuracy may have a large effect on late-time lepton asymmetries. On the other hand, physical fermion-number violating interactions do originate from the helicity-conserving coefficients Q (a−) , R (a−) , S (a−) [19], and they do become appreciable in the broken phase (cf. eq. (4.17)).
As far as the parameter values go, the neutrino Yukawa couplings are fixed as specified in ref. [25], by making use of active neutrino properties from ref. [50] as well as the Casas-Ibarra parametrization from ref. [51] (which can be generalized beyond the see-saw limit [52]). Choosing the right-handed neutrino mass to be M ∼ 1 GeV, and noting that complex phases have effects of O(1), the results depend substantially on just two quantities, the mass splitting ∆M and the Casas-Ibarra parameter Im z. The goal now is to map the viable parameter space in this plane. The viability concerns both the baryon asymmetry, Y B = 0.87(1) × 10 −10 [53], and low-scale lepton asymmetries, which we monitor through max({|Y a |}) evaluated at T = 5 GeV. Specifically we fix, following refs. [17,25], the noncritical parameters to the benchmark point Results obtained from the numerical solution of this system are shown in fig. 1. We observe that largest values of | Im z| are obtained for ∆M/M ∼ 10 −8 . Because of the largest neutrino Yukawa couplings and consequently the largest mixing angle with active neutrinos, this situation, studied in more detail in sec. 8, is ideal for the experimental search for right- handed neutrinos. On the other hand late-time lepton asymmetries can be considerably larger than the baryon asymmetry, but are obtained preferably with small values of Im z and a more extreme degeneracy around ∆M/M ∼ 10 −11 , so that leptogenesis takes place as late as possible. Such a situation is studied in more detail in sec. 9.

Accurate solution for large neutrino Yukawa couplings
Consider the filled circle from fig. 1, corresponding to ∆M ≡ 10 −8 GeV, Im z ≡ −5.3. The magnitudes of the neutrino Yukawa couplings are conveniently characterized by the square roots of the eigenvalues of the matrix h h † , which read 0.7 × 10 −5 and 0.2 × 10 −9 . Baryon asymmetry production peaks at temperatures just above the freeze-out one, T ∼ 130 GeV, so that little washout has time to take place while sphaleron transitions are active, even if the washout rate is large. Lepton asymmetries are, however, efficiently washed out once sphaleron processes have decoupled, at T < 130 GeV. Let us mention that the numerical integration of the basic equations is somewhat demanding in this case. In the language of eqs. (2.9) and (2.10), the dimensionless rate coefficients are a Re(h Ia h * J a ) Q + IJ ∼ a Im(h Ia h * J a ) Q − IJ ∼ 10 3 at T ∼ 160 GeV. Therefore a very fast equilibration process is taking place, and needs to be tracked with high accuracy, in a regime in which the rate coefficients vary rapidly [19]. We have written two independent routines for the integration, utilizing different languages and platforms, and verified that the results agree in general down to the 1...2% level.
The results from the numerical integration are shown in fig. 2, where they are also compared with the momentum-averaged treatment of sec. 7. The basic feature of this benchmark is that baryon asymmetry freezes out close to when lepton asymmetries are maximal. After the freeze-out, lepton asymmetries are rapidly erased. These qualitative features are correctly reproduced by the momentum-averaged approximation, even if momentum averaging is seen to overestimate the correct result by a factor ∼ 2.

Accurate solution for small neutrino Yukawa couplings
Finally we consider the filled square from fig. 1, corresponding to ∆M ≡ 10 −11 GeV, Im z ≡ −0.15. In this case the square roots of the eigenvalues of the matrix h h † are 4.1 × 10 −8 and 3.0 × 10 −8 . Most of the lepton asymmetry generation takes place after sphaleron processes have ceased to be active, i.e. at T < 130 GeV.
Our numerical solution is shown in fig. 3, where we have also compared with the momentum averaged treatment (grey lines). Baryon asymmetry is seen to freeze out already during an early stage of lepton asymmetry generation (left panel). In this particular case the T / GeV grey: mom.average momentum-averaged treatment is seen to underestimate the full result by a factor ∼ 7.
The most remarkable feature of our solution concerns the lepton asymmetries, which are shown in the right panel of fig. 3. We observe that lepton asymmetries obtain a constant value below T ∼ 15 GeV, which is furthermore the same in all flavours. This is the case for low-temperature lepton asymmetries in general. The existence of such a state was proposed in ref. [18], whose eq. (61) can be derived from our eq. (9.3) by summing over both active and sterile flavours, integrating over momenta, and approximating susceptibilities.
The reason for this behaviour can be understood as follows. Consider a state in which the helicity-symmetric density matrix has equilibrated, ρ + = diag(n F , n F ), and the helicityasymmetry is diagonal, ρ − = diag(ρ − 11 , ρ − 22 ). In eq. (2.1), only the first and last term play a role. In eq. (2.5), only the third and fifth term play a role. At low temperatures, the rate coefficients are dominated by the helicity-conserving components, so that Q + ∀a : Numerically, we find that the process towards the stationary state starts with the equilibration of ρ + II . Later on this is followed by ρ − II and the lepton asymmetries, cf. fig. 4. Afterwards the system remains in this state at least as long as πT > ∼ M .

Conclusions
The purpose of this paper has been to carry out a "precise" study of two carefully tuned benchmark points of GeV-scale resonant leptogenesis. By precision we mean that the rate equations and most rate coefficients have been consistently determined to complete leading order in Standard Model couplings in the parametric regime gT ≪ k, πT , where g 2 = 4πα w and k is the right-handed neutrino co-moving momentum. Due to soft thermal effects, a resummation of the naive loop expansion was necessary for achieving this goal. Based on an analysis of lepton number susceptibilities (i.e. relations of chemical potentials and lepton asymmetries), which play a role in our master equations and for which higher-order corrections have been determined [37,38], we expect the theoretical uncertainty to be on the ∼ 20% level. There is one ingredient which was not fully resolved yet, namely the effect of the "chiral" chemical potential µ Z on the "direct" rate coefficients (cf. sec. 6) in the intermediate domain v ∼ T , however we expect the numerical influence from here to be on the ∼ 1% level. In addition there are non-perturbative uncertainties which are difficult to quantify at present, such as that the non-perturbative crossover is at T ∼ 160 GeV whereas within our perturbative treatment the Higgs phenomenon sets in at T ∼ 150 GeV.
As main ingredients of our analysis, we track both helicity states of right-handed neutrinos; consider both the symmetric and broken phase of the electroweak theory; allow for kinetic non-equilibrium; and include a large set of chemical potentials (including gauge field "tadpoles"). To contrast this with extensive recent parameter scans, kinetic non-equilibrium, helicity-conserving rates, a smoothly evolving sphaleron rate, hypercharge chemical potential, as well as all indirect contributions relevant for the broken phase, were omitted in ref. [17]. In ref. [26], kinetic non-equilibrium, the term 2kΓ K in the helicity-flipping indirect contribution (cf. eq. (4.18)), the running of Standard Model couplings, as well as the chemical potential dependences of the rates B ± , D ± (cf. eqs. (2.3), (2.4), (2.9), (2.10)) and of the mass corrections β, κ, δ, were omitted. On the µ-dependence of B ± , D ± we remark that even if such effects are formally of second order in deviations from equilibrium, it may be prudent to include them, given that ρ ± can deviate from equilibrium by O(1). Nevertheless, the results in our fig. 1(left) agree semi-quantitatively with ref. [26].
The first of our benchmarks (cf. the filled circle in fig. 1, and sec. 8) concentrated on large neutrino Yukawa couplings, whereas the second (cf. the filled square in fig. 1, and sec. 9) focussed on small ones. On the methodological side, our main finding was that kinetic equilibrium, even if not justifiable theoretically, is often a reasonable approximation, even if differences of O(10) can be found (cf. figs. 2 and 3). Assuming kinetic equilibrium is attractive in that it accelerates numerics and therefore permits for overall parameter scans.
Apart from kinetic non-equilibrium, another ingredient worth elaborating upon are the mass corrections, parametrized by β ± , κ ± , δ ± in eqs. (2.6) and (2.7). Like the rate coefficients, these can originate either from "indirect" processes (the terms ∝ v 2 ) or from "direct" processes (the terms ∝ T 2 ). We find that implementing precisely the mass corrections has a very important O(10) suppressive effect on late-time lepton asymmetries (less so on Y B ).
On the physics side, our main conclusion concerns the strong interplay between helicity and lepton asymmetries. We have demonstrated that, after undergoing complicated dynamics, the system settles into a stationary state, or "fixed point", at low temperatures (cf. figs. 3,4), in which there is flavour equilibrium both in the active and sterile sectors (cf. eq. (9.3)). The temperature at which this happens lies typically in the range T ∼ 15...50 GeV. The remnant lepton asymmetries can reach values |Y a | > 10 −7 ≫ |Y B |.
The significance of this finding originates from its connection to dark matter physics [10]. Thanks to flavour equilibrium, values |Y a | > ∼ 10 −5 would be large enough to permit for resonant keV-scale sterile neutrino dark matter production at T ∼ 0.1...1.0 GeV [34], proceeding via the Shi-Fuller mechanism [54]. The existence of a stationary state suggests that the leptogenesis and dark matter processes nicely factorize from each other.
Our finding should motivate further work in this direction. Even if our results got to |Y a | > 10 −7 , they fell short of |Y a | ∼ 10 −5 for the parameters in eqs. (7.1), (7.2) (cf. figs. 1  and 3). This justifies broader parameter scans, as well as further refinements of the theoretical framework. For instance, at very low temperatures T ≪ M/π, an additional contribution to lepton asymmetries could originate from the non-equilibrium decays of the right-handed neutrinos [10,12]. Considering such contributions leads to the need to include many new mass effects (such as from m τ , m c ). Finally, it should be clarified whether the sterile neutrino helicity asymmetries that we observed (cf. eq. (9.3)) could constitute "reservoirs", which might facilitate dark matter production, thereby rendering the Shi-Fuller mechanism viable even if Standard Model lepton asymmetries remain somewhat below |Y a | ∼ 10 −5 . Unlike the existence of the stationary state itself, this seems to be a dynamical question, whose resolution depends on the values of the conversion rates |h Ia | 2 Q (a−)II (cf. eqs. (9.1), (9.2)). high temperatures, when v ≪ T and µ A = 0, we find Here µ q = µ B /3, and masses have been retained as reminders of the origins of the contributions. For obtaining n a − n B /3 and n B + a n a , we follow the procedure described in sec. 4.3 of ref. [19], writing µ a =μ a +μ B+L and µ B =μ B+L − aμa /3 (sphaleron equilibrium corresponds toμ B+L = 0, see below). Extremizing with respect to µ Y , and going to the massless limit, when χ F and χ B can be approximated according to eq. (3.18), we obtain Derivatives with respect toμ a ,μ B+L yield n a − n B /3, n B + a n a , respectively, and inverting these relations results in  Subsequently, µ a =μ a +μ B+L and µ B =μ B+L − aμa /3. At low temperatures, when v ≫ T , inserting µ A and µ Y from eq. (3.7) and omitting terms proportional to µ Z leads to  In the numerical solution, we include dependences on top, bottom, W ± , Z 0 and Higgs masses, whereby the equations become a bit more complicated.
In the intermediate regime v ∼ T , we employ the full eq. (3.17) rather than (A.1) or (A.4), and both µ A and µ Y need to be extremized simultaneously [36], which leads to a smooth interpolation between eqs. (A.3) and (A.6). The price to pay is that when neither eq. (3.5) nor eq. (3.6) is satisfied, perturbation theory becomes complicated due to the coupling of gauge and scalar modes (cf. sec. 3). To understand when we find ourselves in this situation, we note that the extremal value of µ Z is given by χ eff ≡ χ F (0) 18 − 36s 2 + 152s 4 3 + χ F (m t ) 3 − 8s 2 + 32s 4 3 (A.8) where n u,c,t ≡ i=u,c,t n i ; s 2 ≡ sin 2 (θ) withθ from eq. (B.4); and Eqs. (A.7)-(A.10) show that the assumption in eq. (3.6) is valid for v 2 ≫ χ eff ∼ T 2 . For v ≪ T , s 2 → 0 and m/T → 0, and eq. (A.7) then implies that µ Z → µ Q . Let us now turn to the implications of sphaleron equilibrium on this discussion. The sphaleron rate falls out of equilibrium in the intermediate domain v ∼ T [2], and in this regime neither eq. (A.3) nor (A.6) is accurate. As long as the sphaleron rate is fast,μ B+L re-adjusts itself to zero on a time scale much shorter than we can resolve, so eqs. (A.3) and (A.6) show that the "would-be" equilibrium state has The corresponding rate equation, viz. Y ′ B+L = a F a − γ (Y B+L − Y eq B+L ), where we employ the notation of eq. (4.2) of ref. [25], contains the coefficient (n G ≡ 3, Γ diff is from ref. [28]) (A.12) The sphaleron rate is in equilibrium when γ ≫ 1. According to eq. (A.11), a sudden switch from one limiting treatment to the other would insert a discontinuity in Y B+L if γ ≫ 1. In order to avoid this, we have derived the analogues of eqs. (A.3), (A.6), (A.11), (A.12) from an extremization of the full eq. (3.17) with respect to both µ A and µ Y . It is straightforward to verify that the resulting expressions interpolate continuously between the limiting values.
Our numerical results make use of this continuous interpolation.

C. Coefficients for the Fermi limit of the active neutrino width
We list here the values of the sums appearing in eq. (5.40) when the integrals of eqs. (5.41) and (5.42) are expanded to linear order in chemical potentials; coefficients are inserted from table 2; and eq. (3.6) is made use of (we write the result in terms ofμ B = N cμq and make use of the symmetry of Ξ s (τ ) in µ α ↔ µ β ): D. Phase space integrals for direct 2 ↔ 2 scatterings We list here the phase space integrals appearing in eq. (6.2). The associated coefficients c i and chemical potentials are listed in table 3. The five cases read , (D. 5) where u, t, s are the Mandelstam variables, the polylogarithmic functions appearing on the right-hand sides have been defined in eqs. (5.25)-(5.26), p 0 = p + + p − , p = p + − p − , and dΩ 2↔2 ≡ p 1 p 2 p 3δ The functionδ is defined such that Pδ (P) = 1. Even if not obvious from the right-hand sides of the expressions, the definitions and numerical values of Ξ t0 (+) and Ξ s0 (+) coincide.