Sterile neutrino rates for general $M$, $T$, $\mu$, $k$: review of a theoretical framework

The temperature of sphaleron freeze-out, $T\approx 130$ GeV, sets a scale for mechanisms for generating lepton asymmetries and converting them to baryon asymmetry (leptogenesis). For Majorana masses $M \gg \pi T$, leptogenesis takes place in the non-relativistic regime, while for $M \ll \pi T$, the dynamics is ultrarelativistic. The intermediate case $M \sim \pi T$ is the most cumbersome, as no effective theory methods are available. We review the definitions of and provide integral representations for all rate coefficients and mass corrections of $\mathcal{O}(h_\nu^2)$ in this regime, such that the expressions extrapolate to known limits and are gauge independent to the order computed, both in the symmetric and the Higgs phase, for any helicity, momentum, and set of chemical potentials. Road signs towards a numerical evaluation are offered.

masses, and the related phenomenon that neutrino flavours mix with each other. Even if this shortcoming can be fixed in a simple way, by adding gauge singlet (sterile) right-handed neutrino fields in the theory, much remains unknown. In the active neutrino sector, major challenges are the determination of the absolute scale of neutrino masses [1] and of a CP-violating parameter δ that manifests itself in their mixings [2]. In the sterile sector, the magnitude of the so-called Majorana mass parameters remains undetermined. In principle, they could vanish (in which case neutrinos are said to be Dirac-like), be as heavy as 10 15 GeV, or be practically anything in between (in the latter cases we talk about Majorana-like neutrinos). It would be a major breakthrough to establish the Majorana mass scale.
In concrete terms, the Lagrangian of the Standard Model completed by right-handed neutrino fields reads A possible indirect handle on physics in the neutrino sector can be offered by cosmology. For instance, already several decades ago, it was realized that cosmological considerations require the existence of three light neutrino species. Today the corresponding parameter, called N eff , is precisely measured [3] and even more precisely computed theoretically [4]. The agreement between these numbers is a great success of modern cosmology.
There are also issues in cosmology that the traditional Standard Model cannot explain. Two famous ones are the apparent existence of particle dark matter, and of an asymmetry in the amounts of matter and antimatter (the latter is referred to as baryon asymmetry). Astronomical observations have determined both parameters with percent level accuracy [3]. If we are lucky, the theoretical explanations of these features could be related to neutrino mass generation, and might allow for us to constrain the Majorana mass parameters.
The possibility that the right-handed neutrino fields, and the corresponding mass eigenstates, sterile neutrinos, could contribute to dark matter and baryon asymmetry, have been proposed long ago. Building on earlier works [5,6], it has been realized that for dark matter, one of the sterile neutrinos should be sufficiently light (M ∼ keV) and thus long-lived [7][8][9]. Even if unconfirmed, concrete demonstrations that the astronomical verification of this scenario could be possible, have been put forward [10,11].
For baryon asymmetry, processes originated from sterile neutrinos are referred to as leptogenesis [12]. Here, the mass scale could be anything from M ∼ 0.1 GeV [13,14] to M ∼ 10 15 GeV [12], with the lower bound again related to the well-constrained N eff (cf., e.g., refs. [15][16][17][18] and references therein). Close to the lower bound, it is conceivable that sterile neutrinos could be experimentally detectable in the future [19]. In the full mass range, active neutrino masses are generated through the see-saw mechanism [20][21][22], which implies that the neutrino Yukawa couplings are small for small Majorana masses, spanning the range |h ν | ∼ 10 −8...0 . If M < ∼ 10 6 GeV [23,24], some sterile neutrinos should be degenerate in mass, in order to permit for resonant enhancement of CP violation that is necessary for generating the observed baryon asymmetry [25][26][27]. The case M ∼ 100 GeV, i.e. of the order of the weak scale, has attracted detailed attention only recently (cf., e.g., refs. [28][29][30][31] and references therein), even though it could arguably be the most "natural" one.
The goal of this paper is to assemble up-to-date tools for addressing the cosmology of sterile neutrinos in broad mass and temperature ranges. In particular we define and derive rate coefficients and mass corrections that characterize their interactions, oscillations and decoherence within a cosmological (or astrophysical) setting. The exposition is somewhat theoretical in nature, as the numerical evaluation poses challenges of its own [32,33]; we return to an outlook at the end.
The presentation is organized as follows. After recalling the general form of the rate equations satisfied by right-handed neutrino density matrices and lepton asymmetries (cf. sec. 2), we review what has been known about the associated rate coefficients (cf. sec. 3). To proceed beyond this level, we discuss the Feynman rules pertinent to the problem (cf. sec. 4), given that the presence of chemical potentials leads to non-standard features. The 2 ↔ 2 and 1 ↔ 3 contributions to the rate coefficients are analyzed in sec. 5, whereas 1 + n ↔ 2 + n processes are treated in sec. 6. Thermal corrections to dispersion relations, such as thermal masses, are the subject of sec. 7. We conclude with a summary and outlook in sec. 8. Matrix elements squared relevant for the problem are listed in appendix A.

