Effective field theories for dark matter pairs in the early universe: cross sections and widths

In order to predict the cosmological abundance of dark matter, an estimation of particle rates in an expanding thermal environment is needed. For thermal dark matter, the non-relativistic regime sets the stage for the freeze-out of the dark matter energy density. We compute transition widths and annihilation, bound-state formation, and dissociation cross sections of dark matter fermion pairs in the unifying framework of non-relativistic effective field theories at finite temperature, with the thermal bath modeling the thermodynamical behaviour of the early universe. We reproduce and extend some known results for the paradigmatic case of a dark fermion species coupled to dark gauge bosons. The effective field theory framework allows to highlight their range of validity and consistency, and to identify some possible improvements.


Introduction
One of the major challenges in cosmology and particle physics is to understand the matter content of the universe. Notably, visible ordinary matter appears to be only a small fraction of the matter in the cosmos, whereas the bulk seems to come in the form of non-luminous and non-baryonic particles, dubbed dark matter (DM). Complementary measurements of large scale structures, galaxy formation, gravitational lensing and the cosmic microwave background (CMB) suggest that more than 80% of the matter in the cosmos consists of DM.
The most accurate determination for the DM energy density is provided by anisotropies in the CMB and amounts to Ω DM h 2 = 0.1200 ± 0.0012 [1], where h is the reduced Hubble constant.
Although this is not the only viable option, the interpretation of DM as due to new, yet undiscovered, particles has been put forward in a plethora of models, see e.g. [2,3] for extensive reviews. One can categorize DM particles according to their production mechanism in the early universe.
The thermal freeze-out has been widely adopted to infer the present-day abundance of a DM candidate. It allows to effectively link the particle model parameters, such as couplings and DM mass, with the observed relic density. The freeze-out mechanism has been extensively used for Weakly Interacting Massive Particles (WIMPs), however it equally applies to cases where interactions are stronger. One assumes an initial thermal abundance for the dark species, whose evolution is governed by the interplay between the thermally averaged annihilation cross section σ ann v rel and the expansion rate of the universe, H, namely the Hubble rate. The standard rate equation is a Boltzmann equation of the form [4-6] 1 (∂ t + 3H)n = − 1 2 σ ann v rel (n 2 − n 2 eq ) , (1.1) where n is the total number density of DM particles (n eq is that in equilibrium 2 ) and v rel = |v 1 − v 2 | is the relative velocity of the annihilating pair. For a DM particle of mass M , the chemical freeze-out occurs at a temperature T ≈ M/25. Therefore, at freeze-out the dark matter particles are non-relativistic. 3 It is crucial to calculate σ ann v rel accurately because the present-day DM energy density, as predicted by a given model, depends on it through the solution of eq. (1.1). The DM mass is in turn fixed as a function of the other model parameters to reproduce Ω DM h 2 . A solid prediction for DM mass benchmarks compatible with the relic density is needed to establish the viability of models, guide the experimental searches and put DM phenomenology on a sound theoretical ground. However, determining σ ann v rel by including the full features of each model, and the thermal environment, is not straightforward.
In a variety of theories, DM interacts with gauge bosons or scalars that induce longrange interactions because of repeated soft exchanges. While a successful WIMP-like DM candidate should be weakly interacting with the Standard Model (SM), we cannot say much on the interactions between the DM particles themselves and/or among additional degrees of freedom in the dark sector. In general, there may be non-negligible forces between DM particles mediated by light particles leading to the formation of bound states of genuine WIMPs, sometimes referred to as wimponium [7,8]. The main motivation that makes self-interacting DM welcome are some compelling inconsistencies between the predictions of collisionless cold DM and the observed large-scale structure of the universe [9], in the numbers of the galactic satellite halos, and in the DM density profiles in the galaxies [10][11][12][13][14][15][16][17][18][19][20].
Bound-state effects cannot be avoided in the context of next-to-minimal DM models, where a mediator between the visible and dark sector is charged under some of the SM gauge groups. In coannihilation scenarios, the presence of additional massive states can drastically affect the thermal freeze-out when the coannihilating partner is close in mass with the DM particle [6,21]. Consequently, one has to track also its (co)annihilations. For example, when the partner particle is charged under QCD, long-range interactions mediated by gluons affect severely the annihilation process and the formation of (many) bound states has to be included [22][23][24][25][26].
Accordingly, in recent years there has been some effort in revisiting dark-matter pair annihilation by encompassing near-threshold effects induced by repeated soft exchange. First, the so-called Sommerfeld enhancement has been included in DM freeze-out calculations for several models, resulting in mass benchmarks that give DM energy densities compatible with observation rather different from the case with no threshold effects [27][28][29][30][31][32][33][34][35][36]. Typically, the Sommerfeld enhancement increases the annihilation cross section, for a pair in an attractive channel, and leads to a larger dark matter mass compatible with the observed relic density. More recently, bound-state effects have shown a relevant impact on DM annihilations as well, as pointed out originally in refs. [30,37]. Indeed, DM bound states contribute as an efficient annihilation channel and dilute further the overall DM population. The presence of meta-stable bound states is simply another manifestation of repeated soft exchanges: in the spectrum of a two-particle system there is an above threshold continuum of states along with bound states below threshold.
Treating interacting non-relativistic particle pairs in a thermal environment is complicated by the presence of several energy scales. The energy scales dynamically generated by the relative motion are the momentum transfer, which is also proportional to the inverse of the typical size of the pair, and the kinetic/binding energy of the pair. Such scales are hierarchically ordered as M M v rel M v 2 rel for near threshold particles moving with relative velocity v rel . For Coulombic bound states, the relative velocity of the pair is fixed by the virial theorem to be v rel ∼ α, hence the corresponding hierarchy is M M α M α 2 , α = g 2 /(4π) being the fine structure constant in terms of the coupling g between the DM particle and the force mediator. The thermodynamical scales include the plasma temperature T and the Debye mass m D , which is the inverse of the chromoelectric screening length; in a weakly-coupled plasma m D ∼ gT . For kinetically equilibrated particles and pairs the typical momentum is of order √ M T and the typical relative velocity is then v rel ∼ T /M ; clearly, after freeze-out v rel 1. Depending on the plasma temperature, the scale √ M T may be larger or smaller than the scales M α and M α 2 .
In this work, we compute production cross sections and decay widths of near threshold weakly-coupled dark matter particle-antiparticle states in a thermal bath of SM particles describing the early universe. We exploit the hierarchy of energy scales typical of near threshold states and a weakly-coupled plasma by replacing the fundamental DM theory with a sequence of non-relativistic effective field theories (EFTs). In practice, in the main body of the paper, we deal with the simple case that consists in having a DM sector made of a new species of fermions, dark fermions, and photons, dark photons, that interact like in QED, see section 2. To remark the new particle content of the theory, we call it QED DM . The EFT that follows from QED DM by integrating out modes associated with the energy scale M is non-relativistic QED (NRQED DM ) [38]. The EFT that follows from NRQED DM by integrating out dark photons of energy or momentum of order M v rel is potential nonrelativistic QED (pNRQED DM ) [39][40][41]. We sketch the tower of energy scales and EFTs considered in this work in figure 1. Potential NRQED DM is made of low-energy particles only and well suited to compute near threshold observables. Many of them have been computed also elsewhere, hence the main motivation of this work is to present a systematic and coherent framework for these computations highlighting at the same time their regime of validity and possible improvements. We consider the case of a DM sector made of dark fermions coupled to SU(N ) dark gauge bosons in appendix D. Near threshold effects not only affect directly the annihilation cross section in eq. (1.1), but also add to it new production and decay mechanisms that may eventually modify the abundance of DM particles significantly. We summarize them in figure 2. The diagram labeled a) shows a transition from a scattering state to a bound state via a gauge boson emission [23,25,37,[42][43][44][45][46][47][48]. The diagram labeled b) shows an inelastic collision of a scattering state with a constituent of the thermal bath that turns it into a bound state [24,[49][50][51][52][53][54][55][56][57][58]. Finally, diagram c) shows the decay of the emitted gauge boson mediator, if this couples to some light degrees of freedom in the thermal bath of the early universe [57,58]. The latter two processes entail further model-dependent details on the interactions among the mediator and the additional light degrees of freedom. The relative importance of the various mechanisms depends on the temperature. In this work, we focus on the boundstate formation/dissociation as induced by the radiative emission of a gauge vector boson, namely process a). In the hierarchy of scales that we assume in our work (cf. eq. (2.2)), this is expected to be the dominant process [59]. Moreover, we do not include in our model couplings of the dark photons with light degrees of freedom, this excludes explicitly the possibility of having the processes b) and c). It also excludes the generation of a Debye mass. Recent investigations of near-threshold effects induced by scalar mediators and their applications to dark matter models can be found in refs. [60][61][62][63][64][65]. The outline of the paper is as follows. In section 2, we introduce the DM model, and in sections 2.1 and 2.2 its low-energy effective field theories upon integrating out modes carrying energies and momenta of order M and M v rel . In section 3, we consider annihilations, and derive the corresponding observables in section 3.1 and 3.2. Electric transitions are covered in section 4. Using the optical theorem, we derive the bound-state formation cross section in section 4.1 and the bound-state dissociation thermal width in section 4.2. Details of the derivations and the treatment of excitations, de-excitations and bremsstrahlung effects are collected in the appendices A, B and C. Appendix D discusses a non-abelian SU(N ) model. In section 5, we scrutinize the interplay between the coupling and velocity expansion, and the inclusion of excited states on the dark matter energy density. Moreover, we present contour plots for the dark matter mass and coupling. Finally, conclusions and outlook are in section 6.

