Relic density computations at NLO: infrared finiteness and thermal correction

There is an increasing interest in accurate dark matter relic density predictions, which requires next-to-leading order (NLO) calculations. The method applied up to now uses zero-temperature NLO calculations of annihilation cross sections in the standard Boltzmann equation for freeze-out, and is conceptually problematic, since it ignores the finite-temperature infrared (IR) divergences from soft and collinear radiation and virtual effects. We address this problem systematically by starting from non-equilibrium quantum field theory, and demonstrate on a realistic model that soft and collinear temperature-dependent divergences cancel in the collision term. Our analysis provides justification for the use of the freeze-out equation in its conventional form and determines the leading finite-temperature correction to the annihilation cross section. This turns out to have a remarkably simple structure.


Introduction
The most widely studied and arguably most natural mechanism of generating the presentday abundance of dark matter (DM) is its thermal production in the early Universe followed by chemical decoupling (freeze-out) of the DM particles from the background plasma. For temperatures much higher than the mass of the DM particles the dark matter component remains in both chemical and kinetic equilibrium. When the temperature of the plasma drops, the interactions are not strong enough to keep the dark matter component in chemical equilibrium and the DM particle number freezes out.
The precise moment when chemical decoupling happens is determined by two different physical processes: the expansion of the Universe governed by the Hubble rate H and the annihilation rate Γ of the DM particles. The precise description of the decoupling process is possible in kinetic theory, where the evolution of the system is given by the transport equations for the phase-space distribution functions f (p). If one assumes that i) the Compton wavelength of DM particles is small with respect to inhomogeneity scale and ii) one can adopt the quasi-particle approximation, one arrives at a semi-classical description. In this case the transport is governed by the Boltzmann equation and its solution can be used for the determination of the DM relic density for a given particle physics model. The DM particle number density then follows the simple equation dn χ dt + 3Hn χ = σ χχ→ij v rel n eq χ n eq χ − n χ nχ , (1.1) where H denotes the Hubble rate, and σ χχ→ij v rel the thermal average of the sum over all annihilation cross sections to two-particle final states ij.
In recent years there has been an increasing interest in higher-order corrections to scattering and annihilation processes involving DM particles. The main phenomenological importance of such corrections is related to the modification of the annihilation spectra relevant for the indirect searches. Quite generally, the increasing precision of dark matter observations will require more accurate computations of the scattering and annihilation processes, in some cases at full next-to-leading order (NLO) in the coupling constant. In particular, it has also been noted recently that corrections to the annihilation rate can affect non-negligibly the relic density computation [1][2][3][4][5][6][7][8]. With this in mind the first numerical codes including the higher-order corrections are being developed, SloopS [9][10][11] and DM@NLO [12,13]. What is usually done is to compute the virtual and real radiation corrections to the two-particle processes χχ → ij using standard quantum field theory methods at zero temperature.
This procedure raises a number of questions, especially for relic density computations, since freeze-out occurs when the temperature of the Universe is small, but non-negligible compared to the DM particle mass.
• Why should the time evolution of n χ be described by inclusive two-particle cross sections and a Boltzmann equation of the form applicable to 2 → 2 reactions? The real radiation amplitude involves three-particle final states, typically containing an additional photon or gluon, which are themselves abundant in the plasma. Moreover, absorption processes exist, but are neglected in the computation.
• How do the soft and collinear infrared (IR) divergences cancel at finite temperature? It is well-known that IR divergences are more severe at finite temperature due to the enhancement from the Bose distribution at small momenta. Moreover, virtual and real scattering matrix elements, which are separately divergent, appear to be multiplied by different statistical factors in the full Boltzmann equation.
• Assuming IR finiteness can be shown, what are the leading finite-temperature effects on the annihilation cross sections and the relic density?
• Does the transport equation itself receive quantum corrections when it is derived from general principles of non-equilibrium quantum field theory (QFT) to NLO accuracy?
In this paper we address these questions. We demonstrate that the IR cancellation happens in the sum over "cuts" of individual self-energy diagrams similar to the situation at zero temperature, but involving the additional processes that occur in the plasma. The form of (1.1) remains valid under the typical conditions of DM freeze-out, but the annihilation cross section is modified by a small and calculable finite-temperature correction. Remarkably, the final NLO finite-temperature correction has a very simple structure and can be computed directly from the zero-temperature tree level cross section. The outline of the paper is as follows. In section 2 we discuss the IR problems that arise at NLO at finite temperature. Section 3 reviews the well-known derivation of the Boltzmann equation from non-equilibrium QFT, with emphasis on the application to the freeze-out process. Next we discuss the computation of the collision term in section 4 and demonstrate the general procedure for a bino-like DM model. In section 5 we present and discuss the IR divergences cancellation and the result for the finite correction from thermal effects. We conclude in section 6.

