Efficient numerical integration of thermal interaction rates

In many problems in particle cosmology, interaction rates are dominated by ${2}\leftrightarrow{2}$ scatterings, or get a substantial contribution from them, given that ${1}\leftrightarrow{2}$ and ${1}\leftrightarrow{3}$ reactions are phase-space suppressed. We describe an algorithm to represent, regularize, and evaluate a class of thermal ${2}\leftrightarrow{2}$ and ${1}\leftrightarrow{3}$ interaction rates for general momenta, masses, chemical potentials, and helicity projections. A key ingredient is an automated inclusion of virtual corrections to ${1}\leftrightarrow{2}$ scatterings, which eliminate logarithmic and double-logarithmic IR divergences from the real ${2}\leftrightarrow{2}$ and ${1}\leftrightarrow{3}$ processes. We also review thermal and chemical potential induced contributions that require resummation if plasma particles are ultrarelativistic.


Introduction
In an interacting system, the simplest processes leading to thermal reaction (or damping, or interaction) rates are 1 → 2 decays as well as 2 → 1 "inverse decays". However, the phase space for such processes is kinematically constrained, as exemplified by the impossibility of elements squared must be cancelled by the identification and inclusion of the pertinent virtual corrections to 1 ↔ 2 processes. For the subsequent numerical evaluation, we generalize the parametrizations that were introduced in ref. [4], to a situation with an arbitrary mass and chemical potential spectrum of the participating particles. The numerical and algebraic scripts implementing these ingredients to a particular example are made publicly available, with the hope that they can be adapted to other problems with modest effort.
The presentation is organized as follows. We start with an exposition of how to concisely represent 2 ↔ 2 and 1 ↔ 3 scatterings, and how to cancel the poles that appear in them, in sec. 2. A regularization permitting to evaluate various contributions before the cancellation of poles is introduced in sec. 3. Subsequently parametrizations for the remaining integrations are introduced, in sec. 4. Even if all poles have been cancelled, the results still suffer from divergences in certain limits, as reviewed in sec. 5. Sec. 6 summarizes the algorithms implementing the ingredients of secs. 2-4. Conclusions are offered in sec. 7, deferring formal proofs of thermal crossing relations to appendix A, integration prescriptions for specific channels to appendix B, and angular averages for virtual corrections to appendix C.
It is apparent from eqs. (2.2)-(2.10) that a minus sign in front of a label belonging to an initial state, e.g. −a, corresponds to an inversion of the sign of the corresponding fourmomentum and chemical potential. This pattern continues to hold once matrix elements squared are added. While the chemical potentials have been treated within the functions N , it is convenient to interpret scat n → m as an "operator", acting on the momenta that appear 3 The overall factor 1 2 can be viewed as a convention, but it also guarantees a straightforward connection to Boltzmann equations on one hand (so that Θ in eq. (2.12) corresponds to a "matrix element squared"), and to the retarded correlators considered in appendix A on the other (originating then as Im 1 i ǫ i −i0 + = πδ( i ǫi) = 1 2 (2π)δ( i ǫi), with 2π included as part of the usual energy-momentum conservation constraint). show the 1 → 2 process leading to Γ Born 1→2 = scat 1 → 2 (ℓ, φ) 4E · P ℓ , and the 1 → 3 processes leading to eq. (2.13). Double lines correspond to right-handed neutrinos; arrowed lined to leptons (ℓ); wiggly lines to gauge fields (γ); dashed lines to scalars (φ).
in the matrix element squared, with the rule that We now posit that for considering 2 ↔ 2 and 1 ↔ 3 scatterings, it is sufficient to define the integrand needed for 1 → 3 decays, with the other processes following from it through crossing symmetries. The reason for selecting the 1 → 3 decays is that they display maximal symmetries, with all thermalized particles appearing in the final state (of course, in phase space integrations, inverse processes are included as well, cf. eq. (2.4)). A formal proof of this statement can be found in appendix A.1.
As an example, consider right-handed neutrino interactions in the symmetric phase of a Standard Model plasma [3,4]. Then 1 → 3 decays, illustrated in fig. 1, can go to a final state with a lepton, a gauge boson, and a scalar particle (a, b, c = ℓ, γ, φ). 5 The matrix element squared |M| 2 for this process, summed over final-state helicities and with the substitution τ u kτūkτ → / E for the initial-state helicity sum, but for convenience factoring out Yukawa 4 We are somewhat lax with terminology here: depending on how Θ is chosen, the "rate proper" may still require a division of Γ by ω, as appears for instance in the collision term of a Boltzmann equation. 5 In the full physical computation, there are also channels in which γ, φ are replaced by a fermion-antifermion pair, originating from a scalar decay, however the algebraic structure of the corresponding matrix element squared is rather simple, and brings in nothing new to the case that we consider.
For the following, it is helpful to rewrite the physical cross section in eq. (2.13), by denoting the locations of the propagator poles with specific symbols, which may be called auxiliary masses. This is helpful not only as a labelling of pole locations, but also serves as an intermediate IR regulator (cf. sec. 3), and in addition permits to take parametric derivatives with respect to specific propagators, as is needed later on. With this motivation, we replace Γ Born 1→3 2(g 2 1 + 3g 2 2 ) → scat 1 → 3 (ℓ, γ, φ) (2.14) The (would-be) masses of the decay products are denoted by m 2 ℓ ≡ P 2 ℓ , m 2 γ ≡ P 2 γ , m 2 φ ≡ P 2 φ .

Corresponding virtual corrections
If the energy entering one of the propagators in eq. (2.14) equals the on-shell energy of that line -that is, if three energies are related to each other like in a 1 ↔ 2 process -then the integrand diverges, and the value of the integral is ambiguous. In accordance with the KLN theorem [23,24], such divergences are cancelled by virtual corrections to 1 ↔ 2 processes. Conversely, requiring the cancellation of divergences, we can reconstruct virtual corrections.
A key observation for cancelling the divergences is that the thermal distribution functions appearing in eqs. (2.4), (2.7), (2.10) can always be factorized, in fact in several ways: In the cases involving energy differences, another useful representation can be obtained through the identity n σ (−x) = −1 − n σ (x). Eq. (2.15) can be put in two alternative forms by the exchanges a ↔ c and b ↔ c, whereas eqs. (2.16) and (2.18) give additional identities via the exchanges b ↔ c and a ↔ b, respectively.
In order to make the statement concrete, let us denote thermally averaged phase space integrals for 1 ↔ 2 processes, interpreted as operators acting on the momenta P a and P b with the rule introduced in eq. (2.11), by Furthermore we denote

Single pole
Let us start by considering a first-order pole in a variable s ab . The corresponding "residue" is denoted by Θ where Θ(P a , P b , P c ) is the matrix element squared appearing in eq. (2.12).
There are various ways to establish what kind of 1 ↔ 2 processes are needed in order to cancel this pole. Probably the most economical is to consider a master sum-integral, as explained in appendix A.2. In this case, the 2 ↔ 2 and 1 ↔ 3 processes and the virtual corrections to 1 ↔ 2 processes that cancel their poles are generated simultaneously from a single structure that respects the KLN theorem. Alternatively, one could take the cancellation of divergences as the construction principle, and then make use of the factorization formulae in eqs. (2.15)- (2.19), however the implementation of the latter strategy is tedious.
We adopt a notation where Γ Born 1↔2 denotes the interaction rate corresponding to a Bornlevel 1 ↔ 2 process, and ∆Γ Born 1↔2 the virtual corrections to this process. Our algorithm reconstructs, strictly speaking, only the IR sensitive part of the full ∆Γ Born 1↔2 . Specifically, the virtual correction associated with a residue like in eq. (2.27) reads (2.28) Here we have defined a "bubble" operator, B, acting on the momenta P a and P b , and having P d as an input parameter, as Let us stress again that, taken literally, eq. (2.29) is not well-defined on its own, due to the poles that the integrands contain, however it is well-defined in connection with eq. (2.12), because then the poles cancel.