EFTs for non-relativistic dark matter pairs
We consider a simple model where the dark sector consists of a dark Dirac fermion X that is charged under an abelian gauge group [66][67][68][69][70]. We denote the corresponding dark photon with γ. The Lagrangian density reads where the covariant derivative is D µ = ∂ µ + igA µ , with A µ the dark photon field and we define the fine structure constant as α ≡ g 2 /(4π). The term L portal comprises additional interactions coupling the dark photon with the SM degrees of freedom. A popular interaction is a mixing with the neutral components of the SM gauge fields [71,72]. Such interactions are responsible for the eventual decay of the dark photons, in this way avoiding that their number density dominates the universe at later stages. Moreover, additional fermionic or scalar degrees of freedom may be coupled to the dark photon. In an abelian theory, and at temperatures T M , they are responsible for quantum corrections to the dark photon propagator, whose pole develops a real and an imaginary part [73][74][75]. 4 The real part introduces a screening (Debye) mass m D of order gT for the temporal dark photon, whereas the imaginary part of the pole originates from 2 → 2 scatterings with plasma constituents, also referred to as Landau damping. It is beyond the scope of this work to elaborate further either on L portal or on the physical effects of additional fermionic or scalar degrees of freedom. From now on we set L portal = 0.
The Lagrangian (2.1) describes also processes involving two dark fermions close to threshold, i.e. processes where the fermions are non-relativistic and move with relative velocity v rel 1.
For v rel ∼ α, (ladder) photons exchanged between the pair contribute with a relative factor of order α/v rel ∼ 1 and need to be resummed. The resummation generates bound-state poles of order M α 2 at negative energies and a continuous scattering spectrum at positive energies. The typical momentum exchanged between the pair when v rel ∼ α is M α, which is of the order of the inverse Bohr radius of the bound state. The dynamically generated scales M α and M α 2 are the more separated the smaller α is: M M α M α 2 . We call them soft and ultrasoft scales, respectively, to distinguish them from the hard scale associated with the mass M . These energy scales affect significantly various processes in the near threshold momentum region, like dark fermion pair annihilation, formation and transition via emission or absorption of photons. The use of the full Lagrangian (2.1) to compute near threshold observables is in general unpractical as all energy scales get entangled in the amplitudes. It is more convenient, instead, to take advantage of the fact that the energy scales are hierarchically ordered and replace systematically (2.1) with a hierarchy of (non-relativistic) effective field theories along what has been done for near threshold fermion pairs in QED and QCD [38,40,41].
Another relevant energy scale is the inverse correlation length characterizing the medium made of dark fermions, dark photons and SM particles. We assume the medium to be thermalized and identify the inverse of its correlation length with the temperature T . If the dark fermions are also thermalized then v rel ∼ T /M .
In the following, we compute near threshold observables affecting the evolution of the dark matter density in the early universe. In particular, we compute annihilation, dissociation and formation cross sections of dark matter fermion-antifermion pairs. We compute these quantities by means of the tower of non-relativistic effective field theories depicted in figure 1. We include thermal effects due to the medium assuming that the temperature T is about the ultrasoft scale M α 2 or smaller. This guarantees that thermal effects do not enter the potential, which may be taken as Coulombic. The typical momentum of the thermalized dark fermions is then M v rel ∼ √ M T M α, which implies v rel α. Moreover, in the temperature ranges considered in this work it also holds that √ M T M α 2 (more precisely √ M T M α 2 /4). These conditions qualify M v rel as a soft scale and M v 2 rel ∼ T as an ultrasoft scale. Our ensemble of thermalized dark fermions and antifermions realizes, therefore, the hierarchy of energy scales shown in figure 1: The hierarchy (2.2) is of phenomenological interest for the study of near threshold effects in the minimal dark matter model under consideration [37,46]. Indeed, since the decoupling from chemical equilibrium happens at around T /M ≈ 1/25, the condition (2.2) may be fulfilled for most of the time after the decoupling. 5 In the model (2.1) with L portal = 0, bound-state formation happens via radiative emission. For complementary temperature regimes and bound-state formation processes, the latter triggered by additional degrees of freedom added to the model (2.1), see also refs. [49,50,52,56,57].

NRQED DM
At energies much smaller than M , the effective degrees of freedom are non-relativistic dark fermions and antifermions, low energy dark photons and SM particles. The effective field theory that follows from (2.1) by integrating out dark photons and fermions of energy or momentum of order M has the form of NRQED [38]. It is organized as an expansion in 1/M and α and its Lagrangian density up to O(1/M 2 ) reads Here, ψ is the two-component Pauli spinor that annihilates a dark matter fermion, χ † is the Pauli spinor that annihilates an antifermion, E is the (dark) electric field, E i = F i0 , and B is the (dark) magnetic field, B i = − ijk F jk /2. The first two lines in eq. (2.3) describe how non-relativistic dark fermions and antifermions propagate and interact with low-energy dark photons of energy smaller than M . The third line describes the propagation and effective self interaction of the photons; it also contains two four-fermion operators.
To keep track of the thermalization of the physical fields, we do not redefine the fermion and antifermion fields ψ and χ to reabsorb their mass terms, which we leave explicit. In the matching, the thermal scales and any other energy scale below M can be set to zero. Hence, upon our assumption M T , no finite temperature effect enters the EFT Lagrangian (2.3). The one-loop expressions of the matching coefficients c F , c D , c S and d 2 in the MS scheme can be found in ref. [79] taking the abelian limit. The coefficients of the kinetic operators are fixed to be one at all orders in the coupling by reparametrization (Poincaré) invariance [80,81]. As for the four-fermion dimension-six operators, the matching coefficients read at order α 2 [82] where µ is the renormalization scale associated with factoring hard from low energy scale contributions. The renormalization of the coupling is discussed at the end of section 3.
The imaginary parts in the four-fermion operator matching coefficients originate from the particle-antiparticle annihilation diagrams. Annihilation processes happen at the scale 2M and are therefore integrated out in the non-relativistic EFT. The four-fermion operators shown in eq. (2.3) encode the annihilation of S-wave fermion-antifermion pairs. Higherdimensional four-fermion operators encode the annihilation of fermion-antifermion pairs with non-vanishing orbital angular momentum. For instance, dimension eight four-fermion operators encode the annihilation of P-wave fermion-antifermion pairs. The leading order contribution to the imaginary part of the dimension-six operators comes from the twophoton annihilation process XX → γγ, when the XX pair is in an S-wave, see figure 3; the imaginary part in eq. (2.4) may be obtained by cutting the loop diagrams along the photon propagators [83,84]. For the computation of the annihilation cross section, we refer to section 3.1. In observables, the µ dependence of the matching coefficients cancels against low-energy matrix elements.
Matching between annihilation diagrams in the relativistic theory at one loop (left diagrams) and the corresponding four-fermion interactions in the non-relativistic EFT (right diagrams). For initial and final pairs in an S-wave, the latter correspond to the spin-singlet and spin-triplet four-fermion operators in (2.3). The associated matching coefficients are given in (2.4) and (2.5), where only d s has an imaginary part at order α 2 . Instead we have Im(d v ) = O(α 3 ), reflecting the fact that spin singlets decay into two dark photons, whereas spin triplets decay into three dark photons. Explicit expressions at order α 3 are given in eqs. (3.4) and (3.5).