IR divergences and the Boltzmann equation
In general the one-loop scattering amplitudes contain soft and collinear IR divergent terms. At zero temperature the Bloch-Nordsieck cancellations and Kinoshita-Lee-Nauenberg (KLN) theorem [14][15][16] ensure that physical observables are free of both of these divergences, as they involve summation over initial and final degenerate states, in the sense of inclusiveness or experimental resolution in energy and angles. At finite temperature no general proof of such a theorem is known, nevertheless the cancellation was observed in all the particular cases studied in the literature (see e.g. [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]). In both situations the main prerequisite for the cancellation is the inclusion of a soft or collinear gauge boson in the NLO computation.
Whether and how the cancellation happens in NLO relic density computations has not yet been investigated. 1 To formulate the problem, let us consider a typical freeze-out scenario of a weakly-interacting massive dark matter particle (WIMP), for which thermal and chemical equilibrium at temperatures larger than the freeze-out temperature T f ≈ m χ /20 is maintained by 2 → 2 scattering processes at leading order. In the Friedmann-Robertson-Walker background the semi-classical Boltzmann equation for the evolution of the phase-space distribution function reads (2.1) It can be rewritten as an equation for the number density where h i is the number of internal degrees of freedom of particle i and we assume that the thermal plasma is "unpolarized" with respect to the internal degrees of freedom. The integrated collision term for the leading-order (LO) annihilation (production) process χχ ↔ ij is given by with p f (p i ) denoting the sum of the final (initial) state momenta of the process, and moreover assumed CP invariance, which implies |M χχ→ij | 2 = |M ij→χχ | 2 . The ± signs are chosen according to whether the particle is a boson (+) or fermion (−). Note that here the squared matrix elements are defined to be summed over the internal (spin) degrees of freedom of both the initial-state DM and final-state Standard Model (SM) particles. We shall make the standard assumption that the SM particles are kept in thermal equilibrium by frequent scatterings and that asymmetries are negligible. Therefore f i,j = f eq i,j , and also the photon phase-space distribution f γ = f eq γ introduced below, are either Bose-Einstein or Fermi-Dirac distributions with vanishing chemical potentials. For the DM particles we assume that kinetic equilibrium is maintained by frequent elastic scatterings with the particles from the thermal bath, resulting in f χ ∝ f eq χ , where the factor of proportionality depends on temperature but not on the energy and momentum. Chemical equilibrium of the DM particles, however, is lost when the temperature falls below the freeze-out temperature T f , and DM particle-number changing processes occur at insufficient rates. Since in the 2 → 2 annihilation reaction all energies are of order O(m χ ), and since T f ≪ m χ , all distribution functions are exponentially suppressed, and we can approximate 1 ± f ≈ 1. Under these assumptions the integrated collision term (2.3) takes the standard form denotes the thermally averaged cross section times velocity and we used f χ = n χ /n eq χ × f eq χ . As long as the amplitudes are computed at tree level, the cross section σ χχ→ij v rel and hence the collision term is evidently IR finite. When the relic density computation described above is extended to NLO, what has been done up to now is to compute the zerotemperature annihilation cross section to NLO, while keeping the form of the remaining equations. This involves the one-loop correction to the 2 → 2 annihilation processes χχ → ij, and the tree-level radiation process χχ → ijγ, assuming electromagnetic radiation for definiteness. The sum is IR finite by the usual zero-temperature IR cancellations.
This procedure is conceptually problematic. To see this, consider the NLO collision term where again we used 1 ± f ≈ 1 except for the photon distribution function. The collision term (2.8) contains both annihilation and production contributions, which are however symmetric and can be described by the same thermally averaged cross section, as long as the theory is CP invariant and the DM particles are in kinetic equilibrium. It can be most easily seen by making use of the detailed balance relation for the photon distribution function the Maxwell approximation for the remaining ones and the energy conservation. It follows, that the collision term has the form analogous to (2.5) but with the thermally averaged cross section replaced by with the interpretation dΠ χχij dΠ γ = dΠ χχijγ . For this reason also when discussing the NLO corrections we will consider only annihilation processes. The problematic approximation corresponds to setting f γ → 0, which amounts to computing the thermal average of the zero-temperature cross section. This step is not justified, since there are relevant regions of photon phase space dΠ γ , where the photon energy is small, in which case f γ ∼ E −1 γ is arbitrarily large. However, if one simply keeps f γ in the expression for the collision term, the virtual one-loop and real terms, |M NLO χχ→ij | 2 and dΠ γ |M χχγ→ij | 2 , respectively, are multiplied by different factors, and the standard IR cancellation no longer occurs. Moreover, since f γ ∼ E −1 γ , an additional IR divergence is generated, which is more severe than the zero-temperature, logarithmic divergences.
It is now important to realize that the photons in the plasma contribute not only to the 2 → 3 emission and 3 → 2 absorption processes, but also to the virtual, one-loop two-body amplitude. Indeed, it has been shown in the special cases of muon decay [26] and the right-handed neutrino production rate [27][28][29][30][31] relevant to leptogenesis, that when finite-temperature Feynman rules are used in the computation of the decay or production rate, the additional IR divergence cancels. In particular, leptogenesis also involves a nonequilibrium situation. The proof of cancellation of all divergences in the general case does not seem to exist, though some partial results can be found in [25,[33][34][35]. Let us therefore add and make explicit the finite-temperature correction to the virtual correction by replacing likewise for the inverse process. 2 Since the SM particles may have masses smaller or of order of T f , we also abandon the assumption that f i,j are exponentially small in the 2 → 3 and 3 → 2 processes, where particles i, j need not have energy of order m χ . We can then extend and reorganize the NLO thermally averaged cross section into the expression Note that we have neglected terms with more than three distribution functions, as they are necessarily exponentially suppressed relative to those given, since the kinematics of 2 ↔ 3 processes allows only one particle to be soft. The NLO collision term also includes the processes χχj ↔ iγ, χχi ↔ jγ, which appear first at this order.
In the cross section above there are both T -independent and T -dependent IR divergences. The former are present in the second line on the right-hand side of (2.12). However, the expression in the parentheses is IR finite by the standard T = 0 KLN cancellations, and we will not discuss it further in this work. Our main interest is in the remaining two lines which contain the finite-temperature correction to the one-loop virtual amplitude and emission and absorption processes multiplied by additional phase-space distribution functions. Our aim is to show that these terms are IR finite and to evaluate the leading correction. Indeed, our main result will be that the relic density can be obtained by solving the equation analogous to (2.5) with collision term C NLO = σ NLO v rel T =0 n eq χ n eq χ − n χ nχ , (2.13) and the NLO thermally averaged cross section replaced by (2.12), which now depends also on T through a finite-temperature correction.
Owing to the presence of the M NLO T =0 χχ↔ij term, in order to obtain a meaningful result at NLO, one needs to perform the computation of the amplitudes in the thermal field theory formalism. The starting point for a systematic treatment is non-equilibrium quantum field theory and the closed time-path (CTP) formalism [36,37]. In the next section we  Figure 1. The contour C in the complex time plane. The value t max can be taken to be +∞ for practical computations. review the derivation of the Boltzmann equation from the Kadanoff-Baym equations [38] with application to relic density computations. The general strategy of this derivation is well-known and gives a prescription for the computation of the collision term, which consistently takes into account all the thermal corrections. We recapitulate it here to set up the notation for the concrete calculations to follow. These are performed in an example model for DM annihilation, where we can directly observe the cancellation of both soft and collinear divergences. As we will show, the IR finiteness of the collision term is related to the finiteness of DM particle self-energy diagrams in the thermal background. The formalism allows us to compute the finite-temperature correction and we find that the naive zero-temperature NLO relic density computations are accurate up to corrections of order O(ατ 2 ), where τ ≡ T /m χ ≪ 1 and α is the fine structure constant. The correction has a remarkably simple form.

Derivation of the Boltzmann equation
In this section we briefly review the derivation of the kinetic equation for non-equilibrium propagators, and from this the Boltzmann equation for the phase space density functions by performing the Wigner transformation and gradient expansion (see, for example, [39][40][41][42][43]). We start with the closed time-path (CTP) formulation of non-equilibrium quantum field theory (reviewed, e.g., in the book [44]), where all correlation functions are defined on a complex time plane along the contour C, see fig.1. The contour Green function for a fermion is defined as where T C denotes the time ordering operation along the contour. It corresponds to four Green functions with real-time arguments: where T c (T a ) denotes chronological (anti-chronological) time ordering in real time. 3 The brackets . . . imply averaging over an ensemble at time t min . The free-field momentum-space propagators are summarized in appendix A. The formalism describes a general nonequilibrium system, where the physical macroscopic observables are averages over an ensemble. The CTP formulation originates from the need to describe the time evolution of an operator expectation value ("in-in formalism") rather than a scattering matrix element ("in-out formalism").
In an interacting non-equilibrium system the two-point Green functions depend on both space-time coordinates, which may be chosen as relative coordinate r = x − y, and averaged (macroscopic) coordinate X = x+y 2 . In equilibrium the system can depend only on the relative coordinate due to translation invariance. Therefore, for systems not far from equilibrium it is useful to perform the Wigner transform and define the Green functions (and similarly the self-energies). The dependence on p describes the fluctuations on the microscopic scale of particle interactions, while the coordinate X describes the macroscopic space-time variations. In equilibrium, the Wigner-space Green functions depend only on the momentum p.
The contour Green functions obey the Dyson-Schwinger equation where the superscript '0' denotes the free propagators, and Σ is the self-energy. The Dyson-Schwinger equations lead to the Kadanoff-Baym equations [38] ( where the collision is defined as Here the subscript h denotes the hermitian part, Σ h = Σ c − 1 2 (Σ > + Σ < ) and analogously for the Green functions. 4 These equations for the Green functions are exact functional equations, but too difficult to solve. At this point we apply the approximations described in the introduction. First, we transform to Wigner space and take t min = −∞. Then we perform the gradient expansion up to the first order in gradients, upon which (3.6) becomes where the Green functions now depend on X and p, ∂ µ ≡ ∂ ∂X µ , and denotes an analogue of the Poisson bracket with respect to the coordinates X and p. The expanded collision term in Wigner space reads By performing the gradient expansion we assume that the variations of physical quantities in the coordinate X are small with respect to the typical inverse momenta of the plasma excitations. The latter are of the order of plasma temperature T . As concerns the former, in the homogeneous Universe, the macroscopic variation of dark matter particle number density is set by the expansion rate H and the annihilation rate e −mχ/T , both of which are of the same order, when the number density freezes out at m χ /T ∼ 20. Thus gradients ∇ ∼ Γ are exponentially suppressed and it is a good approximation to keep only the zeroth order in the gradient expansion, which corresponds to neglecting all terms with Poisson brackets in (3.8). In addition there is also an expansion in the coupling constants of the interactions, as long as they are weak. As we show below, when the collision term is evaluated at lowest non-vanishing order (and in zeroth order of the gradient expansion), one recovers the standard freeze-out equation for the DM number density. Since in the cases of interest the next order in the coupling expansion is much more important than higher-order gradient terms, but still allows for a perturbative expansion in α. Thus, in the following, when we consider relic density computations at NLO, we mean next-to-leading order in the coupling constants in the collision term, but leading order in the gradient expansion. The relevant equation is then Separating the hermitian and anti-hermitian parts leads to constraint and kinetic equations where here {·, ·} and [·, ·] denote the anti-commutator and commutator, respectively. The constraint equation (3.13) to zeroth order takes the simple form It describes the spectral properties of the quasi-particles and in particular puts constraints on the structure of the Green function. Inserting the most general parameterization of the Dirac matrix structure compatible with spatial isotropy, the constraint equation (3.15) leads to the conditions Hence the Green function must be of the form The general solution of the constraint equation (3.15) can also be written in the form of the Kadanoff-Baym ansatz where A 0 (X, p) = ψ(x),ψ(y) 0 describes the spectral properties of the quasi-particles in zeroth order, and is given by Finally, the Boltzmann equation follows from combining the kinetic equation (3.14) in zeroth order of the gradient expansion with the quasi-particle approximation and the solution (3.18), (3.22-3.23) of the zeroth-order constraint equation. We first note that the term p · γγ 0 + m χ γ 0 , iγ 0 S < > vanishes with the above ansatz for S < > . Next we examine the terms containing commutators with self-energies. We assume that the deviation of the DM particle distribution from thermal equilibrium is sufficiently small that the self-energy can be computed with propagators (3.19-3.20). Then at one-loop we can use parametrizations where α,β, σ, a < > , b < > and c < > are scalar functions of the momentum. With this ansatz one can check that both iΣ < > γ 0 , γ 0 S h and Σ h γ 0 , iγ 0 S < > are proportional to p · γ and for this reason, after taking the trace over spinor indices, will not contribute to the Boltzmann equation.
Then, multiplying (3.14) by 2Θ(p 0 ), taking the trace, and integrating over p 0 , we obtain, after using (3.18) with (3.22-3.23): In the FRW background, the time derivative needs to be replaced by the covariant one, from which we recover (2.1) with the collision term C[f ] given by the right-hand side of (3.25). The Boltzmann equation for the number density (2.2) has the integrated collision term given by In the next section we will demonstrate the consistency between the above equation and (2.3), and compute the collision term to NLO.

The collision term
In this section we present the computation of the collision term at zeroth order in the gradient expansion in a "bino-like" DM model. In this model the DM Majorana fermion χ annihilates at tree-level into SM fermions via a 2 → 2 process, mediated by t-and uchannel exchange of a heavy scalar particle φ. The collision term is computed including the thermal NLO contributions to the annihilation process. After introducing the model, we illustrate explicitly the calculation of the right-hand side of (3.12) at tree-level in the CTP formalism to show how the collision term can be expressed in terms of annihilation cross sections and phase-space distributions, as in the standard expression (2.3). We then proceed to the main part of the paper and compute the thermal NLO corrections. As we are primarily interested in the infrared divergence cancellation at finite temperature and the leading finite-temperature correction, we drop the terms that can be associated with the T = 0 NLO correction to the annihilation cross section, even though its finite part is parametrically larger than the finite-T correction. The computation of the zero-temperature cross section is already well understood and could be straightforwardly included in the formalism. We note that while we focus on a particular model to perform explicit calculations, the procedure itself is much more general and can be applied to a variety of different DM scenarios.

The model
We consider the extension of the Standard Model by an SU (2) × U (1) singlet Majorana fermion and a scalar doublet φ = (φ + , φ 0 ) T . The relevant terms in the Lagrangian read In this model the only interaction involving the DM particle χ is the Yukawa interaction with the "sfermion" φ and SM (light) fermion doublet f , of which we include only the charged component. The neutral component would affect the inclusive tree-level cross section through the λχP L f 0 φ 0 interaction, which allows χχ → f 0f 0 , however this process receives no radiative corrections since it contains only electrically neutral particles. The scenario we have in mind, realized in the minimal supersymmetric SM (MSSM) if the dark matter is the bino, is an electroweak or TeV scale DM particle, and a scalar (sfermion) with mass m φ > m χ ≈ O(0.1−1 TeV). In this situation the freeze-out occurs after the electroweak phase transition. In the covariant derivative D µ = ∂ µ − ieA µ we therefore keep only the electromagnetic term.
The motivation for studying this particular scenario follows from its relevance to the dark matter phenomenology of the MSSM, and from its relative simplicity. Moreover, in such a model the zero-temperature NLO corrections have been shown to be significant [45], since they lift the helicity suppression of the LO annihilation process.

Calculation of the collision term at LO
In the CTP formalism the fermion collision term (3.10) to leading order in gradient expansion is given by In the calculation of the self-energies the phase-space distribution functions of all the interacting particle species appear in the finite-temperature propagators, see appendix A. The two terms Σ < and Σ > account for all possible processes involving the interacting species, which includes annihilation, production and scattering processes for χ, as well as absorption processes characteristic of the finite-temperature plasma. In the kinetic equation for the particle number density, the contributions from particle-number preserving scattering processes χf → χf cancel out after summing over the two terms on the righthand side of (4.2) and after taking the trace and performing the integral over the particle four-momentum in (3.26). These terms will therefore be omitted right away. We start from the calculation at leading order in the coupling (loop) expansion to show the correspondence between the self-energy diagrams and annihilation processes. The oneloop self energy, shown in fig. 2, describes 1 ↔ 2 processes, which are not relevant for the relic-density computation, because they are kinematically forbidden or exponentially suppressed. 5 Therefore, the LO annihilation process χχ ↔ ff must be contained in the two-loop self-energy diagrams of fig. 3. The self energies Σ <,> are computed from the diagrams discussed above by applying the Feynman rules of the CTP formalism with the propagators of appendix A, and the proper treatment of the fermion-number violating interactions for Majorana fermions [46]. We denote the propagator of the charge-conjugate fermion field as S ′ ab (p) ≡ C S ab (p) T C −1 , where C is the charge-conjugation matrix and the transpose is with respect to the spinor indices only. Figure 3. The DM self-energy at two loops. The same diagram topologies with reversed arrows are not shown for simplicity. In the following they are into account and denoted by a superscript rev. Let us consider the contribution to Σ > (q) from diagram A in fig. 3. Since Σ > = Σ 21 , the left vertex is of the type '1' and the right of type '2', while one has to sum over both types of internal vertices. This leads to the sum of the four diagrams in fig. 4, where uncircled and circled vertices denote type '1' and type '2' vertices, respectively. Fixing the fermion flow and assigning the momenta as in fig. 3 the whole expression appearing in the collision term reads At this point we note that the thermal part of the sfermion propagator is exponentially suppressed, since m φ > m χ ≫ T f . Dropping this part implies that only the 11 and 22 components of ∆ ab are non-vanishing, so the ends of a scalar (dashed) line must either both be circled or not. Hence the only diagram in fig. 4 that we have to compute is A III . Taking the trace over the spinor indices, which accounts for the polarization sum in the number density equation, the previous equation simplifies to Since in the scalar part S we need only the T = 0 part of the propagators, we have (4.5) -13 -In the fermion part F both the T = 0 and the thermal parts contribute, in principle. However, the expression involves only the purely thermal off-diagonal CTP propagator, leaving The last two lines of (4.6) lead to 16 distinct terms describing different processes in the thermal plasma. Half of them vanish after multiplying by Θ(q 0 ) as needed for (3.26).
Out of the remaining 8 terms 5 are kinematically forbidden, since they refer to 4 ↔ 0 and 1 ↔ 3 processes. One is left with two terms corresponding to scatterings χf → χf and χf → χf , which do not contribute to the number-changing processes and cancel out after including the Σ < S > contribution, and one term describing the annihilation process χχ → ff. Only this last term contributes to the integrated collision term. As explained in the introduction, we assume the background plasma to be in thermal equilibrium with zero chemical potential and therefore the SM fermion distribution function is the Fermi-Dirac one for both particle and antiparticle, ff = f f = f eq f and Finally we get where all the momenta are on-shell. Adding the hermitian conjugate and integrating this expression with −h χ d 4 q/(2π) 4 1 2 Θ(q 0 ) as appropriate to the collision term for the χ number density (3.26) and accounting for the factor 1/2 in (4.2), the structure of the result is now manifestly as in (2.3), with a zero-temperature annihilation cross section times velocity multiplied by the statistical factors corresponding to the process The matrix element squared can be recognized as the interference term between the two tree-level diagrams for the annihilation process χχ → ff , as shown in fig. 5. Specifically, where the trace refers to the first line of (4.6). The same procedure applied to the diagram B in fig. 3 and to the corresponding diagrams with reversed arrows leads to the identifications (4.10) Diagram C of fig. 3 does not contribute, since as discussed above we can ignore any contribution with an off-diagonal sfermion CTP propagator. The calculation of Σ < S > is analogous and reproduces the first term in (2.3), which corresponds to the production process ff → χχ. We therefore conclude that -as anticipated -at LO in the CTP formalism, that is, inserting the DM self-energy at two loops into (4.2) and (3.26) for the integrated collision term, leads to the standard Boltzmann equation (2.2). This is true under the assumptions of the gradient expansion and quasiparticle approximation, which are well-satisfied for the standard scenario of freeze-out of an initially thermal DM particle population. At LO in the coupling expansion the integrated collision term is, provided the tree level 2 → 2 processes are χχ ↔ ff , as in (2.3).

Tree-level annihilation cross section
For later reference, we give the tree-level χχ → ff cross section in our model. More precisely, we give the cross section times velocity expanded in the small velocity (partial waves) in the non-relativistic regime, 6 where v is the CM velocity of one DM particle and we extracted the flux factor 4E 2 χ so that a tree and b tree are dimensionless, and correspond to the s-and p-wave contributions, respectively. It proves useful to adopt variables rescaled by the DM mass, i.e. Then (4.14) Note that always ξ ≥ 1 and ǫ ≤ 1 2 , but typically ǫ ≪ 1. In the first term the appearance of the ǫ 2 factor implies the well-known helicity suppression of s-wave annihilation of a Majorana fermion into SM fermions.

Collision term at NLO
Now that we understand how to match the collision term in the CTP formalism to the form of the freeze-out equation with the standard computation of annihilation cross sections, we are ready to consider the NLO thermal corrections. They are encoded in the threeloop DM self-energy diagrams obtained by adding a photon line to the diagrams A and B in fig. 3 in all possible ways. 7 From the annihilation amplitude point of view they can be arranged into three classes: i) processes corresponding to thermal emission and absorption, ii) thermal internal virtual corrections and iii) thermal corrections to mass and wavefunction renormalization on the external legs. We use this classification for organizing the discussion of the computation, even though it is somewhat artificial from the self-energy diagram point of view. The reason is that we want to show a clear connection between the usual way of doing calculations and the quantities appearing in the collision term as derived from CTP formalism. When showing the results for IR divergence cancellation and leading thermal correction we revert to the more natural classification based on different self-energy diagrams.
At NLO there are 20 three-loop self-energy diagrams contributing to Σ > of a Majorana fermion. 8 They are given in tables 1 and 2, together with the corresponding processes they describe after associating the terms in the CTP sums with matrix elements squared. Since the thermal part of the propagators always contains the on-shell delta function δ(p 2 − m 2 ) we refer to these contributions as "cuts" of the self-energy diagrams.
In the remainder of this section we describe the method of performing the calculations emphasizing the differences with respect to the T = 0 case. The results and their discussion will follow in section 5.
As an example we consider the self-energy diagram in fig. 6, which corresponds to the last diagram in table 2. Following the rules described in section 4.2, we note that every three-loop CTP self-energy contribution to Σ <,> contains 2 4 = 16 different terms from the CTP sum over circled and uncircled vertices. Most of them, however, do not contribute as they involve the exponentially suppressed thermal part of the sfermion. In case of the example shown in fig. 6 only four terms remain, which we may associate with virtual and real photon NLO corrections. To confirm this interpretation, we consider the second cut of Figure 6. An example three-loop self-energy contribution to iΣ > decomposed into a sum over "cuts" and the interpreation of the cuts as scattering processes. iΣ > is obtained by taking the sum over all possible diagrams in which the vertex attached to the external line on the left (right) is of type '1' ('2'). The correspondence between reversing the charge flow arrows and crossing the external fermion legs is the same as displayed in fig. 5. For simplicity, from this figure on, we will denote with a single diagram with no arrows the sum of the two diagrams with and without reversed arrows.
the diagram in fig. 6, labelled CA, and show that it corresponds to the interference term of the two real photon emission amplitudes from the different final state legs multiplied by the associated Bose enhancement factors. Proceeding as in section 4.2 we obtain for this contribution the expression In the scalar part S it is again sufficient to keep only the T = 0 part of the propagators, while the photon propagator V contains only the thermal part. Omitting for brevity the traces over the numerator Dirac matrices we get From the above expressions we see that the contributions from the thermal parts of the Real Virtual S 11 (k 1 +s) and S 22 (−k 2 −s) are vanishing, since the combination of δ-functions multiplying those terms has no support. We are then left with 32 terms, from which again half vanishes after multiplying by Θ(q 0 ). Out of the remaining terms 8 describe emission and 8 absorption of a photon attached to the tree-level diagram of type B. Among those terms there are 6 that correspond to processes which are kinematically forbidden and from the remaining ones 6 describe scatterings and 4 annihilation. Only the annihilation terms eventually contribute to the Boltzmann equation for χ particle number, hence (4.15) simplifies to where the equilibrium distribution function for the photon is given by the Bose-Einstein statistics The factors |M CA | 2 collect the traces contained in the definition of F 1 and F 2 , coupling constants, as well as the non-thermal propagator denominators. The first one (the third line of (4.19)) can be identified with the interference of zero-temperature emission amplitude, namely |M CA (q, k 1 , k 2 )| 2 = M C (M A ) * . By using the crossing symmetry one can identify the remaining ones with the parts of the amplitudes for absorption processes: where the minus sign comes from interchanging the fermions between initial and final states. The example shows that the surviving terms from the three-loop CTP self-energy correspond precisely to the collision term in the form of (2.12). In the following we discuss the computation of the IR divergent and leading IR finite thermal correction, separately for the real and virtual cuts. The results are summarized and discussed in the following section 5.

