Dark Matter bound-state formation at higher order: a non-equilibrium quantum field theory approach

The formation of meta-stable dark matter bound states in coannihilating scenarios could efficiently occur through the scattering with a variety of Standard Model bath particles, where light bosons during the electroweak cross over or even massless photons and gluons are exchanged in the t-channel. The amplitudes for those higher-order processes, however, are divergent in the collinear direction of the in- and out-going bath particles if the mediator is massless. To address the issue of collinear divergences, we derive the bound-state formation collision term in the framework of non-equilibrium quantum field theory. The main result is an expression for a more general cross section, which allows to compute higher-order bound-state formation processes inside the primordial plasma background in a comprehensive manner. Based on this result, we show that next-to-leading order contributions, including the bath-particle scattering, are i) collinear finite and ii) generically dominate over the on-shell emission for temperatures larger than the absolute value of the binding energy. Based on a simplified model, we demonstrate that the impact of these new effects on the thermal relic abundance is significant enough to make it worthwhile to study more realistic coannihilation scenarios.


Introduction
One of the leading dark matter (DM) candidates are Weakly Interacting Massive Particles (WIMPs) [1][2][3], which can account for all the present DM energy fraction [4] through their thermal production mechanism in the early Universe. While current upper bounds on the coupling strength to Standard Model (SM) particles put a variety of thermal dark matter candidates around the electroweak scale under tension, the TeV mass region and above remains less constrained and an attractive possibility. In such a high mass region, however, even the heaviest gauge bosons of the SM start to act as a long-range force between annihilating WIMP pairs, leading to a variety of quantum mechanical phenomena. These are important to include for predicting i) the thermal relic abundance precisely and ii) the flux of final state SM particles produced from DM annihilation in, e.g., the galactic center.
One famous example in the Minimal Supersymmetric extension of the SM (MSSM) is the traditional case of wino-like neutralino, for which the important role of quantum mechanical effects was pointed out in seminal works [5][6][7][8]. Already at the TeV mass region, the ladder exchange of the most massive SM gauge bosons enhances the probability that a slowly moving wino pair annihilates. This quantum mechanical effect is called Sommerfeld enhancement (SE) [9] or Sakharov enhancement [10] and lowers the predicted abundance by about 50 % compared to a tree-level computation in this case [8]. In turn, the enhancement in the cross section allows for a 30% larger wino mass to compensate for the effect. The flux of SM particles from the present neutralino annihilation in, e.g., the galactic center is sensitive in particular to the predicted mass, since small variations can lead to Sommerfeld resonances [6]. This example shows that it is required to predict the DM mass precisely once quantum mechanical effects can drastically change the observational outcome, in order to constrain the model or to estimate the required exposure time for fully testing the thermal case. Various follow-up works extended the studies of the Sommerfeld effect in the MSSM to more general cases, see, e.g., refs. [11][12][13][14][15][16][17][18][19][20][21][22][23] and [24][25][26][27][28][29][30] for formal aspects. For recently refined calculations of the present day SM flux in wino and higgsino models, see refs. [31][32][33].
Additional quantum mechanical effects caused by attractive long-range interactions are the existence of meta-stable bound states in the two-particle spectrum of WIMPs. Their formation and subsequent decay into SM particles can be seen as an additional channel depleting the relic density [34] and therefore allows for even larger dark matter masses. On the one hand, bound-state effects turn out to be negligibly small [35] in the traditional wino case within the current description of bound-state formation (BSF) via the emission of an on-shell mediator. On the other hand, large corrections to the predicted mass were found in other electroweak coannihilation scenarios, e.g., the quintuplet case in the context of minimal dark matter [35]. Similar strong effects were identified for coannihilation scenarios with colored charged particles [36][37][38][39]. In addition to the electroweak gauge boson, photon, and gluon induced bound states, also the Higgs boson [40,41] can attractively contribute to confine DM into a meta-stable bound state. Self-interacting DM [42] with new light mediators [43][44][45][46][47][48][49][50][51][52][53][54], motivated from a bottom-up approach to alleviate the diversity problem in galactic rotation curves [55,56] or in certain cases even simultaneously ameliorate the Hubble tension [57,58], are further examples where long-range interactions can affect the dark matter relic abundance.
So far, the formation of dark matter bound states via the emission of an on-shell mediator was considered as the dominant process. However, it was pointed out recently that this does not necessarily have to be the dominant BSF channel during the thermal dark matter freeze-out [59]. The conversion between scattering and bound states in the radiation dominated epoch can also efficiently take place via bath-particle scattering, where a mediator is exchanged virtually between a DM two-body pair and a relativistic primordial plasma particle.
While in ref. [59], a mediator with a mass larger than the binding energy was investigated, we consider massless or lighter mediators in the present work. The masses of SM force-carries is temperature dependent during the electroweak cross over, which motivates to investigate these cases more carefully. The massless case is also relevant in coannihila- Graphic shows examples of bound-state formation processes via bath-particle scattering. The amplitudes for these processes are divergent in the forward scattering direction of the bath particles. Parallel black solid lines represent the bound state, while open lines correspond to the initial two-particle scattering state.
tion scenarios with electroweakly charged or colored particles. Examples with photons or gluons are illustrated diagrammatically in fig. 1, where bound states are formed via SM bath-particle scattering. Since they all come with a large multiplicity, one may ask an interesting question by how much these additional processes contribute to the previously studied case of BSF via the emission of an on-shell mediator only.
The problem is that the bound-state formation amplitudes diverge for massless mediators in the forward scattering direction of the in-and out-going bath particles in fig. 1 (see also ref. [59]). Since these processes are temperature dependent, the Kinoshita-Lee-Nauenberg (KLN) theorem [60,61] is not directly applicable. Therefore this problem can not be addressed through the computation of higher-order amplitudes in the collision term of the Boltzmann equation in zero-temperature quantum field theory. Throughout this paper, we refer to this approach as the "conventional Boltzmann formalism". Thus, the question by how much these processes could additionally deplete the dark matter relic abundance can not be answered given the conventional methods. Naively, one might insert a Debye screening mass as a regulator. However, even this procedure is not justified for all temperatures. The energy flowing inside the mediator propagator from the inelastic BSF process is at least as large as the binding energy. Thus, hard thermal loop (HTL) effective field theory [62] would break down for temperatures about the binding energy or below.
The main purpose of this work is to develop a more general method of how to calculate dark matter bound-state formation processes inside the Early Universe plasma. One main focus is to address the problem of collinear divergences once massless mediators are involved and to refine the description for light mediators during the electroweak cross over. We manage to derive the BSF collision term in a more general framework of non-equilibrium quantum field theory, explicitly demonstrated for a vector mediator model in section 2. A more general BSF cross section is defined from the collision term, expressed in terms of thermal correlation functions. Based on this novel result, we demonstrate in section 3 that our BSF cross section correctly contains the on-shell mediator emission at the perturbation order of free correlation functions. In section 4, a full next-to-leading order computation for the direct capture into the ground state is presented, cancellation of forward scattering divergences for bath-particle scattering and the off-shell decay is proven in a model independent way, as well as the renormalization of the ultraviolet vacuum divergence are discussed. The impact on the thermal relic density from all these new next-to-leading order effects is shown in section 5 and our findings discussed in section 6. We concluded in section 7.

