Neutrino dynamics below the electroweak crossover

We estimate the thermal masses and damping rates of active (m<eV) and sterile (M ~ GeV) neutrinos with thermal momenta k ~ 3T at temperatures below the electroweak crossover (5 GeV<T<160 GeV). These quantities fix the equilibration or"washout"rates of Standard Model lepton number densities. Sterile neutrinos interact via direct scatterings mediated by Yukawa couplings, and via their overlap with active neutrinos. Including all leading-order reactions we find that the washout rate generally exceeds the Hubble rate for 5 GeV<T<30 GeV. Therefore it is challenging to generate a large lepton asymmetry facilitating dark matter computations operating at T<5 GeV, whereas the generation of a baryon asymmetry at T>130 GeV remains an option. Our differential rates are tabulated in a form suitable for studies of specific scenarios with given neutrino Yukawa matrices.


Introduction
There has been recent interest in the cosmological role that right-handed (sterile) neutrinos with masses in the GeV range could play. The dynamics of such particles may lead to a generation of a lepton asymmetry, which could explain the observed baryon asymmetry of the universe [1,2]. The lepton asymmetry generation could continue after it cannot be converted into a baryon asymmetry any more (T < ∼ 130 GeV); if the resulting lepton asymmetry is several orders of magnitude larger than the baryon asymmetry, it might be resonantly converted into keV scale sterile neutrino dark matter at T < ∼ 1 GeV [3]. Computations of lepton asymmetry generation have been carried out for two almost degenerate generations of sterile neutrinos [4][5][6][7][8][9][10], as well as for three possibly less degenerate generations [11,12]. If the masses fall below the 0.1 GeV range (but above the 0.1 MeV range), the production is so efficient that such particles should already have been observed through the total energy density that they carry [13]. If the masses are higher but still relatively close to this lower bound, roughly 0.5...5 GeV (m K ...m B ), these particles can be searched for with the planned SHiP experiment at CERN [14]. Therefore, there is a need to refine the theoretical understanding concerning the behaviour of such particles within the ultrarelativistic plasma that filled the early universe.
Computing reliably the lepton asymmetry generated in a specific scenario is extremely challenging, because both CP violation and complicated plasma physics play a role. At the same time it is relatively straightforward to establish a constraint on whether any lepton asymmetry would be washed out after its generation (cf. e.g. refs. [15,16] for reviews). This is in analogy with attempts to explain baryon asymmetry through an electroweak phase transition: the washout constraint ∆ φ † φ > ∼ T 2 /2, where φ denotes a Higgs doublet and ∆ a discontinuity across the transition temperature, already rules out a large class of theories, including the Standard Model which has no genuine transition. The goal of the present paper is to establish similar washout constraints for generating a large lepton asymmetry.
More precisely, we discuss a number of quantities that characterize the behaviour of GeV scale right-handed neutrinos within a plasma. One is called their "production rate": it tells how fast the particles are being produced if their initial density is much below the equilibrium value. Another is their "equilibration rate": it tells how fast the particles can adjust their number density if they are initially close to equilibrium but the temperature of the equilibrium ensemble evolves, as is the case in the early universe. Despite being conceptually different, it turns out that these quantities can be related to each other [17]. The main quantity that we consider is the lepton number "washout rate", which tells how fast any initial lepton asymmetry is being depleted in the presence of right-handed neutrinos. This rate can also be related to the rates mentioned above [18]. In the following, we frequently refer to the "production rate", with the relations to the other quantities specified in sec. 2. mentioning. For small masses (M ≪ πT ) the production rate has been computed at low temperatures T < 5 GeV both for vanishing [19] and non-vanishing [20,21] lepton asymmetries. At high temperatures T > 160 GeV it has been computed in the non-relativistic M ≫ πT [22][23][24], relativistic M ∼ πT [25,26], and ultrarelativistic M ≪ πT regimes [27,28]. Computations in the ultrarelativistic regime are challenging, because they require a nested resummation of the loop expansion in order to generate a consistent weak-coupling series. An interpolation applicable for any M/(πT ) has also been suggested for T > 160 GeV [29].
The status as outlined above means that there is a gap in our understanding in the range 5 GeV < ∼ T < ∼ 160 GeV. A rough estimate was presented in appendix A of ref. [4], however this was not based on a controlled computation but just included Born level 1 ↔ 2 processes with vacuum masses evolved through a changing Higgs expectation value. The main goal of the present paper is to fill the gap 5 GeV < ∼ T < ∼ 160 GeV.
The paper is organized as follows. The basic observables considered are defined and the structures of the corresponding results are outlined in sec. 2. In sec. 3 we discuss righthanded neutrino production through direct 1 ↔ 2 processes, as well as the so-called LPM resummation of the 1 + n ↔ 2 + n reactions that contribute at the same order in the soft regime M, m W ≪ πT . In sec. 4 the production through direct 2 ↔ 2 scatterings is considered. Sec. 5 is devoted to "indirect" production, via an overlap with left-handed neutrinos, pointing out the importance of 2 ↔ 2 scatterings mediated by soft gauge boson exchange at high temperatures and of 1 → 2 gauge boson decays at low temperatures. Numerical results are collected in sec. 6, and we conclude in sec. 7. A number of technical details and remarks concerning NLO effects are relegated to five appendices.