Thermal emission and absorption
The computation of the emission and absorption processes at non-zero temperature follows the same procedure as is well known from the T = 0 case, simply because at NLO all the contributions from the thermal part of internal propagators are either exponentially suppressed or kinematically forbidden at this order, as explained for the example above. The only difference comes from the fact that the external particles can be thermal, in which case the corresponding external leg is multiplied by the phase-space distribution function.
From (2.12), the thermal emission and absorption contributions to the annihilation process are given by Let us focus first on photon emission and absorption in χχ annihilation as given by the first line of (4.22). In the freeze-out situation the phase-space distribution of the photons is always close to equilibrium, and therefore emission and absorption of hard photons with energies ω of order of m χ ≫ T ∼ T f is exponentially suppressed by the distribution function f γ . The scattering matrix elements can therefore be evaluated in an expansion in ω ∼ T f ≪ m χ , that is, in the soft-photon regime. In particular, the leading IR divergence could be obtained from the amplitudes in the eikonal approximation. However, since we are interested also in the leading finite thermal correction, which turns out to be of order O(τ 2 ), we compute the full amplitude. After performing the integration over all phasespace variables except the energy ω of the emitted or absorbed particle, we are left with an expression of the form (4.23) where the range of the integration for ω is determined by the phase space delta function contained in dΠ ijγ . For absorption ω max = ∞, while for emission there is a kinematic limit ω max . Since ω max = O(m χ ), the upper limit is not relevant, however, because f γ is already exponentially suppressed, and the limit may be extended to infinity. At this point we can perform an expansion of S(ω) retaining terms up to linear order in ω, 9 and the final integral over ω can be expressed in terms of The divergence for n = 0 follows from the Bose enhancement f B (ω) ∼ 1/ω of soft photons and implies a stronger divergence than at zero temperature, where the soft IR divergence is only logarithmic. There is no such enhancement for fermions, hence when the same considerations are applied to the SM fermion emission and absorption terms in the last line of (4.22), the relevant integrals are Of particular relevance will be the integrals In order to obtain the above analytic expression I 1 for the leading thermal correction in the case of thermal fermion emission and absorption, we have to assume that the fermions are massless. We briefly discuss the size of the corrections from finite fermion masses in section 5.3. Returning to photon emission, the amplitude for the annihilation process χ(p A )χ(p B ) → f (k 1 )f (k 2 )γ(q) can be written as 10 where the letters A, B, C refer to the amplitudes as given in table 1. After the Fierz rearrangement the three terms are (4.28) For the absorption process, due the crossing symmetry, the amplitude squared summed over polarization can be obtained from the emission process by changing the sign of the four momentum of the particle emitted and absorbed from the thermal bath, as in (4.21).
Although the emission and absorption contributions have different phase-space integration limits, we have already seen that this is irrelevant up to exponentially small terms in m χ /T . Thus, when the emission contribution is expanded in the form eq. (4.21) implies for the corresponding absorption process. Since (4.22) always involves the sum of emission and absorption, the even terms in the expansion in ω cancel. Eqs. (4.24), (4.25) then imply that the leading finite-temperature correction is of order The contributions from thermal photon emission and absoprtion, though divergent, can be computed without regularization in four dimensions, since the cancellation of the linear IR divergence proportional to J −1 with the thermal virtual correction will be shown algebraically below before integration over the photon energy. The integration over the remaining phase-space variables that was already done to arrive at the function S(ω) is finite, since the non-vanishing fermion mass plays the role of the regulator for collinear divergences. This is no longer the case when the thermal fermion emission and aborption processes are considered, since the integral over photon energy contained in S(ω) has to be regularized. In this case we perform the phase-space integration in dimensional regularization with D = 4 − 2η and η < 0.