Generalized bound-state formation cross section
In order to study the infrared divergence structure of higher order bound-state formation processes within the primordial plasma environment for the case of a massless mediator, we consider for simplicity a QED-like model: where χ is a heavy Dirac fermion, (e.g. coannihilating chargino partners) and the Lagrange density of the plasma environment contains all SM fields, e.g. the interaction with light fermions L env. ⊃ −gēγ µ eA µ . The final form of the derived cross section in this section will be independent of the underlying particle content in the environment and therefore L env can be chosen later. Temperatures around the typical chemical decoupling and below are considered, where the χ fields can be assumed to be non-relativistic to a good approximation. We commence from the potential non-relativistic (pNR) effective theory [63], which allows for naturally writing the two-body bound and scattering states in terms of wave function fields. In this framework, the transition between different two-body states are described by an "electric" dipole operator The dipole interaction Hamiltonian contains all possible types of two-body conversion processes, i.e., scattering-scattering, bound-bound, and scattering-bound state transitions. The general two-body field operator O sr (x, r, t) consists of the whole energy spectrum of χ-field pairs, whose individual components are converted into each other through dipole interactions. Spin configurations are conserved in the transitions, represented by s and r for particle and anti-particle, respectively, and summed over. The components with negative energy eigenvalues E B (binding energy) are created by fundamental bound-state operatorŝ c † with quantum number B = {nlm}. The positive energy spectrum with kinetic energy E k = k 2 /(2µ) = µv 2 rel /2 and reduced mass µ is created by the particle operatorâ † and anti-particle operatorb † . The expansion coefficients involve the bound state ψ B (r) and scattering state ψ k (r) wave functions, which are the solution of the time-independent Schrödinger equation with a static Coulomb potential. The wave function dependence on the spin is neglected, and r, x denote relative and center-of-mass coordinates, respectively.
1 At this stage we would already like to remark that scattering and bound states in Eq. (2.3) are assumed to be separable, which is however not valid once the bound-state energy spectrum starts to overlap with the continuum in the high temperature regime for large thermal widths. This phenomena depends on what is contained in Lenv, see discussion in Sec. 6 and Fig. 10 for more details about the limitation of our formalism.
In the following, the time evolution of the particle number density is derived in the density matrix formalism, relating the collision term to the dipole Hamiltonian. For the sake of completeness, we first clarify the underlying assumptions. We take an unperturbed Hamiltonian so that the number densities for the scattering and bound states are conserved and compute the BSF and dissociation rates perturbatively. LetĤ,Ĥ 0 , andĤ int being the full, free, and interaction Hamiltonian, namelyĤ =Ĥ 0 +Ĥ int . The unperturbed Hamiltonian is given by the summation of that for the scattering state, the bound state, and the thermal environment,Ĥ 0 =Ĥ S +Ĥ B +Ĥ env . We take the initial density matrix which is factorized into a tensor product of the scattering state, the bound state, and the thermal environment:ρ(t = 0) =ρ S ⊗ρ B ⊗ρ env . By definition, it commutes with the unperturbed Hamiltonian: Our starting point is the von Neumann equation in the interaction picture: iρ I (t) = [Ĥ I (t),ρ I (t)], where the density matrix and interaction Hamiltonian in the interaction picture are given byρ I = e iĤ 0 tρ e −iĤ 0 t andĤ I = e iĤ 0 tĤ int e −iĤ 0 t respectively. The time evolution of the particle number density is obtained from the expectation value ofn kχ ≡ sâ s † kχâ s kχ as follows: In the second line, we perform the time dependent perturbation with respect H I and take the leading order. 2 We also take t → ∞ by assuming that a typical time scale of our interest is much slow compared to that in the thermal plasma. The time evolution of the anti-particle and bound-state number density can be expressed in the same way. Notice that the expectation value in the right-hand side of Eq. (2.4) should be taken by the initial factorized density matrix.
To derive the change of the particle number density under dipole transitions, we have to evaluate the double commutator in Eq. (2.4) for the dipole interaction Hamiltonian in Eq. (2.2). While the dipole interaction Hamiltonian contains all possible types of two-body conversion processes, only the conversion of a scattering state into a bound state and the reverse process can change the particle number density n χ . Thus, it is sufficient to consider only the mixed contributions containing the product c † ab and its hermitian conjugate part in Eq. (2.2). These two contribution lead to bound-state formation and the reverse process called dissociation. The computational details of the double commutator are discussed in Appendix A.
Here we highlight the most important part in the computation in Appendix A. Since the dipole Hamiltonian contains the electric field, the computation of the double commutator in Eq. (2.4) involves photon two-point functions. As already emphasized, this commutator should be evaluated by means of the factorized initial density matrix, i.e.,ρ(t = 0) = ρ S ⊗ρ B ⊗ρ env withρ env ∝ e −Ĥenv/T . Therefore the photon two-point functions appearing in the computation are nothing but the thermal propagators: As well known, thermal propagators fulfill the Kubo-Martin-Schwinger relation [64,65] (see also chapter "Real-time formalism prerequisites" in [66]): with f eq γ being the equilibrium phase-space distribution obeying Bose-Einstein statistics. The KMS condition relates both two-point functions to the spectral function, defined in coordinate space by D ρ . This photon spectral function encodes all interactions of the primordial plasma environment. We will outline how to estimate the effects from the thermal plasma perturbatively in Eq. (2.12).
With these relations, the change of the particle number density in Eq. (2.4) due to the dipole interactions in Eq. (2.2) is found to be given by: with ∆E = E k − E B , k = (k χ − kχ)/2, K = k χ + kχ, and p being the three-momentum of the mediator with fixed P 0 = ∆E. The transition matrix elements T of scattering and bound states are proportional to the dipole overlap integrals: T i k,B (P 0 , p) ≡ g iδ ss δ rr √ g χ gχ P 0 d 3 rψ B (r)r i ψ k (r). (2.8) The transition element fulfills current conservation P µ T µ k,B (P ) = 0, as a consequence the global symmetry of the model in Eq. (2.1).
The χ fields and the bound states are assumed to be in kinetic equilibrium and dilute. Their phase-space densities take the classical Maxwell-Boltzmann statistics and can be written as: f X = f eq X n X /n eq X (T ). By noticing the standard non-relativistic thermal average in Eq. (2.6), a bound-state formation cross section can be defined as: Replacing in the l.h.s. of Eq. (2.6), the time with cosmic time derivative and adding the usual DM annihilation part in the r.h.s., the particle number density equation can be brought into the standard form: n eq χ n eq χ n eq B − σ an v rel n χ nχ − n eq χ n eq χ . (2.10) Finally, we can identify the generalized bound-state formation cross section as This cross section is one of the central results of this work and can be diagrammatically expressed in terms of the self-energy diagram depicted in fig. 2. The vertices represent the leading order dipole transition contained in T k,B . The black blob indicates the interacting spectral correlation function D ρ of the mediator, which encodes all interactions with the primordial plasma environment and is the key difference compared to previous literature. It can be seen as a probability density of having a mediator excitation at a certain p for fixed ∆E. The excitation can be on-shell but also virtually induced through the thermal environment via, e.g., bath-particle scattering. It is convenient to compute the spectral function from the retarded correlator via D ρ µν = 2 Im iD R µν , since the latter obeys also in thermal field theory the Dyson-Schwinger equation, given in momentum space by In the next section 3, we demonstrate that BSF via on-shell mediator emission is reproduced from the free retarded propagator. In section 4, the first interaction term containing the retarded self-energy Π R is analyzed.
Although the cross section was derived for a particular model, we expect the factorization into the interacting mediator spectral function and the transition elements at the Born level to be a rather model-independent feature. In ref. [59], a direct relation between the transition matrix elements T k,B and relativistic amplitudes has been given. The advantage of this relation is that results of the previous literature, computing the on-shell emission only, can be used to directly determine T k,B . For example, ref. [39] provides already expressions for amplitudes for non-abelian mediators at the leading order dipole approximation, ready to be used with Eq. (2.11). In case of a mediating gluon, this would allow to study BSF via gluon scattering (triple vertex), see the fourth diagram in fig. 1, which is expected to occur at the first order of interactions in the gluon spectral function. For Yukawa interactions involving scalar mediators, one can drop the Lorentz indices in Eq. (2.11). The corresponding transition elements, as well as many other cases, can be found in ref. [67,68]. 3 Recovering on-shell emission at the leading order Consider the dominant direct capture into the ground state nlm = 100 for the QED-like model in Eq. (2.1). The overlap integral in Eq. (2.8) has been analytically computed in, e.g., ref [67]. From this result one can easily get an expression for our transition matrix elements as: with ζ ≡ α/v rel and ∆E = m χ α 2 (1 + ζ −2 )/4 for the ground-state capture. The nice feature now is that Eq. (3.1) can be reused in Eq. (2.11) at any order in the two-point correlation function of the photon field to determine leading and higher-order bound-state formation cross sections inside the primordial plasma environment. Focusing in this section on the leading order, the free retarded propagator is given by The imaginary part of this expression determines the free photon spectral function, which is nothing but the on-shell contribution: The cross section for capture into the ground state via the emission of an on-shell vector mediator is recovered (cf., e.g., [67]) by inserting the free spectral function in Eq. (3.3) together with the transition matrix elements in Eq. (3.1) into the main formula Eq. (2.11): At leading order of the mediator spectral function, the generalized BSF cross section reduces to the result one would obtain in the Boltzmann formalism based on vacuum amplitudes. In this example, the result is the same as for SM neutral hydrogen recombination, with the emission of a photon. The correspondence between the leading order self-energy diagram and the amplitudes in the vacuum field theory is shown in fig. 3. In the regime where T ∆E, ζ 1, and using lim ζ→∞ ζacotζ = 1, one obtains that the leading order BSF cross section Eq. (3.4) is by approximately a factor three larger compared to the s-wave Sommerfeld enhanced annihilation cross section into two photons, see ref. [67] for details. We turn now to BSF processes arising at the first order of interactions with the primordial plasma environment. The imaginary part of the first interaction term in the Dyson Eq. (2.12) of the retarded mediator correlation function defines the next-to-leading order contribution to the spectral function and encodes the first interactions with the environment:

Next-to-leading order
Hereby, the self-energy Π αβ R (P ) depends on the specific model. Under the assumption of a SM photon as mediator, its self-energy is a sum over various contributions with Wbosons, quarks and charged leptons running in the thermal loop. This is equivalent with considering all processes interacting with the plasma at next-to-leading order, for instance BSF via bath-particle scattering and off-shell mediator decay into a bath-particle pair is intrinsically taken into account in the thermal self-energy.
In the following, we concentrate on the interactions with ultra-relativistic Dirac fermions ψ in the primordial plasma such that we identify L env = −gψγ µ ψA µ , which resembles the interaction with SM quarks and charged leptons. Hence, the retarded self-energy contains a fermion loop as illustrated by the left diagram in fig. 4. We derive the retarded self-energy directly from the Wightman functions for massless fermions as described Appendix B. By inserting the obtained expression for the retarded photon self energy into the spectral function Eq. (4.1), and the latter into the general expression for the BSF cross section Eq. (2.11), we obtain the following expression: When integrating the above equation, the imaginary part of the integral contains a product of possible double and single poles. The former originate from the squared photon propagator D 2 and the latter from the self-energy Π R in Eq. (4.1).
Taking the imaginary part of single poles corresponds to putting the fermions in the loop as on-shell, while the photons remain virtual. From top to bottom, the four single poles belong to: BSF via particle 3 and anti-particle scattering, as well as BSF via off-shell decay of the photon into a pair of bath particles and the reverse process. The reverse process of off-shell decay is an exception and does not contribute to the imaginary part for kinematical reason (Note that ∆E is always positive). All others are graphically represented by the top and middle diagram on the right-hand side of fig. 4.
The contributions arising from the imaginary part of the double pole correspond to the on-shell emission at next-to-leading order with a temperature dependent fermion loop. We call these contributions in the following interference terms, collectively represented by the bottom right diagram in fig. 4. The total BSF cross section amounts to a sum over all imaginary parts of the aforementioned single and double poles. The bath-particle scattering and off-shell decay, which arise from the imaginary part of single poles, diverge in the case of a massless photon for zero opening angle of the bath particles (k 1 ·k 2 → 1). We will demonstrate that this divergence is exactly canceled by the collinear divergence appearing in the interference terms. In the conventional Boltzmann approach (i.e. considering only bath-particle scattering without 3 To identify the bath-particle scattering case, the identity 1 + f eq Similarly, all three other possibilities can be identified. We checked that the resulting cross sections from the imaginary part of the single poles only are identical to the cross sections obtained in the conventional Boltzmann formalism, and that they are collinear divergent. on-shell emission at NLO and finite temperature fermion loop), there is no other term regulating this divergence, while in our approach the cancellation will be shown to happen automatically. The loop contains also UV divergent vacuum parts, i.e. f ψ independent terms in Eq. (4.3), which can be renormalized with the standard counter terms.
In the following, we show in section 4.1 the main steps for computing BSF at nextto-leading order for the special case of capture into the ground state. In section 4.2, we provide a general proof for the cancellation of the appearing collinear divergences in our formalism. The result for the next-to-leading order BSF cross section is compared to the leading order contribution (on-shell mediator emission) in section 4.3. We set our results into context to the massive mediator case in section 4.4.