Summary of the setup and main results
One of the physical quantities that we are interested in is the production rate of right-handed neutrinos. The produced right-handed neutrinos have a momentum k ≡ |k| and a mass M . The corresponding on-shell four-momentum is denoted by K = (k 0 , k), where k 0 = √ k 2 + M 2 . In light of the scenario relevant for SHiP [14], we consider here small masses, M ≤ 16 GeV, and high temperatures, T ≥ 5 GeV. Then the right-handed neutrinos can be considered to be "ultrarelativistic", with momenta k ∼ πT and masses M ≪ πT .
Let h Ia be a neutrino Yukawa coupling, a ∈ {1, 2, 3} a left-handed lepton generation index, and I ∈ {1, 2, 3} a right-handed neutrino generation index, defined in a basis in which the Majorana mass matrix is real and diagonal: N I h Iaφ † a L ℓ a +l a a Rφ h * Ia N I .
(2.1) Hereφ = iσ 2 φ * is a Higgs doublet; a L , a R are chiral projectors; and ℓ a = (ν e) T a is a left-handed lepton doublet. For notational simplicity we normally suppress the generation index in a Majorana mass, i.e. M I → M . We consider time scales large enough that all Standard Model (SM) degrees of freedom are in thermal equilibrium. However quantities whose interactions involve the h Ia , notably right-handed neutrino phase space distributions (≡ f Ik ) and Standard Model lepton densities (≡ n a ), can be out of equilibrium. The task is to determine the equilibration rates of these observables. We note that n a can be carried both by neutral and by charged leptons.
For T ≥ 5 GeV, all Standard Model leptons can be considered degenerate and massless (πT ≫ m τ ≈ 1.8 GeV). Then the dynamical information concerning the rates of interest is contained in the 2-point function of the operator to which the right-handed neutrinos couple [17]. For computational convenience we first define the corresponding imaginary-time correlator, where K ≡ (k n , k) and k n is a fermionic Matsubara frequency. Moreover, X ≡ (τ, x) denotes a Euclidean space-time coordinate and ... T an equilibrium expectation value. The retarded correlator Π R can be expressed as an analytic continuation of Π E as and the rate observables are proportional to the spectral function ρ ≡ Im Π R . We choose a normalization for the phase distribution function f Ik such that the total number density of right-handed neutrinos, summed over the two spin states, reads Denoting k 0 ≡ E I for the mass M I , it can be shown that [17] f where the right-hand side was expanded to leading order in small lepton densities,ḟ refers to a covariant time derivative in an expanding phase space background, and n F denotes the Fermi distribution (similarly, n B denotes a Bose distribution). 1 The coefficient γ Ik can be called the (spin-averaged) "equilibration rate", and is given by This relation applies to all orders in Standard Model couplings. Normally, when referring to the right-handed neutrino "production rate", it is assumed that their number density is small, f Ik ≪ n F . For this case eqs. (2.4) and (2.5) imply thaṫ The same processes by which right-handed neutrinos equilibrate or are produced also violate lepton densities carried by Standard Model particles. Because lepton numbers are violated, their equilibrium values vanish. Close to equilibrium, the lepton densities evolve aṡ where the matrix of decay coefficients, or "washout rates", can be written as [18] Here Ξ ab = ∂n a /∂µ b | µ b =0 ∼ T 2 is a susceptibility matrix related to lepton densities. It was determined up to next-to-next-to-leading order (NNLO) in Standard Model couplings at T > ∼ 160 GeV in ref. [30], and leading-order results valid for T < ∼ 160 GeV are given in appendix A. We note that Ξ is non-diagonal, because the plasma as a whole is charge neutral, so that changes in the number densities of different lepton flavours are correlated. As is clear from eqs. (2.6), (2.7) and (2.9), the dynamical information that we need is contained in the function Im Π R , obtained from eq. (2.3). We now turn to its determination.
In order to carry out a theoretically consistent computation, power-counting rules need to be established for the various scales appearing in the problem. We denote by h t the renormalized top Yukawa coupling; by N c ≡ 3 the number of colours; by g 1 , g 2 the hypercharge and weak gauge couplings; and by λ the Higgs self-coupling. The notation g 2 refers generically to the couplings g 2 1 , g 2 2 , h 2 t , λ which are taken to be parametrically of the same order of magnitude, and "small" in the sense that g 2 ≪ π 2 .
Suppose that we are at a temperature T < 160 GeV so that, in gauge-fixed perturbation theory, the neutral component of the Higgs field has an expectation value. The expectation value is denoted by v; at T = 0, v ≃ 246 GeV. We mainly consider a regime in which v < ∼ T , even though the case m W > ∼ πT , i.e. v > ∼ πT /g, is considered as well. For v < ∼ T vacuum masses ∼ gv are of the same order as thermal masses ∼ gT but much smaller than typical momenta k ∼ πT . In other words, all particles can be considered to be ultrarelativistic. Based on various numerical tests, this regime is numerically applicable in a rather broad temperature range, 30 GeV < ∼ T < ∼ 160 GeV . Examples of 2 → 2 processes for the generation of left-handed neutrinos which subsequently oscillate into right-handed ones. The notation is as in fig. 1. The complete set for case (a) is shown in fig. 1 of ref. [29] and for case (b) in fig. 7 below.
At lower temperatures, Higgs and gauge bosons become non-relativistic and need to be decoupled from the computation (the top quark becomes non-relativistic already at a somewhat higher temperature).
In the regime of eq. (2.10), there are two types of contributions to Im Π R . First, the Higgs fieldφ in eq. (2.2) can represent a propagating mode (Goldstone or Higgs). This leads to the same processes as have previously been considered in the symmetric phase [27,28]; examples of 1 + n ↔ 2 + n processes are shown in fig. 1(a) and of 2 ↔ 2 processes in fig. 2(a). Second, the Higgs field could be replaced by its expectation value,φ ≃ (v 0) T / √ 2. Then we are left to consider processes experienced by an active (left-handed) neutrino. Examples of amplitudes are illustrated in figs. 1(b) and 2(b). We refer to first type as a "direct" contribution and to the second as an "indirect" one.
When amplitudes such as those in figs. 1 and 2 are squared, there are no interference terms between the direct and indirect sets, provided that we adopt a class of gauges (such as the R ξ gauge) in which scalar and gauge fields do not transform to each other. Then the rate can be written as where the "direct" processes are like in sets (a) of figs. 1 and 2. Like in the symmetric phase [27,28], the direct term has the parametric magnitude Im Π R | direct ∼ g 2 T 2 (recalling that |h Ia | 2 has been factored out). In contrast the indirect term has a more complicated structure (cf. sec. 5.1), where m ℓ is the active neutrino thermal mass in the ultrarelativistic regime (cf. eq. (3.2)), 2 and Γ is its thermal width (cf. eq. (5.7)). The term in eq. (2.12) is proportional to v 2 because it originates from processes induced by electroweak symmetry breaking, and to M 2 because in the massless limit helicity and fermion number conservation would forbid transitions between left and right-handed states. We note in passing that, in an alternative language, the combination |h } originating from eq. (2.12) can be interpreted as a medium-modified mixing angle squared. This weights the part of the interaction rate k 0 Γ of the weak eigenstates that is transmitted to the sterile mass eigenstates. The interaction rate k 0 Γ also appears in the denominator of the effective mixing angle, thereby contributing towards the "unitarity" of the conversion process. Nevertheless, because M 2 and m 2 ℓ can cancel against each other, the medium-modified mixing angle can be much larger than the vacuum one. Now, an essential ingredient in our analysis is the determination of the active neutrino interaction rate Γ. We find that, because of strong infrared enhancement, Γ ∼ g 2 T /π in the regime of eq. (2.10), cf. eq. (5.33). For k 0 ∼ k ∼ πT we thus get (2.13) This implies that for M ∼ gT we get Im Π R | indirect ∼ v 2 ∼ T 2 , i.e. the indirect production dominates over the direct one. The direct production dominates only if we go to the symmetric phase In order to consolidate these findings, we proceed to discuss Im Π R | direct and Im Π R | indirect . We start from the former, considering the 1 + n ↔ 2 + n and 2 → 2 contributions in turn, and return to the indirect processes in sec. 5.