Rate equations
In order to study the cosmology of sterile neutrinos, and particularly their contribution to dark matter and/or baryon asymmetry, a suitable theoretical framework is needed. The early universe represents a multiparticle system, with a maximal temperature likely much higher than any of the known particle masses. In this situation multiple interactions take place. The traditional tool of particle physics, based on Feynman diagrams, describes individual scatterings of a few particles, and is not trivially suited to such a situation. Rather, it has to be combined with tools of statistical physics.
An important notion in statistical physics is that of thermal equilibrium. A system in thermal equilibrium is simple: it carries a minimal amount of information (entropy is maximal), which can be characterized by a handful of parameters, namely the temperature and the chemical potentials associated with conserved charges. Of course, the universe as we observe it is not in thermal equilibrium. For instance, for dark matter, equilibrium would imply that density is exponentially suppressed by ∼ (M T /π) 3/2 e −M/T , which would lead to an energy density much smaller than the observed one, if M ∼ keV and T ∼ 10 −4 eV. Baryon asymmetry cannot be generated in thermal equilibrium either, as famously stated by Sakharov [34].
Even if the complete universe cannot be in thermal equilibrium, most of the Standard Model particles were in equilibrium for most of the time [35]. The key question for explaining dark matter and baryon asymmetry is to find particle types which decoupled from equilibrium at some point, or never entered it in the first place.
A suitable framework for describing such dynamics is to classify degrees of freedom (that is, the densities or density matrices of various particle species), as either "fast" or "slow". We call degrees of freedom fast, if they experience reactions often enough to stay in thermal equilibrium. In contrast, slow variables can fall out of equilibrium.
For fast variables, it is easy to determine their density (or energy density), as it is completely characterized by the temperature and chemical potentials, with the distribution functions taking the Fermi-Dirac or Bose-Einstein form. In contrast, for slow variables we do not know the density in advance. We need to carry out a theoretical computation to find it out. The computation is based on a set of rate equations, which in turn are characterized by the expansion rate of the universe (Hubble rate) as well as the most important quantities concerning us here, thermal rate coefficients. In general, this theoretical framework can be referred to as the Landau theory of thermodynamic fluctuations [36].
Focussing on our problem, the most important slow variables are: • lepton and baryon asymmetries carried by Standard Model particles, which in the following are denoted by n a and n B , with a ∈ {e, µ, τ }; • two 3 × 3 density matrices for the sterile neutrinos, which in the following are denoted by ρ (±) , or their linear combinations, defined according to eq. (2.1). The subscripts (±) refer to two helicity states, as carried by massive spin-1/2 particles, and the requirement for 3 × 3 matrices originates from the three different generations.
The evolutions of these variables are coupled to each other. We note in passing that at very high temperatures further slow variables emerge [37], but we restrict to the simplest situation in the following, viable for T < ∼ 10 5 GeV [35]. The evolution equations for the slow variables are parametrized by a number of coefficients. One of the coefficients is a special one, namely the rate of anomalous baryon plus lepton number violation that is present in the Standard Model [38,39] (conventionally this is referred to as the "sphaleron rate"). In many leptogenesis scenarios, lepton number production continues down to temperatures where the anomalous rate rapidly switches off. Given that the anomalous rate plays an essential role in converting lepton asymmetries to the physically observable baryon asymmetry, the precise features of the switch-off do need to be incorporated [40]. Fortunately, this rate and its switch-off have been precisely determined with the help of numerical simulations [41]. This rate coefficient is also "universal", in the sense that it is independent of the properties of the sterile neutrinos. In the following, we concentrate on the other rate coefficients.
Unless the right-handed neutrinos are extremely heavy (M ≫ 10 15 GeV), the see-saw formula suggests that the neutrino Yukawa couplings are small (|h ν | 2 ≪ 1). The smallness of the couplings is the very reason why the sterile neutrinos are slow variables. In this situation, we can carry out a perturbative expansion, with the first non-trivial order being O(h 2 ν ). There are important new effects at O(h 4 ν ), notably those related to CP-violation, however such effects often factorize into a product of contributions of O(h 2 ν ). At O(h 2 ν ), the sterile neutrino rate coefficients are related to the 2-point "retarded" correlation function of the composite operator to which the sterile neutrinos couple. This correlation function is a matrix in Dirac space, and its matrix elements determine the helicity components of the rate coefficients. The general framework by which one arrives at this 2-point correlator is similar to a classic one considering active neutrino oscillations in a statistical background [42], however a more general derivation was needed, in order to properly include vacuum masses, both helicity states, as well as the fact that electroweak gauge bosons become light at high temperatures, and eventually the electroweak symmetry gets restored. Building on a number of previous steps, general derivations of the rate equations and definitions of the rate coefficients have been provided in refs. [43,44].
In order to specify the rate equations in a compact form, we consider helicity-symmetrized or antisymmetrized density matrices, If the sterile neutrinos are not degenerate, the off-diagonal components of the density matrices average out, up to effects suppressed by 1/M . Now, the evolution equations for lepton asymmetries take the forṁ where k ≡ d 3 k (2π) 3 , n F (ω) denotes a matrix with Fermi distributions on the diagonal, and ω I ≡ k 2 + M 2 I are the corresponding frequencies. The density matrices evolve aṡ where again ω is a diagonal matrix. The coefficients A and C, appearing in the terms without density matrices, are proportional to chemical potentials, cf. eqs. (2.5) and (2.7). After separating the neutrino Yukawa couplings into C-even and C-odd parts, via A + (a)II =μ a φ + (a)II Q + (a) , where we have denoted leptonic chemical potentials byμ a ≡ µ a /T . The Q's and U 's are absorptive (i.e. rate) and dispersive (i.e. mass) coefficients, respectively, capturing the dynamics of the heat bath. Specifically, like in eq. (2.1), Q ± (a) = [Q (a+) ± Q (a−) ]/2 denote symmetrization and antisymmetrization with respect to helicity, and similarly for U ± (a) . The values Q (a±) , U (a±) , and the corresponding C-odd partsQ (a±) ,Ū (a±) , are obtained from the imaginary and real parts of a retarded correlator, as defined in eqs. (4.15) and (4.16).