Computational results for direct capture into the ground state
The discussion of the cancellation of the forward scattering divergences and the renormalization of the UV divergences can be separated by noticing that the self-energy can be written as a sum over finite temperature and zero temperature parts. According to this separation, it is useful to define cross sections respectively: where the angular averaged parts are defined as These two cross sections are computed separately in the following. The transition elements T for the direct capture into the ground state are given in Eq. (3.1). Out of computational reasons it is already reasonable to perform here the angular average over the orientation of the DM relative momentum and it naturally arises in the thermal average in Eq. (2.9). The following identities for the direct capture into the ground state will be helpful to simplify appearing expressions:

Vacuum contribution
The renormalized retarded self-energy at zero temperature can be obtained from the conventional Euclidean time-ordered self-energy through analytic continuation. In the MS-scheme, the retarded photon self-energy for vacuum polarization with approximately massless fermions running in the loop is given by: The logarithm has a complex contribution for time-like P 2 , which originates from the offshell decay into a massless bath-particle pair. Computing the integral over the imaginary parts in Eq. (4.5), we realize that the NLO cross section factorizes into the LO contribution as given in Eq. (3.4) and the NLO contribution with (4.11) Performing the integral and taking → 0, we obtain a finite result (see appendix C.1 for the details on the contour integration): These vacuum corrections are shown in fig. 12, shared and discussed in more detail in appendix C.2. From now on, the renormalization scale is fixed to the ground state binding energy µ 0 = E 100 = µα 2 /2. As naively expected at T = 0 from the KLN theorem, the collinear divergences in real and virtual corrections cancel each other. Hence, for a finite, physical cross section, both contributions, off-shell decays as well as the virtual correction to the on-shell emission have to be taken into account. This similarly suggests that when considering off-shell decays in a thermal plasma that we have to consider similarly on-shell emission at next-to-leading order at finite temperatures.
While a consistent treatment of BSF at NLO is non-trivial in the conventional approach (e.g. with respect to double counting of real intermediate states), our advocated formalism will take care of all crucial subtleties such as finite temperature effects, double counting of real intermediate states, and the "automatic" cancellation of all appearing infrared divergences, as we will discuss in the following.
Finite temperature contributions Similar to Eq. (4.10), the finite temperature part of the NLO cross section can be written as where the remaining part to compute is the dimensionless function: 4 (4.14) In this compact notation, the summation over σ i ∈ {+, −} in Eq. (4.13) takes into account all the finite temperature contributions contained in the self-energy in Eq. (4.3). To arrive here, integration over k 2 in Eq. (4.3) was performed over the momentum conserving delta function and k 1 was relabeled by k. The angular integration variable is τ ≡p ·k. The function G σ 1 σ 2 contains a product of the single and double poles, as given by where for the numerator we get In principle, it does not matter which integral is chosen to perform the integration over the poles.
However, we would like to share some insights why we find this particular order especially suited. Here, we have chosen the integral over k2 in Eq. (4.3) to first perform the simple integration over the momentum conserving delta functions, and then p to perform integration over the poles. In this order, the collinear divergence occurs at τ ≡p ·k1 → 1, which is physically not anymore the zero opening angle between the bath particles but related to that divergence. In a different order where p integration over the delta functions is performed first, the collinear divergences occur when the opening angle of the bath particles approaches zero, i.e.,k1 ·k2 → 1. This order introduces coordinate singularities in addition, which one can get rid off by transforming into certain elliptical coordinates. The final result in these coordinates would be the same as in the order of integration chosen in this section, but more difficult to physically interpret. Finally, we would like to remark that one is not allowed to choose different integration orders for the single and double poles. This is because only the simultaneous integral over both contributions is in general finite, while the individual terms can diverge. When separating the divergent integrals over single and double poles, we find that one can still show that the collinear divergences cancel at the end, but one misses important boundary terms and the final result would be different.
For the integration over |p| in Eq. (4.14), we consider the analytic continuation |p| → z and perform the integration over the single and double poles of the function in Eq. (4.15). In appendix C.3, we show that due to the imaginary part in Eq. (4.14), F σ 1 σ 2 holomorphic on the real line, and integration range from 0 to infinity, only the real and positive poles contribute. All the real and positive poles are listed in table 1, together with their existence criteria. We denote the double pole as z 0 and the possible two single poles as z p/m for each σ 1 σ 2 configuration in Eq. (4.15). The table can be used to simply write down the R σ 1 σ 2 functions expressed in terms of the residues as: The different signs in front of the residues originate from the i in Eq. , is subject to further discussion. BSF via bath-particle scattering is contained in the single pole contribution Res(G ++ 0 , z p ) and the finite temperature part of the off-shell decay of the vector mediator into a bath-particle pair is contained in the z p residues in R −+ . The individual residue would diverge in the limit τ → 1, reflecting the collinear divergence. However, the double pole contribution with z 0 diverges in the limit as well, but with opposite sign, resulting in the fact that the sum over double and single pole contributions remains finite in collinear direction. A rather general proof for the cancellation of collinear divergences is presented in the next section. Therein, it is also explained why the terms where no single poles exist are finite. The remaining k and τ integrals are all finite and we show their numerically obtained values in fig. 14, shared in the appendix C.4. Worthwhile to note is that the most dominant contributions are R ++ and R −− , which contain BSF via particle and anti-particle, respectively.