LPM resummation in the symmetric phase
An essential ingredient in the physics of the processes illustrated in fig. 1(a) is the proper inclusion of the so-called Landau-Pomeranchuk-Migdal (LPM) resummation. Because of phase space suppression, the rate of the 1 ↔ 2 process (after factoring out |h Ia | 2 ) is ∼ m 2 ∼ g 2 T 2 , where m denotes vacuum or thermal masses. On the other hand, adding gauge scatterings to the 1 ↔ 2 result leads to no further suppression, because the exchanged gauge boson is soft, with a virtuality ∼ g 2 T 2 . Therefore all the gauge scatterings need to be resummed in order to obtain the correct leading-order result.
Starting with the symmetric phase, the basic equations for the LPM resummation can be summarized as follows [27]. We definê where ∇ ⊥ is a two-dimensional gradient operating in directions orthogonal to k, and the thermal masses of hard particles (with k ≫ m) read where m H ≈ 125 GeV is the physical Higgs mass. Soft gauge scatterings are represented by a thermal width which reads where d 1 ≡ 1, d 2 ≡ 3, and K 0 is a modified Bessel function. The Debye masses associated with the hypercharge and SU(2) gauge fields are defined as Here n S ≡ 1 is the number of Higgs doublets and n G ≡ 3 the number of fermion generations. The Debye masses appear frequently in the remainder of this paper. The Hamiltonian plays a role in the inhomogeneous equations From the solutions of these equations, the LPM-resummed contribution to Im Π R reads Im ∇ ⊥ · f (y) .

LPM resummation in the broken phase
In the broken phase, the scalar sector splits up into Higgs and Goldstone modes. The contribution of the Goldstone modes depends strongly on the gauge choice; at tree-level, it is straightforward to verify that both the "direct" and "indirect" contributions are gaugedependent, but their sum is gauge-independent. Once LPM resummation is incorporated, it is complicated to carry out computations in a general gauge, because this implies the presence of many different masses and correspondingly a large matrix of gauge and scalar states mixed by gauge interactions. In the following we restrict ourselves to the Feynman R ξ gauge, which minimizes the number of different states and masses. In this gauge, the Goldstone modes correspond to the physical W ± and Z 0 bosons, and we denote With non-degenerate scalar masses, the Green's functions in eq. (3.5) split up into several components, g 0 , ..., g 3 , and similarly for f . The LPM-resummed 1 ↔ 2 contribution can be expressed as a generalization of eq. (3.6), Here the contributions from µ = 1 and µ = 2 are equal, given that the charged scalar fields φ 1 and φ 2 are degenerate. The task now is to determine the HamiltonianĤ for this situation. As a first step, we introduce notation for defining the gauge field propagators in the broken phase. Because the temporal gauge field components get thermal masses, given by eq. (3.4), the temporal and spatial gauge fields mix differently. In fact, the self-energies contain more structure than just thermal masses; in general the mixing is momentum-dependent. However, because of a sum rule derived in ref. [31] and a more general argument presented in ref. [32], the quantities of our interest (see below) can be reduced to (static) Matsubara zero mode propagators. Within the regime of validity of the Hard Thermal Loop (HTL) description [34,35], the static self-energies are momentum independent. Therefore mixing can be described by constant angles, separate for spatial and temporal gauge fields.
The standard weak mixing angle can be defined as (3.9) Denoting s ≡ sin θ, c ≡ cos θ and adopting a convention according to which a covariant derivative acting on the Higgs doublet reads where Tr (T a T b ) = δ ab /2, the spatial gauge field components can be diagonalized as where Q µ is the massless photon field and B µ is the hypercharge field. The mixing angle of the static temporal components is denoted byθ, and is given by Denotings ≡ sinθ,c ≡ cosθ, the zero components are diagonalized by All diagonal fields have non-zero masses because of the thermal corrections in eq. (3.4): 14) The gauge field combinations to which neutral and charged left-handed leptons couple, respectively, are In the diagonal basis the corresponding propagators become These structures appear frequently below. Now, we return to the thermal width, denoted by Γ in eq. (3.1). Choosing k to point in the x 3 -direction, nearly light-like particles couple to the gauge field components A 0 and A 3 . In a thermal plasma, the soft scatterings mediated by these components are not identical, and the final result originates from the difference of the two contributions. Because of a sum rule [31,32], the result can most simply be expressed in terms of the static Matsubara zeromode sector (q n = 0) related to these gauge potentials. In the static limit the propagators of temporal components can be expressed as in eqs. (3.15) and (3.16). We define the widths related to W ± , Z, and Z ′ exchanges as where q ⊥ ≡ d 2−2ǫ q ⊥ (2π) 2−2ǫ and q ⊥ ≡ |q ⊥ |. Dimensional regularization has been used for defining the value of an infrared divergent integral in eq. (3.19), related to soft photon exchange, even though this divergence soon drops out (cf. the discussion below eq. (3.21)). The full width matrix, in the space of neutral and charged scalars and leptons that participate in the production of right-handed neutrinos, ordered as νφ 0 , νφ 3 , eφ 1 , eφ 2 , reads (3.20) The arguments 0 and y correspond to self-energy and exchange contributions, respectively. The combination 2Γ W (0)+Γ Z (0) corresponds to the active neutrino width or interaction rate, re-derived in some more detail in sec. 5.5 (cf. eq. (5.33)).
Given that the pairs eφ 1 and eφ 2 are degenerate, we can choose one of them as a representative. Then, the matrix in eq. (3.20) can be reduced into a 3 × 3 form, A nice consequence of this reduction is that infrared divergences related to photon exchange cancel in the combination Γ Z ′ (0) − Γ Z ′ (y). To be explicit, whereas the difference in the bottom-right component takes the form With the width determined, let us generalize the Hamiltonian of eq. (3.1) to contain a diagonal mass matrix, The Green's functions g and f are generalized to 3-component vectors. With the 3 × 3 width Γ 3×3 , we can then solve eq. (3.5), and insert the result into eq. (3.8).
As a crosscheck, we note that in the symmetric phase, the parameters appearing in eqs.
Moreover all the pairs νφ 0 , νφ 3 , eφ 1 and eφ 2 become degenerate, so we can reduce the 3 × 3 matrix into a single function, Noting that lim m→0 [K 0 (my) + ln my 2 + γ E ] = 0, this agrees with eq. (3.3). For a numerical solution, we make use of the general approach of ref. [33], adapted to the problem at hand in ref. [29]. The idea is to express the solutions of the inhomogeneous equations, eq. (3.5), in terms of the solutions of the homogeneous equation which are regular at origin. Choosing the normalizations of the regular solutions as where ρ ≡ m E2 y and ℓ is an angular quantum number, we find Im Π LPM,broken Here, as before, k 0 ≡ √ k 2 + M 2 and the kinematic range M ≪ k is assumed. The numerical solution is straightforward, with a result as illustrated in fig. 3 (the solid lines at high temperatures). 3