pNRQED DM
Consistently with the hierarchy of energy scales (2.2), the next degrees of freedom to integrate out to describe threshold phenomena at the ultrasoft scale are dark photons of energy or momentum of order M v rel , which encompasses the scales M α and √ M T . At energies much smaller than M v rel the dynamical degrees of freedom are dark fermions and antifermions with energy of order M v 2 rel , and ultrasoft dark photons with energy and momentum of order M v 2 rel , which encompasses the scales M α 2 and T . The effective field theory that describes them has the form of potential NRQED (pNRQED) [39,85] and we denote it as pNRQED DM . The case of pNRQED at finite temperature has been studied in refs. [74,75]. We can proceed as in section 2.1 and integrate out the soft scale by setting to zero all lower (ultrasoft) energy scales, which include the temperature characterizing the thermal distribution of the dark photons. The matching is done at weak coupling, i.e. order by order in α, although the EFT is suited to accommodate a non-perturbative framework as well [41,86].
Threshold phenomena affect fermion-antifermion pairs, hence it is convenient to project the EFT on the fermion-antifermion space and express it in terms of gauge singlet fermionantifermion bilocal fields φ(t, r, R), where r ≡ x 1 − x 2 is the distance between a fermion located at x 1 and an antifermion located at x 2 and R ≡ (x 1 + x 2 )/2 is the center of mass coordinate. In order to ensure that the photons are ultrasoft, photon fields are multipole expanded in r. Hence the pNRQED DM Lagrangian density for the dark matter theory (2.1) is organized as an expansion in 1/M , inherited from NRQED DM , r and α (at weak coupling): where H(r, p, P , S 1 , and S 1 = σ 1 /2 and S 2 = σ 2 /2 are the spin operators acting on the fermion and antifermion, respectively. The static potential is the Coulomb potential: Because T M α 2 , the potential does not get, by construction, thermal contributions at any order. The power counting of the EFT goes as follows: the inverse of the relative coordinate r scales like M v rel , whereas the inverse of the center-of-mass coordinate R can at most change by M α 2 or T , if the DM fermion-antifermion pair recoils against ultrasoft dark photons. The fact that the variation in R is larger than r guarantees the validity of the multipole expansion. The dots in eq. (2.6) stand for irrelevant operators of higher order in the 1/M and multipole expansion. The relative momentum p = −i∇ r and the center of mass momentum P = −i∇ R are the conjugate variables of r and R, respectively.
The expression of the potential V (r, p, P , S 1 , S 2 ) in the center of mass frame including V (0) , V (1) and V (2) can be found in ref. [40]. Here, besides in the static potential V (0) , we are interested in the contributions to V (r, p, P , S 1 , S 2 ) that are responsible for the annihilation or creation of dark fermion-antifermion pairs. Because annihilation happens at the energy scale 2M , such contributions are inherited from the imaginary parts of the four-fermion operators in L NRQED DM and enter V , starting from the V (2) /M 2 term, as where S = S 1 + S 2 is the total spin of the pair and the dots stand for terms of higher order in the EFT power counting. At order r 0 in the multipole expansion, the equation of motion of the fermion-antifermion pair is a Schrödinger equation with potential V (r, p, P , S 1 , S 2 ). Hence the leading order fermion-antifermion propagator in pNRQED DM automatically accounts for bound-state effects and multiple Coulomb rescatterings (Sommerfeld enhancement) in physical observables. We discuss some of the higher-order corrections to eq. (2.10) and multiple Coulomb exchanges in section 3. At order r, the term φ † (t, r, R) r · E(t, R) φ(t, r, R) in the pNRQED DM Lagrangian describes the electric dipole interaction of the dark fermion-antifermion pair with ultrasoft dark photons that comprise thermal photons. The matching coefficient of the electric dipole interaction has been taken equal to one. Fermion-antifermion pairs above threshold form scattering states of positive energy and fermion-antifermion pairs below threshold form bound states of negative energy. It may be therefore convenient to decompose the fermionantifermion field φ(t, r, R) into a scattering component and a bound-state component [87], where φ † n (P ) creates a bound state, |n, P = φ † n (P )|0 , with center of mass momentum P , quantum numbers n, energy E n and wavefunction Ψ n (r) S ij , whereas φ † p (P ) creates a scattering state, |p, P = φ † p (P )|0 , with center of mass momentum P , relative momentum p, energy E p and wavefunction Ψ p (r) S ij . The indices i, j are spin indices. In particular, S-wave dark fermion-antifermion pairs may be either in a spin-singlet state, in which case S ij = δ ij / √ 2, or in a spin-triplet state, in which case S ij = (σ · ) ij / √ 2, where σ are the Pauli matrices and is the polarization vector of the spin-triplet pair. The sum over spin in the second line of eq. (2.11) is a sum over all spin configurations; in the first line, this sum is included in the sum over the quantum numbers n. If the dark fermion-antifermion pair is bound we may call it darkonium, which, in the S-wave case, we may further distinguish into a spin-singlet paradarkonium state, and a spin-triplet orthodarkonium state. The various transitions between scattering and bound states induced by the electric dipole vertex are shown in figure 4. We discuss these transitions in section 4.

Annihilation cross section
The dynamical quantity entering the rate equation (1.1) is the thermal average of the annihilation (ann) cross section, σ ann , times the relative velocity. The annihilation processes that we consider are XX → γγ or XX → γγγ .
The thermal average is defined as the normalized integral of the cross section over the incoming momenta weighted by the Boltzmann distribution [5,32]. This assumes kinetic equilibrium of the dark matter particles at and after the chemical freeze-out, and that the mass M of the dark matter particle is much larger than the temperature, so that we can describe the statistical distribution of the dark matter particles by means of the Maxwell-Boltzmann distribution. Thermal effects contribute to the fermion-antifermion annihilation as ultrasoft corrections to the fermion-antifermion wavefunction. These corrections are subleading with respect to the dynamics governed by the pNRQED DM Hamiltonian in (2.7). If we neglect thermal effects, σ ann v rel is an in-vacuum quantity that is Lorentz invariant and independent of the center of mass momentum P . Hence, in the thermal average the integral over P factorizes and cancels against the normalization; the same happens to the statistical factor e −2M/T . We are left with where v rel is related to the relative momentum p by M v rel /2 = p. Dark fermion-antifermion annihilation happens at the energy scale 2M . At this energy scale the process is best described by NRQED DM , as we see in section 3.1. Multiple soft dark photon exchanges between the fermion-antifermion pair may however significantly distort the fermion-antifermion wavefunction either above threshold, when they form a scattering state, or below threshold, when they form a bound state. Multiple soft photon exchanges are best described by pNRQED DM , as we see in section 3.2.

Annihilation cross section in NRQED DM
Fermion-antifermion annihilation is encoded in NRQED DM in the imaginary part of the matching coefficients multiplying the four-fermion operators. Accounting only for dimension-six four-fermion operators and describing the fermion and antifermion as free noninteracting spinors we get from the optical theorem the cross section 6 where M NR (ψχ → ψχ) is the non-relativistic fermion-antifermion scattering amplitude, and Im(d s ) and Im(d v ) are the imaginary parts of the dimension-six NRQED DM operators introduced in section 2.1. At order α 3 they read [88,89] Im(d s ) = πα 2 The quantity σ NR ann v rel does not depend on the momentum of the fermion-antifermion pair. At leading order (LO), when inserting the explicit expression Im(d s ) = πα 2 , one recovers the well known result [90] At order α 2 the annihilation occurs when the fermion-antifermion pair forms an Swave in a spin-singlet configuration, while the annihilation of the fermion-antifermion pair in a spin-triplet S-wave may only occur via an odd number of photons, i.e. starting from order α 3 . Velocity suppressed contributions to the cross section can be systematically included by adding four-fermion operators of dimension higher than six to the Lagrangian (2.3). Higher dimensional four-fermion operators also account for the annihilation of fermion-antifermion pairs with higher orbital angular momentum.
The EFT provides a straightforward framework to compute higher-order corrections to the annihilation cross section, either in terms of α corrections to the matching coefficients of the four-fermion operators, or in terms of subleading terms in the velocity expansion by including higher dimensional four-fermion operators to the NRQED DM Lagrangian (in this respect see also refs. [91,92] for neutralinos DM in the context of the MSSM and [83] for the case of heavy quarkonium in QCD). Computing next-to-leading order (NLO) annihilation cross sections has been for long pursued in a variety of models, e.g. [93][94][95][96][97][98][99][100]. Neglecting the effect of velocity suppressed operators and multiple soft photon rescatterings, which is important near threshold and we discuss in the next section, the expression at NLO in α of the annihilation cross section in the abelian dark matter model (2.1) reads The corrections are negative and make the cross section smaller. Taking α ≈ 0.4, the NLO cross section is reduced by about 17% with respect to the LO cross section (even larger couplings have been considered in the literature, see e.g. [37]).

Annihilation cross section in pNRQED DM
Multiple soft dark photon exchanges between the annihilating dark fermion-antifermion pair modify the fermion-antifermion wavefunction near threshold and lead to a significant change in the annihilation cross section. The annihilation process of S-wave fermionantifermion pairs is described in pNRQED DM by the imaginary local potential (2.10) directly inherited from the dimension six four-fermion operators of NRQED DM . Soft photon exchanges are encoded in pNRQED DM in the potential (2.8). The fermion-antifermion wavefunction in pNRQED DM , which at leading order in the multipole expansion is the solution of the Schrödinger equation with the potential (2.8), accounts by construction for the effect of multiple soft photon rescattering. We consider first the effect of multiple soft photon exchanges on fermion-antifermion scattering states. To compute the annihilation cross section we use the optical theorem as in the first equality of eq. (3.3), but now we compute the scattering amplitude from an initial to a final scattering state |p, 0 . The annihilation cross section reads where we have set from the beginning the center of mass momentum P = 0, δV ann (r) has been defined in eq. (2.10), and Im(d s ) and Im(d v ) have been given at order α 3 in eqs. (3.4) and (3.5). Differently from eq. (3.3), the quantity (σ ann v rel )(p) depends on the relative momentum p. The scattering wavefunction, Ψ p , has been expanded in partial waves, Ψ p , where is the orbital angular momentum quantum number; only the S-wave component of Ψ p contributes at leading order in the non-relativistic expansion, since dimension six fourfermion operators in NRQED DM and the ensuing annihilation potential (2.10) only project on S-waves. The factor S ann (ζ) is called Sommerfeld factor [101] and for an attractive Coulombic potential, like in our case, it reads (see e.g. [102,103]) Equation (3.8) shows manifestly the factorization of the different energy scales: the hard dynamics is contained in the NRQED DM matching coefficients Im(d s ) and Im(d v ), whereas the soft dynamics is contained in the wavefunction squared |Ψ p0 (0)| 2 . The Sommerfeld factor modifies the cross section from σ NR ann v rel , which is the cross section computed at the level of NRQED DM without Coulomb rescattering effects, see eq. (3.3), into (σ NR ann v rel )S ann (ζ). For v rel α, the annihilation cross section is significantly enhanced by the Sommerfeld factor and the prediction from (1.1) changes accordingly. 7 We consider now the effect of multiple soft photon exchanges on fermion-antifermion pairs below threshold. This leads to the formation of Coulombic bound states, whose masses up to order α 2 are where a 0 is the Bohr radius. To compute the annihilation width, we compute the amplitude in the first equality of eq. (3.3) from an initial to a final bound state |n, 0 . 8 At the order we are working, as we have seen in the case of scattering states, only S-wave bound states contribute. S-wave bound states may be in a spin-singlet paradarkonium state (para) or in a spin-triplet orthodarkonium state (ortho). At leading order paradarkonium and orthodarkonium are described by the same radial wavefunction R n0 (r). The paradarkonium and orthodarkonium annihilation widths read Γ n,para A derivation of the Sommerfeld enhancement for S-wave pair annihilation that includes the regime of very small momenta (velocities) for the unbound pair has been presented in [104][105][106]. The main result is a saturation of the Sommerfeld factor and a regular behaviour for v rel → 0. Diagrammatically this amounts to resum the annihilation term, namely the local four-fermion interactions shown in figure 3 (right). In this work, we assume to be away from such regime for unbound states. For bound states this resummation is never needed because the momentum of the particle in the pair is constrained to be of order M α. Finally, it is worth noticing that the thermally averaged cross section (3.2) removes the singularity at vanishing v rel . 8 The different normalization of the states leads automatically to a cross section for scattering states and a decay width for bound states.
Γ n,ortho For 1S states, R 10 (r) = 2 e −r/a 0 /a 3/2 0 , and we can write, keeping order α 3 terms in Im(d s ) and Im(d v ), (3.14) We conclude this section with two observations. First, we remark that the inclusion of the orthodarkonium decay width in the DM evolution equations requires for consistency the inclusion of the order α 3 corrections to Im(d s ) in the paradarkonium decay width and (σ NR ann v rel ) NLO in the computation of the annihilation cross section for scattering states, since all these terms are of the same order and originate from the same set of four-fermion matching coefficients in NRQED DM . This has not always been the case in some of the original literature [37].
The second observation is on the coupling constant. In the fundamental theory (2.1), dark photons couple only with the dark fermions X, hence the coupling runs with one flavor at scales greater than M and is frozen at the value α ≡ α(M ) at scales below M . The coupling would still run below the mass scale if the gauge bosons would be coupled to themselves like in non-abelian gauge theories (see appendix D) or to other light degrees of freedom. In NRQED DM and pNRQED DM , the dark matter fermions have been integrated out, hence the coupling appearing there is α, and it does not run. Because in NRQED DM and pNRQED DM the coupling remains the same at the mass scale and at the scale of the binding, it is parametrically consistent to include subleading corrections in α to the pair annihilation of the scattering states and to the decay widths of the bound states originating from higher-order corrections to the matching coefficients of the fourfermion operators, while at the same time ignoring higher-order corrections in α in the near threshold wavefunctions, which are suppressed by α 2 , as we have done. The situation is different in QCD-like non-abelian theories, where corrections to the wavefunction are typically as relevant as or more relevant than corrections to the four-fermion matching coefficients, on one hand because there, due to asymptotic freedom, the coupling at the mass scale is smaller than the coupling at the near threshold scales, and on the other hand because corrections to the wavefunction are suppressed by just α [107].

Bound-state formation and dissociation cross sections
The pNRQED DM Lagrangian (2.6) contains electric dipole terms that are responsible for dark fermion-antifermion pair transitions between scattering states, (XX) p , and bound states, (XX) n , as shown in figure 4. In section 4.1 we address transitions that lead to bound-state formations and in section 4.2 transitions that lead to bound-state dissociations. We complete the analysis in appendix B by providing results for the bound-state to bound-state transitions, which amount to excitation and de-excitation processes, and in appendix C by considering bremsstrahlung and thermal absorption processes in scatteringstate to scattering-state transitions.
Electric dipole transitions involve the emission or absorption of an ultrasoft dark photon. The medium breaks explicitly Lorentz invariance and we have to choose a reference frame. A convenient reference frame choice is the one where the medium is at rest. This choice made, the cross section depends on the center of mass momentum P . The center of mass momentum in the thermal average of the cross section times v rel scales like √ M T , which is the momentum scale in the Boltzmann distribution. Although √ M T is a soft scale, its effect on electric dipole transitions turns out to be subleading. The reason is that the center of mass momentum dependence is carried by the kinetic energy of the fermionantifermion pair, which changes by an amount is subleading with respect to T or M α 2 even in the momentum region P ∼ √ M T . We may therefore systematically expand in the center of mass momentum P . If we retain the leading order term, this amounts to set P = 0 in the cross section, 9 which is our choice in the following sections.

Bound-state formation
In this section, we compute at leading order in the non-relativistic and coupling expansion the cross section for the process where a bound state (XX) n is formed from a scattering state (XX) p via the emission of an ultrasoft dark photon. The transition is an electric dipole transition. In pNRQED DM , it corresponds to the inverse process of the one depicted in the first diagram of figure 4. The energy of the incoming pair is E p = 2M + p 2 /M + · · · = 2M + M v 2 rel /4 + . . . , and the energy of the resulting bound state is E n , defined in eq. (3.10). Since E p > E n , the process (4.1) is allowed in vacuum. Scattering-state to scattering-state transitions of the type (XX) p → γ + (XX) p are discussed in appendix C.
The cross section can be computed at leading order from the imaginary part of the oneloop self energy in pNRQED DM shown in the left diagram of figure 5. The photon could be a thermal photon and therefore the computation needs to be performed in the thermal field theory version of pNRQED DM . We use the real-time Schwinger-Keldysh formalism [109][110][111]. The real-time formalism necessarily leads to a doubling of the degrees of freedom called of type 1 and 2. The type 1 fields are the physical ones, i.e. those that appear in the initial and final states. Propagators are represented by 2 × 2 matrices, as they may 9 We come to the same conclusion if we choose the reference frame of the center of mass of the dark fermion-antifermion pair. In this case, the velocity of the medium is about T /M , which is much smaller than one, the velocity of light. Expanding in it, the thermal distribution of the dark photons reduces at leading order to the thermal distribution of the medium at rest [108]. involve fields of both types. The thermal dark photon propagator in Coulomb gauge reads The choice of the gauge is irrelevant for our computation of the electric dipole transition, since the electric field is gauge invariant. The thermal propagator of the fermion-antifermion field φ in pNRQED DM reads where p 0 is the energy of the dark fermion-antifermion pair. Note that the fermionantifermion pair behaves like a boson. In the heavy mass limit, recalling that H = 2M +. . . , n B (H) is exponentially suppressed as e −2M/T , so that it holds It has been remarked in [59] that, because the 12 component of a heavy-field propagator vanishes in the heavy-mass limit, the physical heavy fields do not propagate into type 2 fields. Therefore, to a large extent the type 2 fermion-antifermion fields decouple and may be ignored in the heavy-mass limit, which makes the real-time formalism convenient when dealing with heavy fields.
From eqs.
where p 0 is the energy of the incoming pair. We have regularized in d dimensions to get rid of scaleless integrals. In eq. (4.5) one can distinguish the in-vacuum and thermal contributions originating from the photon propagator, while the heavy fermion-antifermion pair has been taken to propagate in vacuum. In order to extract the imaginary part of (4.5), we can use the identity P stands for the principal value prescription. Following refs. [59,77], we obtain where ∆E ≡ p 0 − H and we have set d = 4 since the imaginary part of Σ 11 is finite. Finally, the imaginary part of the self energy is related to the cross section of the process (XX) p → γ + XX by the optical theorem We have projected on the scattering state r|p = Ψ p (r); p 0 = E p is the energy of the incoming fermion-antifermion pair. The difference ∆E is positive on any intermediate state that is a bound state, and in that case the process may happen both in vacuum and in the medium. If ∆E < 0 the process may only happen through the absorption of a dark photon from the medium. We can now proceed and compute the bound-state formation (bsf) cross section by projecting eq. (4.8) on intermediate bound states r|n = Ψ n (r), We have defined which holds at LO; σ n bsf is a shorthand notation for the cross section σ (XX)p→γ+(XX)n . An alternative way of computing the bound-state formation cross section is from the 21 component of the self energy, Once projected on intermediate bound states, this leads to the bound-state formation cross section given in eq. (4.9). If we would keep in our derivations the thermal distribution of the heavy fermionantifermion pair in the self-energy loop, which amounts to keep the temperature dependent part of (4.3), then the bound-state formation cross section would modify into Clearly, the thermal distribution of the heavy fermion-antifermion pair vanishes exponentially for T M ∼ E n /2, in which case the above expression reduces to (4.9). In order to compare with existing literature and to provide an application of the above expression, we consider the formation of the lowest-lying 1S bound state, whose wavefunction is r|1S = R 10 (r)/(4π). In this case, only scattering states in the partial wave = 1 contribute, whose wavefunction is r|p1 = Ψ p1 (r). The bound-state formation cross section reads with p = M v rel /2 and If we select the spin of the final state, the bound-state formation cross section for paradarkonium is σ 1S,para bsf = σ 1S bsf /4 and for orthodarkonium is σ 1S,ortho bsf = 3 σ 1S bsf /4. The evaluation of the dipole matrix element squared | n|r|p | 2 can be found in appendix A for a generic scattering-state to bound-state transition (see also refs. [26,42,112]).
By rewriting the result in (4.14) in terms of ζ = α/v rel , we recover in the zero temperature limit, i.e. by setting to zero the Bose-Einstein distribution in eq. (4.14), the expression derived in refs. [37,42], and also the abelian limit of the QCD expression derived in ref. [77]. We also agree with the finite temperature expression presented in ref. [57].
From the Maxwell-Boltzmann distribution we can define a most-probable velocity v p = 4T /M , which is the velocity that maximizes Both velocities depend on the ratio T /M and are much smaller than one, consistently with the nonrelativistic assumption. In figure 6, we give the particle-antiparticle annihilation cross section for scattering states and the bound-state formation cross section normalized to that of free pairs at LO, πα 2 /M 2 , for α = 0.1 and α = 0.05. Solid lines stand for cross sections computed at v rel = v p , dashed lines for cross sections computed at v rel = v m , and the shadowing for cross sections at v m < v rel < v p . Both cross sections increase as the temperature (and the typical particle velocity) decreases. For the same M/T , smaller coupling constants imply less enhancement and the bound-state formation takes longer to overcome the pair-annihilation cross section. We show the bound-state formation cross section for the 1S, 2S, 2P and 3S state; explicit expressions for the cross section of states above the ground state can be found in appendix A. Our results agree with ref. [112].
We can also define, similarly to what done for the annihilation cross section, cf. eq. (3.2), an average bound-state formation cross section σ bsf v rel . In figure 7, we give the thermally averaged annihilation cross section (solid black curve) and bound-state formation cross section (solid red curve). In order to show the impact of the thermal average on the different quantities, we also include as dashed lines the corresponding cross sections with a fixed velocity, namely v p . We notice that for the annihilation cross section the difference is very small, whereas it is larger for the bound-state formation cross section. Finally, the effect of the Bose-Einstein distribution for the emitted photon is way more important in the thermally averaged cross section than in the cross section at fixed-velocity (cf. orange curves). This is due to the fact that for

Bound-state dissociation
Dark photons from the plasma may trigger bound-state dissociation (bsd) through the reaction The process is also called photo-dissociation. As the plasma temperature decreases, fewer darkonium states can be ionized. In the literature, the dissociation width has been often derived from the in-vacuum bound-state formation cross section through the Milne relation [25,37], while the thermal character of the process has been recovered by thermally averaging over the Bose-Einstein distribution of the photon at a later stage. In the following, we compute the dissociation width of a bound state starting from pNRQED DM at finite temperature. In this way, the thermal photon distribution appears from the beginning in the self energy. Bound-state to bound-state transitions of the type γ + (XX) n → (XX) n and (XX) n → γ + (XX) n are discussed in appendix B. At one loop, we can obtain the photo-dissociation width of a bound state with quantum numbers n from the self-energy diagram shown in the right panel of figure 5. The optical theorem relates the bound-state width to the imaginary part of the self energy either in the form Γ n = −2 n|Im[Σ 11 (E n )]|n or Γ n = n|[−iΣ 21 (E n )]|n with E n being the energy of the incoming bound state. The expressions of Im Σ 11 and Σ 21 have been given in eq. (4.5) and eq. (4.11), respectively. Projecting on intermediate unbound fermionantifermion pairs of relative momentum p selects the part of Im The Bose-Einstein distribution n B (∆E p n ) vanishes in the T = 0 limit, which reflects the fact that the decay of a bound state into an unbound pair is kinematically forbidden in vacuum. The photo-dissociation width is therefore a purely thermal width, which reads where we have expanded the scattering states into partial waves of orbital angular momentum ; note that the state |p is an eigenstate of the orbital angular momentum, but no more of the momentum p. Keeping the thermal distribution of the heavy fermionantifermion pair in the self-energy diagram would modify the photo-dissociation width into Taking into account that the thermal distribution of the heavy fermion-antifermion pair vanishes exponentially for T M ∼ E p /2, the above expression reduces to the one in eq. (4.18).
For photo-dissociation of the lowest-lying darkonium, eq. (4.18) becomes 10 is the binding energy of the 1S state. Note that the photon needs to have a threshold momentum to trigger the breaking of the bound state. The result agrees in the abelian limit with the gluo-dissociation width of a color singlet quark-antiquark bound state in the temperature regime T ∼ M α 2 s [77,113]. The gluo-dissociation width of a heavy quarkonium in the static limit was obtained in ref. [59] and for a hydrogen atom in QED in ref. [74].
The thermal width can be also understood as the convolution of the in-vacuum ionization (ion) cross section of the bound-state for the process γ + (XX) n → (XX) p , which we denote σ n ion , with the thermal distribution of the incoming photon, where 2 is the number of final state photon polarizations and the relative velocity between the darkonium and the photon from the bath has been set equal to one. Comparing the above equation with eq. (4.20) and using the expression of the dipole matrix element given in appendix A, we obtain σ 1S ion (|k|) = α with w 1 (|k|) ≡ |k|/|E b 1 | − 1. The result in eq. (4.22) agrees with the ionization cross section given in ref. [37], where it was obtained through the Milne relation. 11 Here, we did not rely on an explicit use of the Milne relation, but on thermal field theory alone. Thermal field theory provides, by construction, the dissociation width with the correct temperature dependence. It is straightforward to compute the dissociation width of excited bound states. We consider 2S, 2P and 3S states, with binding energies E b 2 = −M α 2 /16 and E b 3 = −M α 2 /36, respectively, and compute the corresponding in-vacuum photo-dissociation cross sections in appendix A.1. In the left panel of figure 8, we plot the corresponding photo-dissociation widths, accordingly to eq. (4.21). At very small temperatures the thermal width for the 1S state vanishes faster than the one for the 2S state, at higher temperatures the 1S thermal width is larger than the 2S one.

Milne relation at finite temperature
The processes of bound-state dissociation and formation have a rather different behaviour depending on the temperature of the plasma. Typically, dissociation gets suppressed at low temperatures, whereas formation becomes larger. Nevertheless dissociation and formation cross section are related by the Milne relation. Let us consider, as an example, the 1S bound-state case. In terms of v rel = 2p/M , the ionization cross section (4.22) reads (4.23) The ratio of σ 1S ion and σ 1S bsf from eq. (4.14) is then .

(4.24)
In the T = 0 limit at fixed v rel one recovers the in-vacuum Milne relation as given in ref. [37], (4.25) If the initial state dark photons in the process (4.17) are thermally distributed then it could make sense to define a temperature-dependent bound-state dissociation cross section It is precisely this quantity that enters the width (4.21) for the ground state. Differently

Dark matter energy density
In this section, we collect numerical results for the dark matter energy density. The main motivation for a precise calculation of the DM energy density relies on the accurate determination of such observable, Ω DM h 2 = 0.1200±0.0012 [1], which features a 1% uncertainty. As previously mentioned, in this work we rely on Boltzmann equations as far as the rate equations are concerned. Upon including bound states in the network of Boltzmann equations, their solution can be cumbersome. We use an approximation, first introduced in ref. [114] and commonly adopted in the literature, that is based on an effective treatment of dark matter bound states. In a typical cosmological setting, the annihilation rate of bound states is pretty efficient, Γ ann H, so that they quickly adjust to their equilibrium number densities. Upon neglecting the (de-)excitations between bound states, one obtains a single Boltzmann equation that depends only on the density of scattering states n, and is governed by an effective cross section (see footnote 16). The latter comprises the effects of DM annihilation via unbound pairs and bound states, and bound-state formation cross sections and dissociation widths. Recent studies have further investigated the validity and generalization of this approach when transitions between different bound states are kept in the network of Boltzmann equations [26,115].
The Boltzmann equation that we use in the following numerical analyses reads The dissociation width and different transition widths contributing to the effective cross section σ eff entering eq. (5.1) are shown as a function of the temperature in figure 9 for the case of the 1S state. When bound-state to bound-state transitions are much smaller than Γ n ann and Γ n bsd , the thermally averaged effective cross section can be written as The first term stems from the annihilation of a pair in a scattering state, yielding eq. (1.1), whereas the second term encodes the reprocessing of an unbound pair into a bound state. The sum over n extends over all bound states: the 1S spin singlet paradarkonium state, the 1S spin triplet orthodarkonium states, and so on. The combination Γ n ann /(Γ n ann + Γ n bsd ) includes the information about the dissociation width of the bound state. It is small at large temperatures, but it is close to one at temperatures of the order of the binding energy or smaller, when bound state dissociation becomes less efficient than annihilation, i.e. Γ n bsd Γ n ann . At sufficiently small temperatures, the bound-state to bound-state transitions become as large as or larger than the dissociation widths and cannot be neglected with respect to them (see figures 9 and 16). In this situation, we replace eq. (5.2) with the expression that can be found in ref. [26], which takes into account transition rates. Finally, we remark that late in the evolution of the universe, i.e. at temperatures such that the DM particles are very diluted, DM particle interactions become negligible with respect to the universe expansion rate that dominates the evolution equation (5.1).
It is convenient to recast the Boltzmann equation (5.1) in terms of the yield Y ≡ n/s, with s being the entropy density. The present-day DM relic density is then Ω DM = M s 0 Y 0 /ρ crit,0 , where Y 0 , s 0 and ρ crit,0 denote respectively the present yield, entropy density and critical density. The values for s 0 and ρ crit,0 can be taken from e.g. [116], and one obtains Ω DM h 2 = (M/GeV) Y 0 /(3.645 × 10 −9 ), where h is the reduced Hubble constant. The temperature-dependent relativistic degrees of freedom entering the Hubble rate in eq. (5.1) are assumed to be those of the SM with the addition of the dark photon [37].
The main advantage of the EFT framework is to allow for a rigorous derivation and a systematic inclusion of corrections to the relevant observables. In the problem at hand, we are interested in the cross sections and widths entering the thermally averaged effective cross section. We refer to the scheme in figure 10 to pictorially illustrate the different corrections. As for the scattering states, which are integrated over all momenta, we can improve cross sections and widths by adding relativistic corrections. This amounts at including higher-dimensional operators in both NRQED DM and pNRQED DM . Furthermore, we can add radiative corrections to the matching coefficients of the effective field theories (see, for instance, section 3.1). As for the bound states, the situation is similar. Higherdimensional operators account for higher-order relativistic corrections, which at the same time may open new decay channels. Radiative corrections improve the matching coefficients. Moreover, we add a third dimension to the scheme that accounts for the number of bound states to be included in the analyses consistently with the other corrections. The inclusion of excited states can be seen as a further improvement towards an accurate description of the actual physical system, and hence a correction to the simplest possible situation, when only the ground state is considered. We provide below some examples that may help to assess the relative importance of the different corrections.
In the first case that we examine, we want to investigate the relative importance of NLO radiative corrections to the hard matching coefficients, which are corrections of higher order in α, with respect to relativistic corrections, which are corrections of higher order in the relative velocity of the dark fermion-antifermion pair. For illustration, we consider a particular subset of higher-order corrections in v rel , those stemming from four-fermion operators of dimension 8, which show up at order 1/M 4 . 13 Those operators can be read off from ref. [83], and the corresponding matching coefficients at O(α 2 ) from refs. [83,117]. The resulting annihilation cross section is where S =1 ann (ζ) = (1 + ζ 2 )S ann (ζ) is the Sommerfeld enhancement for a scattering state in a P-wave [102,103], and (σ NR ann v rel ) LO is given in eq. (3.6). We solve the Boltzmann equation (5.1) without bound-state effects, i.e. we neglect the second term in the right-hand side of eq. (5.2). The result is shown in figure 11. The black solid line is the ratio, R, between the energy density Ω DM h 2 as obtained with (σ NR ann v rel ) NLO S(ζ), where (σ NR ann v rel ) NLO is given in eq. (3.7), and (σ NR ann v rel ) LO S(ζ) as input for σ ann v rel , respectively. As elaborated in section 3.1, the NLO corrections to the matching coefficients make the cross section smaller, and hence a more abundant dark matter population is found for each value of α. Accordingly, we find R > 1, as shown in the plot. The dashed orange line stems for the ratio of the energy density as obtained from the cross section in eq. (5.3) and (σ NR ann v rel ) LO S(ζ). The trend here is different. The corrections to the cross section make it larger for each α, and accordingly we find a smaller DM energy density that results in R < 1. The P-wave contribution overcomes the negative correction of the velocity dependent S-wave correction, especially at large values of α. The result has a rather mild dependence on the specific value of the DM mass (set to M = 1 TeV in the plot). In the second case that we consider, we include bound-state effects, i.e. we reinstate the second term in the effective cross section in eq. (5.2). We work at leading order in the 1/M expansion. Therefore, we do not include dimension-8 operators neither in the annihilation cross section, nor in the annihilation widths of the bound states that can only be S-wave states, since P-wave states decay through dimension-8 operators at order 1/M 4 . We solve the Boltzmann equation with both contributions from scattering and bound states at order α 3 and α 2 . We restrict the sum over n just to 1S states, which means that at order α 2 we consider only paradarkonium and at order α 3 we include also the formation and decay of orthodarkonium. In the left panel of figure 12, we show the ratio of DM energy densities for the α 3 case over the α 2 case. As one may see, the ratio is R < 1 because the effective cross section (5.2) is larger whenever we add the spin-triplet bound state. The formation and decay of orthodarkonium induces an effect that is much larger than the few-per-cent decrease in the annihilation cross section due to O(α 3 ) terms. The overall effect is larger than the uncertainty on the relic density, which is about 1%, already for α = 0.05, and is larger than 13% for α = 0.1. We remark that the conditions required in order to use eq. (5.2) are satisfied in the orthodarkonium case. Indeed, the annihilation rate of the orthodarkonium (3.14) is much larger than the Hubble rate, and the transition rate from the spin-triplet to the spin-singlet 1S state satisfies Γ 1S,ortho→1S,para /Γ 1S,ortho In our last case, we compare the effect of adding excited states (n > 1) to the second term in the right-hand side of eq. (5.2) at order α 2 with having only σ ann v rel plus the effect of the 1S state at order α 3 . As for the excited states, we consider four additional states, namely the states 2S and 2P m with m = −1, 0, 1. Moreover, for n > 1 we also compare the no-transition case, i.e. the energy density obtained from the effective cross section given in eq. (5.2), with the energy density obtained from the effective cross section derived in ref. [26] that accounts for transitions between bound states. For these three approximations, we compute the DM energy densities normalized, like in the previous case, with respect to the DM energy computed from the effective cross section at order α 2 including the 1S bound state and scattering states. In the right panel of figure 12, we show the result for α = 0.1. The formation and decay of orthodarkonium for the 1S state (black-solid line) gives a larger contribution than including excited states with n = 2 in the spin-singlet configuration, when neglecting transitions among them (orange-dotted line). The reason for this is that the bound-state formation cross section for orthodarkonium is significantly larger than the bound-state formation cross sections for excited states. Upon including the bound-state to bound-state transitions (see appendix B for the corresponding widths), 2P m=0,±1 states may decay into the 1S ground state providing an additional contribution to the depletion of DM. This results in the orange-dashed line lying below the orange-dotted line, where transitions are neglected, and the black line.
In figure 13, we show the energy density contours for Ω DM h 2 = 0.1200 in the model parameter space (M, α). The orange line corresponds to the inclusion of paradarkonium 14 The transition between ortho-and paradarkonium is a magnetic transition, i.e. it is triggered by a magnetic-dipole operator. The width is proportional to α times the ratio of the hyperfine splitting to the third power over M 2 , which gives Γ 1S,ortho→1S,para ∼ M α 13 ; see, e.g. ref. [118]. and the annihilation cross section at order α 2 , whilst the brown line also includes effects at order α 3 . The brown line is systematically below the orange curve because the additional annihilation channel via orthodarkonium makes the effective cross section in eq. (5.2) larger, so that the same energy density is reproduced for smaller α. The black line includes the effects of order α 3 corrections as well as the excited states 2S and 2P m=0,±1 with transitions among them. This setting allows even smaller values of α to reproduce the observed energy density because of the additional channels for DM annihilations. In this section, we have solved the Boltzmann equation over a wide range of temperatures; as mentioned in the introduction, the freeze-out temperature can be estimated to be about M/25. The condition T M α 2 in eq. (2.2), or more conservatively T ≤ |E b 1 |, is satisfied for all times after freeze-out only if α 0.4, see figure 14. For smaller couplings, there exists a period of time after freeze-out when T is larger than M α 2 . In general, when T > M α 2 , thermal effects modify the dark fermion-antifermion potential and therefore the wavefunction. In the model (2.1) with L portal = 0, it has been shown that these modifications are subleading with respect to the Coulomb potential as long as T M [74]. 15 This is why we have ignored them here and computed widths and cross sections on Coulombic wavefunctions. Another important observation is that when T approaches M α, the multi- 15 In Coulomb gauge, the photon propagator is given by eqs. (4.2). Temporal photon propagators do not depend on the temperature at leading order and do not get screened by light fermion loops, which are absent, but get thermal corrections starting from two insertions of the F µν D 2 Fµν /M 2 operator appearing in the third line of eq. (2.3), whose contribution is T 4 /M 4 suppressed. Spatial photons, on the other hand, couple to the heavy DM fermions and antifermions through 1/M suppressed couplings, cf. the first two lines of eq. (2.3). Hence thermal corrections do not modify the Coulomb potential (2.9), but only 1/M suppressed potentials, whose contribution to the wavefunction is subleading as long as T M .
pole expansion for thermal dark photons breaks down (and one has to treat them at the level of NRQED DM ). However, for all the couplings considered in this section to compute the DM energy density, it holds that T M α after freeze-out. The numerical results for the DM energy density presented in this section are based on the single effective Boltzmann equation (5.1) or its modification that includes transitions among bound states. If instead we solve the coupled Boltzmann equations for the ground state and the unbound pairs numerically either by taking Im(d s ) at order α 2 and Im(d v ) = 0 or Im(d s ) and Im(d v ) at order α 3 , in which case an additional Boltzmann equation due to orthodarkonium has to be added, we get results that at most differ by 1% from the ones presented here, and so are within the uncertainty of the measured relic density. 16

Conclusions and outlook
In recent years, there has been a considerable effort in transferring and extending known techniques adopted in atomic and heavy-quarkonium physics to dark matter models, both at zero and finite temperature. Indeed, whenever the dark particles experience selfinteractions via vector (or scalar) mediators, and the DM abundance is fixed through the freeze-out mechanism, the dark particle dynamics shares some similarities with the one of heavy quarkonium in a quark gluon plasma. 16 The coupled Boltzmann equations for the 1S bound-state number densities n para 1S , n ortho 1S and the sum of the dark matter particle and antiparticle number densities n = nX + nX = 2nX are [5,37] (∂t + 3H)n = − where n ortho 1S,eq = 3n para 1S,eq , and n para 1S,eq = (M T /π) 3/2 e −E 1 /T and neq/2 = nX,eq = 2(M T /2π) 3/2 e −M/T (the factor of 2 for nX,eq is due to the spin polarization of the fermion). The coupled equations reduce to equation (5.1) with the effective cross section given by (5.2) through the following steps and approximations. i) From the detailed balance condition at equilibrium between bound-state formation and dissociation it follows that bsf v rel n 2 eq = Γ 1S bsd n ortho 1S,eq .
The detailed balance condition is equivalent to requiring (∂t + 3H)neq ≈ 0, i.e. that at equilibrium the number density changes only because of the universe expansion. ii) By using the detailed balance condition in the right-hand sides of the equations for n para In the context of the freeze-out mechanism, the relevant processes that determine the non-equilibrium dynamics of dark matter single particles and pairs, namely annihilations and decays, bound-state formation and dissociation, occur in a non-relativistic regime, i.e. when T M . Dark matter particles may therefore form non-relativistic bound states, whose energy scales together with the thermal scales characterize the dark matter dynamics. The computation of physical observables becomes cumbersome and not very transparent when all the scales are intertwined as in the full relativistic field theory. Instead, one can take advantage of the fact that the energy scales are hierarchically ordered by systematically replacing the relativistic field theory by a tower of simpler non-relativistic effective field theories. The advantage in using non-relativistic effective field theories is that one deals with the degrees of freedom that describe the physics of interest at the Lagrangian level, and that long-and short-range contributions are factorized for any observable. Moreover, it is straightforward to systematically improve the determination of physical quantities by including higher-order corrections in the coupling and in the relative velocity, guided by the power counting of the effective theory. In this paper, we have scrutinized the implementation of non-relativistic effective field theories to address the dynamics of heavy dark matter particles in the thermal environment provided by the early universe.
The main aim of this work is to provide a detailed step-by-step construction of nonrelativistic effective field theories for a simple dark matter model, a dark version of QED given in eq. (2.1) (its SU(N ) extension is in appendix D). We complement the existing literature with some theoretical and technical aspects that, in our perception, have not been sufficiently highlighted. The hierarchy of scales that we assume has been given in eq. (2.2). We show how the relevant reactions, which comprise pair annihilations of scattering states, decays of the bound states, their formation, dissociation and transition rates in a thermal environment, can be accounted for in the effective field theories NRQED DM and pNRQED DM , which are NRQED and pNRQED applied to dark fermions and photons, respectively. In particular, pNRQED DM is well suited to describe particle-antiparticle pairs in scattering states and in bound states, as well as transitions among them.
For pairs in scattering states, we derive the annihilation cross section in pNRQED DM , see eq. (3.8). It accounts for the effect of multiple soft photon exchanges that ultimately leads to the Sommerfeld enhancement. Our expression makes manifest the factorization between contributions from hard modes, which are encoded in the four-fermion operator matching coefficients of NRQED DM , and the soft dynamics captured by the wavefunction of the pair. For pairs in bound states, we write the paradarkonium and orthodarkonium decay widths in eqs. (3.13) and (3.14), respectively. We stress that if the orthodarkonium decay width, which is an effect at order α 3 , is included in the network of Boltzmann equations for the evolution of the dark matter pairs, then one has to include at the same order the relevant matching coefficient of NRQED DM also in the scattering state annihilation cross section, which has not always been the case in the literature. Indeed the two contributions originate from the same four-fermion operator, the only difference consists in the two-particle states onto which they are projected.
For transitions between heavy pairs mediated by ultrasoft/thermal photons, the leading effect, as made manifest by pNRQED DM , is due to an electric-dipole operator. In the main body of the paper, we focus on the transition between a scattering and a bound state, and its reverse process, and compute the bound-state formation cross section, eq. (4.9) for the general case and eq. (4.14) for the 1S case, and the bound-state dissociation width, eq. (4.18) for the general case and eq. (4.20) for the 1S case. We show how they can be derived from the self energy of the dark matter fermion pair in the thermal field theory version of pNRQED DM . The formation cross section and dissociation width turn out to depend on the temperature through the (unsuppressed) statistical distributions of the involved particles, which appear naturally in the thermal field theory expression of the propagators. In the appendices of the paper, we collect expressions for the scattering-state to bound-state and bound-state to bound-state transitions for generic states, both for the abelian model presented in the main body of the paper as well as for a non-abelian dark matter model. Among the differences between abelian and non-abelian models, we remark here that, while in an abelian model a small value of the coupling at the hard scale is enough to guarantee a weak-coupling treatment for threshold observables, in a non-abelian dark matter model a weak-coupling treatment requires that the coupling remains small also at the ultrasoft scale. If we take the ultrasoft scale to be of the order of the temperature, then the smallest scale considered in this work is T ≈ 10 −5 M . At one-loop running, α(2M ) needs to be quite small to keep α(T ) < 1. For instance, in the SU(3) non-abelian model considered in appendix D, the weak-coupling condition requires α(2M ) 0.04. This should caution about computing at weak coupling the bound-state formation cross section and dissociation width entering the network of Boltzmann equations for the extraction of the DM energy density, whenever dealing with QCD-like charged particles in coannihilation dark matter scenarios. Nevertheless, the EFT framework holds also if the ultrasoft scale is strongly coupled, which is a situation familiar to QCD [41,119].
Bound-state formation and dissociation rates are routinely used in the network of Boltzmann equations in order to extract the dark matter energy density. The derivation of these rates relies on the evaluation of matrix elements of the electric-dipole operator following from the multipole expansion in pNRQED DM . The multipole expansion holds for thermal photons as long as the typical distance of the fermion-antifermion pair is smaller than 1/T . At large temperatures, T ∼ M α, the multipole expansion breaks down. Estimating DM formation and dissociation in this situation requires computing thermal effects in NRQED DM and matching to a version of pNRQED DM that does not contain thermal photons as degrees of freedom, but encodes them in a temperature dependent potential. 17 A quantitative assessment of such scenarios may be needed in order to solve the Boltzmann equations over a range of couplings that include values making M α smaller than the freeze-out temperature, i.e. α 0.04. This range encompasses coupling values typical of the electroweak SM sector, which may be relevant for genuine WIMP dark matter particles, with or without coannihilating partners, e.g. supersymmetric model realizations [27,28,35,47,91,104] and the inert doublet model [120][121][122][123].
As a further outlook for future developments of this work, a richer dark sector, as well as a minimal non-abelian model, would modify, also significantly, the thermal dynamics, for instance by generating another thermal scale in the form of a Debye mass for the gauge boson. The Debye mass modifies the thermal propagator of the gauge boson. Moreover, light degrees of freedom coupled to the gauge boson induce new processes for bound-state formation and dissociation. A careful investigation of the Debye mass scale within the framework of non-relativistic effective field theories, in particular its role in the boundstate formation via gauge boson emission, is being actively pursued [58,124].
Finally, in the framework of the Boltzmann equations, thermal rates are just ingredients to be computed independently to fix the dynamics of the evolution equations. It would be desirable, however, to derive out-of-equilibrium evolution equations for dark matter particles from non-relativistic effective field theories along the lines worked out for similar systems in the SM and QCD eventually leading to Lindblad-like equations [87,[125][126][127][128][129][130]. It would then be the solution of these equations that provides the energy densities of the dark matter particles.

A Bound-state formation and dissociation: general expressions
In this appendix, we give the bound-state formation (bsf) cross section and the bound-state dissociation (bsd) width for a bound state with generic quantum numbers n, and m: n is the principal quantum number, the orbital momentum quantum number and m the magnetic quantum number. The expression for the bsf cross section of a generic bound state follows from eq. (4.9) and reads (σ n m bsf v rel )(p) = where the dark photon energy is given in eq. (4.10). The electric-dipole matrix element enters also in the reverse process, namely the dissociation of a bound pair into unbound DM particles by absorption of a dark photon from the plasma. The thermal break-up width of a bound state with arbitrary quantum numbers n, and m can be written as a convolution integral (see eq. (4.21)) with the ionization cross section given by .
Both the bsf and ionization cross sections depend on the same electric-dipole matrix element. We derive it in the following. We build our derivation of the electric-dipole matrix element on refs. [131,132]. 18 The necessary ingredients are the wavefunctions of the Coulombic scattering and bound states. The Coulomb wavefunction for a DM scattering state of positive energy p 2 /M reads where p ≡ |p| ≡ M v rel /2 is the momentum in terms of the relative velocity of the pair and 1 (a, b, c) is the confluent hypergeometric function. Choosing the coordinate system in such a way that p points to the z-direction, the scattering wavefunction can be expanded into partial waves Ψ p (r) = r|p as (2pr) (2 + 1)! (2 + 1)P (cos θ)e ipr where P (x) are Legendre polynomials. The Coulomb wavefunction for a DM bound state of quantum numbers n, and m, Bohr radius a 0 = 2/(M α) and negative binding energy (A.7) The matrix element n m|r|p is then where e x , e y and e z are the unit vectors in Cartesian coordinates. Due to the selection rule for electric-dipole transitions, only the partial waves = ± 1 give a non-vanishing contribution. The constants N , X and Y are whereas G 1 and G 2 are combinations of hypergeometric functions, 2 F 1 (a, b, c, d), (A.13) The squared matrix element (A.8) is | n m|r|p | 2 = |N | 2 2 ( + 1)(δ m,1 + δ m,−1 ) + 4( + 1) 2 δ m,0 |X| 2 |G 1 | 2 In particular, the squared matrix element for the 1S state appearing in the bound-state formation cross section (4.14) and in the dissociation width (4.20) is We notice that the relative momentum of the unbound pair, p = M v rel /2 , and the Bohr radius of the bound state, a 0 = 2/(M α), appear in the matrix element in the combination a 0 p = v rel /α.

A.1 Bound-state formation and dissociation for 2S, 2P and 3S states
In this section, we give the explicit expressions of the bound state formation cross section and dissociation width for 2S, 2P and 3S states as an application of the general formulas (A.1)-(A.3) and (A.14).

2S state
A 2S bound state has quantum numbers n = 2, = 0 and m = 0, and binding energy E b 2 = −M α 2 /16. The formation cross section reads In the second equality, we have inserted the expression of the electric-dipole matrix element squared | 2S|r|p1 | 2 that we compute from (A.14) and find to be Thermal break-up of 2S darkonium can be triggered by the absorption of a photon from the background plasma. The momentum threshold of the photon is smaller than for the 1S state because the 2S state is less tightly bound. The bound-state dissociation width for the 2S state following from eqs. (A.2) and (A.3) is .

2P states
We consider the formation of 2P bound states, whose quantum numbers are n = 2, = 1 and m = 0, ±1. Applying eq. (A.14) for this special case, only the two partial waves with = 0 and = 2 contribute; we obtain (A.21) In order to account for the degeneracy in the magnetic quantum number m and to compare explicitly the bsf of a 2P state with that one of a 2S state, we sum up the results (A. 20) and (A.21) Since at O(α 2 ) the binding energies depend only on the principal quantum number n, the states 2S and 2P m=0,±1 have the same energy. The thermal dissociation width, averaged over m, is . (A.25)

B Bound-state to bound-state transitions
In this section, we derive the transition rates between two DM bound states. In pNR-QED DM the relevant diagram is the middle diagram of figure 4, which shows the transition from a bound state with quantum numbers n to a bound state with quantum numbers n mediated by an electric-dipole vertex. Excitation refers to the process where a DM bound state absorbs an ultrasoft photon from the thermal bath. 19 Due to the energy gain, the state jumps to a higher energy configuration, like in the transition γ + (XX) 1S → (XX) 2P . De-excitation by emission of a photon is the reverse process, it involves an excited state that looses energy to a more tightly bound state. Since these processes are mediated by the electric-dipole operator, the angular momentum of the bound state must change by one unit, |∆ | = 1, whereas the spin is left unchanged. In pNRQED DM , both excitation and de-excitation processes can be described at leading order by cutting the self-energy diagram shown in figure 15. Here, the matter states inside and outside the loop are bound states, at variance with the diagrams in figure 5 that show transitions between bound states and scattering states. The computation of the transition width closely follows the one carried out in section 4. De-excitation transitions require the energy conservation condition E n = E n + E γ , with E n the energy of the incoming bound state, E n the energy of the outcoming bound state and E γ > 0 the energy of the on-shell dark photon. For Coulombic bound states, 19 The energy of the incoming photon is not large enough to break the bound state into an unbound DM pair, thus the excitation process can be distinguished from the thermal break-up process by requiring the energy Eγ + En to be negative. where the distribution of the bound states may be approximated by the Maxwell-Boltzmann distribution, as M T , so that n B (E n )/n B (E n ) ≈ e ∆E n n /T . Radiative transitions have been considered also in ref. [26], however, there the photon distribution has been set to 0 in the de-excitation width, which amounts at ignoring stimulated emission, and to e −|∆E n n |/T in the excitation width. Both approximations are valid only at very small temperatures, T |∆E n n |. In our numerical computations, we employ the transition widths in eqs. where we have now averaged over the magnetic quantum numbers of the decaying states. We plot the ratio of the excitation and de-excitation widths over the paradarkonium annihilation width in figure 16. One recognizes that in the region M α 2 T bound-state to bound-state widths are suppressed with respect to the leading-order annihilation width Γ 1S,para ann = M α 5 /2. Moreover, one can notice the qualitative different dependence on the temperature of the excitation and de-excitation processes. The former process becomes exponentially suppressed for smaller temperatures (solid lines), since the photons that allow for the excitation are progressively less abundant in the plasma, whereas the deexcitation has an in-vacuum contribution that dominates at small temperatures (dotted lines).

C Scattering-state to scattering-state transitions
Similarly to what happens for bound-state to bound-state transitions, scattering states can undergo processes of emission (or bremsstrahlung) and absorption of an ultrasoft/thermal photon. In pNRQED DM the relevant diagram is the most right diagram of figure 4, which shows the transition from a scattering state of relative momentum p to a scattering state of relative momentum p mediated by an electric-dipole vertex. Choosing the reference frame such that the relative momentum p = M v rel /2 of the incoming scattering state is oriented along the z-direction, the outcoming relative momentum p can be written as p = (cos (φ p ) sin (θ p ), sin (φ p ) sin (θ p ), cos (θ p ))M v rel /2. The matrix element p|r|p reads (C.7)

D Non-abelian SU(N ) model
In this section, we consider the DM fermion to be in the fundamental representation of a non-abelian SU(N ) gauge group, N ≥ 2. The corresponding Lagrangian reads where D µ = ∂ µ + igA a µ T a , T a are the group generators, A a µ the SU(N ) dark gauge fields and G a µν the corresponding field strength tensor. As in the abelian case discussed in the main body of the paper, we set L portal = 0.
We are interested in the energy regime M M α √ M T M α 2 T Λ, where Λ denotes the non perturbative scale where a weak-coupling expansion in α = g 2 /(4π) breaks down. 20 This choice of hierarchy for the energy scales allows us, like in the abelian case, to integrate out modes associated with the hard and soft scale by setting to zero the temperature characterizing the thermal distribution of the dark gauge fields. Moreover, computations at the hard, soft and ultrasoft energy scale may be done in weak-coupling perturbation theory. The ultimate EFT, made just of dark fermion-antifermion pairs of energy M v 2 rel and dark gauge bosons of energy and momentum M v 2 rel , which encompasses the scales M α 2 and T , or Λ, is a pNRQCD-like EFT [39,86] for DM fields and an SU(N ) gauge group. We dub it pNREFT DM . The Lagrangian at order r in the multipole expansion reads (D. 2) The field S = S 1 N ×N / √ N is an SU(N ) singlet field made of a dark fermion and antifermion and O = O a T a / √ T F is the corresponding SU(N ) adjoint field; they depend on time, the relative coordinate r and the center of mass coordinate R. The Hamiltonians, H s and H o , can be written as in (2.7) and (2.8) with the leading order potentials given respectively by 3) The Casimir of the fundamental representation is C F = (N 2 − 1)/(2N ) and the Dynkin index is T F = 1/2. The adjoint field O is a color octet field for N = 3, which is the QCD case.
The electric dipole terms in the second line of eq. (D.2), with V A (r) = V B (r) = 1 at leading order, allow for transitions between an unbound adjoint dark fermion-antifermion pair and a bound or unbound singlet pair. Such transitions are responsible for bound-state formation and dissociation. Moreover, they allow for scattering-state to scattering-state transitions among dark fermion-antifermion pairs in the SU(N ) adjoint representation. At variance with the abelian case, SU(N ) singlet transitions, either involving bound states or scattering states, cannot happen in this model because of the SU(N ) charge conservation.
Another crucial difference with the abelian model (2.1) is in the running of the coupling. For the non-abelian model (D.1), despite the fact that we do not consider additional light fermions or scalars, the interactions among gauge bosons are enough to make the renormalized coupling run also at scales below M , the first coefficient of the beta function being β 0 = 11N/3. Therefore, the coupling constant assumes different values when computing processes occurring at the hard, soft and ultrasoft scale. We typically renormalize α at the scale 2M in annihilation processes, at a soft scale of order M α in the wavefunction 21