Proof for the cancellation of collinear divergences
From the mathematical point of view, collinear divergences occur since in the collinear limit τ → 1 the single pole z p approaches the double pole z 0 as To make that clear, let us take a closer look on the residues of a function with a double pole at z 0 and a single pole at z p , approaching each other in the collinear limit as, e.g., in Eq. (4.21). Every holomorphic function having such kind of pole structure can be written as where H is holomorphic at z = z 0 and z = z p . Using this general form, the residues of G are given by: In the last line we expanded H(z p ) around z p = z 0 . Clearly, each individual residue is divergent in the collinear limit, when z p → z 0 . However, the colliner divergent terms from Eq. (4.23) and Eq. (4.24) occur with opposite sign and hence cancel each other in their sum. Let us comment on the remaining term 1 2 H (z 0 ) whose appearance one can intuitive understand. Since the double pole merges with the single pole in the collinear limit, the function G| τ =1 = H(z) (z−z 0 ) 3 has a triple pole. And indeed, its residue gives the result expected from Eq. (4.23) and Eq. (4.24) as without any collinear divergence. This proves the collinear finiteness for any H holomorphic at z = z 0 and z = z p and is in particular fulfilled for any F σ 1 σ 2 in Eq. While the residue of the single pole Res(G, z p ) and the residue of the double pole Res(G, z 0 ) diverge individually in the collinear limit, we have shown in summary that the sum of both terms remains finite. In the Boltzmann formalism only the single pole Res(G, z p ) would occur, and hence the collision term would be ill-defined. This makes the thermal field theory approach necessary for studying BSF at higher order, at least if massless gauge bosons or light mediators (for how light see section 4.4) are involved.
It remains to discuss the residue at z m , only present in R −+ 0 . Similarly, Res(G 0 , z p ) and Res(G 0 , z m ) have poles if z p = z m . This situation occurs if |τ | = 2∆E|k|−∆E 2 |k| 2 =: τ * and therefore it is only present in the second line of Eq. (4.19). Because of the relative signs these singularities do not cancel. However, due to the choice of the order of the integration (see discussion in Sec. 4.1, paragraph Finite temperature contributions) the singularities only grow like (τ − τ * ) − 1 2 . For that reason the integrating over τ does not cause any problems. Finally, the last exception occurs if z 0 = z p = z m . This is the case for τ = 1, |k| = ∆E. In this case the function F −+ given in Eq. (4.16) goes to zero and the integral remains finite.  fig. 5. Overall one can recognize that the NLO starts to dominate over the LO cross section for temperatures larger than the absolute value of the ground state binding energy. While the thermally averaged BSF cross section via the on-shell mediator emission increases for decreasing temperature, the NLO contribution becomes larger for higher temperature.

Comparison of ground state capture at leading and next-to-leading order
The reason for this behavior can be understood from the lower panel of fig. 5, where the non-averaged quantities as a function of inverse relative-velocity are shown. In the high temperature regime around the Bohr-momentum κ, the total cross section differs in shape and amplitude compared to the LO. The overall increase in the amplitude towards higher temperatures is mainly caused by the number density of relativistic bath particles, scaling as T 3 . Additionally, the NLO cross section has a stronger dependence on the relative-velocity for v rel α, which turns into the milder scaling v −1 rel for v rel α. The temperature and velocity enhancement leads to the fact that the overall BSF probability distribution in the thermal average is sharply peaked at the DM relative momentum ∼ κ and contributes less at the conventional one ∼ (m χ T ) 1/2 . These features are less pronounced in the on-shell emission case, qualitatively explaining the differences of the thermally averaged quantities.

When is a thermal field theory approach required?
In the Boltzmann formalism, the bound-state formation cross section for bath-particle  Figure 6. Graphic compares the thermally averaged bound-state formation cross section for bathparticle scattering in the conventional Boltzmann and the thermal field theory approach.
scattering is always finite for massive mediators. However in the limit mass to zero, the cross section has a logarithmic divergence, which originates from the forward scattering divergence of the bath-particles, see also ref. [59]. In this section, we would like to answer the following question for practical purposes: what is the critical mediator mass, above which the conventional Boltzmann formalism is a sufficient description and below which a more sophisticated thermal field theory analysis as presented in our work is required?
To estimate the critical mediator mass, only the dominant contributions of R ++ and R −− are considered, which contain the BSF via particle and anti-particle scattering, respectively. Additionally, these functions also contain interference terms originating from the double poles of the photon propagator. For a vector mediator with a mass m V , one can simply replace the photon propagator by a massive one. This changes only the double pole into: where the single pole and F ++ remain unaffected. One can immediately recognize that the double pole at |p| =z 0 = (∆E 2 − m 2 V ) 1/2 only exists for ∆E > m V . This is always true for mediator masses smaller than the absolute value of the binding energy. In general, for massive mediators one can replace R ++ by: As a reminder, the term Res(G ++ 0 , z p ) corresponds to BSF via bath-particle scattering and reproduces ref. [59]. From this equation one can see that the double pole contribution, which causes the difference to the Boltzmann formalism, can be neglected for mediator masses much larger than the binding energy due to θ(z 2 0 ). This statement is confirmed numerically as shown in fig. 6, where we compute the BSF cross section based on the full Eq. (4.27) and without the double pole contribution. The Boltzmann formalism overestimates the effect of the bath-particle scattering for m V E 100 . For m V E 100 , the Boltzmann formalism is reliable.

Impact on thermal relic abundance
The results of the previous sections allow us to explore the impact of higher-order BSF processes on the evolution of the thermal relic abundance for the following simplified model with N bath particle species: An s-wave χ-pair can annihilate into 2γ and N effectively massless ψ pairs. The total s-wave Sommerfeld enhanced annihilation cross section, averaged over the initial spin and summed over final degrees of freedom, can be written as where the first two factors account for the tree-level part and |ψ k,l=0 (0)| 2 = 2πζ(1−e −2πζ ) −1 is the Sommerfeld enhancement factor [9,10]. For the lowest s-wave χ bound states, we consider the spin singlet (S) decay into 2γ [69,70], and the spin triplet decay into 3γ [71] and N ψ-pairs [72]. The decay widths can be written as where |ψ 100 (0)| 2 = (µα) 3 /π, and c 3γ = 4(π 2 − 9)α/(9π) [71]. We compute the BSF cross section in Eq. (2.11) for the capture into the ground state up to NLO in the dark photon spectral function: The effective cross section W (T ) comprises all information of the bound states. It consists out of a sum over the individual singlet and triplet BSF cross sections (here spinindependent), weighted by branching ratios containing the individual decay and dissociation rates. The latter corresponds to the inverse process of the bound-state formation and is therefore related via detailed balance as n eq χ n eq χ n eq 100 , (5.8) where the equilibrium number density of the bound state contains one spin d.o.f. The effective cross section W features two asymptotic regimes: Γ dis i (out-off io. eq.).