Limit of low temperatures
Once we go deep in the broken phase, the masses m φµ defined in eq. (3.7) eventually become large, m 2 φµ /k 0 ≫ g 2 T /(8π). Then Γ 3×3 represents a small correction compared with the mass terms in eq. (3.5), and can be omitted. However, the collinear approximation m 2 φµ /k 0 ≪ k 0 that is employed in the formalism of the LPM resummation also breaks down in the same regime. In this situation the rate is given just by the 1 ↔ 2 processes, without any resummation nor kinematic approximation. The hard thermal lepton mass m ℓ can also be 3 Numerics can be sped up by realizing that the off-diagonal elements in eq. (3.21) fall off exponentially For large enough ρ one can then switch to three separate solvers for the three independent u r ℓ,µ (ρ). This is particularly advantageous for T < ∼ 60 GeV, where a large tree-level term is present, which requires integration to large values of ρ to reach the required accuracy. omitted at low temperatures. Then the result can be given in a closed form, where we have defined Actually, M ≪ m φµ so that eq. (3.32) could be simplified by setting M → 0 (cf. eq. (5.15)).
In our numerical solution we switch from the LPM resummed result of eq. (3.30) to the Born term of eq. (3.31) when the two results cross at low T , cf. fig. 3.

Ultrarelativistic regime
We now move on to discuss direct 2 → 2 scatterings, illustrated in fig. 2(a). As long as we are in the ultrarelativistic regime, m W ≪ πT , or v < ∼ T , the masses of the "real" particles participating in these processes play no practical role, because all the scatterers have hard momenta ∼ πT . Therefore, to leading order, the computation can be directly taken over from the symmetric phase [28]. The techniques we employ here are similar to those in ref. [28], except for the treatment of soft momentum transfer, where a mild modification is adopted.
In this section we describe the computation of the direct 2 → 2 scatterings in some detail, in order to prepare the ground for the generalization to the indirect case in secs. 5.4 and 5.5.
In the 2 → 2 scatterings of fig. 2(a), the particles mediating t-channel exchange can have soft momenta. However, an explicit computation shows that only the lepton exchange is so infrared (IR) sensitive that the thermal lepton mass plays a role. The computation is organized by first determining the contribution from hard momentum transfer by using massless propagators, and subsequently treating the case of soft momentum transfer more carefully.
In naive massless perturbation theory, the hard part can be written as Here dΩ n→m denotes the usual phase space integration measure with 4-momentum conserva- The three-momenta of incoming particles are denoted by p i , with p i ≡ |p i |; those of outgoing particles are k i , with k m ≡ k the right-handed neutrino momentum. The matrix elements squared read The phase space integrals can be reduced into 2-dimensional integrals as explained in ref. [28]. Different parametrizations are introduced for s and t-channel exchange (the uchannel can be transformed into the t-channel by the exchange p 1 ↔ p 2 ). Defining the Here Φ s1 refers to bosonic and Φ s2 to fermionic s-channel exchange, and Φ t2 to fermionic t-channel exchange; the notation Φ t1 is reserved for bosonic t-channel exchange that does not appear here (or rather, appears as a diagram but does not lead to non-trivial kinematic dependence). The functions appearing in eq. (4.10) are The s-channel functions Φ s1 and Φ s2 remain finite in the whole integration range, whereas the t-channel function Φ t2 has a divergence in the vicinity of the origin (q, q 0 ≪ k 0 ): This divergence is not integrable and its proper treatment requires HTL resummation [34,35], as we now explain.
Suppose that we compute Im Π R within HTL resummed perturbation theory. The HTL scalar propagator has no cut, so the HTL result has no part which would correspond to 2 → 2 scatterings with soft Higgs exchange. Therefore the Higgs can be taken to be a "hard" external particle, and its thermal mass can be omitted. This leads to 15) where the lepton spectral function is given in eqs. (B.1)-(B.3). The lepton spectral function is parametrized by the mass given in eq. (3.2), which is purely of thermal origin, so that all left-handed leptons are degenerate. Setting k = k 0 and restricting to q, q 0 ∼ m ℓ ≪ k 0 where the HTL structures play a role, the constraint δ(k 0 − q 0 − |k − q|) in eq. (4.15) leads to Simplifying the integrand with this approximation but keeping the full integration range (the reason should become clear in a moment), we get whereρ s andρ 0 are from eqs. (B.2) and (B.3) and terms of O(m ℓ /k 0 ) have been omitted. Now, the contribution from hard momentum transfer, eq. (4.10), is IR divergent because of the term in eq. (4.14). The reason for this divergence is that the computation leading to eq. (4.10) did not incorporate HTL resummation. Fortunately we can correct for this mistake a posteriori. In order to do so, we need to subtract from eq. (4.10) the "would-be" HTL contribution, which appears there in a naive perturbative form. This is obtained from eq. (4.17) by formally expanding in a weak coupling, i.e. by assuming q, q 0 ≫ m ℓ . According to eqs. (B.2) and (B.3) this yieldŝ where terms of O(m 4 ℓ ) have been omitted. Within this approximation eq. (4.17) becomes Taking note of eq. (3.2) and of the changes of integration variables 20) this agrees exactly with eq. (4.14). The philosophy thus is to subtract eq. (4.19) from the "naive" computation of eq. (4.10). Subsequently the "soft" contribution from eq. (4.17) is added in its proper form.
Let us now compute eq. (4.17) properly. The integral contains two scales, k 0 and m ℓ , and we evaluate it in the approximation m ℓ ≪ k 0 . The leading contribution originates from q, q 0 ∼ m ℓ . In order to evaluate this contribution, it is advantageous to change integration variables from q, q 0 to q ⊥ , q 0 , where where we made use of eq. (4.16). Then where we employed a sum rule derived in ref. [28]. 4 Putting everything together, we obtain This expression is IR finite and agrees with ref. [28]. Parametrically, Im Π R | direct ∼ g 2 T 2 . A numerical evaluation is shown in fig. 4 with a dashed line ("ultrarel. 2 ↔ 2"). 4 Within O m 4 ℓ /k 2 0 accuracy the argument of the logarithm can be simplified, cf. eq. (4.23). If however the result is evaluated numerically for small k < ∼ gT where it is not leading-order correct but represents an extrapolation, it is advantageous to employ eq. (4.22) in order to avoid spurious negative expressions. We have adopted this recipe for our numerics.