Multiple poles in a single variable
If some process can be mediated by different particle types (e.g. Higgs and Z 0 bosons), then it is possible that interference terms contain two different poles in a single kinematic variable. In this case, we can partial fraction the dependence on that variable, and then make use of eq. (2.28).
A related situation is met if there is a second order pole in a single variable. The corresponding residue is denoted by The simplest procedure in this case is to make use of the result for a single pole, just taking a mass derivative thereof: Here we have assumed a notation like that introduced in eq. (2.14), whereby m 2 d appears as a pole location. We note that the mass derivative acts both on scat 1 ↔ 2 (d, c), and on the coefficient function B(P d ; a, b) Θ (2) ab , which is a function of m 2 d via its dependence on P d .

Single poles in separate variables
The most complicated case is when there are poles in two separate kinematic variables, as can happen in interference terms. Say, if the variables s ab and s bc can go on-shell, then the corresponding residue is defined as where Θ is the function in eq. (2.12).
It is important to realize that if a residue like in eq. (2.33) exists in a 1 → 3 matrix element, then there is necessarily also another decay channel, namely one in which the particles d, b, e are in the final state, whereas the propagators of a and c appear in the matrix element squared. Concretely, noting that the direction of the line b gets inverted in the second "cut", poles of this type must appear in the combination . (2.34) Once again, a convenient way to prove that interference terms necessarily appear in this form goes through the consideration of master sum-integrals, as shown in appendix A.1. We note that if the particle b is neutral (i.e. carries no chemical potential), and a = d, c = e, then eq. (2.34) can be reduced to a single structure. The virtual corrections cancelling the poles of eq. (2.34) can now be expressed as where we have defined a "cubic" or "triangle" operator, C, which is similar to the bubble operator in eq. (2.29), but acts now on three momenta, with two incoming momenta (those related to scat 1 ↔ 2 ) appearing as parameters. The incoming momenta are constrained by their sum equalling K (cf. eq. (4.19)), which permits to simplify some of the structures appearing: .

(2.36)
A rather tedious analysis shows that eq. (2.35) indeed cancels all divergences of eq. (2.34); alternatively, eq. (2.35) can be obtained by considering the cuts of a corresponding master sum-integral, as discussed in appendix A.2.
It is important to stress that the two terms in eq. (2.35) do not correspond one-to-one to the two terms in eq. (2.34). Rather, to cancel the poles of the first term in eq. (2.34), parts of both terms of eq. (2.35) are needed. Nevertheless, since both terms of eq. (2.34) are guaranteed to appear in a correct computation of Γ Born 1→3 , we can effectively adopt an algorithm which derives the first term of eq. (2.35) from the first term of eq. (2.34), and similarly for the second terms. This does produce the correct ∆Γ Born 1↔2 for eq. (2.35).

Example
For the case in eq. (2.14), there are poles corresponding to the variables s ℓγ and s γφ , located at m 2 ℓ and m 2 φ , respectively. We also note that the interference term respects the symmetry described below eq. (2.34), so that only one of the structures in eq. (2.35) needs to be included. This then leads to