(5.9)
In the first asymptotic regime, the bound and scattering states are in ionization equilibrium (io. eq.). Here, W is i) independent of the BSF cross section and dissociation rate [66], and ii) maximum for a given bound-state decay rate and temperature. In fig. 7, the maximum value of W is indicated by the green line for electroweak (left panel) and strong couplings (right panel). In the second asymptotic regime at later times, the number density depletion depends only on the total BSF cross section, since the bound states immediately decay without being dissociated back into the scattering states (BSF without "back reaction"). The transition between these two regimes occurs when the dissociation rate is comparable to the decay rates. We estimate the maximum relative size of W with respect to the SE annihilation cross section 7 , assuming io. eq.: In this model, the ratio is independent of the number N of bath-particle species. In io. eq., W first grows rapidly for decreasing temperature with a power law (m χ /T ) 3/2 , becomes comparable to the SE for temperature around twice the binding energy, and has an exponentially growing factor for lower temperature.
The relevance of NLO BSF contributions depends, for instance, on for how long the LO BSF can keep the system in ionization equilibrium. In fig. 7, we show that for only one bathparticle species in the plasma, the on-shell emission of a vector mediator is still effective enough to keep W close to its maximum value for sufficiently long time. In this case, the effect of the NLO contributions is only marginally relevant. For larger N , however, the relevance of the NLO contributions is more important. This is because the SE annihilation, NLO BSF cross section, and triplet decay rate are proportional to N , while the LO BSF is independent. BSF via bath-particle scattering as the dominant NLO process, can keep the system in io. eq. until temperatures close to the binding energy even for large N .
Moreover, the relevance of NLO BSF contributions depends on the relative size of W compared to the SE annihilation around the decoupling time from ionization equilibrium. Close to the binding energy, the depletion of the DM number density is exponentially sensitive to the ionization decoupling temperature according to Eq. (5.10). While in the conventional Boltzmann framework one can artificially push the decoupling temperature below the binding energy by lowering the vector mediator mass, our thermal field theory result suggests even for a massless vector mediator a decoupling before but close to the strong exponentially enhanced regime. Based on this result, we may expect a significant impact on the predicted relic abundance in the model under consideration mainly for strong couplings α ∼ 0.1.
To demonstrate this more clearly, we choose a strong coupling value and more conservatively N = 3 ("one quark") in fig. 8. Here, the effect of NLO corrections can significantly change the relic abundance due to a delayed decoupling from ionization equilibrium. The Overclosure bound Figure 9. Under the assumption of a strong coupling α = 0.1 and one quark in the plasma, we show the predicted relic abundance when considering a tree-level cross section only (blue), including the Sommerfeld effect (red), considering BSF via LO on-shell emission (pink), as well as BSF at NLO (black). resulting corrections to the upper limit of the DM mass, above which the Universe would contain too much DM (overclosure bound), are shown in fig. 9. For choices of larger N , the differences between LO and NLO would further increase.

Discussion
From the results of previous sections, we learned that DM bound-state formation inside a relativistic thermal bath can be dominated by higher order processes for a certain temperature range for massless mediators. It was demonstrated that next-to-leading order processes contained in the mediator spectral function have the potential to change the thermal relic abundance significantly. Although we analyzed a simplified model, this message should be recognized in the broader context of coannihilation scenarios, where the number of SM bath particles can be much larger than considered here in this work.
The formalism of this work also allows to study the impact of the ambient primordial plasma on excited DM state transitions. One can simply change the initial scattering state into a bound state, which transforms the BSF cross section Eq. (2.11) into a leveltransition rate. The factorization between the spectral function and the level-transition matrix element holds in this case as well. Contributions from higher DM states were however dropped in the previous section and remain a major uncertainty in the prediction of the final relic abundance.
The bound-state formation and the reverse process have also been investigated in the context of heavy annihilating quarkonia inside the quark-gluon plasma produced, e.g., at the Large-Hadron-Collider, which has many similarities compared with heavy annihilating dark matter in the primordial plasma. For fairness, we compare our approach to the exist- Figure 10. Graphic illustrates goodness of different descriptions, based on the assumption that the mediator is coupled to the primordial plasma background. LO and LO+NLO stand for the boundstate formation description via the emission of an on-shell mediator and higher order couplings of the mediator to the plasma, respectively. The bottom bar line represents HTL resummed corrections to the SE annihilation and bound-state decay rate under the assumption that ionization equilibrium is maintained.
ing methods in the quarkonia literature in the following. A similar formalism is developed in ref. [74], where the LO cross section for quarkonium formation and dissociation is derived from the Lindblad equation. In ref. [75], a particular part of cancellation of collinear divergences is suggested while an ad hoc insertion of a finite temperature self-energy into vacuum amplitudes was performed. Our work can be regarded as a first complete analysis at NLO in the mediator spectral function. Starting from first principles, we have derived Eq. (2.11) and demonstrated that the mediator spectral function evaluated at NLO automatically leads to a collinear safe cross section and includes all the possible processes simultaneously. For example, we have shown that even the finite part of the interference terms is important for predicting the relic abundance precisely (see Fig. 14), which cannot be addressed without the full NLO treatment.
The limitation of the formalism presented in this work is set by the assumptions of potential non-relativistic effective field theory, including the leading order dipole approximation. For temperatures larger than the Bohr-momentum κ, the plasma can entirely probe the typical size of the bound state, and the validity of our generalized BSF cross section Eq. (2.11) breaks down. At the critical value T ∼ κ, the bound and scattering states become non-separable because of the rapid transition between them, which calls for another formalism capable of coping with this situation. At the same time, it was shown that BSF and the back reaction is such efficient, that the system is robustly in ionization equilibrium. We emphasize that the collision term is independent of the BSF cross section in such a state (see Eq. (5.9)). By using another formalism valid at high temperature but limited to ionization equilibrium, environmental corrections to the upper formula of (5.9) have been studied in the literature [76][77][78][79][80][81]. The method is based on an effective in-medium potential, which shows next to the expected Debye screening mass, a temperature dependent energy shift (Salpeter correction), as well as an imaginary thermal width [82,83]. A conservative estimate shows that the thermal width in the effective in-medium potential can lead to an entire melting of all bound states for T κ [84,85]. Such melting phenomena have been experimentally observed, e.g., in the decay spectra of bottomonium systems in a quark-gluon plasma, see fig.1 in ref. [86]. In ref. [66], it has been emphasized that the formalism for computing the relic abundance in the effective in-medium potential description is limited to the assumption of ionization equilibrium (see also ref. [87]), illustrated by the bottom bar line in fig. 10.
One of the main formal achievements of this work was to develop a complementary method, which overlaps with the validity region of the previous effective in-medium potential approach. Combined, we can now describe the evolution of the DM system for a broader range of temperatures, ranging from the melting of bound states in the high temperature regime down to arbitrarily below the decoupling from ionization equilibrium. In particular, we have now addressed the temperature regime below the Bohr-momentum, where the thermal width can be treated as a perturbation and bound states can be characterized by their usual quantum numbers. Note again that the expansion given in Eq. (2.3) implicitly assumes a clear separation between the scattering and bound states, which is a crucial assumption when the continuum starts to overlap with the discrete energy spectrum in the high temperature regime. For lower temperatures around the binding energy, these corrections are expected to be negligible compared to our improved description of the decoupling from ionization equilibrium. In our model, we found that the depletion of the dark matter density is exponentially sensitive to the precise decoupling temperature from ionization equilibrium. The presented formalism in this work now allows to accurately resolve the number density evolution especially during this regime.