Limit of low temperatures
All the 2 ↔ 2 scattering reactions depicted in fig. 2(a), leading to eqs. (4.2)-(4.4), involve a particle in the initial state whose contribution becomes exponentially suppressed when m W > ∼ πT . Therefore, the contribution of eq. (4.23) rapidly switches off once we exit the regime of eq. (2.10). Because of the resummations that were needed for obtaining eq. (4.23) it is non-trivial to obtain a general expression which has the correct high and low-temperature limits and is a smooth function in between. However, we can easily determine the lowtemperature limit. The formally dominant contribution originates from Higgs mediated bottom quark scatterings. Accounting for these through a Fermi type computation and making use of the same notation as in eqs. (4.5)-(4.10), we obtain Here 03 is the bottom quark Yukawa coupling. The result has been illustrated with a dotted line in fig. 4 ("Fermi 2 ↔ 2"). However, in practice this contribution is so small within its range of validity (T < ∼ 30 GeV) that it can be omitted.
At the same time, it is important to switch off the massless 2 ↔ 2 scatterings in this regime. We have done this by multiplying Im Π R | direct, 2 → 2 of eq. (4.23) by a phenomenological factor κ(m W ). For this we choose a "susceptibility" related to W ± bosons (cf. eq. (A.5)), normalized to the massless limit: ( 4.27) This has been included in the numerical "ultrarel. 2 ↔ 2" results shown in fig. 4. In principle it would be interesting to carry out a consistent computation for this regime, however in practice this is not necessary because, as we will see, the indirect contribution dominates the full result by many orders of magnitude at T < ∼ 30 GeV.

General structure
We now proceed to the indirect contributions, illustrated in figs. 1(b) and 2(b). As a first step, let us justify the form of eq. (2.12). According to eq. (2.2), the indirect contribution reads It is advantageous to resum the neutrino propagator to all orders, so that the production rate remains finite even for small virtualities. The inverse neutrino propagator is of the form If we make use of the property / Σ (−K) = − / Σ (K), valid in a CP-symmetric plasma, 6 then These results apply in the Feynman R ξ gauge. In a general gauge the Goldstone mode part changes; the gauge dependence cancels against similar 2 ↔ 2 indirect contributions, of the type discussed in sec. 5.6. 6 We reiterature that, like in eq. (2.5), we work close to equilibrium, with vanishing lepton asymmetries.
After the analytic continuation in eq. (2.3), we write Σ → Re Σ + i Im Σ. Then a few steps lead to 7 It is clear from the eq. (5.4) that an essential role in the indirect production is played by the real and imaginary parts of the (retarded) active neutrino self-energy. In the regime in which weak gauge bosons are ultrarelativistic (m Z ≪ πT ), the self-energy has a Hard Thermal Loop form corresponding to eqs. (B.2) and (B.3), where m 2 ℓ is given by eq. (3.2) and L ≡ 1 2k ln k 0 +k k 0 −k . Then we get for the real part When m Z > ∼ πT , the result changes; up to a normalization, 2K · Re Σ is then referred to as a finite-temperature matter potential, whose structure is reviewed in sec. 5.2. As far as the imaginary part of the active neutrino self-energy goes, we can write where a Lorentz-violating term proportional to the four-velocity of the heat bath u ≡ (1, 0) has been singled out. Then lim Subsequently the final expression for Im Π R , valid for M ≪ k, takes the form in eq. (2.12).

Finite-temperature matter potential
An important role in the indirect contribution discussed above is played by the real part of the active neutrino self-energy, 2K · Re Σ, cf. eq. (5.4), which up to a normalization is also called the finite-temperature matter potential. We review here its general form in the phase in which the Higgs mechanism is operative. We work in a regime in which the right-handed neutrinos and all active leptons are ultrarelativistic, M, m τ ≪ πT . Then 2K · Re Σ is a function of two dimensionless ratios, k/(πT ) and m G /(πT ), where m G refers to weak gauge boson masses, i.e. G ∈ {W, Z}. Defining a straightforward computation yields where n F and n B are the Fermi and Bose distributions, respectively. This is a limit of the results in ref. [36]. At high temperatures the potential can be approximated as which directly leads to eq. (5.6) [37]. At low temperatures we get which corresponds to the result in ref. [38] (cf. also refs. [39,40]). For a numerical evaluation, the q-integral in eq. (5.10) can be divided into two ranges as Then numerical evaluation poses no problems; the result is illustrated in fig. 5.

Interaction rate from 1 + n ↔ 2 + n scatterings
We now move on to Γ as defined by eq. (5.8), frequently called the active neutrino width or damping or interaction rate. The 1 + n ↔ 2 + n contributions to Γ are illustrated in fig. 1(b), and correspond physically to the decays W ± , Z 0 , γ → ℓN . From the kinematics point of view these processes are similar to those appearing in the direct contribution of sec. 3.2, if we simply replace the Goldstone modes by gauge fields; the difference is that the cubic coupling is now g rather than h. Thus, in the Feynman R ξ gauge in which gauge field propagators are similar to scalar propagators and no additional structures appear, 8 the parametric magnitude of these processes is δk 0 Γ ∼ g 2 m 2 ∼ g 4 T 2 . This turns out to be of NNLO compared with the contribution from 2 → 2 scatterings, and is therefore negligible at high temperatures. At low temperatures, when m W > ∼ k 0 ∼ πT , there is no need for resummation, cf. sec. 3.3. 9 Then the relevant 1 ↔ 2 processes are the decays of the W ± and Z 0 gauge bosons. We can 8 In other gauges the longitudinal combination QµQν/m 2 W appears in the gauge field propagator. This leads to large O(π 2 /g 2 ) effects for q ∼ πT, mW ∼ gT and violates our power counting setup. 9 Resummation becomes important when the ultrarelativistic 1 ↔ 2 and the full 1 + n ↔ 2 + n LPM lines depart from each other in fig. 3, i.e. T > ∼ 60 GeV. the active neutrino self-energy, or "finite-temperature matter potential", 2K · Re Σ/T 2 , from eq. (5.9). The high-temperature limit is given by eq. (5.6), and the low-temperature limit corresponds to the Fermi model, cf. eq. (5.12). Right: minus the sterile neutrino mass squared in the same units. The crossing of the two curves implies a "resonant" conversion from active to sterile neutrinos, however the resonance is parametrically fairly broad because of a large width k 0 Γ, cf. eq. (2.13) and fig. 6.
write a Born rate like in eq. (3.31), where F is from eq. (3.32). It is appropriate to remark that Γ is gauge independent only on the mass-shell of active neutrinos, i.e. M → 0, in accordance with eq. (5.8). Thereby we obtain .
The contribution of eq. (5.14) in this limit is illustrated in fig. 6 ("Born 1 ↔ 2"), and it represents the dominant process for T < ∼ 30 GeV.