IR regularization
Having found a representation for 2 ↔ 2 and 1 ↔ 3 scatterings (cf. eq. (2.12)) as well as the virtual corrections that render these expressions finite (cf. eq. (2.37)), the next task is to carry out the integrals in practice. However, the two sets are ill-defined separately. Therefore we introduce generic masses as intermediate regulators that render both sets finite, and check in the end that the results are stable if the masses are taken to their physical values. For the following, it is helpful to assume that only first order poles are present. Second order poles, like the one appearing in eq. (2.14), are to be viewed as parametric derivatives, , as was already done in eq. (2.37). If we partial fraction the energy dependence of the 2 ↔ 2 and 1 ↔ 3 matrix elements squared, it appears in products of fractions like Let now ǫ 3 represent an energy variable that becomes singular in the sense of eq. (2.27). In the 2 ↔ 2 and 1 ↔ 3 scatterings, these fractions always appear pairwise, such that when put together we end up having combinations like where the numerator cancels against a similar term from the integration measure. This reflects the fact that ǫ 1,2 are energies of external state particles, whereas ǫ 3 is that of a virtual state. In the virtual corrections, the same would-be divergence appears as (3.3) Now ǫ 1 or ǫ 2 does not appear as an integration variable, but rather represents the energy of a virtual particle, ǫ 2 1 = (p 2 ± p 3 ) 2 + m 2 1 or ǫ 2 2 = (p 1 ± p 3 ) 2 + m 2 2 .
We define integrals over eqs. (3.2) and (3.3) as principal values. Proceeding first with the virtual corrections, the principal value concerns the integration over the directions of a momentum like p 2 or p 1 . As shown in sec. 4.2 and appendix C, these angular averages can be carried out analytically, and the results are non-singular as long as masses are finite, apart from integrable (logarithmic) singularities that influence final energy integrations.
Returning then to 2 ↔ 2 and 1 ↔ 3 scatterings, the situation is more complicated. In many cases we can take s ± 12 = (ǫ 1 ± ǫ 2 ) 2 − (p 1 ± p 2 ) 2 as an integration variable in eq. (3.2), and then the singularity is located at the position where s ± 12 = m 2 3 . However, in interference terms, two separate poles appear. Even though one of the singularities can still be easily localized, the other one is more challenging. The principal value integration per se can be dealt with in connection with azimuthal averaging (cf. eq. (4.10)), however a remnant divergence can manifest itself in energy variables.
We note that, if the singularities of the energy variables cannot be localized analytically, it is possible to take care of them numerically. This goes by implementing principal value integration as a limit, viz.
In general δ has to be set to a small value for reliable results, e.g. δ ∼ (10 −4 T ) 2 . Finally, it is appropriate to stress that actual poles do not always appear. To this end, consider a general structure appearing in eq. (2.34), . (3.5) The allowed ranges of the kinematic invariants for the 1 → 3 process ("Dalitz plots") and its crossings can be established as usual (cf. sec The 3 → 1 channels can be realized in three different ways, and poles are met if For the example in eq. (2.14), with a → ℓ, b → γ, c → φ, d →l, e →φ, and setting subsequently mφ → m φ and ml → m ℓ , and assuming that m φ and M are macroscopic whereas m ℓ and m γ are small, poles can only exist in the 2 ↔ 2 channels, and then only if In the physical world, both of these masses vanish in an unresummed computation. The results depend continuously on them, so that the massless limit is well-defined. If we take m γ → 0 first, eq. (3.12) indicates that there are no poles in any of the real scatterings. Even if unproblematic in principle, the case in which poles do appear is discussed from a different perspective in sec. 5.4.

Phase space integrals for t-channel 2 ↔ 2 reactions
Among the reactions in eq. (2.12), we start by considering 2 ↔ 2 processes, which in many cases are physically the most important ones. A corresponding analysis of 1 → 3 and 3 → 1 processes can be found in appendix B.
For purposes of the practical integration, it is convenient to label the initial and final-state four-momenta, energies, masses and chemical potentials with different symbols (K a vs. P a , E a vs. ǫ a , M a vs. m a , and ν a vs. µ a , respectively). Here on-shell energies are denoted by so that K a = (E a(ka) , k a ) and P a = (ǫ a(pa) , p a ). The momentum label can often be left out from the energies without a danger of confusion. With the given momenta and channels, Mandelstam invariants are defined as usual, For 2 ↔ 2 processes, we introduce two parametrizations, the t-and s-channel ones [4]. In principle both of them can be used for any matrix element squared. However, if only one propagator appears in the integrand, the kinematic variable should be chosen to correspond to that variable, because then we can easily identify the pole location (if one appears), and implement the corresponding principal value integration (cf. sec. 3). In the interference terms, where two propagators appear, it is convenient to choose the parametrization according to which is the "most singular" kinematic invariant, in order to understand the behaviour in the massless limit (cf. sec. 5.1). The most singular variable is identified by putting auxiliary masses such as those introduced in eq. (2.14) to zero. Furthermore, if u appears as a singular variable, we can make use of the substitution P 1 ↔ P 2 , in order to rename it into t. Thus, the most singular variable can always be chosen as t or s.
With this prescription, and numbering the momenta and chemical potentials corresponding , so that the sign flips associated with initial-stage particles have already been included when showing K 1 and ν 1 , the crossings of eq. (2.12) yield the following t-channel rate for the example of (2.14) (the case of s-channel scatterings is discussed in appendix B.2): Here E can be K, or the medium four-velocity, U ≡ (1, 0), or a linear combination thereof.
The key idea now is to have t as an integration variable, so that most of the poles are easily resolved. To achieve this we introduce a four-momentum Q such that t = Q 2 , and rephrase the integration measure from eq. (2.6) as where we have integrated over p 2 and k 1 . The Dirac-δ's fix two angles as where p i ≡ |p i |. The other Mandelstam variables can be expressed as If products like 1/(u m s n ) appear, with m, n > 0, the dependence on k · p 1 should be partial fractioned to appear as powers of inverse first-order polynomials. The azimuthal angle between k and p 1 is unconstrained as of now. We note that this angle does not appear inside the thermal distribution functions in this parametrization (cf. eq. (4.18)). Choosing q as the z-axis, we may write k · p 1 = kp 1 cos θ q,k cos θ q,p 1 ≡a + kp 1 sin θ q,k sin θ q,p 1 ≡b cos ϕ , (4.9) where cos θ q,k and cos θ q,p 1 are fixed through eq. (4.7). Let us define a generating function, incorporating the possibility of a numerical principal value regularizationà la eq. (3.4), as This average vanishes if |z − a| < |b| (for δ → 0), i.e. when poles are actually crossed. Nevertheless the function is singular when approaching this domain from the outside, i.e. |z − a| → |b| + , so that a regularization is still needed, either numerically through δ, or analytically, by resolving the singular domains of the outer integrations. Through a series in 1/z, the azimuthal averages of positive powers of k · p 1 are readily extracted from eq. (4.10). In this case we can put sign( z − a ) → 1, yielding Derivatives with respect to z yield averages of negative powers, however they are not needed for eq. (4.5), unless we use a t-channel parametrization for an s-channel process (cf. eq. (B.9)). With all angles resolved, the remaining integration measure can be worked out. As a first step, we may consider separately the two "vertices" (or energy conservation constraints) appearing in eq. (4.6), finding where ǫ ± 1 are from eq. (4.16). Subsequently it is advantageous to replace the variables q 0 , q through t = q 2 0 − q 2 , q 0 , yielding a Jacobian The integration range of t can be established as (−∞, (m 2 − M ) 2 ), and for the part t > 0 we find that the sign of q 0 depends on M − m 2 , such that θ(t)θ((M − m 2 )q 0 ) may be inserted in the integrand. Altogether, eq. (4.12) can thus be converted into where q = q 2 0 − t and Here the Källén function is defined as Finally, the thermal distributions associated with 2 ↔ 2 scatterings, cf. eq. (2.7), can be factorized as in eq. (2.16), which after the insertion of q 0 = ǫ 1 − E 1 = ω − ǫ 2 from eq. (4.6) and renaming chemical potentials as mentioned above yields The first factor implies that the q 0 -integral is exponentially localized around modest |q 0 |, even if the boundaries from eq. (4.15) can obtain large values; likewise, the integration over ǫ 1 is exponentially convergent at large ǫ 1 .
We note that if τ 1 σ 1 = +, the first factor in eq. (4.18) has a pole at q 0 − µ 1 + ν 1 = 0, but the second factor vanishes at the same point, lifting it. In order to avoid numerical issues with this, it is helpful to replace eq. (4.18) by an alternative representation if

Phase space integrals for virtual corrections
The virtual corrections, discussed in sec. 2.2, have the structure of a 1 ↔ 2 phase space average, which we call an "outer integral", convoluted with a bubble or triangle function, which we call an "inner integral" (cf. eqs. (2.28), (2.32), (2.35)). Let us discuss these in turn.

Outer integration
An outer integral like scat 1 ↔ 2 (d, c) in eq. (2.28) contains three different channels, corresponding to 1 → 2 decays and two different 2 → 1 inverse decays, cf. eq. (2.26). Remarkably, the three channels can be combined into a single expression, once we make use of the identity n σ (−x) = −n σ (x) in the 2 → 1 channels and substitute ǫ i → −ǫ i in one of them. For instance, if we choose ǫ d as the integration variable, then where √ λ is from eq. (4.17), we have expressed the spacetime dimension as D = 4 − 2ǫ, and the integration bounds are analogous to eq. (4.15), viz.
Furthermore the direction of p d is fixed similarly to the second constraint in eq. (4.7), The sign function in eq. (4.19) is negative in the 2 → 1 channels, in one because of the inverted sign of ǫ d and in the other because of the negative sign of ω − ǫ d . 6 For completeness we have shown the part of O(ǫ) in eq. (4.19), because the inner integrals normally contain 1/ǫ divergences;μ denotes the scale parameter of the MS scheme. In this context is also appropriate to note that in D dimensions, some of the coefficients in eq. (2.13) contain additional parts proportional to D−4, which had been omitted. If we worry about the renormalization of the coupling associated with the vertex in fig. 1(left), these terms should be included, but normally such renormalization effects are genuinely small in a weakly coupled system, unlike the IR effects that we are mostly interested in.

Inner integrations
For the inner integrals, we need to consider B from eq. (2.29) and C from eq. (2.36). For generality we start with C, which contains one more angular variable. Let us recall that all propagators within the thermal averages are assumed regularized as principal values.
We may now express the integration measure in spherical coordinates as usual. Assuming that the inner integral is UV finite (we return to this below), we write and similarly for p b and pc in eq. (2.36).
For the integrand, we consider two possible structures, defined in eqs. (C.1) and (C.2). There is freedom in choosing the axis with respect to which θ is measured. This way, the angular integrals can be carried out explicitly, as detailed in appendix C.
The two middle terms of eq. (2.36), associated with ǫ b , require a more careful look, as the dependence on the angles needs to be partial fractioned, in order to bring the results in a form in which we can make use of eq. (C.2). For this we may write where we made use of P d + P e = K.
Collecting together the contributions, and making use of H from eq. (C.2), we obtain

(4.24)
We have used the sign ≃ because the integrals are, in general, UV divergent; below eq. (4.25) we turn to how the divergences can be subtracted.
For the integral B from eq. (2.29), where only one propagator appears, the results can be obtained from eq. (4.24) by appropriately dropping structures and renaming variables: Here G is from eq. (C.1).

UV subtraction
The nature of the integrals in eqs. (4.24), (4.25) depends on the functions Φ which, in turn, depend on the spin structure of the non-equilibrium particle considered. For the example in eq. (2.14), which leads to eq. (2.37), Φ contains the helicity projectors E · P i . In this case, or if Φ is independent of the integration momenta, only the function B needs a subtraction. In case of a quadratic momentum dependence, C would need a subtraction as well. However, even if C were finite, we endorse adopting a subtraction procedure, specified as follows.
The most straightforward implementation of a subtraction is to omit the vacuum contributions from eq. (4.25), and to compute them separately. That is, we may write where B T is obtained by dropping the temperature-independent factors 1 2 from eq. (4.25): (4.28) The part B T vanishes at zero temperature. By returning back to eq. (2.29), it is possible to verify that the vacuum factors, in turn, amount to a usual vacuum integral, 7 . (4.30) Now, the vacuum integral can be worked out with standard methods. If the function Φ in eq. (4.30) depends on momenta, the result can be reduced to integrals with trivial numerators by making use of Passarino-Veltman relations, however the details depend on the helicity structure considered. For the case of eq. (2.37), where the momentum dependence is The integrals on the second line of eq. (4.32) can be evaluated, where terms of O(ǫ) have been omitted. The integral over x could be carried out [25], however the integral representation is a practical tool. The are mild singularities within the integration . The presence of 1/ǫ-divergences in eqs. (4.33) and (4.34) reflects the fact that we are considering a bare expression. In a full computation, the 1/ǫ-divergences of virtual corrections 7 What is meant here is that the factors in eq. (2.29) combine into where the integral has been defined after Wick rotation to Euclidean signature. Rotating back, the principal value (cf. sec. 3) corresponds to the real part of a usual Feynman integral, as indicated in eq. (4.30).
cancel against the renormalization of the parameters that appear in the leading-order 1 ↔ 2 process. For simplicity, we assume in the following that this renormalization has been taken care of, specifically that the 1 ↔ 2 process is expressed in terms of MS scheme couplings and masses, evaluated at the scaleμ = 2πT . Under this assumption, the 1/ǫ divergences can be dropped from eqs. (4.33) and (4.34) and the scale can be set toμ = 2πT . It is appropriate to stress, however, that were we considering the bare rate, before renormalization, then we should recall that the outer integral of eq. (4.19) contains a part of O(ǫ), as well as possible overall coefficients proportional to D − 4, which would give a finite contribution once multiplied by the 1/ǫ-divergence.
A splitup into thermal and vacuum parts is also possible for the triangle integral C, even if it does not always contain UV divergences: (4.36) The vacuum limit from eq. (2.36) can be expressed as 8 . (4.37) In the case of a linear dependence like in eq. (2.37), we can write The first structure points in a direction spanned by the latter two vectors, where , (4.41) . (4.42) The Feynman representation for the first rows of eqs. (4.41) and (4.42) reads where we have followed ref. [25], denoting a ≡ m 2 , y 0 ≡ −(d + eα)/(c + 2bα), y 1 ≡ y 0 + α, y 2 ≡ y 0 /(1 − α), y 3 ≡ −y 0 /α. The kinematics of the outer integration guarantees that α is real (cf. eq. (2.35)); the sign in front of √ λ is best chosen so that no large cancellation takes place in the numerator of α. If the arguments of the logarithms have real zeros, which happens for λ(m 2 a , m 2 b , m 2 d ) ≥ 0, λ(m 2 a , m 2 c , M 2 ) ≥ 0, λ(m 2 b , m 2 c , m 2 e ) ≥ 0, respectively, we may express the corresponding integral in terms of (real parts of) four dilogarithms [25], otherwise the integral representations are quite efficient. Such divergences imply that the naive perturbative series needs to be re-organized, or "resummed". Unfortunately it is difficult to automate this treatment, as it requires a carefully tuned subtraction-addition step, in order to avoid double counting when implementing the resummation that eliminates the divergence. Here we illustrate the procedure for the example of fig. 1, and briefly mention other widely discussed applications.

2 ↔ 2 scatterings via soft t-channel exchange
We start by considering what may be referred to as the high-energy small-angle limit, parametrically λ 2 ≪ −t ≪ s, where λ 2 ∼ m 2 i , M 2 and s ∼ (πT ) 2 . In vacuum the structure of scattering amplitudes has been studied in great detail in this "Regge" domain. If λ 2 ≪ (πT ) 2 , the thermally averaged scattering rates in general diverge, as we now review.
If we expand the integrand in eq. (4.5) in a Laurent series in q 0 , q, the leading term may come with a negative power of q 0 , q and a positive power of ǫ 1 . The latter integral is weighted by the thermal distribution functions, and therefore turns into a positive power of the temperature. To compensate for the dimensions, the domain of small q 0 , q then leads to an IR divergence.
In order to isolate this divergence, we may work out the angles from eq. (4.7), the azimuthal average from eq. (4.10), as well as all kinematic variables, in the massless limit. To facilitate this task, it is essential that the parametrization through Q has been so chosen that the divergences are associated with small values of q 0 , q. Denoting the massless value of ǫ 1 by p 1 , we may approximate the integrand by taking q 0 , q ≪ p 1 . (5.1) In this limit eq. (4.7) implies that We remark that k is often approximated as ∼ p 1 ≫ q 0 , q, however we have included the correction (q 2 − q 2 0 )/(2qk) in cos θ q,k , given that the leading term cos θ q,k ≃ q 0 /q cancels against cos θ q,p 1 , when estimating the magnitude of inverse powers of u and s in eqs. (5.6) and (5.7). In the massless limit, the integration measure from eq. (4.12) becomes implying that t < 0, i.e. |q 0 | < q, so that the angles in eq. (5.2) are well-defined.
It should be stressed that, as long as masses are non-zero, these IR sensitivities are not divergences in a literal sense. Their presence simply implies that if we consider the limit m i , M ≪ πT , the thermally averaged cross sections become anomalously large. The way to handle the large terms goes through Hard Thermal Loop (HTL) resummation [26][27][28][29], which gives a thermal mass ∼ g 2 T 2 to the would-be massless modes. Subsequently the Born-level result can be rectified through a subtraction-addition step, with the addition part involving the HTL-resummed result and the subtraction part its unresummed form, eliminating double counting. Adapting a clever trick developed in ref. [30], whereby thermal light-cone observables can be analytically continued to static ones, HTL-resummed results have been worked up to the NLO level for certain interaction rates characterizing massless or nearly massless probe particles [31,32].
For the concrete example of fig. 1, HTL resummation affects the lepton propagator, whereas the scalar propagator only experiences a mass correction. 9 Let us denote the lepton spectral function by where the temporal and spatial parts read (q 0 ≡ Re q 0 + i0 + ) [33,34] ρ µ (q 0 , q) = Im 1 − In general, HTL computations include vertex corrections as well. There is none for our example with a Yukawa vertex, and more generally they are unimportant when considering "hard" momenta, k ∼ πT .
For |q 0 | < q, the overall imaginary part originates from Im L = −π/(2q), whereas for |q 0 | > q, it originates from a zero of the denominator. The thermal lepton mass reads The interaction rate Γ HTL is determined by making use of eq. (5.8). Integrating over angles, the remaining measure can be cast in the form of eq. (4.12), with q · k fixed according to eq. (4.7) and the integral over ǫ 1 represented by eq. (5.10). In the HTL spectral function, it is the spacelike domain, with θ(−t), which corresponds to 2 ↔ 2 scatterings.
Subsequently, we have to subtract the part already included in eq. (4.5). This subtraction can be obtained by "re-expanding" eq. (5.9) to O(g 2 T 2 ) ∼ O(m 2 ℓT ). With a slight abuse of notation, we however keep a lepton mass in the denominator, renaming it into m 2 ℓ , whereby it serves as an intermediate IR regulator in the sense of sec. 3. This then leads to , where E ≡ (ε 0 , ε ), and the integration bounds are a special case of eq. (4.15), viz.
Let us mention that if we choose ml = m ℓT in eq. (5.11), then the numerical value of eq. (5.11) is quite small. This implies that using the thermal value m ℓT in the unresummed computation of sec. 4.1 would yield a reasonable approximation to the resummed result. That said, the error would not be parametrically of higher order; eq. (5.11) is of higher order only as far as the domain q 0 , q ≫ m ℓT is concerned, where the subtractions are effective.
In eq. (5.11), we have introduced a function Λ(q 0 , ω) that necessitates further elaboration. In order to correctly implement HTL resummation for q 0 , q ∼ m ℓT ≪ πT , we require lim q 0 →0 Λ(q 0 , ω) = 1. What we do outside of this domain, where the effect is formally of higher order, is a matter of choice. A frequent logic is to set Λ → Λ ⋆ , where A benefit of this choice is that the remaining integrals can be solved analytically in the socalled ultrarelativistic regime (m φ , M ≪ k ∼ πT ). On the other hand, a numerical evaluation with general masses and momenta is more straightforward by setting Λ → 1. With the latter choice, the magnitude of eq. (5.11) will be illustrated in fig. 3.

2 ↔ 2 scatterings off soft bosons
Another possible source for an IR divergence in 2 ↔ 2 scatterings is when one of the external scatterers becomes soft. Technically, this is the case when the argument of a bosonic thermal distribution vanishes in eq. (4.18). Let us recall that chemical potentials associated with bosons should be smaller than the particle masses, otherwise we are driven to Bose-Einstein condensation. In cosmology, chemical potentials are in any case small compared with the temperature. In most computations it is then sufficient to expand to first order in chemical potentials. In this case, the domain of a small bosonic energy exhibits a quadratic pole, n ′ B (q 0 ) = −βn B (q 0 )[1 + n B (q 0 )] ≈ −T /q 2 0 for q 0 ≪ T , leading to a logarithmic divergence if the mass has been omitted. Such a divergence can arise for instance from the domain q 0 ∼ ω inside the thermal distribution n σ 2 in eq. (4.18) [35]. It has also been pointed out that a similar effect can arise even without chemical potentials, if we consider virtual corrections to 1 ↔ 2 scatterings and pick up an enhancement ∼ (T /q 0 ) 2 from two Bose distributions [36]. Both of these effects are related to kinematic endpoints for scattering off soft Higgs bosons, and can be regularized simply by keeping the thermal Higgs mass finite.
In the presence of chemical potentials, some care is needed with virtual lines as well. For instance, as discussed below eq. (4.18), q 0 may cross zero despite the particle having a mass, and then n τ 1 σ 1 has a would-be pole at q 0 = µ 1 − ν 1 if τ 1 σ 1 = +, which is lifted by the second factor. If we expand in chemical potentials, the pole is of second order, but it is still lifted. In order to account for this, it is helpful to make use of the representation given below eq. (4.18).

Small-angle 1 + n ↔ 2 + n reactions
Even if at first sight simpler than 2 ↔ 2 and 1 ↔ 3 scatterings, it turns out that the IR structure of 1 ↔ 2 reactions is more complicated, as the phase space is strongly constrained by masses. This implies that a subclass of 1 + n ↔ 2 + n reactions, with n ≥ 1, can be as important as the 1 ↔ 2 reactions, and needs to be incorporated through Landau-Pomeranchuk-Migdal (LPM) resummation [1,3,37]. The procedure requires the inclusion of "asymptotic", i.e. large-momentum thermal mass corrections [33]. The full information cannot be easily deduced from a matrix element squared like eq. (2.14), and therefore requires an effort specific to the application in question. Profitting from ideas put forward in ref. [30], LPM-resummed results have been extended up to NLO in some cases [38][39][40].
Similarly to sec. 5.1, LPM resummation requires a subtraction-addition step, manifestly avoiding the danger of double counting. It is important here that as LPM resummation is normally implemented after making use of kinematic simplifications pertinent to the ultrarelativistic regime, the subtraction should adhere to the same simplifications, guaranteeing that the resummation only has a negligible (higher-order) effect in domains where it is not justified. It is an open problem, deserving further study, how LPM resummation could be smoothly connected to kinematic domains beyond the ultrarelativistic one.
For the concrete example of fig. 1, even though we do not discuss LPM resummation itself here, for the reason just mentioned, the subtraction term can be deduced from the same HTL computation that led to eq. (5.11). The difference to eq. (5.11) is that in 1 ↔ 2 reactions, the lepton is on-shell and timelike, and therefore we need to pick up the pole part from eq. (5.9). 10 Subsequently we again need to re-expandρ µ up to O(g 2 T 2 ) ∼ O(m 2 ℓT ) to obtain the terms to be subtracted. However, great care is needed here: in eq. (2.37), we have reconstructed only a part of the virtual corrections to 1 ↔ 2 scatterings, namely those that have a counterpart on the side of 2 ↔ 2 and 1 ↔ 3 scatterings. This subset does not include the corrections that modify the lepton mass but do nothing else. Therefore, only terms from the numerator of eq. (5.9) play a role in the re-expansion. Renaming the mass in the denominator to m 2 ℓ , which is viewed as an IR regulator like in eq. (5.11), this leads to (5.14) Here the integration boundaries and angle are from eqs. (4.20) and (4.21), and the expression in curly brackets represents an angular average of the type in eq. (C.1). Eq. (5.14) is closely related to and indeed subtracts most of the first term on the first line of eq. (2.37). The numerical magnitude of eq. (5.14) will be illustrated in fig. 3. We note that the IR divergences of eqs. (5.11) and (5.14) cancel against each other when ml → 0 [36]; the IR divergence of eq. (5.14) originates from the region of small ǫl, where the phase space distributions in scat 1 ↔ 2 (l, φ) take the same form as on the right-hand side of eq. (5.13).

How about real intermediate states?
One issue frequently discussed in the Boltzmann equation literature is that of "real intermediate states" (cf., e.g., ref. [41]). Consider a reaction in which the Mandelstam variable entering a line can coincide with the mass-squared of that particle. This case emerges, for instance, if we consider the process in fig. 1(left), and let the scalar decay into a tt pair, assuming that with thermal modifications it could be possible to have M − m ℓ > m φ > 2m t . The general conditions for the appearance of real intermediate states were listed in eqs. (3.6)- (3.11).
If the propagator in question appears quadratically in the cross section, the thermal phase space average looks potentially divergent. In our procedure, as specified in sec. 3, this is regularized by treating the quadratic propagator as a mass derivative of a principle value, rendering the average integrable. There is necessarily also a 1 ↔ 2 reaction in which the resonant intermediate state appears as a real particle. Virtual corrections to these 1 ↔ 2 processes (a closed tt loop in the above example) reflect the same singularity. Specifically, after angular averaging, the remaining energy integrals, from eqs. (4.19) and (4.25), cross a hypersurface where the argument of a logarithm vanishes. Even though this is again integrable, an efficient numerical procedure normally requires dividing the integration domain into subregions, so that the singularities appear at their boundaries.
To summarize, real intermediate states (or, in the language of sec. 3, actual poles) in the matrix elements squared lead to no IR singularities in our framework. However, their appearance makes it more tedious to render the practical integrations efficient.

Description of algebraic and numerical procedures
Let us summarize in a procedural fashion an algorithm implementing the ingredients introduced in secs. 2 and 4.
The general philosophy is to start from a 1 → 3 amplitude, like in eq. (2.14), describing how a "non-equilibrium" particle decays into "plasma" particles (the process does not need to be kinematically allowed in practice). Let us stress again that even if we refer to a decay, inverse processes are always included as well. The non-equilibrium particle represents a slow variable which may fall out of equilibrium or never enter it in the first place, whereas the plasma particles are fast ones, with their density matrices fully characterized by a temperature and a handful of chemical potentials. There are really two algorithms, an algebraic [a] and a numerical one [n]. Their ingredients are: • input parameters describing plasma particles [n]: These parameters include the temperature; various chemical potentials (lepton and baryon asymmetries, gauge field backgrounds); effective couplings and masses like g 1 , g 2 , ml, mφ in eq. (2.14), which may incorporate "hard loop corrections" from fast processes, whereby these are in general functions of the temperature; masses m a , chemical potentials µ a , and statistics σ a , associated with the final state of the 1 → 3 process; a specification for negative index choices, according to eq. (A.2).
• input parameters describing the non-equilibrium particle [n]: Momentum k, mass M , and energy ω = √ k 2 + M 2 ; polarization state, e.g. through the four-momentum E in eq. (2.14); Standard Model quantum numbers of the vertex to which the non-equilibrium particle couples, e.g. through its possible dependence on the lepton flavour in eq. (2.14).
• reflection of 1 → 3 decay rate into other channels [a]: The full (Born-level) interaction rate, incorporating the other channels and inverse processes, is obtained from the 1 → 3 decay rate according to eq. (2.12). Subsequently, we can make use of permutations of momenta, to label the four-momenta according to the conventions that we have adopted. For the 1 → 3 process there are 3! = 6 permutations, whereas for the 2 ↔ 2 and 3 → 1 cases the identification of one initialor one final-state particle, respectively, leaves over a two-fold freedom that can be used to rename variables. An automated implementation of this step requires a software with wildcard pattern matching capabilities, such as form [42].

• identification of virtual corrections [a]:
The 1 → 3 rate permits for the automatic identification of IR sensitive virtual corrections to 1 ↔ 2 processes, by searching for poles in the matrix elements squared, determining the corresponding residues, and attaching these to the proper thermal average (cf. sec. 2.2).
• integration measure and thermal distributions for given channel [a]: Proceeding to the integration, two angles are fixed by energy-momentum conservation, e.g. from eq. (4.7); all kinematic invariants are fixed in terms of the chosen integration variables, e.g. from eq. (4.8); thermal distributions can be represented in a factorized form, e.g. from eq. (4.18); integration measure can be inserted, e.g. from eq. (4.14).

• azimuthal average [a]:
A key ingredient of the parametrization is that one the remaining integrals, over the azimuthal angle, can be carried out analytically. The dependence on the azimuthal angle originates through a single scalar product, for instance k · p 1 in eq. (4.10), which does not appear in thermal distributions. After the full expression has been partial fractioned into terms containing positive powers of this variable, or inverse powers of first-order polynomials, the azimuthal average can be inserted, from eqs. (4.10), (4.11). For virtual corrections, the full angular averages can be resolved, from appendix C.
• phase space integrals [n]: The remaining at most three-dimensional integration can be carried out numerically, e.g. from eq. (4.14). In virtual corrections the integration is normally two-dimensional.

• tests [n]:
The procedure contains some redundancies, which permit to crosscheck for its correct implementation. For 2 ↔ 2 and 3 → 1 processes, if we override the permutations used to select the "optimal" parametrization, we should get the same result for both t-and s-channel parametrizations. For 1 ↔ 3 processes, all 3! relabellings of the final-state momenta should lead to the same result.
Let us conclude by showing example plots for the case of eq. (2.14), for the plasma input parameters g 1 = 1/3, g 2 = 2/3, µ L = 10 −3 T , µ Y = 2 × 10 −2 T , where µ Y refers to the hypercharge chemical potential. The properties of the decay products are set as The auxiliary masses are set to ml = m ℓ and mφ = m φ . For the non-equilibrium particle we set M = 0.3T or M = 3T , choosing the helicity projection E = K or E = U . Recalling our convention of referring to the side of the non-equilibrium particle as the initial state, 1 → 3 channels are open for the mass M = 3T , and 3 → 1 channels for M = 0.3T . The 2 ↔ 2 scatterings are allowed in any case. The results are shown in fig. 2, including for comparison also the results originating from 1 ↔ 2 processes. For further illustration, in fig. 3 we show the influence of thermal resummations from eqs. (5.11) and (5.14), whose effect is to replace the mass ml by the proper thermal lepton mass m ℓT . Our final results are obtained by summing together figs. 2 and 3; the outcome is shown in fig. 4, this time as a function of M/T , for various fixed momenta k/T . The conclusions that we draw from these plots are discussed in sec. 7.
Finally, even if we are able to produce accurate results for several test cases, it is appropriate to acknowledge that numerical integrations become less efficient in certain limits. Roughly speaking, the integrations are simple if all masses and momenta are of order T , whereas large scale hierarchies are challenging to handle. Examples are the non-relativistic limit M ≫ T , where a huge cancellation between real and virtual corrections demands exquisite numerical precision, and the ultrarelativistic regime M ≪ T , particularly with k ≫ T , where the literal integration ranges can be broad but the integrands are strongly localized. Moreover, particle spectra leading to real poles can be costly, if principal value integrations are regularized numerically rather than analytically (cf. secs. 3 and 5.4). We cannot exclude numerical inaccuracies if parameters are pushed to domains which happen to have eluded our tests.

Conclusions and outlook
We have described a method to represent and evaluate thermal 2 ↔ 2 and 1 ↔ 3 scattering rates, including a way to regularize and subtract the poles that appear in the matrix elements squared. Choosing a language in which the side of the non-equilibrium particle is called the 2↔2,1↔3 , with the matrix element from eq. (2.14), and ∆Γ Born 1↔2 , from eq. (2.37). The sum of these two is compared with Γ Born 1↔2 , with the matrix element Θ(P ℓ , P φ ) = 4E · P ℓ . The mass is M = 0.3T (left) and M = 3T (right), with the other parameters explained at the end of sec. 6. The top row shows results for E = K (normalizing to T 2 ), the bottom row for E = U (normalizing to T ). The results are for m γ = 0.01T , noting that real and virtual corrections separately depend strongly on this IR regulator, however the full result (solid line) is independent of it. 2↔2 , from eq. (5.11), and ∆Γ HTL 1↔2 , from eq. (5.14), for M = 0.3T (left) and M = 3T (right). The top row shows results for E = K (normalizing to T 2 ), the bottom row for E = U (normalizing to T ), with the other parameters explained at the end of sec. 6. The purpose of these corrections is to replace the auxiliary mass ml = 0.1T through the physical thermal lepton mass m ℓT ≈ 0.3T , which parametrizes HTL propagators. The 2 ↔ 2 and 1 ↔ 2 corrections depend strongly on the IR regulator ml, however the full result (solid line) is almost independent of it, as long as we stay away from the constrained domain M ∈ (m φ − ml, m φ + ml).  2↔2,1↔3 + ∆Γ Born 1↔2 and the HTL subtraction-addition contribution ∆Γ HTL 2↔2,1↔2 , as the latter cannot remove threshold singularities that originate from non-HTL structures, particularly the virtual doublepole corrections (last line of eq. (2.37)). The spikes could be eliminated by sending ml → 0, whereby the physical thermal lepton mass m ℓT ≈ 0.3T is reinstated by the HTL contribution. initial state (even if both processes and inverse processes are always included), the idea is to give as input a 1 → 3 matrix element squared, which displays maximal symmetries, as all thermalized particles are in the final state (the process does not need to be kinematically allowed). The 2 ↔ 2 and 3 → 1 scattering rates are obtained by crossing relations, and the virtual corrections to 1 ↔ 2 rates that cancel their poles are automatically identified (cf. sec. 2.2). Vacuum contributions can be pulled apart, so that the final step is to carry out an exponentially convergent three-dimensional numerical integral. The results have been worked out for general chemical potentials, and can therefore be applied not only to cosmology but, potentially, to dense astrophysical environments as well.
In our framework, the full interaction rate can be represented as The first part, Γ Born 2↔2,1↔3 , captures 2 ↔ 2 and 1 ↔ 3 rates, and the second, ∆Γ Born 1↔2 , the virtual corrections that cancel their poles. The third part, ∆Γ HTL 2↔2,1↔2 , represents a subtractionaddition step implementing HTL resummation (cf. secs. 5.1 and 5.3), which replaces auxiliary masses, used as an intermediate IR regulator, through physical thermal masses. The last part, Γ LPM 1+n↔2+n , with n ≥ 0, sums together 1 + n ↔ 2 + n processes. In the present paper, Γ LPM 1+n↔2+n has been approximated through its lowest-order term (n = 0), Γ Born 1↔2 , given that the corresponding formalism has not been generalized to the case that some particles cease to be ultrarelativistic (in figs. 2-4 we have chosen m φ = T to be relatively "heavy", whereby restricting to n = 0 should be a fair approximation).
The numerical importance of the virtual corrections can be appreciated from fig. 2(right). Due to the presence of logarithmic and double-logarithmic IR divergences [17], the inclusion of only 2 ↔ 2 scatterings would overestimate the correct result by a factor ∼ 10 3 for these parameters. Moreover, 1 ↔ 3 rates, which are often overlooked, play an equally important role as 2 ↔ 2 scatterings. When virtual corrections are included, the result becomes much smaller than that from 1 ↔ 2 processes, and can be safely omitted. This conclusion is not changed by HTL resummation, as illustrated in fig. 3(right).
The situation is very different for the parameters in fig. 2(left), where the non-equilibrium particle is ultrarelativistic (M = 0.3T ). Even though there is still a substantial cancellation between 2 ↔ 2 and 3 ↔ 1 rates and virtual corrections, the remainder is now larger than that from the 1 ↔ 2 processes, provided that k > ∼ (a few) × T . This conclusion is not changed by HTL resummation, whose influence is smaller at similar momenta, as is illustrated in fig. 3(left). Furthermore the conclusion turns out to be strengthened if m φ is reduced towards its physical high-temperature value, m φ ≃ 0.4T .
To summarize, if the non-equilibrium particle is ultrarelativistic (M ≪ T ), it is essential for quantitative investigations to include 2 ↔ 2 and 1 ↔ 3 rates as well as virtual corrections to 1 ↔ 2 ones. If the non-equilibrium particle is non-relativistic (M ≫ T ), it would be dangerous to incorporate 2 ↔ 2 or 1 ↔ 3 rates, without a full account of the virtual corrections to 1 ↔ 2 processes that may cancel most of the result. 11 A c code implementing the numerical parts of our procedure, and a form code implementing the algebraic ones, are attached as ancillary files to this paper. Even if their details are specific to the example in eq. (2.14), relevant for leptogenesis scenarios, the structures and main steps are quite general. Therefore we hope that they can be applied to other problems as well, for instance in the context of freeze-in dark matter production, where the need for determining 2 ↔ 2 and 1 ↔ 3 scattering rates has been underlined recently [12].

A. Proof of thermal crossing relations
A.1. 2 ↔ 2 and 1 ↔ 3 real corrections In order to handle any 2 ↔ 2 or 1 ↔ 3 reaction, we start by introducing the concept of a "master" sum-integral. The master sum-integral originates by viewing the interaction rate as an imaginary part ("cut") of a retarded correlator, in analogy with the optical theorem. The retarded correlator can in turn be represented as an analytic continuation of an imaginarytime (Euclidean) correlator. This representation, even if sounding formal, is quite helpful, as the imaginary-time expression automatically encodes many independent reactions as well as the crossing symmetries between them. In particular, each such master structure is IR finite by itself, containing no poles in accordance with the KLN theorem [23,24].
As can readily be verified pictorially, 2 ↔ 2 and 1 ↔ 3 processes correspond to cuts of 2-loop diagrams. Any 2-loop contribution can, in turn, be represented as a linear combination of master sum-integrals. Inspired by ref. [43], we define a 2-loop master sum-integral as . (A.1) 11 The latter statement applies as such to a mass spectrum similar to that in our example, where the nonequilibrium particle can be heavier than the plasma particles and experiences only inelastic reactions. If it can participate in elastic scatterings, or if one of the plasma particles has a mass close to that of the non-equilibrium one, 2 ↔ 2 reactions can dominate even in the non-relativistic regime.
Here Σ P ≡ T pn p is a Matsubara sum-integral, with p n referring to a bosonic or fermionic Matsubara frequency, and P ≡ (p n , p). The imaginary-time external four-momentum is denoted by K, the corresponding Minkowskian four-momentum by K. The indices i x , j x are integers, whereas the a x label particle species (a x ∈ {h, W, Z, Q, ν a , e a , u, d, ...}). Masses and chemical potentials appear through where µ ax is the chemical potential, m ax is the mass, and p ≡ |p|. This notation implies that ∆ −P ;ax = ∆ P ;−ax , where −a x labels an antiparticle. Momenta denoted byP ≡ (p n + iµ ax , p) include a shift by the chemical potential. The H factors in the numerator denote helicity projections, for instance for spin-1 2 particles where E is some external four-momentum, e.g. K or the medium four-velocity U . For spin-1 fields the projectors are quadratic in momenta, for instance H iP = p 2 − (p · k) 2 /k 2 for the sum over transverse polarization of on-shell photons. As the helicity projection is a linear operation, we have assumed a linear dependence on H in eq. (A.1); this is a simplification, even if we believe that the result holds more generally.
Once the Matsubara sums are carried out, eq. (A.1) contains spatial momentum integrals, weighted by Bose-Einstein and Fermi-Dirac distribution functions. We now analytically continuek n to a Minkowskian frequency, and take the cut. The complete cut involves virtual corrections to 1 ↔ 2 scatterings, as well as real processes, namely 2 ↔ 2 and 1 ↔ 3 scatterings. We first focus on the real processes, deferring virtual corrections to appendix A.2.
Consider, for instance, a cut illustrated as .
Finally, consider the other diagonal cut, going through the lines a 1 , a 3 , a 5 . By substituting P ↔ Q in eq. (A.1), the result can be directly obtained from eq. (A.4), just by inverting a 3 → −a 3 and setting P 3 → −P 3 in the numerator. This confirms eq. (2.34).

A.2. 1 ↔ 2 virtual corrections
In analogy with the 2-loop master in eq. (A.1), we define a 1-loop master sum-integral as where the propagator structures are defined according to eq. (A.2). As originally demonstrated with scalar field theory [44], the imaginary part of eq. (A.5) can be expressed in terms of phase space integrals, in particular Here we have set i 1 = i 2 = 1; results for higher powers can be obtained by taking derivatives with respect to the masses m 2 a 1 and m 2 a 2 . This motivates eq. (2.26). In order to verify eq. (2.28), we may start from eq. (A.4) but choose, for instance, i 5 = 0, eliminating one of the poles. Simplifying other index choices as well, the part corresponding to 2 ↔ 2 and 1 ↔ 3 processes then amounts to Im I 0j 2 j 3 j 4 0j 6 1 1 1 1 0 0 (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 with the other channels following by crossings according to eq. (A.4). We note that this channel has a pole of the type in eq. (2.27), with a residue −1. It is then sufficient to work out the corresponding virtual contributions, i.e. the ones in which the lines a 1 and a 4 are cut. This yields precisely the structure in eq. (2.28), with the same overall −1.
Finally, for eq. (2.35), we need to consider two possibilities, cutting the lines a 1 , a 4 on one hand, and a 2 , a 5 on the other. According to eq. (A.1), the results can be related to each other through the substitution P ↔ Q. In any case, eq. (2.35) can be confirmed.

B. Further thermal averages
In this appendix we supplement the procedure described in sec. 4.1 for t-channel 2 ↔ 2 scatterings, by working out thermal averages for the other channels appearing in eq. (2.12).

B.1. 1 → 3 reactions
For 1 → 3 reactions, all particles (apart from the "external" one, carrying the momentum K) appear on one side and can thus be interchanged as far as momentum labellings are concerned. Let us use this freedom to choose s 12 = (P 1 + P 2 ) 2 as a potentially IR sensitive Mandelstam variable. If another Mandelstam variable appears, it can be chosen as s 23 through further permutations. For the example in eq. (2.14), with the momenta for scat 1 → 3 (a 1 , a 2 , a 3 ) labelled as P i ≡ P a i , this gives The goal now is to have s 12 as an integration variable. To this end we introduce a fourmomentum Q such that s 12 = Q 2 , and write the phase space integration measure from eq. (2.3) as where we have integrated over p 1 and p 3 . The Dirac-δ's fix two angles as The other Mandelstam variables can be expressed as As discussed below eq. (4.8), the dependence on k · p 2 should be partial fractioned. The azimuthal average of the angle between k and p 2 can be worked out like in eqs. (4.10)-(4.11), with the exchange p 1 ↔ p 2 . Resolving the energy-conservation constraints in eq. (B.2), and replacing subsequently q through s 12 = q 2 0 − q 2 , the integration measure for the 1 → 3 channel becomes where q = q 2 0 − s 12 and, making use of √ λ from eq. (4.17), The phase space distribution associated with 1 → 3 scatterings, cf. eq. (2.4), is conveniently factorized as in eq. (2.15), which after the insertion of q 0 = ǫ 1 + ǫ 2 = ω − ǫ 3 from eq. (B.2), as well as the employment of n σ (−ǫ) = −1 − n σ (ǫ), yields Integrals over ǫ 2 are exponentially convergent at large ǫ 2 ; those over q 0 are localized close to the lower bound q − 0 . In the massless limit, the integration range in eq. (B.5) collapses to a point. Therefore thermal IR divergences of the type discussed in sec. 5.1 are absent in 1 → 3 decays.

B.2. s-channel 2 ↔ 2 reactions
The part of 2 ↔ 2 reactions that cannot be put in the form treated in sec. 4.1 are s-channels reactions. With the labelling for scat 2 → 2 (−a 1 ; b 1 , b 2 ) chosen as , so that the sign flips necessary for initial-state momenta and chemical potentials are already explicit when using K 1 and ν 1 , the s-channel part originating from eq. (2.14) becomes Γ Born 2↔2(s) The goal now is to have s as an integration variable. For this aim we introduce a fourmomentum Q such that s = Q 2 , and define s-channel parametrization of the integration measure from eq. (2.6) as where we have integrated over p 1 and k 1 . The Dirac-δ's fix two angles as The other Mandelstam variables can be expressed as The azimuthal average of the angle between k and p 2 can be worked out like in eqs. (4.10)-(4.11), with the exchange p 1 ↔ p 2 . Resolving the energy-conservation constraints in eq. (B.10), the integration measure becomes where ǫ ± 2 are from eq. (B.16). Subsequently we replace q through s = q 2 0 − q 2 , resulting in where q = q 2 0 − s and When using the s-channel parametrization, the phase space distribution associated with 2 ↔ 2 scatterings from eq. (2.7) is conveniently factorized as in eq. (2.17), which after the insertion of q 0 = ǫ 1 + ǫ 2 = ω + E 1 from eq. (B.10) yields The latter factor guarantees that integrals over ǫ 2 are exponentially convergent at large ǫ 2 ; those over q 0 are localized close to the lower bound q − 0 . In the massless limit, the integration domain in eq. (B.13) becomes Therefore q 0 is never small, and q could be small only in the vicinity of q 0 = 2k, however around that point s ≈ 4k 2 is large. Therefore s-channel 2 ↔ 2 scatterings do not lead to IR divergences from small q 0 , q, of the type that were discussed in sec. 5.1.
In contrast, if we expand n τ 1 from eq. (B.17) to first order in ν 1 , a second order pole emerges. Thus a logarithmic divergence can originate from the domain around q 0 ≈ ω, as outlined in sec. 5.2, which is regularized by finite masses [35]. Likewise the final-state distributions n σ 1 , n σ 2 can become singular if the particles are massless bosons and they carry finite chemical potentials.

B.3. t-channel 3 → 1 reactions
Part of the 3 → 1 reactions from eq. (2.12) can be called t-channel ones, with the Mandelstam variables defined according to eq. (4.4). For eq. (2.14), with momenta and chemical potentials labelled for scat 1 , so that the sign flips necessary for initial-state momenta and chemical potentials are already included when using K i and ν i , the t-channel processes amount to Γ Born Some of these channels are not allowed kinematically, but this is taken care of automatically once we work out the integration measure, cf. eq. (B.23) below.
The goal now is to have t as an integration variable. To achieve this we introduce a fourmomentum Q such that t = Q 2 , and write the phase space integration measure for 3 → 1 reactions from eq. (2.9) as (B.20) The Dirac-δ's fix two angles as The other Mandelstam variables can be expressed as u = M 2 + M 2 1 − 2(k · k 1 − ωE 1(k 1 ) ) , s = m 2 1 + M 2 2 − t + 2(k · k 1 − ωE 1(k 1 ) ) .

(B.26)
Integrals over E 1 are exponentially convergent at large E 1 ; those over q 0 are localized close to the upper bound q + 0 . In the massless limit, the integration domain of eq. (B.23) shrinks to a point. Thus there are no thermal IR problems of the type discussed in sec. 5.1.

(B.34)
Integrals over E 2 are exponentially convergent at large E 2 ; those over q 0 are localized close to the lower bound q − 0 . In the massless limit, the integration domain of eq. (B.31) shrinks to a point. Thus there are no thermal IR problems of the type discussed in sec. 5.1.

C. Angular averages for virtual corrections
We define here angular averages that appear in the virtual corrections discussed in sec. 4.2.
Let (θ, ϕ) be the spherical angles associated with a loop momentum, which in the following we denote by p a . The axis with respect to which θ is measured will be specified later on, but the choice plays no role, as all angles are integrated over. The angular average is taken in the presence of two further vectors, denoted by p d and k, and sometimes it is also convenient to employ p e ≡ k − p d . We are concerned with two types of averages, G pa;p d ,k;z P n (p a · k) ≡ where P n , Q n are polynomials of degree n, and È denotes the principal value.
Starting with the latter average, the first step is to write Q n (p a · k) = Q n (z 2 ) + Q n (p a · k) − Q n (z 2 ) ≡ Q n (z 2 ) + (z 2 − p a · k) Q n−1 (p a · k) .

(C.3)
If we express the original polynomial as Q n (p a · k) = n i=0 a i (p a · k) i , then Q n−1 (p a · k) = n−1 j=0 b j (p a · k) j , where the coefficients read b j = − n i=j+1 a i z i−j−1 2 . Thereby eq. (C.2) becomes H pa;p d ,k;z 1 ,z 2 Q n (p a · k) = H pa;p d ,k;z 1 ,z 2 Q n (z 2 ) + G pa;p d ,k;z 1 Q n−1 (p a · k) . (C.4) For the first term, we combine the denominators with a Feynman parameter, 1 (z 1 − p a · p d )(z 2 − p a · k) = 1 0 ds 1 [s z 1 + (1 − s) z 2 − p a · (k − s p e )] 2 , (C.5) where we made use of p e = k − p d . The angle is now chosen as θ ≡ θ pa,k−s p e . There is no dependence on ϕ, so that both integrals are readily carried out, (C.6) The denominator is a second order polynomial in s, and the integral over s yields where the value of A α,β,γ is given by (α, β, γ ∈ Ê) Let us then turn to the average G from eq. (C.1). The idea here is to choose the angle θ pa,p d to play the role of θ. The other scalar product can be expressed as p a · k = p a k cos θ pa,p d cos θ k,p d + sin θ pa,p d sin θ k,p d cos ϕ , (C.10) where θ k,p d is fixed by the outer integral (cf. eq. (4.21)). Inserting this as the argument of P n yields a polynomial in cos ϕ. The azimuthal averages can now be carried out, . (C.11) As only even powers contribute, the dependence on sin θ pa,p d is quadratic, and can be expressed in terms of cos 2 θ pa,p d . It is convenient to write cos θ pa,p d = p a · p d /(p a p d ). Thereby P n (p a · k) ϕ = n k=0 c k (p a · p d ) k ≡ R n (p a · p d ) , (C.12) and eq. (C.1) takes the form G pa;p d ,k;z P n (p a · k) = +1 −1 d cos θ pa,p d 2 To carry out the remaining integral, we repeat the logic of eqs. (C.3) and (C.4). We write R n (p a · p d ) = R n (z) + (z − p a · p d ) R n−1 (p a · p d ) , (C.14) where R n−1 (p a · p d ) = n−1 l=0 d l (p a · p d ) l , with the coefficients d l = − n k=l+1 c k z k−l−1 . Then G pa;p d ,k;z P n (p a · k) = +1 −1 d cos θ pa,p d 2 È R n (z) z − p a p d cos θ pa,p d + n−1 l=0 d l p l a p l d cos l θ pa,p d .
(C.15) Both parts are easily integrated, yielding G pa;p d ,k;z P n (p a · k) = R n (z) 2p a p d ln z + p a p d z − p a p d +