Summary and conclusion
In the conventional Boltzmann formalism, the amplitude for higher-order bound-state formation processes can become collinear divergent in the case of massless mediators. A large number of possible bound-state formation channels through SM particle scattering, via the virtual exchange of, e.g., photons or gluons in coannihilation scenarios, can not be investigated in this context. Based on non-equilibrium quantum field theory techniques, we derived a more general cross section, Eq. (2.11), which addressed this issue without the need of an ad hoc screening mass regulator. We presented for the first time a full computation and analysis of the thermal one-loop correction for the dominant capture into the ground state in the case of ultra-relativistic fermions in the primordial plasma environment, resembling the interactions with light SM leptons and quarks. For temperatures larger than the absolute value of the binding energy, we found that bound-state formation via bath-particle scattering dominates also for massless vector mediators over the so far considered on-shell mediator emission. Our results are complementary to the one in ref. [59], where massive mediators where investigated.
The key quantity in the generalized bound-state formation cross section is the spectral two-point correlation function of the interacting mediator. It was demonstrated that a perturbative expansion in the coupling parameter of this spectral function successively generates on-shell and virtual mediator contributions in a proper thermal field theoretical framework. The known result for the capture into the ground state via the emission of an on-shell mediator was reproduced at the lowest order in the perturbative expansion. The higher-order BSF processes via bath-particle scattering and off-shell mediator decay were identified at the first interaction level. These processes are collinear divergent for massless mediators in the conventional Boltzmann approach. It was shown that other terms automatically occur in the spectral function in addition, canceling the collinear divergences and resulting in a finite collision term. A rather general analytic proof for the collinear finiteness, applying in particular to our full first interaction term, was presented in section 4.2.
Based on our extended analysis in section 4.4, we conclude that a thermal field theory approach is required for models where the mediator mass is smaller than the absolute value of the binding energy. This implies that our approach can become important also for massive SM gauge bosons or the Higgs field during the electroweak cross over, featuring a large number of possible BSF channels via SM particle scattering as well.
In the case of our simplified model, we demonstrated in section 5 that the impact of the new higher-order bound-state formation effects on the thermal relic abundance can be especially significant if the on-shell emission is not enough to keep the system in ionization equilibrium until temperatures close to the binding energy. Regarding more realistic coannihilation scenarios, the size of the corrections to the upper bound on the DM mass from higher-order bound-state formation processes still remains an open question from the perspective of this work, since the details of the underlying model might play a crucial role. This is mainly due to the fact that the impact on the relic abundance is exponentially sensitive to the precise decoupling time from ionization equilibrium if it is maintained until temperatures around the binding energy. Model-dependent order one variations in the ionization equilibrium decoupling temperature are crucial for the impact of the higher-order effects. With an accurate formalism by hand, together with the insights from simplified DM models, makes it nevertheless worthwhile to analyze more realistic scenarios in future work. and bound states: where the overlap integral is given by Here we use a shorthand notation for integrals: k ≡ d 3 k (2π) 3 and P ≡ d 4 P (2π) 4 . The positive quantity ∆E ≡ E k − E B is the total energy emitted in the inelastic conversion, i.e. relative kinetic energy plus the absolute value of the negative binding energy.
All one needs to do is to evaluate the double commutator in the right-hand side of Eq. (2.4) by using Eq. (A.1) in order to obtain the collision term for the BSF and dissociation rates. Below we summarize commutators and expectation values relevant for our computations: for commutators, and for expectation values of creation and annihilation operators. Here note that these expectation values are taken by the initial factorized density matrix: • = Tr[•ρ(t = 0)]. We also need two-point functions of the electric field defined by For later convenience, we rewrite this electric correlator in terms of its gauge field as: . Since the expectation value is taken by the factorized initial density matrix and the environment is assumed to be in thermal equilibrium, these correlators coincide with the thermal correlator, fulfilling the KMS relation: where the photon spectral function is given by D ρ µν (x−y) ≡ [A µ (x), A ν (y)] . f eq γ represents the Bose-Einstein distribution for the gauge boson.
Using these equations, one finally obtain the collision term associated with the BSF and dissociation: with ∆E = k 2 /(2µ) + |E B |, k = (k χ − kχ)/2, and K = k χ + kχ. Here we have dropped the Bose enhancement and Fermi suppression factors for the scattering and bound states because they are dilute. To obtain this compact form, we have used the following relation between the overlap integral and tensors defined in Eqs. (2.7) and (2.8): where g χ and gχ represent the spin degrees of freedom for particle and anti-particle respectively. The collision term can be expressed as the following simple form One might wonder why this cross section contains 1 + f eq γ (∆E) even for the inverse process. To avoid this confusion, we provide an explicit proof for the inverse process starting from Eq. (A.8): In the first line we insert the definition of the generalized cross section. In the second line we use the relation f eq γ (∆E)/[1 + f eq γ (∆E)] = e −∆E/T and the approximation of kinetic equilibrium f B = f eq B n B /n eq B . In the third line we use the diluteness of the bound state f eq B = e −(E B +K 2 /8µ)/T and the energy conservation e −∆E/T e −(E B +K 2 /8µ)/T = e −k 2 χ /4µ e −k 2 χ /4µ . This completes the proof of the second term in Eq. (A.10).