Interaction rate from 2 ↔ 2 scatterings with hard momentum transfer
We now turn to the 2 → 2 contribution to Γ. Proceeding first with Feynman diagrams, the result can be written in a form analogous to the direct contribution in eq. (4.1): eq. (5.14) ("Born 1 ↔ 2"), the Fermi model result for 2 ↔ 2 scatterings from eq. (5.34) ("Fermi 2 ↔ 2"), and the soft 2 ↔ 2 scattering contribution from eq. (5.33) ("soft 2 ↔ 2"). The total result has been obtained by taking the smaller between the Fermi and the soft 2 ↔ 2 scattering results, which limits both to their ranges of applicability, and adding to it the Born 1 ↔ 2 rate. On the right, the total rate is shown for a number of momenta.
The corresponding diagrams are shown in fig. 7. In the massless limit (this will be rectified below), we obtain In order to simplify the last equation we have symmetrized the integrand in p 1 ↔ p 2 and made use of the identity u 2 /(st) + t 2 /(su) + s 2 /(ut) = 3. If the phase space integrals were In analogy with eq. (4.10), the phase space can be reduced into a 2-dimensional one: (5.20) The most important case Ξ t1 corresponds to bosonic t-channel exchange. The functions appearing in eq. (5.20) are given in appendix C. The s-channel functions Ξ s1 and Ξ s2 remain finite in the whole integration range. In contrast, the t-channel functions Ξ t1 and Ξ t2 have non-integrable divergences at q, q 0 ≪ k 0 : Here all terms that need to be subtracted in order for the integrals to be finite have been shown; in eq. (5.21) one more order is needed, because the multiplier of Ξ t1 in eq. (5.20) contains the divergent factor n B (q 0 ). It may be noted that eq. (5.22) has precisely the same type of logarithmic divergence as eq. (4.14). In contrast, eq. (5.21) leads to a power-divergent integral. This indicates that our naive estimate concerning the magnitude of k 0 Γ is incorrect; in fact soft gauge scatterings boost the width, so that its correct magnitude is k 0 Γ ∼ g 4 T 4 /(g 2 T 2 ) ∼ g 2 T 2 . We now turn to the determination of this IR contribution. 5.5. Interaction rate from 2 ↔ 2 scatterings with soft momentum transfer The particle mediating soft t-channel exchange in fig. 7 can be either a gauge boson or a lepton. However only one particle can be soft at a time. We start by briefly discussing the simpler case that the exchanged particle is a lepton, because the analysis is then analogous to that in eqs. (4.15)-(4.22), however in the end this contribution will turn out to be parametrically subdominant (it amounts to an NNLO contribution).
When the lepton is soft, the gauge boson is hard. A hard gauge boson can be treated like a free massless particle in the symmetric phase. Then the HTL contribution looks much like in eqs. (4.15) and (4.17), We again treat this in two different ways. A subtraction term is obtained by expanding like in eq. (4.18), and yields k 0 Γ| HTL-lepton, expanded This agrees with the divergence originating from eq. (5.22). The philosophy is to subtract eq. (5.25) from eq. (5.20), and add the corresponding "soft" contribution from eq. (5.24) in its proper form. Following eq. (4.22), we readily obtain where terms of O(m ℓ /k 0 ) were omitted. With eq. (5.25) subtracted and eq. (5.26) added, the indirect contribution from the function Ξ t2 to k 0 Γ| 2 → 2 is finite and of O(g 4 T 2 ).
If, in contrast, the lepton is hard, it can be treated like a free massless particle in the symmetric phase. The exchanged gauge boson needs now to be HTL resummed. We get 10 where (for k = k 0 ) The HTL spectral functions ρ T2 , ρ E2 , ρ TZ , ρ EZ are given in appendix B. We again treat the soft contribution in eq. (5.27) in two different ways. A subtraction term is obtained by going to the symmetric phase, like in the computation based on eqs. (5.17)- (5.19), and expanding like in eq. (4.18). For q, q 0 ≫ gT this yields (5.29) Thereby k 0 Γ| HTL-gauge, expanded Inserting the masses from eq. (3.4) this agrees exactly with the leading 1/q 4 divergence as shown in eq. (5.21). The philosophy is thus to subtract eq. (5.30) from eq. (5.20), and add the corresponding "soft" contribution from eq. (5.27) in its proper form. The soft contribution originates from q, q 0 ∼ m Ei ≪ k 0 . Changing variables from q to q ⊥ (cf. eq. (5.28)) the leading contribution from eq. (5.27) (in an expansion in O(m 2 Ei /k 2 0 )) becomes 10 Vertex corrections are of higher order and can be omitted.
At this point we can change the order of integrations like in eq. (4.22) and make use of a sum rule derived in refs. [31,32], (5.32) This structure corresponds to that in the Matsubara zero-mode sector, and equals the integrand in eq. (3.17). Similarly, the Z channel case leads to the integrand in eq. (3.18). The integral over q ⊥ can now be carried out, yielding We note that eq. (5.33) is of O(g 2 T 2 ) and is finite in the broken phase. In the notation of eq. (3.20), Γ| HTL−gauge, soft 2→2 = 2Γ W (0) + Γ Z (0). Eq. (5.33) is among our main results. Given that eq. (5.33) represents an IR sensitive result, being dominated by momentum transfer of O(gT ), it could experience large radiative corrections. In fact, as is typical of observables determined by the Debye scale, these are only suppressed by O(g/π). We have not computed these NLO corrections, but a way to do this is outlined in appendix D. Let us stress again that, in contrast, the contribution from hard momentum transfer is of O(g 4 T 2 ), i.e. NNLO, once the IR sensitive parts have been subtracted and treated properly.