What is currently known and not known about rate coefficients
The dynamical parts of the rate coefficients, denoted by Q and U above, originate from matrix elements of a 2-point correlator. Their definitions will be given in eqs. (4.15) and (4.16), once the notation needed has been introduced. Given the definitions, the correlator and its matrix elements need to be evaluated. This task has only been completed in a number of special cases. The status strongly depends on the "regime" considered, i.e. the relation of the mass of a sterile neutrino (M ) and the temperature (T ). In thermal field theory, the temperature naturally appears in the combination πT . The following regimes can be identified: • non-relativistic regime, M ≫ πT . In this situation thermal corrections are small.
Actually, there has been a long discussion about how small they are. Based on leadingorder computations, one might naively think that thermal corrections are Boltzmannsuppressed, i.e. exponentially small. It turns out, however, that in general thermal corrections are only power-suppressed [45].
• relativistic regime, M ∼ πT . In this situation thermal corrections are of O(1). Then they must be determined without approximation, leading to cumbersome computations.
• ultrarelativistic regime, M ≪ πT . In this situation thermal corrections dominate the dynamics, and physics looks very different from that in vacuum. Fortunately, a lot of work has been carried out for this situation in the context of QCD, and effective field theory methods and resummations have been developed, notably an ingenious simplification relevant for light-cone correlators, originating from considerations of causality in a thermal setting [46].
Beyond order-of-magnitude estimates or leading-order computations in the non-relativistic regime (cf., e.g., refs. [47][48][49][50][51] for reviews), the state of the art of the determination of the rate coefficients in the different domains can be summarized as follows: • in the non-relativistic regime, the simplest CP-even rate coefficient, which is of O(h 2 ν ) in neutrino Yukawa couplings, has been determined up to next-to-leading order (NLO) in Standard Model couplings [52][53][54][55][56], partly by making use of the methodology of ref. [45]. Subsequently, the simplest CP-odd rate coefficient, which is of O(h 4 ν ) in neutrino Yukawa couplings, has also been determined up to NLO in Standard Model couplings in the non-relativistic regime [57][58][59][60].
• in the relativistic regime, the NLO level has been reached for the simplest CP-even rate coefficient in the symmetric phase [61,62].
• in the ultrarelativistic regime, the simplest CP-even rate coefficient has been determined to full leading order (LO) in the symmetric phase [63,64], which in this case requires a systematic resummation of the loop expansion.
• a smooth interpolation between NLO results in the non-relativistic and relativistic regimes and resummed LO results in the ultrarelativistic regime has been been worked out for the simplest CP-even rate coefficient in the symmetric phase [65].
• making use of the methods of ref. [46], the NLO level has been reached for the first rate coefficient even in the ultrarelativistic regime in the Higgs phase [71].
Despite this progress, open goals remain. A few shortcomings that need to be overcome to permit for phenomenological studies in the full parameter space, are as follows: (i) in the relativistic regime, only the sum over helicity states has been addressed at NLO; the difference between helicity-conserving and helicity-flipping rates has only been determined in the ultrarelativistic regime. This is a problem, given that helicity asymmetry has the same quantum numbers as the lepton asymmetry, and therefore plays an important role in the rate equations. (Extrapolations towards the relativistic regime, however without a full computation, have been presented in ref. [29].) (ii) in the relativistic regime, the chemical potential dependence of the rate coefficients has not been resolved at NLO. This needs to be done, given that chemical potentials are proportional to lepton asymmetries, and therefore give again a formally order unity contribution to the rate equations.
(iii) in the relativistic regime, rate coefficients have been systematically determined only in the symmetric phase at NLO. Their determination becomes much more complicated in the Higgs phase, where Standard Model masses play a role. This is a shortcoming for leptogenesis computations, given that the electroweak crossover is at T ≈ 160 GeV [72,73] whereas the sphaleron rate switches off at T ≈ 130 GeV [41], i.e. on the side of the Higgs phase. Moreover, lepton asymmetry generation can continue at T < 130 GeV [43,[74][75][76][77], and this may play an important role for sterile neutrino dark matter, given that large lepton asymmetries need to be present at T ∼ 0.2 GeV when dark matter production peaks [8,[78][79][80][81].
(iv) apart from a few test cases [70,82,83], rate coefficients have been evaluated after momentum averaging (cf., e.g., refs. [84][85][86][87][88][89][90][91] for recent works). This assumes that density matrices depend on the momentum k through the shape of the Fermi distribution, with only their overall amplitude appearing as a dynamical variable. Even if often accurate up to a factor of order unity, discrepancies of an order of magnitude are also possible [70], so the full momentum dependence of the rate coefficients needs to be established. In the relativistic regime, this has so far only been achieved for the simplest CP-even rate coefficient in the symmetric phase at NLO.
(v) the gauge independence of the NLO rate coefficients has not been demonstrated in the Higgs phase, and the possibility to separate their origin into "direct" and "indirect" contributions has been asserted [69] but not properly proven.
At this point, it is helpful to be more specific about the order that needs to be reached for phenomenologically relevant computations. Above, NLO results were mentioned for the nonrelativistic and relativistic regimes. In these regimes, the LO reactions are 1 ↔ 2 processes, such as the decays of the sterile neutrinos into Standard Model particles or, if the sterile neutrinos are light, the decays of Higgs or gauge bosons into sterile neutrinos and Standard Model leptons. NLO reactions then comprise of virtual corrections to 1 ↔ 2 processes, supplemented by new classes of real processes, notably 2 ↔ 2 and 1 ↔ 3 scatterings. Now, when we go to the ultrarelativistic regime (m i , M ≪ πT ), the picture changes. The 1 ↔ 2 processes become phase-space suppressed compared with the scaling dimension T , whereas 2 ↔ 2 reactions experience no suppression. The phase-space suppression of 1 ↔ 2 processes also implies that they are sensitive to small corrections, necessitating the summation of soft 1 + n ↔ 2 + n scatterings, with n ≥ 0. The upshot is that a full LO computation in the ultrarelativistic regime must include all 2 ↔ 2 reactions [64] as well as 1+n ↔ 2+n scatterings in a resummed form [63]. The latter goes under the name of Landau-Pomeranchuk-Migdal (LPM) resummation [92][93][94].
What is of interest to us is to obtain a result which is LO accurate in all domains, and interpolates smoothly between them [65]. It then needs to include 1 ↔ 3, 2 ↔ 2 and resummed 1 + n ↔ 2 + n scatterings, in a way that that thermal mass corrections are included where necessary, but double countings are avoided. General numerical methods for determining the 1 ↔ 3 and 2 ↔ 2 and the virtual corrections to 1 ↔ 2 processes that cancel mass singularities, have been developed in ref. [32], whereas a systematic approach to 1 + n ↔ 2 + n scatterings and the interpolation between the ultrarelativistic and relativistic regimes, has been put forward in ref. [33]. In the remainder of this paper, we turn to how the ingredients necessary for refs. [32,33] can be determined, for general M, T, µ, k.

Feynman rules for Standard Model in a statistical background
In order to compute all rates, we need to establish the Feynman rules relevant to the problem. The way to include chemical potentials cannot easily be located in literature so, for completeness, the ingredients are outlined in sec. 4.1, following ref. [70] but generalizing to an arbitrary gauge parameter. Subsequently, definitions of the dynamical coefficients Q and U , alluded to in sec. 2, are motivated in sec. 4.2.

R ξ gauge in the presence of chemical potentials
The presence of chemical potentials induces non-zero expectation values for temporal gauge field components, so we need to include these backgrounds when deriving Feynman rules. Within the imaginary-time formalism, we denote the gauge backgrounds by where B µ is a hypercharge field (cf. eq. (4.2) for sign conventions) and g 1 , g 2 are the U Y (1) and SU L (2) gauge couplings, respectively.
If we find ourselves on the side of the Higgs phase, it is useful to "diagonalize" the chemical potentials used, via a linear transformation. To this end we may write  Here µ B is the baryon chemical potential, µ Y and µ A are defined in eq. (4.1), and µ a is the chemical potential associated with the active lepton flavour a ∈ {e, µ, τ }.
with temporal gauge field components [70]. 1 After the transformation, µ Z couples to the neutral weak current, whereas µ Q is the chemical potential associated with electromagnetism.
The values of µ Y and µ A are to be found dynamically, by extremizing an effective potential V (v, µ Y , µ A ) [95]. The solution implies that when the symmetry gets restored (v ≪ T ), µ A → 0, s → 0, and µ Z → µ Q . In contrast, deep in the Higgs phase (v ≫ T ), µ Z ∼ µ Q T 2 /v 2 ≪ µ Q and µ Y + µ A ≈ 0. As discussed in sec. 3 of ref. [70], it is practical in the Higgs phase to choose µ Z = 0 in the tree-level Feynman rules, and let its small non-zero value be generated by 1loop tadpole corrections, at the same stage when 1-loop corrections to dispersion relations are derived (cf. eq. (7.6)). The tree-level chemical potentials induced by the gauge backgrounds to different Standard Model particles are listed in table 1.
As a powerful crosscheck of the computations, it is useful to carry them out in a general R ξ gauge, which needs to include gauge field backgrounds. With the conventions where σ a are the Pauli matrices, the gauge fixing term appearing in the imaginary-time action When Feynman rules and ghost properties are derived from these gauge conditions, the usual Feynman rules are recovered, except that charged Goldstone modes and ghost fields carry the same chemical potentials as the W ± gauge bosons, and the momenta appearing in vertices are shifted by chemical potentials. Defining where p n is a Matsubara frequency and µ ax a chemical potential associated with a particle of type a x , the upshot is that Feynman rules are like without chemical potentials, but with all momenta appearing in the shifted formP ax . Chemical equilibrium requires that xPax = 0 at each vertex, so momentum conservation can be used as before. It should be noted, however, that the substitution P → −P implies that the sign of a chemical potential gets inverted.

Sterile neutrino correlator and its helicity projections
At O(h 2 ν ), the coefficients parametrizing the rate equations of sec. 2 are all related to the real and imaginary parts of a particular retarded correlation function, We represent the retarded correlator as an analytic continuation of an imaginary-time one, Here a ∈ {e, µ, τ } labels an active lepton flavour, µ a is the corresponding chemical potential, For future reference, let us rewrite eq. (4.11) more explicitly, after going to the Higgs phase (v ≃ 246 GeV) and to momentum space. In a compact notation, implying a sum-integral over the four-momenta P, Q in products where they appear, this yields where Ω = X is the space-time volume. The first structure, proportional to v 2 , represents an "indirect" contribution in the language of ref. [69]. The last structure, containing no v, is a "direct" contribution, as the operator couples directly to propagating modes. The middle structure in eq. (4.12) represents a cross term between these two possibilities; the corresponding diagrams can be found in fig. 1. Helicity matrix elements are obtained from the correlator Π R a as where a L , a R are chiral projectors, and we inserted on-shell spinors in the form (4.14) Coefficients denoted by Q andQ in sec. 2 parametrize C-even and C-odd parts, respectively, of the "absorptive" part of the matrix element, normalized by the energy, whereas U andŪ parametrize its "dispersive" part, The C-even and C-odd parts can be extracted by symmetrizing and antisymmetrizing in chemical potentials, respectively. (Note that neutrino Yukawa couplings have been factored out from the Standard Model correlator.) In order to simplify eq. (4.13), it is helpful to consider the sum and the difference of the helicity states. Resolving η τ in terms of 2-component spinors, defined as eigenstates of k · σ, the sum is seen to amount to Making use of the Clifford algebra we then get where a mass term was eliminated by the chiral projectors.
It is a bit less standard to work out the helicity difference. Making use of an explicit (Dirac or Weyl) representation of the Dirac matrices, it can be verified that Inserting into eq. (4.13), and dropping again terms eliminated by chiral projectors, yields As a final step, the sum in eq. (4.18) and the difference in eq. (4.20) can be employed for determining the positive and negative helicity contributions. To summarize, if we define and work in the medium rest frame, denoting its four-velocity by U ≡ (1, 0), then the helicity matrix elements are given bȳ As an example, at leading order in the Standard Model, Π R a ⊃ 2a L / P a R C (cf. sec. 6.1). It then follows that Θ E ⊃ 4E · P C, and correspondingly the positive helicity component experiences the rateū k+ a L Π R a a R u k+ = 2(ω + k)(ǫ − p ) C, where we wrote P ≡ (ǫ, p) and p ≡ p · k/k for the component parallel to k. Similarly, the negative helicity component

Rates from 2 ↔ 2 and 1 ↔ 3 scatterings
Moving on to the evaluation of the expectation value in eq. (4.12), we start by considering the contribution given by 2 ↔ 2 and 1 ↔ 3 scatterings. This amounts to a Boltzmann-like computation. For Boltzmann equations, matrix elements squared are needed, and we first discuss how they can be obtained, following ref. [32].

Thermal crossing relations for real corrections
Let Θ E (P 1 , P 2 , P 3 ) denote the matrix element squared for a 1 → 3 decay of a sterile neutrino of four-momentum K into final-state particles with momenta P 1 , P 2 and P 3 . Here E defines a generic polarization sum, in the sense of eq. (4.21). Thermal averaging is denoted by where we have defined in order to permit for a simultaneous treatment of Bose-Einstein (σ = +) and Fermi-Dirac distribution functions (σ = −). Furthermore on-shell energies are denoted by where p a ≡ |p a | and k ≡ |k|, so that four-momenta read P a = (ǫ a , p a ) and K ≡ (ω, k). At the Born level, the full rate originating from 2 ↔ 2 and 1 ↔ 3 reactions can be expressed as Γ Born 2↔2,1↔3 = scat 1 → 3 (a 1 , a 2 , a 3 ) + scat 2 → 2 (−a 1 ; a 2 , a 3 ) + scat 2 → 2 (−a 2 ; a 3 , a 1 ) where the negative index labels indicate that the signs of the corresponding momenta are to be inverted. Crossed phase space integrals read The upshot from eq. (5.6) is that we only need to compute Θ E (P 1 , P 2 , P 3 ) explicitly, with the six other cases obtained "for free" from crossing symmetries [32]. It is essential here that the 1 → 3 decays do not need to be kinematically allowed; only the analytic structure of Θ E (P 1 , P 2 , P 3 ) is needed.

Virtual corrections and mass singularities
The poles in the matrix element squared, Θ E (P 1 , P 2 , P 3 ), and their residues, permit to determine IR-sensitive virtual corrections to 1 ↔ 2 scatterings. Generalizing on the KLN theorem [96,97], their role is to cancel logarithmic and double-logarithmic mass singularities from the real 2 ↔ 2 and 1 ↔ 3 scatterings. We do not repeat the rules for this recipe here, noting just that a computer-algebraic implementation can be found in ref. [32].

How to obtain matrix elements squared
According to secs. 5.1 and 5.2, the challenge is to extract the matrix element squared Θ E (P 1 , P 2 , P 3 ) associated with 1 → 3 decays, with other ingredients following from it through established relations. In this section we describe how to achieve this.
One way to obtain Θ E is to determine the retarded correlator of eq. (4.9), and then to extract its cut. Viewing the operators as in eq. (4.12), the Feynman diagrams contributing to the retarded correlator are those in fig. 1.
This type of computations are conveniently implemented with FORM [98]. There are, however, challenges. The result obtained after Wick contractions contains many ambiguities, related to choices of sum-integration variables. Even just to demonstrate the gauge independence of the retarded correlator, we need to make use of renamings of variables and of integration-by-parts identities, in order to eliminate the ambiguities. Even if this can be done and gauge independence verified, there is a price to pay for these reductions, namely that when the cut is extracted, the remaining integrands do not display their natural symmetries, as all symmetries have been eliminated in favour of an unambiguous representation.
In order to illustrate the situation, let us define "master" sum-integrals, general enough such that all 2-loop contributions can be determined as their linear combinations: Here Σ P ≡ T pn p is a Matsubara sum-integral, with p n referring either to a bosonic or fermionic Matsubara frequency, and P ≡ (p n , p). The inverse propagators read ∆ P ;ax ≡ (p n + iµ ax ) 2 + p 2 + m 2 ax , (5.14) where µ ax is the chemical potential associated with the particle species a x , and m ax is its mass, whereas H iP denotes a helicity projection, for instance H iP = E · iP for our problem. In this language, symmetries depend on how many of the indices i 1 , ..., i 5 are positive and whether some of the corresponding propagators carry identical particles. For instance, if i 4 ≤ 0 and a 1 = a 3 , we can make use of the substitution P → Q − P to interchange the 1st and 3rd propagator, replacing the original master by 1/2 times the sum of the original and reflected master. However, the momenta come with opposite signs after the substitution, so this is a symmetry only if the particle in question carries no chemical potential. Furthermore the numerator structures weighted by j 1 , j 3 and j 4 change, and need to be projected back to the basis. Finally, if i 4 < 0, this inverse propagator needs to be expressed in terms of the new inverse propagators, via This assumes that the chemical potentials of the different species are related to each other in a way guaranteeing chemical equilibrium (cf. the paragraph following eq. (4.8)). Now, if the numerator structures are invariant in the reflection, the appearance of a minus sign in front of ∆ −K−P ;a 4 on the right-hand side ensures that this part is eliminated when the original and reflected term are averaged over. Had we not considered this symmetrization, a redundant term, with a possibly gauge-dependent coefficient, would have remained.
The most extensive symmetrizations appear when there are just three propagators, which is the minimal case leading to 2 ↔ 2 and 1 ↔ 3 cuts. Prime examples are I j 1 ···j 6 1i 2 1i 4 1i 6 , with i 2 , i 4 ≤ 0, and I j 1 ···j 6 i 1 111 i 5 i 6 , with i 1 , i 5 ≤ 0. In these cases, 3! = 6 momentum permutations need to be implemented, in order to eliminate ambiguities and establish gauge independence.
As a bottom line for this approach, we note that it is ideal if a minimal set of master integrals are evaluated afterwards. If, however, a nice-looking intermediate expression is needed, such as our Θ E (P 1 , P 2 , P 3 ), which should display maximal symmetries, then the elimination of all redundancies is somewhat counterproductive.
Another strategy is to determine directly the matrix elements squared appearing in the decay rate, by considering the amplitudes depicted in fig. 2. For this kind of computations, many automated tools are available, though we find it simplest to again implement the algebra with FORM [98].
As a remark on the latter approach, we note that even if it works well for simple amplitudes, the more complicated cases, with many interference terms, are not identified in a nice form. Indeed, as kinematic variables are not independent, several equivalent representations can be given to any matrix element squared. In general, human inspection is required for judging which variables yield the "nicest" expression. Our results after some undoubtedly incomplete efforts in this direction are collected in appendix A.

Splitup into direct and indirect contributions
As indicated by eq. (4.12) and illustrated by figs. 1 and 2, the result for 1 → 3 decays does not get trivially separated into indirect and direct contributions, but there are interference terms between these two sets. Nevertheless, such a splitup is still possible, as we now show. 2 The idea is to first compute all the contributions in fig. 1 or 2. At this point, we have verified the gauge independence of the full expression. Then, we view the result as an expansion in the virtuality of the sterile neutrinos, M 2 . The expansion contains two singular orders, 1/M 4 and 1/M 2 . It turns out that the coefficients of these terms are not independent, but are related to each other in a specific way. This permits for us to "resum" these singular terms into an indirect contribution. More concretely, let us define where E is the external four-vector from eq. (4.21). As demonstrated by the explicit expressions in appendix A, all appearances of inverse powers of M 2 can be represented through scalar products containing I. Subsequently, we can write 3 The self-energy Σ represents the indirect contribution. The parts of Θ E containing scalar products with E (amounting to direct contributions) and I (amounting to indirect contributions) are listed in eqs. (A.10) and (A.12) for hadronic and leptonic effects, respectively.

HTL resummation
If the temperature is large compared with all masses (m i ≪ πT ), then the thermal masses generated by Standard Model interactions (δ T m i ∼ gT ) can be as large as the vacuum masses, and needed to be included for a consistent computation. The systematic way to do this goes through Hard Thermal Loop resummation [99][100][101][102].
Considering first either the side of the symmetric phase (cf. eqs. (A.9) and (A.11)), or direct contributions in the Higgs phase, i.e. effects proportional to E in eqs. (A.10) and (A.12), the terms are integrable after regularizing lepton propagators by a small mass [32]. As reviewed in sec. 5 of ref. [32], the regulator could subsequently be replaced by a physical thermal lepton mass, via an addition-subtraction step analogous to that in eq. (5.19). However, if M ∼ πT , this is not necessary at leading order, because the large energy in M guarantees that the soft phase space region gives a subdominant contribution. 4 Turning to the indirect contribution, it can be anticipated from eq. (5.18) that Im Σ can be singular if v → 0, i.e. if masses are sent to zero. Put another way, the division by v 2 implies that masses ∼ g 2 v 2 turn into couplings ∼ g 2 . Therefore, as can also be deduced from fig. 2, ω Im Σ is of O(g 4 ) in terms of couplings originating from vertices. At the same time, the phase space average is quadratically divergent at small momenta if masses are sent to zero [69]. Given that the v-independent scales characterizing the thermal average are ∼ πT and M , the divergence starts to play a role if v ≪ max{M, πT }.
If we are in the regime M ≪ πT , the would-be divergence can be cured by HTL resummation, which gives the gauge bosons thermal masses ∼ gT . The thermal masses "regulate" the propagators in the regime gv ≪ gT . For the problem at hand, leading-order HTL resummation has been worked out in sec. 5.2 of ref. [70], and the NLO level has been reached in ref. [71]. The magnitude of ω Im Σ is O(g 2 T 2 ) after the resummation.
In a practical computation, we normally need to interpolate between domains in which HTL resummation is important or not. In order to avoid double counting, HTL resummation must then be implemented through an addition-subtraction step, where the subtraction removes terms that would be part both of / Σ Born and / Σ HTL,full . Let us stress again, however, that if M > ∼ πT , HTL resummation is not necessary at leading order for the observables that we are concerned with. Figure 3: Left: 1-loop self-energy diagrams possessing a non-vanishing cut, with the same notation as in fig. 1. Right: the corresponding amplitudes, with the same notation as in fig. 2. 6. Rates from 1 + n ↔ 2 + n scatterings Finally we consider sterile neutrino rate coefficients originating from 1 + n ↔ 2 + n reactions.
Again we obtain both direct and indirect contributions. Following ref. [33], we carry out the basic discussion on the level of 1 ↔ 2 reactions (cf. sec. 6.1), noting that LPM resummation over n ≥ 0 (cf. sec. 6.2) can be implemented through a subsequent addition-subtraction step, analogous to eq. (5.19).

Splitup into direct and indirect contributions
To set the stage, we start by considering Feynman diagrams in the Higgs phase, shown in fig. 3 (the diagrams on the left are analogous to the diagrams in fig. 1 but at 1-loop level).
As was the case with the 2-loop diagrams in fig. 1, the two contributions are not gaugeindependent by themselves. Thus, gauge independence needs to be verified first, before splitting the contributions into the direct and indirect ones. Let us describe this procedure, which is similar to that introduced around eq. (5.18), in some more detail. 5 Omitting tadpole-like terms that are related to the renormalization of the tree-level contribution from the first row of eq. (4.12), the diagrams in fig. 3 where / ∆ −1 and ∆ −1 are fermion and scalar propagators, respectively; a sum-integral over P is implied; GZ, GW stand for neutral and charged Goldstone modes; and Z ′ , W ′ refer to the longitudinal gauge modes, with m 2 Z' ≡ ξm 2 Z and m 2 W ' ≡ ξm 2 W where ξ is a gauge parameter.
It is obvious from eqs. (6.1) and (6.2) that the two parts are not gauge independent separately, as the former depends on the Goldstone masses and the latter on m 2 Z' and m 2 W ' . Adding together and noting that in the tree-level vacuum, m 2 GZ = m 2 Z' and m 2 GW = m 2 W ' , gauge dependence drops out, because of opposite signs on the first rows.
How the cancellation operates is that the first row of eq. (6.2) effectively replaces the Goldstone masses in eq. (6.1) through physical gauge boson masses. The same result could be obtained by evaluating eq. (6.1) in Feynman gauge, as in that case m GZ = m Z and m GW = m W . We may call this the direct contribution.
The remaining terms in eq. (6.2) are proportional to inverse powers ofK 2 . They become large if we analytically continueK 2 → −M 2 and then make M 2 small. These terms need to be "resummed" into an indirect contribution, in analogy with eqs. (5.17) and (5.18). In total, the 1 ↔ 2 processes can thus be represented as The zeroth order term in eq. (6.4) originates from a tree-level evaluation of the first row of eq. (4.12). If we go to the symmetric phase, i.e. v → 0, then only the direct contribution is present. That result is given by eq. (6.3), with the substitutions ν a , e a → ℓ a and h, Z, W → φ.

LPM resummation
If the particles participating in 1 ↔ 2 scatterings are ultrarelativistic (m i ≪ ǫ i , M ≪ ω) and if their virtualities are small compared with the rate of thermal scatterings (m 2 i /ǫ i , M 2 /ω ≪ g 2 T /π), then many scatterings of the type 1 + n ↔ 2 + n take place, and need to be summed to all orders, through so-called LPM resummation [63,[92][93][94]. An implementation of LPM resummation guaranteeing that the Born limits of sec. 6.1 are correctly recovered when we exit the ultrarelativistic regime, has been worked out in ref. [33]. However, if M is not small compared with thermal scales, for instance if we find ourselves in the regime M ∼ πT , then no LPM resummation is needed.

Real part of the direct contribution
The retarded correlator of eq. (4.9) contains also a real part, yielding dispersive corrections according to eq. (4.16). In contrast to the imaginary part, where a "cut" leads to phase-space constraints and separate classes of diagrams (cf. secs. 5 and 6), the real part contains an unconstrained integral, whose magnitude is simpler to estimate. Normally, it is sufficient to restrict to 1-loop level in its evaluation. Then the direct contribution originates from eq. (6.3) and the indirect one from eq. (6.5). 6 To be concrete, let us define a master sum-integral where the notation adheres to that in eq. (5.13). Then we may rewrite eq. (6.3) as The real part is, in general, ultraviolet divergent. As the divergences are related to NLO corrections, notably wave function renormalization, we address only the part proportional to the plasma four-velocity, which is finite and plays a qualitatively more prominent role. Following ref. [70], this part is defined by writing Re J − 1 2 00 By projecting onto the part proportional to / U , a straightforward analysis leads to

Full indirect contribution
Consider finally the indirect contribution, given by eq. (6.4). Here the self-energy / Σ includes the 2 ↔ 2 and 1 ↔ 3 contributions from sec. 5, the 1 + n ↔ 2 + n contributions from sec. 6, as well as the mass corrections from the real part of eq. (6.5).
For the real part we only need the contribution proportional to / U , as in eq. (7.3). In addition, as alluded to in sec. 4.1, at this point it is convenient to include the chemical potential representing a Z-tadpole (ref. [70], eq. (A.7)), We now return to eq. (6.4). As the self-energy / Σ includes both a real and an imaginary part, the real and imaginary parts of the inverse read Before the inversion, we should disentangle the imaginary part of the self-energy just like the real part in eq. (7.3) [104], viz.
Matrix elements taken according to eqs. (4.21) and (4.22) then reduce to the corresponding matrix elements of / K and / U , which are straightforward to evaluate.

Conclusions
The purpose of the present review has been to collect together the theoretical background for computing the thermal rate coefficients and mass corrections that parametrize rate equations for sterile neutrinos in the early universe at O(h 2 ν ) in neutrino Yukawa couplings, while keeping the mass scale M and temperature T general. This problem originates from a minimal extension of the Standard Model according to eq. (1.1), and may play a physical role in certain dark matter and leptogenesis scenarios. For dark matter, temperatures much below the electroweak crossover (deep in the Higgs phase) are relevant as well.
In particular, we have verified the gauge independence of the appropriate retarded correlator (cf. eqs. (4.15) and (4.16)), both in the symmetric and in the Higgs phase of the electroweak theory, and after incorporating all 1 ↔ 2, 2 ↔ 2 and 1 ↔ 3 processes. We have shown how in the Higgs phase a subset of the effects can subsequently be resummed into an "indirect contribution", even though this is far from obvious at first sight, given that on the amplitude level there are interference terms between direct and indirect sets (cf. fig. 2). The indirect contribution dominates at low temperatures, and amounts to processes where active neutrinos interact with the plasma and are then converted into sterile ones. The computations have included generic momenta and chemical potentials, and permit for arbitrary helicity projections.
The missing step of the program would be a numerical evaluation of all rate coefficients. Even though tools for this have become available recently [32,33], which have been demonstrated to be adaptable to a broad range of problems [105], their practical implementation to the current theory represents a challenge of its own. One reason is that the equations incorporating LPM-resummation [33] are inhomogeneous differential equations and can be time-consuming in practice. Another is that the number of objects to be considered "proliferates": we need to determine both the real and imaginary part of the retarded correlator; for both helicity states; for both the direct and indirect contributions. Thus there are many observables to determine; all of them as functions of M, T, µ i , k, where µ i stands for the set of chemical potentials and k is a momentum; in broad numerical ranges (say, M = 10 1...3 GeV; T = 10 −4...4 GeV; |µ i | = 10 −10...−5 T ; k = 10 −2...2 T ). The numerical integrations become demanding if some particle mass is much lighter or heavier than the temperature, and therefore cannot be performed "on the fly". Rather, the results are to be tabulated on a dense grid in the aforementioned parameter space, suitable for a subsequent interpolation. We hope to attack this challenge in the future, noting that partial results in a similar spirit have already been collected on the web pages associated with refs. [65,69,83].

A.1. Notation
In this appendix we specify the matrix elements squared that correspond to 1 → 3 decays of right-handed neutrinos of mass M into Standard Model particles. As discussed in secs. 5.1 and 5.2, the same matrix elements also fix, via crossing symmetries, the thermal interaction rates originating from 2 ↔ 2 and 3 → 1 processes as well as those from virtual corrections.
In order to display the results, we make use of Mandelstam variables, viz.
which are related through Here m i denote the masses of the decay products. The U Y (1) and SU L (2) gauge couplings g 1 , g 2 are often represented through In the symmetric phase, the only non-vanishing mass (apart from M ) is that associated with scalar particles, and is denoted by m φ . In the Higgs phase, most particles carry masses, however those of the charged leptons and active neutrinos have been omitted for simplicity (the same concerns lepton Yukawa couplings). By γ we denote a massless neutral gauge boson, specifically photon in the Higgs phase. In the case of charged gauge bosons, the symbol W denotes W + ,W denotes W − . Antiparticles are sometimes denoted with a minus sign, e.g. −W ≡W , as they come with opposite chemical potentials.
For the representation of the matrix elements squared, we make use of eq. (5.16), viz.
Then "direct" contributions, which do not involve the insertion of a Higgs vev, are expressed in terms of E, whereas "indirect contributions", proportional to 1/M 4 or 1/M 2 , are expressed in terms of I. This can be achieved by employing the relation in terms multiplied by 1/M 2 . It turns out that powers of 1/M 2 come with coefficients c n such that c n (s jk − m 2 i ) = O(M 2 ), so that eq. (A.5) eliminates the appearance of 1/M 2 . Going in the opposite direction, I's not multiplied by a Higgs vev can be eliminated by Unfortunately, the splitup into E and I is not unique. This should not be surprising, given that the direct and indirect sets are not gauge independent by themselves. In addition, powers of M 2 cannot be identified uniquely, as kinematic variables are related through eq. (A.2). We have attempted to employ this freedom in order to express the results in a compact fashion, however the procedure has ambiguities. Furthermore, the appearance of kinematic variables in the numerator leads to several alternative representations. In the end, to verify the equivalence of two representations, one needs to choose some lexicographic ordering and an algorithm putting the expression in this form, which leads generically to lengthier expressions, but unique ones in the chosen basis.
We express the result as a contribution to Θ E from eq. (4.21). In terms of matrix elements, the results correspond to |M| 2 , with all final-state spins summed over. For the righthanded states, eq. (4.21) implies that we need to adopt the effective replacement with E fixed according to eq. (4.22). For thermal phase space averaging we make use of the notation in eqs. (5.1)-(5.3).
The expressions for the matrix elements squared enjoy certain symmetries that provide for stringent crosschecks of their correctness. This concerns particular the interference terms, which necessarily appear in the structure of eq. (2.34) of ref. [32], viz. . (A.8) By making use of eq. (A.5), it can be verified that the interference terms in eqs. (A.11) and (A.12) indeed satisfy this property.
Machine-readable versions of eqs. (A.9)-(A.12) are attached to this paper as ancillary files.

A.2. Hadronic effects, symmetric phase
In the symmetric phase, decays into final states containing hadrons amount to Here h t and h b are the top and bottom Yukawa couplings, the others having been omitted.

A.3. Hadronic effects, Higgs phase
When the Higgs mechanism is active, both direct and indirect contributions are present, and many different masses play a role. In this section we display the processes involving hadronic final states. It can be checked that for m i → 0 the Higgs phase results go over into eq. (A.9), when the latter is evaluated with m φ → 0. For simplicity we keep only the top and bottom quark masses finite (m t , m b ), whereas the other quarks are treated as massless, with furthermore the labels set to those of the up and down quark, viz. c → u, s → d. Then the results read (A. 10) The last two terms, originating from charged currents, have been simplified slightly, suppressing small off-diagonals from the CKM matrix elements |V td | 2 , |V ts | 2 , |V ub | 2 , |V cb | 2 as well as parity-odd contributions involving the Levi-Civita tensor.

A.4. Leptonic effects, symmetric phase
In the symmetric phase, the processes only involving leptons in the final state amount to (A.11) Here φ denotes a scalar field and γ a massless gauge boson.