D.1 Pair annihilation
In the non-abelian theory, there are two more four-fermion dimension six operators contributing to fermion-antifermion annihilation than in the abelian case due to the SU(N ) adjoint configurations. The contribution of the four-fermion dimension six operators to the part of the pNREFT DM Lagrangian relevant for annihilation is [134] δL ann (D.4) At order α 3 , the imaginary parts of the singlet matching coefficients read [83,135]: 22 and the imaginary parts of the SU(N ) adjoint matching coefficients are [83,135,136]: 23 The coupling α in the potentials (D.3) typically runs at a scale of order 1/r. 22 The matching coefficients Im(f [1] ( 3 S1)) reported in [83] and [135] agree for N = 3, however they do not for N = 3, as the factor 1/54 in [83] reads 1/(18N ) in [135]. Equation (D.6) is the one that can be found in [135]. 23 The matching coefficient Im(f [adj] ( 3 S1)) vanishes at order α 2 because the gauge bosons in (D.1) do not couple to light fermions. At order α 3 the expression in (D.8) is the one reported in ref. [136], while the expression in ref. [135] is three times larger.
with the free annihilation cross section at NLO given by (σ NR ann v rel ) The Sommerfeld factor for the singlet wavefunction at leading order is 24 The subscript NLO * reminds that our expression for the cross section is complete at NLO only in the hard part, the one encoded in the four-fermion matching coefficients, while it lacks NLO corrections in the wavefunction, i.e. in the Sommerfeld factor. Similarly, for the annihilation cross section of a pair in the adjoint representation we find (σ NR ann v rel ) [adj] with ζ defined as in (D.11). The total cross section reads (D. 15) In contrast to the abelian case, the NLO corrections increase the annihilation cross section both for SU(N ) singlet and adjoint pairs. However, the Sommerfeld factor has a different impact on the cross sections. For the attractive singlet channel, the Sommerfeld factor (D.11) is larger than one and increases the cross section, whereas for the adjoint repulsive channel, the Sommerfeld factor (D.14) is smaller than one and consequently decreases the cross section. At O(α 2 ) our results in eqs. (D.10) and (D.13) agree with those in ref. [137], once the cross sections in ref. [137] are expanded in v rel , and only annihilations into gauge bosons are taken into account (we have no light fermions in the model).
As mentioned above, the expressions (D.11) and (D.14) for the Sommerfeld factors do not include relativistic corrections to the wavefunctions. It is important to realize that in a non-abelian theory these corrections may eventually turn out to be more important than the NLO corrections to the free annihilation cross section, as the first ones are governed by the coupling at the soft scale, whereas the latter ones are governed by the coupling at the hard scale.
Bound states are sustained only by the SU(N )-singlet configuration. The annihilation width of a spin-and SU(N )-singlet pair in an nS wave reads at NLO in Im(f [1] ( 1 S 0 )) Γ nS,para 16) where at leading order we have distinguished between the coupling coming from the fourfermion matching coefficient, which is renormalized at the scale 2M , and the coupling coming from the wavefunction, which is renormalized at a scale µ s of the order of the soft scale. At next-to-leading order the natural scale of α appearing in (D.16) is 2M , since it originates from Im(f [1] ( 1 S 0 )). The decay width (D.16) does not include relativistic corrections to the wavefunction, which, as we have remarked, may be as much important as or more important than the included NLO corrections to the four-fermion matching coefficient. 25 At leading order, the decay width of a spin-triplet SU(N )-singlet pair in an nS wave reads Γ nS,ortho