Thermal virtual corrections
Thermal virtual corrections arise from terms in the CTP sum, to which the thermal parts of the diagonal 11 or 22 photon and fermion propagators contribute. As the sfermion is at least as heavy as the DM particle it has a negligible thermal contribution and we do not consider the corresponding amplitudes. We only need to include the terms when one of the virtual particles is thermal. When two are thermal this gives the imaginary part, which does not contribute to the real part of the interference with the tree diagram (see e.g. [47]), while when three are thermal at least one of them has to have momentum of order m χ and is therefore exponentially suppressed by the phase-space distribution function.
We denote the relevant amplitudes by M i with i = 1, ..., 4, and the contribution from the thermal part of the photon (SM fermion) propagator by M γ i (M f i ). The corresponding diagrams are displayed in table 1. The general form of every virtual contribution is where m γ = 0 and for thermal photons and for thermal fermions, (4.40) Note that F f 3 = 0, since there are no internal fermion lines in this diagram, and that F f 4 has to be counted twice, since any one of the two fermion lines in the loop can be thermal. Given these expressions we first perform the integral over q 0 , which leads to with ω ≡ | q | for photons Changing the integration variable q → − q in the second integral gives F γ,f i (q) + F γ,f i (−q) in the bracket. Then we compute the interference of the resulting expression with the tree-level amplitude and perform the integration over the two-body phase space together with the angles of q. We are left with an integral over the ω similar to (4.23), which can be computed in expansion in T /m χ by expanding the integrand in ω. The result involves the same integrals J n , I n as was the case for the emission and absorption terms.