Limit of low temperatures
At low temperatures, m W > ∼ gT , the mass ratios appearing in eq. (5.33) behave as m 2 W /m 2 W ≈ 1 + m 2 E2 /m 2 W , so that the interaction rate decreases as k 0 Γ ∼ k 0 g 4 T 3 /m 2 W . This is not the correct low-temperature asymptotics, however; the approximations made break down when m W ≫ gT , and the correct form is k 0 Γ ∼ k 2 0 g 4 T 4 /m 4 W . This asymptotics can be computed within the Fermi model. The corresponding results have been tabulated, up to T ∼ 10 GeV, on the web page related to ref. [19]. For completeness we specify here the result for 5 GeV < ∼ T < ∼ 30 GeV in a form which is easily amenable to a numerical evaluation.
Adding the contribution of the bottom quark to the processes listed in ref. [19]; going to the limit M → 0 in which the result is gauge independent; and making use of the same variables as in eqs. (4.5)-(4.10), we obtain The functions appearing read where G F is the Fermi constant, and Here s = sin θ is the weak mixing angle, defined like below eq. (3.9), and |V ub | 2 and |V cb | 2 have been omitted as vanishingly small. Numerically, the integral evaluates to k 0 Γ ≈ 10G 2 F T 4 k 2 0 . The result is illustrated in fig. 6 ("Fermi 2 ↔ 2").

Numerical results
The contributions of secs. 3 (direct 1 + n ↔ 2 + n scatterings), 4 (direct 2 ↔ 2 scatterings), and 5 (indirect contributions), are collected together into fig. 8. It is immediately clear that the indirect contribution dominates at low temperatures by several orders of magnitude. However, the smaller M is, the sooner does the direct contribution take over, because the indirect rate is proportional to M 2 , cf. eq. (2.12). The total rate, obtained by summing together the direct and indirect contributions from fig. 8, is shown in fig. 9. The k-dependence is illustrated in more detail for M = 0.5 GeV and M = 2.0 GeV in fig. 10. 11 Various physical quantities can be obtained by weighting these rates appropriately and integrating over the spectrum (cf. eqs. (2.6), (2.7) and (2.9)). Because of the appearance of the Fermi distribution in these weights, the phase space is dominated by k ∼ 3T . Two examples of physically relevant quantities are discussed in the next section.

Conclusions
For T < 160 GeV so that the Higgs mechanism is operative, the equilibration rate of a right-handed neutrino of mass M (cf. eq. (2.6)) can be split into "direct" and "indirect" contributions (cf. eq. (2.11)). These correspond to different types of scatterings as illustrated in figs. 1 and 2. In the ultrarelativistic regime, where all particle masses are ≪ πT , the indirect contribution can in turn be expressed in terms of the left-handed neutrino "asymptotic thermal mass", m ℓ , and "interaction rate", Γ, as indicated by eq. (2.12). At lower temperatures the general structure remains intact even though m ℓ gets replaced by a more complicated (momentum-dependent) function, as has been reviewed in sec. 5.2. We have shown that in the regime T > ∼ 40 GeV, the active neutrino interaction rate Γ is dominated by t-channel scatterings mediated by soft gauge boson exchange (referred to as the "soft 2 ↔ 2" contribution in fig. 6). In this situation Γ is "large", Γ ∼ g 2 T /π. The explicit expression is fairly simple, cf. eq. (5.33). This large contribution originates from contributions sensitive to momenta ∼ gT which would be quadratically infrared divergent without the appropriate HTL resummation. There is also a subleading (linear) infrared divergence in eq. (5.21) whose origin can also be understood (cf. appendix D).
For the masses M ∼ 0.5...2.0 GeV, relevant for the SHiP experiment [14], the right-handed neutrino equilibration rate peaks at low temperatures, T ∼ 5...30 GeV (cf. fig. 9). In this regime the active neutrino interaction rate Γ is dominated by 1 → 2 decays (cf. fig. 6) and Im Π R is dominated by the indirect contribution (cf. fig. 8). It is again possible to express the dominant contribution to Γ in a simple analytic form, cf. eqs. (5.14) and (5.15).
In order to illustrate the physics significance of these results, let us first compare the righthanded neutrino equilibration rate γ Ik from eq. (2.6) with the Hubble rate H = 8πe/(3m 2 Pl ), where e is the energy density of the universe and m Pl is the Planck mass. For simplicity we consider a seesaw scenario with hierarchical neutrinos, and assume that only one neutrino Yukawa coupling contributes to a given mass difference. Then active neutrino mass differences are of the form |∆m| = |h Ia | 2 v 2 /(2M ), whereby we can eliminate |h Ia | 2 from γ Ik to get Inserting e as tabulated in ref. [41] (cf. also ref. [42]), the result is illustrated in fig. 11(left). eq. (7.1), for k = 3T . Active neutrino masses correspond to |∆m| sol ≈ 8.7 × 10 −3 eV; for the atmospheric neutrino value |∆m| atm ≈ 4.9 × 10 −2 eV the rate is faster by a factor ∼ 5.6. Right: The (diagonal) lepton number washout rate compared with the Hubble rate, cf. eq. (2.8).
We conclude that in the mass range M ∼ 0.5...16 GeV right-handed neutrinos do equilibrate at temperatures above T = 5 GeV. Increasing the mass above 4 GeV decreases the peak equilibration rate but broadens the temperature range in which the rate is substantial. Turning to our main observable, the lepton number washout rate from eq. (2.8), the flavourdiagonal part of the result is shown in fig. 11(right). The flavour non-diagonal rate is an order of magnitude slower because of the smaller inverse susceptibility, cf. fig. 12(right). The flavour-diagonal rate exceeds the Hubble rate for all masses considered. However we note that this equilibration dynamics rapidly switches off in the range T < ∼ 4 GeV in which dark matter computations have been carried out [19][20][21].
The results of fig. 11(right) indicate that leptogenesis based on right-handed neutrinos with few GeV masses remains an interesting possibility, because these degrees of freedom do not equilibrate at T > ∼ 130 GeV when sphaleron processes are active [43]. In contrast it is difficult to generate a large lepton asymmetry for low temperatures, which could boost dark matter production in the scenario of ref. [3], because at T < ∼ 30 GeV lepton number violating reactions are in equilibrium and therefore an efficient washout process takes place. It should be acknowledged, however, that we have not performed a detailed phenomenological scan of the whole parameter space, so the existence of fine-tuned regions where the window may remain open cannot be excluded. The numerical results tabulated as explained in footnote 11 should hopefully permit for further work to be carried out in this direction. and c and d determining its inverse (cf. eq. (A.2)). The apparent discontinuity at T = 130 GeV originates from the fact that the B + L violating rate was assumed to be in equilibrium at T > 130 GeV and out of equilibrium at T < 130 GeV. Treating this regime precisely requires solving a non-equilibrium problem with a finite B +L violating rate [43], however in practice it may be sufficient to solve separate non-equilibrium problems on both sides and connect the solutions continuously.
∂n a /∂µ b takes the form Its inverse, playing a role in eq. (2.9), reads .
The functions a, ..., d are plotted in fig. 12. We now give details concerning the computation. The computation proceeds by assigning chemical potentials µ a to the different lepton densities; the chemical potential µ q to the quark number density; and by denoting the zero components of the gauge potentials by In the regime T > 130 GeV, the quark chemical potential is eliminated through the sphaleron constraint µ q = − 1 9 a µ a , whereas for T < 130 GeV it is expressed in terms of a baryon chemical potential as µ q = µ B /3; subsequently we make a Legendre transform to an ensemble with a fixed baryon density n B and then set n B → 0 in comparison with lepton densities. The gauge potentials µ Y and µ A are eliminated by requiring charge neutrality, ∂p/∂µ Y = ∂p/∂µ A = 0, where p = −Ω/V is the pressure, Ω is the grand canonical potential, and V is the volume. Subsequently, Ξ ab = ∂ 2 p ∂µa∂µ b . The contributions of various particle species to p are parametrized by the susceptibilities where K 2 is a modified Bessel function. At leading order 12 we obtain 12 We stress that the gauge boson contribution is treated consistently only in the regimes mW ≪ πT and mW ≫ πT . For mW ∼ gT the mass dependence amounts to a correction of O(g) which is not correctly represented by this expression. For mW ∼ πT the susceptibilities are parametrically of the same order as the unknown O(g 2 ) corrections and do not constitute any theoretically well-defined subset thereof.
where v 2 ≃ −m 2 φ /λ is the thermal Higgs expectation value and µ denotes generically all chemical potentials. It can be checked that for µ q , µ a → 0 and taking the temperature to be larger than all masses, this expression reproduces the Debye masses in eq. (3.4).
As far as the coefficients in eqs. (A.1) and (A.2) are concerned, for the regime T > 130 GeV a straightforward minimization leads to The other linear combinations have more complicated expressions, for instance where the χ's are combinations of susceptibilities coupling to different chemical potentials: Going to the symmetric phase, it can be checked that the resulting expressions for c and d agree with the leading-order results given in ref. [18].
(A. 18) Eq. (A.7) and the form of eq. (A.8) remain unchanged. At low temperatures the expressions agree with those given in ref. [20]. Numerical results are shown in fig. 12. The numerical uncertainties of the leading-order susceptibilities have been discussed in refs. [18,30]. Because of effects of the QCD gauge coupling on quark number susceptibilities, and because of infrared sensitive bosonic effects only suppressed by O(g), uncertainties are expected to be on the ∼ 20% level.