D.2 Bound-state formation and dissociation
In this section, we start by considering the formation of bound states out of unbound adjoint states.
where we have indicated the natural scale of the coupling at NLO, in this way keeping distinguished NLO corrections coming from Im(f [1] ( 1 S0)) and the wavefunction. theory due to the SU(N ) charge conservation, i.e. (σ bsf v rel ) [1] (p) = 0. The total boundstate formation cross section, defined similarly to eq. (D.15), is then just (N 2 −1)/N 2 times (σ bsf v rel ) [adj] (p). It is the total bound-state formation cross section that can be found in the literature [26,48]. For the reverse process, namely bound-state dissociation, the thermal width is given by 4n 2 ) is the binding energy and the in-vacuum ionization cross section is . (D.20) In the ionization cross section we have averaged over the N 2 − 1 degrees of freedom of the incoming dark gauge field. The remaining computational effort is in the electric dipole matrix element. We give the result in full generality, as in the abelian case (A.8), with p chosen along the z-direction, (D.24) (D. 26) In the dipole matrix element the natural renormalization scale of the coupling α is µ s , which is of the order of the soft scale.
As an example, we provide some results for the formation of the ground state and its dissociation into an unbound SU(N )-adjoint pair. From the general expression (D.21), the matrix element squared for the ground state reads 27) which is in agreement with the result in ref. [77]. The total bsf cross section for the ground state is (D. 28) with p = M v rel /2 the incoming relative momentum of the adjoint pair, and Our result agrees with the outcome of ref. [48]. 26 Similarly to the abelian case, we thermally average the bsf cross section over the incoming momenta of the adjoint pair following the Maxwell-Boltzmann distribution (4.16). In figure 18, we plot the thermally averaged annihilation and bsf cross section, normalized by the thermally averaged free annihilation cross section at LO of an unbound pair, (D.30) as a function of M/T . We show the behaviour of the cross sections for a constant (dashed lines) and running coupling (solid lines) in an SU(2) and SU(3) dark matter model, respectively in the left and right panel of figure 18. In the case of a constant coupling, we take as a benchmark value α = 0.03. For the running coupling, the same value is taken at the hard annihilation scale, α(2M ) = 0.03, and run down to the smaller soft and ultrasoft scales, 27 where one finds α(2M ) < α(µ s ) < α(µ us ) for typical non-relativistic velocities. We take µ s = M v rel /2 and µ us = M v 2 rel /4. This gives more prominent near-threshold effects, as one may see in the annihilations of the singlet pair (orange dashed lines below the solid lines) and in the total annihilation cross section (black dashed lines below the solid lines). Conversely, for the adjoint SU(3) unbound pairs, the running coupling implies a stronger repulsion (brown solid curve below the dashed curve). We remark that the annihilation cross section of the pair in the adjoint representation vanishes for the SU(2) model at LO: (σ NR ann v rel ) [adj],SU(2) LO = 0. Comparing with the abelian case in figure 7, we see that the annihilation processes for the pair in a color singlet show a similar behaviour, namely a Sommerfeld enhancement that increases for smaller temperatures. However, in the SU(3) case, the contribution for the adjoint pair annihilations is suppressed by a Sommerfeld factor smaller than unity. This is due to the repulsive potential experienced by the adjoint pair, which becomes more relevant for smaller T , and thus lower velocities. The effect of a repulsive potential appears also clearly in the bsf process (red solid and dashed curves). At variance with the abelian case, the rising of the bsf cross section is saturated by the repulsive potential at small temperatures, and the bsf process becomes progressively less likely. This is signaled by the fact that for small v rel , whereas the right-hand side of eq. (4.14) goes like 1/v rel , the right-hand side of eq. (D.28) is exponentially suppressed. The running coupling enhances the bsf cross section with respect to a frozen coupling for the range of M/T shown in the plot.
The dissociation width of an SU(N )-singlet ground state into an unbound adjoint state is 26 In ref. [48], the authors distinguish between the coupling coming from the bound-state wavefunction, renormalized at a soft scale of order M α, and the coupling coming from the scattering-state wavefunction, renormalized at a soft scale of order M v rel . Although it may be useful to keep track of the two different couplings, the distinction between them is a higher-order effect that goes beyond the accuracy of the bsf and ionization cross section formulas that we use in this work. In those formulas, we only keep track of the ultrasoft scale in the electric dipole coupling. 27 Through this section, we evolve α at one loop.  Figure 18. Ratios of the thermally averaged cross sections in the non-abelian model SU(2) (left) and SU(3) (right). The orange lines denote the ratio when taking the thermal average of (D.9) but at LO in the free annihilation cross section, the brown lines when taking the thermal average of (D.12) but at LO in the free annihilation cross section, the black lines when taking the thermal average of (D.15) but at LO in the free annihilation cross section, and the red lines when taking the thermal average of (D.28). We take the ratio with respect to the thermally averaged free annihilation cross section (D. with w 1 (|k|) ≡ |k|/|E b 1 | − 1. The expression is in agreement with the result of ref. [77]. The left panel of figure 19 shows the behaviour of the dissociation width, normalized to the paradarkonium decay width at LO for the ground state, Γ 1S,para ann = C 4 F M α(µ s ) 3 α(2M ) 2 /4 with µ s = M α/2 and µ us = M α 2 /4, as a function of M/T for different non-abelian models, with and without the running of the coupling.

D.3 Energy density
At last we solve the coupled Boltzmann equations. Differently from the abelian case, the number density of scattering states comprises now the unbound pairs in both the SU(N ) singlet and adjoint configurations. We consider the simple case where we include the ground state only, see footnote 16. The coupled equations can be traded with a single effective Boltzmann equation in the form of eq. (5.1) for unbound pairs. Since bound-state to bound-state transitions are zero in the SU(N ) models under study, the thermally averaged effective cross section (5.2) does not get modified by them. We use for the annihilation cross section the thermal average of eq. (D.15), for the bsf cross section the thermal average  figure 19, we plot the thermally averaged effective cross section normalized to (D.30) for the SU(2) (black lines), SU(3) (orange lines) and SU(4) (brown lines) theory. Dashed lines correspond to a constant α = 0.03, whereas solid lines stand for a running coupling that has been evolved from the initial value α(2M ) = 0.03. The bell shape originates from the bound-state formation and decay contribution, σ 1S bsf v rel Γ 1S ann /(Γ 1S ann + Γ 1S bsd ), being dominant with respect to the annihilation term, σ ann v rel , for some N -dependent temperature regions. Note that the solid curves are systematically higher than the dashed curves, i.e. accounting for the running of the coupling increases the effective cross section.
In analogy with the numerical results presented for the abelian model, we first solve the effective Boltzmann equation without bound-state effects, while retaining the Sommerfeld factors in the annihilation cross section. In the left panel of figure 20, we show the ratios of the DM energy density Ω DM h 2 as obtained from the cross section in eq. (D.15) and the one extracted by replacing (σ NR ann v rel ) [1] NLO * with (σ NR ann v rel ) [1] LO , and (σ NR ann v rel ) [adj] NLO * with (σ NR ann v rel ) [adj] LO . Differently from the abelian case, NLO corrections increase the free annihilation cross sections (D.10) and (D.13) and, therefore, dark matter is more effectively depleted for increasing values of the coupling (compare with the different behaviour of the black solid line in figure 11 for the abelian model).
We then reinstate 1S bound-state effects in the effective cross section and numerically solve the effective Boltzmann equation the ratio of the DM energy density, Ω DM h 2 , obtained from annihilation cross sections and decay widths at NLO and the energy density obtained from annihilation cross sections and decay widths at LO. We observe that for the SU(3) and SU(4) models the NLO corrections are more significant, as there is an additional annihilation/decay channel due to the orthodarkonium, than for the SU(2) model, where Γ nS,ortho,SU(2) ann = 0 at order α 6 . As a cross check for the validity of the numerical results obtained by solving the single effective Boltzmann equation, we numerically evaluate also the coupled evolution equations. Considering the running of the coupling α and annihilations at NLO, we observe a 1% difference to the obtained energy density in the SU(2), SU(3) and SU(4) model, which is within the uncertainty of the measured relic density.
As a final remark, we comment on the value chosen for the coupling in this section when treating the running coupling case, namely α(2M ) = 0.03. This is a fairly smaller value than the ones that have been used in the abelian case. The reason for this choice is that larger values of the coupling would endanger the weak-coupling expansion in the lowest temperature region, T ≈ 10 −5 M , considered in this work. The temperature sets the magnitude of the ultrasoft scale. If we require that the coupling at the ultrasoft scale is weak, then we need to require that, in particular, α(10 −5 M ) < 1. This condition is fulfilled for the SU(4) theory at one-loop running if we choose α(2M ) 0.03, and by somewhat larger values of α(2M ) for smaller gauge groups: α(2M ) 0.04 for the SU(3) theory and α(2M ) 0.06 for the SU(2) theory. For these values of the couplings, the absolute value of the binding energy, |E b 1 |, is smaller than the freeze-out temperature M/25. Differently from the abelian case (see figure 14), it is not possible to accommodate weak coupling with having |E b 1 | of the order of the freeze-out temperature, i.e. in the weakly-coupled non-abelian theory, the evolution after freeze-out proceeds necessarily for some time at a temperature that is larger than the binding energy of the dark fermion-antifermion bound states. In dependence of the coupling and the theory, also the soft scale may turn out to be smaller than or of the same order as the freeze-out temperature, which signals the break down of the multipole expansion for thermal gauge fields at freeze out. This happens for the SU(2), SU(3) and SU(4) theories considered here at α(2M ) 0.03. We have not accounted for this in figure 20, where we assume that the multipole expansion holds also at freeze-out.