Thermal corrections to external legs
The remaining part of the virtual correction can be interpreted as a thermal correction to the mass and wave-function renormalization of the external SM fermion lines. Due to the universality of the renormalization factor, we can follow the standard procedure (see e.g. [26,47]) of computing the one-loop corrected thermal propagator (4.43) When the result for the self-energy at one loop is written as with quantities c B,F , K µ B,F to be defined shortly, the propagator is expressed as . (4.45) The subscript B refers to the contribution when the photon propagator in the one-loop self-energy diagram is thermal and the SM fermion propagator is not. Vice-versa for the quantities with subscript F . Then for the thermal photon contribution, and The wave-function renormalization factor is derived from the expansion of the propagator around the particle pole. Letp µ = (p 0 , p ) withp 0 = (m 2 f + p 2 ) 1/2 be the on-shell limit of the external momentum p, and letf denote the value f (p) of a function f (p 0 , p ). Then one finds that c B vanishes on-shell by antisymmetry of the integrand under q → −q, i.e.ĉ B = 0, so that its expansion around the on-shell value reads and J −1 the divergent integral defined in (4.24). The explicit calculation of the integral defining c F in (4.47) shows that the thermal fermion contributionĉ F is only vanishing in the m f = 0 limit, so that in general The coefficientsĉ F andĉ ′ F can be obtained in general by solving the integrals numerically. In the massless case they simplify tô with I −1 defined in (4.25). The vector contribution from the photonK µ B in the on-shell limit readŝ and J 1 given by (4.26). The fermion contributionK µ F in the massless limit is divergent on-shell. We use dimensional regularization (D = 4 − 2η and the MS scheme), which giveŝ , (4.52) where is the D-dimensional generalization of (4.25). The quantities 2p · K B , 2p · K F appearing in the denominator of (4.45) can be related to c B , c F by as follows immediately from the defining expressions (4.46), (4.47). Here contribute to the thermal corrections to the fermion mass. Expanding the fermion propagator (4.45) around the corrected mass-shell, we obtain, up to non-singular terms, In the massless case m f = 0, this simplifies to In summary, the thermal plasma affects the external SM fermion lines, and therefore the annihilation cross section computation, in three ways: 1. Modification of the spinor orthogonality relations This contribution to the annihilation cross-section σ χχ→ff v rel is simply obtained by computing the tree-level diagrams with the modified relation (4.58) and taking the O(α) term. We note from (4.51) and (4.52) that this contribution is finite for thermal photon, while it contains a 1/η pole for thermal fermion in the massless limit. This pole cancels when adding the corresponding real correction "cut".
2. Temperature-dependent wave function renormalization The contribution is simply the O(α) term in (Z T 2 ) 2 − 1 (σ tree χχ→ff v rel ). We note that this contribution is divergent only for the thermal photon case, and it vanishes for the thermal fermion in the massless limit.
3. Shift of the fermion pole mass by the thermal contributions which leads to a change in the phase-space integration. This results in a contribution to cross section that can be written as where we used the short notation σ tree ≡ (σ tree χχ→ff v rel ). This contribution is finite for both thermal photon and fermion. In the massless limit this is ensured by the fact that σ tree is analytic in m 2 f .