B Retarded self-energy for massless fermions
The retarded photon self-energy is defined in terms of greater and lesser self-energies as where the two-point functions of the fermionic bath-particles are defined as In the free limit and in thermal equilibrium, their Fourier transform is given by f eq ψ is the Fermi-Dirac equilibrium phase-space distribution. In the free and equilibrium limit, the terms needed for the self-energy can be expressed in Fourier space as: Where α and β are short notation for the phase space density and the Pauli blocking, respectively.
Using these equations, the retarded self-energy can be written in Fourier space as After rearrangement, the latter equation can be brought into a convenient form which is used in section 4 to analyze the various NLO contributions.
C Next-to-leading order contributions in more detail In this appendix, we share some details of the next-to-leading order computation of the mediator spectral function, as well as some more detailed discussion about the individual contributions.

C.1 Contour integration for vacuum part
To evaluate the vacuum part of the cross section in Eq. (4.10), one has to compute where N specifies the number of particle species in the thermal plasma and thus also the multiplicity of particles that are running in the loop in the NLO contribution to the on-shell emission. As expected from the renormalization procedure of the occurring UV divergence, the final cross section will feature a renormalization scale dependence µ 0 that needs to be fixed by a physical motivated scale. As the characteristic scale of the BSF process lies around the binding energy, we choose accordingly µ 0 = E 100 = m χ α 2 /4. In fig. 12, we show the vacuum contribution at LO and NLO in dependence of the inverse relative velocity v rel for coupling strengths up to the order of the unitarity bound for s-wave SE annihilation, α = 0.5 [88]. In order to estimate the uncertainty that arises from this particular choice of scale, we vary the latter between µ 0 = 2E 100 and µ 0 = E 100 /2 (shaded area). The upper panel features the general expected behaviour: With larger velocity the relative cross section decreases as BSF is suppressed due to the small overlap of the wave functions. Moreover, larger couplings α lead to a relatively larger BSF cross section. From Eq. C.9, we can infer that for larger couplings the NLO correction becomes more relevant and hence also the dependence of the scale uncertainty, as visible in the lower panel. Note that we have taken the conservative choice of N = 1. For larger N not only the correction itself, but also the scale uncertainty is expected to be more pronounced.
For the renormalization of the occurring UV divergence, we had chosen the simplest choice for the renormalization scheme, namely the MS scheme. However, a more physical motivated approach would be to use an on-shell scheme for the renormalization of the photon propagator since the photon is on-shell in the corresponding interference terms. However, this unnecessarily complicates the cancellation of the infrared divergences. While it would be academically interesting to study the application of the on-shell scheme and to compare it with the MS scheme in more detail, we leave this investigation for future work.

C.3 Contour integration for the finite temperature part
The function G σ 1 σ 2 (z) = G σ 1 σ 2 (z, τ, |k|) in the finite temperature part of the cross section in Eq. (4.13) has three types of poles: • the double poles of the squared mediator propagator at ±(∆E + i ), • the single poles of the retarded self-energy at w p/m , w p lying in the upper half and w m lying in the lower half of the complex plane, • and the poles of F σ 1 σ 2 , which do not lie on the real line.
When goes to 0, w 0 ≡ ∆E + i goes to z 0 = ∆E which lies on the positive real line. w p and w m go to z p/m = lim 0 w p/m ≡ −τ σ 1 |k| ± sign(∆E + σ 1 |k|) τ 2 |k| 2 + ∆E 2 + 2σ 1 ∆E|k|, (C. 10) which can lie on the positive real line, as well. Figure 13. Contour used to compute the finite temperature part. Black dots: poles of the function G σ1σ2 ; z 0 , z p , z m : location of the poles of G σ1σ2 0 which lie on the real line; w 0 , w p , w m : location of the poles of G σ1σ2 which go to z 0 , z p , z m as 0; α 1 , α 2 , α 3 , α 4 : paths on the real line; β 0 , β p , β m : paths in the complex plane (half circles) avoiding the singularities of G σ1σ2 .
For the moment let us assume that z p/m are real and positive. To avoid integrating through the singularities we deform the path of integration as shown in fig. 13. The deformation of the path does not affect the value of the integral over |p| as F σ 1 σ 2 is holomorphic in some neighborhood of the real line. Since G σ 1 σ 2 0 has no poles on the new path γ, we get lim 0 ∞ 0 d|p| G σ 1 σ 2 (|p|, τ, |k|) = lim 0 γ dz G σ 1 σ 2 (z, τ, |k|) = γ dz G σ 1 σ 2 0 (z, τ, |k|). (C.11) As G σ 1 σ 2 0 (z) is real valued on the real line, the integration along α 1 , . . . , α 4 only gives a real contribution. Moreover, it follows that the Laurent coefficients of G σ 1 σ 2 0 at z = z 0 , z = z p , and z = z m , respectively, are real. Using Laurent series expansion it is straightforward to verify that If one of the poles z p/m is negative or not real then the corresponding term in Eq. (C.14) vanishes as there is no need to deform the path of integration. In particular, by the same argument the poles on the imaginary axis, originating from the equilibrium phase-space distributions in Eq. The R σ 1 σ 2 functions, in Eq. (4.17) to Eq. (4.20), are all dimensionless. Therefore they can only dependent on the ratio ∆E/T for massless mediators. This fact makes the relic abundance computation efficient, since one can tabulate the sum over all these functions in advance depending on a single variable. For a selected range we show the individual contributions in fig. 14. R ++ and R −− contain the BSF via bath-particle and anti-particle scattering and are the most dominant contributions. The double pole in R +− contributes negatively, and R −+ has some oscillating features. The latter contains the finite temperture part of the off-shell mediator decay into bath particles as well as other single and double pole residues. Depending on which part dominates inside R −+ , the function takes positive and negative values.
Most important is the value of the sum over all functions at T /∆E ∼ 1, since there one would roughly enter the exponential depletion phase of the relic abundance if ionization equilibrium is maintained. One can see that for a very precise determination of the relic abundance, the inclusion of all terms are needed at least in this case. Throughout this work we take all of them into account.