Appendix B. Hard Thermal Loop resummed leptons and gauge bosons
For completeness we list in this appendix the Hard Thermal Loop resummed [34,35] spectral functions corresponding to the lepton and gauge field propagators in the regime where the masses of these particles, including the contribution from the Higgs mechanism, are parametrically at most of O(gT ).
The lepton spectral function, defined as a four-vector originating from the imaginary part of the retarded propagator, has the form ρ ℓ (q 0 , q) ≡ q 0ρ0 (q 0 , q), qρ s (q 0 , q) , (B.1) where L ≡ 1 2q ln q 0 +q q 0 −q and q 0 has a small positive imaginary part, so that Im L = −π/(2q) for q > q 0 . The "asymptotic" thermal mass appearing in these equations is given in eq. (3.2).
The gauge field spectral function also contains two independent structures, associated with the (imaginary-time) projectors For W ± the spectral functions corresponding to these projections read ρ T2 (q 0 , q) = Im 1 where again q 0 has a small positive imaginary part, and m 2 E2 is from eq. (3.4). The Zchannel spectral function, appearing in eqs. (5.27) and (5.31), is more complicated because the self-energies lead to a different mixing angle than in vacuum. It can be expressed as , (B.8) and correspondingly for ρ EZ , where the weak mixing angles c, s have been defined around eq. (3.9). The self-energy Π T1 is like Π T2 in eq. (B.6) but with the replacement m E2 → m E1 . In the symmetric phase, i.e. m Z → 0, the Z channel spectral function simplifies into ρ TZ = s 2 ρ T1 + c 2 ρ T2 , (B.9) so that (g 2 1 + g 2 2 )ρ TZ = 2 i=1 g 2 i ρ Ti . Finally, we note that the photon and the mixed photon-Z propagators can be expressed in a form similar to eq. (B.8), where ∆ is the denominator of eq. (B.8).
Appendix C. Integrated matrix elements for indirect 2 ↔ 2 processes We list in this appendix the functions defined in eq. (5.20), obtained after carrying out all but two of the phase space integrals in eq. (5.16). Making use of the notation of eqs. (4.5)-(4.9) we get Ξ s1 = n S 4 g 4 1 + 3g 4 2 q + 2T (ln + b − ln − b ) Appendix D. Towards soft momentum transfer at next-to-leading order In sec. 5.5 we accounted for the leading divergence at q, q 0 ≪ k 0 given in eq. (5.21). Let us now show that the origin of the next-to-leading divergence can be understood as well.
If we change variables from q to q ⊥ , defined in eq. (5.28), and integrate over q 0 , then Γ can be expressed as an integral over q ⊥ . It turns out that the integrand is equivalent to the "transverse collision kernel", C(q ⊥ ), determined up to NLO in QCD in the domain q ⊥ ∼ m E in ref. [32], or the elastic scattering cross section, dΓ el /d 2 q ⊥ , determined up to O(g 4 ) for q ⊥ ≫ m E in ref. [45]. Concretely, eq. (20) of ref. [32] can for q ⊥ ≫ m E be expanded as where C F ≡ (N 2 c − 1)/(2N c ) and C A ≡ N c are Casimir factors related to the fundamental and adjoint representation, respectively. We now show that the NLO divergence in eq. (5.21) amounts precisely to the NLO term in eq. (D.1).
Let us introduce a scale Λ in the range gT ≪ Λ ≪ k 0 and consider the contribution to eq. (5.20) from momenta q ≤ Λ. In this domain q ⊥ can be approximated as in eq. (4.21), and we can change variables according to 2) The leading and NLO divergences originate from a domain where we can approximate Thereby we are left with the integrals Inserting these, the leading and NLO infrared divergences to Γ become Γ hard,expanded The non-Abelian part agrees exactly with eq. (D.1) after inserting C F = 3/4 and C A = 2; the contribution from the Higgs field has a similar structure but a different group theory factor.
We conclude that accounting properly for the subleading divergence in eq. (5.21) would require an NLO computation similar to that performed in ref. [32] but generalized to the broken phase and including the contribution of the Higgs field.