Results
In this section we summarize the results of the calculation of the thermal correction. We first note that we are interested in the situation T ≪ m χ (τ ≪ 1). The thermal correction arises from soft thermal propagators, ω ∼ T up to exponentially small terms, since larger energies are suppressed by the thermal distribution functions. After expansion in ω/m χ , all the thermal real and virtual corrections can be expressed in terms of the integrals J n for photons, and I n for massless fermions, as defined in (4.24), (4.25). Since we are interested in the infrared divergence cancellation and the leading thermal correction, we only keep terms to order O(τ 2 ). The scattering processes depend further on the masses m χ , m f , m φ , and the DM energy E χ , or the corresponding dimensionless variables ǫ, ξ and e χ = E χ /m χ . Freeze-out occurs when the DM particles are non-relativistic, so that we can expand in e χ ≪ 1. We performed the calculation for the first two terms of this expansion, which correspond to the s-and p-wave terms, respectively, keeping the full dependence on the scalar and SM fermion mass parameters ξ and ǫ. All computations were done in Feynman gauge. In addition, we also computed the result without an expansion in the non-relativistic DM kinetic energy, keeping the full dependence on e χ . However, in this case we performed an expansion for large scalar mediator mass ξ ≫ 1, up to the order O(ξ −10 ). The large-mass expansion may be physically motivated, since often (but not necessarily) the scalar particle in the DM model is significantly heavier than the DM particle itself. Going to high order in 1/ξ is motivated by the observation that one needs to retain terms up to O(ξ −8 ) to see the lifting of helicity suppression of the non-thermal NLO contribution.

IR divergence cancellation
As we have seen above, the basic quantities from which the inclusive annihilation cross section is derived are the CTP self-energies. It is well known that at zero temperature, off-shell self-energy diagrams are IR finite. The cancellation of divergences between virtual and real corrections to an inclusive process occurs after summing all possible cuts of the self-energy diagram. Our first important result is that we find that this is also true at finite temperature. That is, in the sum of all additional contributions from the thermal part of the propagators, the IR divergences cancel. Moreover, as at zero temperature, this happens at the level of individual CTP self-energy diagrams. This ensures that the collision term in the Boltzmann equation is IR finite, since it is directly built out of the self-energies Σ <,> .
In order to show how the cancellation takes place, we discuss in detail the correction from thermal photons to the collision term for s-wave annihilation. We also verified the cancellation for the p-wave term, for the contribution from thermal fermions, and without partial wave expansion in the large-scalar mass expansion, as discussed above.
At one-loop the amplitude can have singular terms of the order O(ω −1 ) at small ω, which at T = 0 lead to the usual logarithmic soft divergence. At finite temperature, the enhancement of the Bose distribution function f B (ω) for small energies results in linear and logarithmic divergences, encoded in the singular integrals J −1 and J 0 , respectively. As already pointed out, the latter vanishes when both the emission and absorption of thermal Table 3. Coefficients of the divergent integral J −1 omitting the overall factor α/(πǫ 2 ) × a tree .
Here "Real" includes both, emission and absorption, while "Virtual" comprises vertex and external leg corrections. An empty space means that the corresponding contribution does not exist, while 0 implies that the diagram exists, but is finite. L denotes the logarithm as defined in table 6. photons are included, due to the different sign of these contributions for even orders in ω. The results for the remaining part proportional to J −1 are given separately for all self-energy diagrams in table 3, where the prefactor α/(πǫ 2 ) × a tree involving the tree-level s-wave annihilation cross section (4.13) has been factored out. 11 We immediately note the aforementioned fact that the sum of all contributions cancels for every self-energy diagram separately. The logarithm present in the last row is defined in table 6 and contains a collinear divergence L ǫ→0 → log ǫ in the limit of small SM fermion mass (ǫ = m f /(2m χ )). This collinear divergence also cancels in the sum.
The tree annihilation cross section is helicity-suppressed, a tree ∝ ǫ 2 . The appearance of terms in table 3 and tables 4 and 5 below, which are not O(ǫ 2 ) for small ǫ, implies that individual terms are not helicity-suppressed.
The finite part J 1 Type A Real Virtual Table 4. Coefficients of the finite O(τ 2 ) correction for the type A diagrams. An overall factor πα/(6ǫ 2 ) × a tree is left out. D, D ξ and polynomials p i and f i are defined in table 6.

Finite-temperature correction from thermal photons
Once the divergent J −1 and J 0 contributions are cancelled, the remaining finite correction is necessarily of O(τ 2 ), proportional to the integral J 1 . Again, we first show the diagramby-diagram results for the s-wave contributions, which can be found in tables 4 and 5 for the diagrams of type A and B, respectively. A common factor πα/(6ǫ 2 ) × a tree is left out. We see that the separate contributions are significantly more complex than was the case for the divergent parts. The first simplification occurs when summing over the different cuts of a given self-energy diagram. At this stage all the logarithms L cancel, which is a sign of cancellation of the collinear divergence on a diagram-by-diagram basis. An even more remarkable simplification occurs upon adding separately all diagrams of type A and B, respectively (given as "Total" at the bottom of the tables). Finally, helicity-suppression is recovered after summing over A and B. The thermal correction to the s-wave annihilation cross section (times velocity, see (4.11)) can be written as It is worth noting that the leading thermal correction is suppressed not only by ατ 2 but also one power of ξ 2 , if the mediator mass is large. This is true not only for the s-wave contributions, but also for all partial waves.
The finite part J 1 Type B Real Virtual Table 5. Coefficients of the finite O(τ 2 ) correction for the type B diagrams. An overall factor πα/(6ǫ 2 ) × a tree is left out. D, D ξ and polynomials p i and f i are defined in table 6. Beyond the s-wave case displayed explicitly in the tables, we computed the thermal correction to the p-wave cross section; further without the partial wave expansion in the limit ξ ≫ 1, up to the order O(τ 2 , ξ −10 ), retaining full dependence on e χ and ǫ. We find that all cases are covered by the remarkably simple expression which implies that the total O(τ 2 ) correction from thermal photons can be obtained directly from the tree cross section without any explicit calculation. To appreciate the simplicity of this expression, note the complicated dependence on ξ and ǫ of the tree-level p-wave cross section given in (4.13). The fact that this formula holds in all limits that we investigated leads us to conjecture that it is generally valid beyond the non-relativistic approximation (partial-wave expansion). The structure of (5.2) is certainly not accidental, and it is not the only example of "universality" of a finite-temperature correction. In charged particle decay [26] the finitetemperature correction was found to be related to the tree-level decay width by the simple factor − π 3 ατ 2 , while in neutral Higgs decay to two fermions it vanishes [48]. This suggests that the leading thermal photon correction is related to the coupling to the electric charge multipoles of the initial or final state. In DM pair annihilation the total charge is zero, but higher moments are not, which may be the reason for the ξ suppression. Further work in this direction is in progress.
The leading O(τ 2 ) thermal correction does not lift the helicity suppression of the swave cross section, even though the NLO T = 0 radiative correction does. This is easy to understand, since it is the hard photon emission from internal bremsstrahlung from the scalar mediator that lifts the helicity suppression in the T = 0 case, while here such contributions are strongly suppressed. Helicity-suppression is absent in the first sub-leading O(τ 4 ) temperature correction. Explicitly, in the limit of massless SM fermions ǫ → 0, we find ∆a ǫ=0 τ 4 = 8π 2 λ 4 ατ 4 45 This thermal correction can be larger than the tree-level s-wave cross section a tree when ǫ is very small (e.g., for SM leptons). Nevertheless, it is always parametrically smaller than the thermal correction to the p-wave cross section, because τ ∼ v 2 around freezeout. It is also smaller than the zero-temperature O(α) NLO correction, which has no τ 4 suppression, while both come from internal bremsstrahlung, and therefore have the same ξ −8 suppression. Finally, we note that (5.2) holds also when the DM particle is a Dirac fermion, as can be expected from the structure of the total in table 5. The difference between the Majorana and Dirac cases is that for Dirac fermions the diagrams of type A are absent altogether, which also implies the absence of helicity suppression.

The finite-T correction from thermal fermions
Like photons the light SM fermions are very abundant in the plasma around freeze-out and also contribute to the finite-temperature correction, see (2.12). The method of computation of these contributions follows the same steps as for thermal photons, and has been described in section 4. However, the results differ considerably between these two cases. This is, because at zero temperature soft fermion radiation does not cause IR divergences. Furthermore, the Fermi-Dirac distribution is finite in the soft limit, hence the degree of divergence is not changed at finite temperature. As a consequence the thermal fermion contributions have no IR divergences from soft fermions. However, for massless fermions there is a divergence from hard-collinear photons, which has the same origin as the corresponding T = 0 divergence. When working in the massless limit, we use dimensional regularization to regulate this divergence. The poles in 1/η cancel in the sum over all cuts for a given CTP self-energy diagram.
We first discuss the leading finite-temperature correction for the case when DM is a Dirac fermion, in which case only diagrams of type B are present. Since the light SM fermion masses satisfy the condition m f ≪ eT , they can effectively be treated as massless. 12 We therefore consider the correction for ǫ = 0. Once again the total result turns out to be simple once all the cuts from a given CTP self-energy diagram are summed up. The s-wave contribution vanishes for each self-energy diagram separately, due to an exact cancellation between real and virtual corrections. The s-and p-wave corrections (for ǫ = 0) that need to be added to the tree cross section are ∆a f = 0 ∆b f = 16 9 ατ 2 λ 4 (1 + ξ 2 ) 3 . (5.4) Turning now to the Majorana case, we find the the s-wave and p-wave O(τ 2 ) contributions from the type A diagrams actually vanish for ǫ = 0. This means that (5.4) also holds for Majorana DM. However, since a tree ∝ ǫ 2 the limit ǫ = 0 is not sufficient for the s-wave contribution, in which case the above only shows that there is no lifting of helicity suppression from the thermal fermion correction. To analyse the result for finite fermion mass we computed the thermal correction numerically in the region ǫ, τ ≪ 1. We find that there is no ǫ 2 τ 2 term in a, which means that indeed there is no leading finite thermal correction to the s-wave cross section, different from the case of thermal photons. On the other hand, for the p-wave cross section the thermal correction from fermions is of the same order as the one from photons.

Conclusions
The present work was motivated by the observation that existing approaches to calculating the DM relic density at NLO are based on zero-temperature calculations of the annihilation cross section in the standard freeze-out equation. This procedure has never been truly justified. In particular, it ignores the potential presence of IR divergences, which, in individual terms, are more severe than at zero temperature. In this paper we showed, using a realistic example model and the CTP formulation of non-equilibrium quantum field theory, that for relic density computations at NLO it is indeed sufficient to treat the annihilation process at the thermal level, leaving the Boltzmann equation intact. That is, under the usual assumptions, the freeze-out equation dn χ dt + 3Hn χ = σ NLO v rel T =0 n eq χ n eq χ − n χ nχ , (6.1) retains its form, and only the annihilation cross section receives a finite-temperature correction. By computing the thermal contributions to the NLO collision term, which includes emission and absorption as well as thermal virtual corrections, we showed that all soft and collinear divergences cancel, a prerequisite for σ NLO v rel T =0 to be well-defined. The cancellation was demonstrated in the non-relativistic expansion for the s-and p-wave cross sections and, additionally, for the full velocity-dependent process but in an expansion in the mass of the t-channel mediator up to order 1/ξ 10 . We find that the cancellations occur at the level of individual CTP dark-matter self-energy diagrams. To our knowledge, this is the first time that the IR finiteness of relic density computations at NLO has been demonstrated.
Finiteness assured, we investigated the leading finite-temperature correction to the annihilation cross section from thermal photons and light SM fermions in the plasma. The result can be summarized as follows: • The leading correction is of order α × (T /m χ ) 2 . Since T ≪ m χ near freeze-out this is parameterically smaller than the zero-temperature NLO correction, thus justifying the naive zero-temperature radiative correction calculations.
• The structure of the O(τ 2 ) correction is remarkably simple. The contribution from thermal photons can be inferred directly from the tree level annihilation cross section, and is given by 2) The simplicity of this result and the amount of cancellations required to arrive at it call for a deeper explanation. Work in this direction is in progress.