Radiative bound-state formation in unbroken perturbative non-Abelian theories and implications for dark matter

We compute the cross-sections for the radiative capture of non-relativistic particles into bound states, in unbroken perturbative non-Abelian theories. We find that the formation of bound states via emission of a gauge boson can be significant for a variety of dark matter models that feature non-Abelian long-range interactions, including multi-TeV scale WIMPs, dark matter co-annihilating with coloured partners and hidden-sector models. Our results disagree with previous computations, on the relative sign of the Abelian and non-Abelian contributions. In particular, in the case of capture of a particle-antiparticle pair into its tightest bound state, we find that these contributions add up, rather than partially canceling each other. We apply our results to dark matter co-annihilating with particles transforming in the (anti)fundamental of SU(3)c, as is the case in degenerate stop-neutralino scenarios in the MSSM. We show that the radiative formation and decay of particle-antiparticle bound states can deplete the dark matter density by (40-240)%, for dark matter heavier than 500 GeV. This implies a larger mass difference between the co-annihilating particles, and allows for the dark matter to be as heavy as 3.3 TeV.


Introduction
In theories with light force mediators, the formation of bound states may have significant implications. Early work pointed out aspects of the effect of bound states on the expected dark matter (DM) detection signatures [1][2][3]. More recently, a variety of phenomenological implications have been identified or explored in greater depth. It has been shown that the formation and decay of unstable bound states can affect the density of thermal-relic DM [4], and is essential in determining the unitarity limit on the DM mass [5]. Moreover, boundstate formation (BSF) processes enhance the expected annihilation signals of symmetric [1,[6][7][8][9][10][11] and asymmetric DM [5,12]. Since BSF cross-sections are typically dominated by different partial waves than the direct annihilation processes, they exhibit different velocity dependence and resonance structure, and can give rise to complementary signatures [7,9]. In the case of asymmetric DM [13], the formation of stable bound states may produce novel direct [14,15] and indirect detection signals [16][17][18][19], as well as affect the DM selfscattering inside haloes [20]. In confining theories, bound states may set the DM mass scale [21,22], and relate it to that of ordinary matter [23]. Finally, DM bound states may have implications for collider experiments [3,[24][25][26]. 1 Here, we consider BSF processes in non-Abelian theories. Non-Abelian interactions are particularly important in scenarios where DM is coannihilating with coloured partners, as well as in models where DM consists of TeV-scale Weakly interacting massive particles (WIMPs). Scenarios that feature DM coannihilation with coloured particles are encountered within the minimal supersymmetric standard model (MSSM) [32][33][34][35][36][37][38][39], and are in part motivated by the measurement of the Higgs mass [40,41]. These scenarios are being probed in high-precision collider experiments; the accurate prediction of the DM density within their parameter space is therefore necessary in order for the cosmological constraints to meaningfully complement those from colliders. On the other hand, TeV-scale WIMP models are being probed mostly via indirect searches. Their annihilation signals exhibit sharp resonances that depend on the DM mass [42][43][44][45][46]. It follows that the precise value of the DM mass, which, under minimal assumptions, is predicted by the observed DM density, determines the viability of these models. Clearly, in both cases, an accurate computation of the DM freeze-out is necessary. It is then essential that the depletion of the DM abundance via BSF in the early universe -which, as we show, can be a leading order effect -is properly accounted for.
In this paper, we compute the cross-sections for the radiative capture into bound states of non-relativistic particles transforming under a non-Abelian gauge group. We consider, in particular, unbroken non-Abelian theories in the regime where the gauge coupling is perturbative. We do not specify the gauge group or the representation of the interacting particles, such that our results are applicable in a variety of models. Rather than an effective field theory approach [47] that has been employed in other studies [6,8,[48][49][50][51][52][53][54][55][56], we use the method described in ref. [57], where the non-relativistic approximation is carried out directly on the relativistic amplitude.
We apply our results to a simplified model where DM is co-annihilating with scalar particles transforming in the fundamental of SU p3q c , and show that BSF has a very important effect on the DM relic density. In the MSSM incarnations of this scenario, the coloured particles typically possess also a sizeable coupling to the Higgs boson, which has been shown to mediate a sufficiently long-range interaction that enhances the annihilation cross-section [58]. While we do not consider it here, the attractive force mediated by the Higgs is expected to make the BSF effect more pronounced [59].
The formation of bound states in non-Abelian theories and their implications for DM have been considered in previous works. Reference [8] considered DM transforming un-der the adjoint of SU p2q L (Wino-like DM), and performed computations in the broken electroweak phase, with the purpose of estimating the DM indirect detection signals. References [48][49][50][51] considered BSF via non-radiative scattering processes that can take place at a high rate in a thermal bath, and computed the effect of these processes on the DM density. Various rearrangement processes in Abelian and non-Abelian theories have been discussed in Ref. [60]. Finally, refs. [55,56] considered radiative BSF in unbroken non-Abelian theories. Our computations disagree with those of [55,56] on the relative sign of the Abelian and non-Abelian diagrams contributing to these processes, but are in agreement with the dissociation rate of heavy quarkonium via gluon absorption computed in earlier work [61]. For the capture of a particle-antiparticle pair into its tightest bound state -which is typically the most significant capture process -our results imply that the leading order contributions add up, rather than partially canceling each other. This has very significant phenomenological implications, as we showcase in section 3.
The paper is organised as follows. In section 2, we compute the radiative BSF crosssections. In section 3, we calculate the DM freeze-out including BSF in the scenario of DM co-annihilation with coloured partners. We conclude in section 4. Several important calculations are included in the appendices. In particular, in appendix B, we compute the overlap integrals that enter the BSF cross-sections. In appendix C, we adopt an effective field theory standpoint, we derive the non-relativistic Hamiltonian of the interactions that determine the formation of bound states, and point out the disagreement with previous works [55,56]. This provides an independent check of the validity of our results.

Radiative bound-state formation in non-Abelian theories
We consider two complex scalar fields X 1 and X 2 , transforming in the representations R 1 and R 2 of a non-Abelian gauge group G. The Lagrangian is where D µ " B µ`i g s G a µ T a is the covariant derivative, with G a µ being the gluon fields (we shall denote the corresponding particles with g, as usual) and T a " T a 1 or T a 2 being the generators of X 1 and X 2 respectively. The fine structure constant is In the following, we shall compute the radiative BSF processes While we express our results in terms of the capture processes (2.3), it is straightforward to generalise them to transitions between scattering states (bremsstrahlung) and bound-state excitation or de-excitation processes, by simply substituting the appropriate wavefunctions for the initial and final states [57]. Since we shall only compute the leading order terms to these transition processes, our results apply also to fermionic X 1 and/or X 2 . Spin-orbit coupling arises only in higher orders in the non-relativistic regime.
We begin in section 2.1, with a summary of various formulae we use in our computations. In section 2.2, we give the interaction potential, and discuss the running of α s . In section 2.3, we compute the amplitude for the transitions (2.3), for general representations and masses of X 1 and X 2 . In section 2.4 we apply the result to particles in conjugate representations, but with arbitrary masses. For transitions involving a colour-singlet scattering or bound state, we compute explicitly the projected amplitude, and in section 2.5, we calculate the corresponding cross-sections.

Definitions and useful formulae
For easy reference, we first summarise various formualae that will be used in the following sections.
We define the total and the reduced mass of the interacting particles, 4) and the dimensionless factors For the momenta p 1 and p 2 of X 1 and X 2 , we shall often use the following momentum transformation, which allows to separate the center-of-momentum (CM) from the relative motion [57,62], Let S 1 pp 1 q and S 2 pp 2 q be the propagators of the X 1 and X 2 , For convenience, we also define Spp; P q " S 1 pη 1 P`pq S 2 pη 2 P´pq , (2.8) To leading order in the non-relativistic regime [57, appendix C], and [57, appendix E] ż dq 0 2π dp 0 2π Spq; Kq S 1 pη 1 P`pq S 0 pq; Kq S 0 pp; P q p2πq 4 δ 4 pq´p´η 2 P g q » 2m 2 p2πq 3 δ 3 pq´p´η 2 P g q, Spq; Kq S 2 pη 1 P´pq S 0 pq; Kq S 0 pp; P q p2πq 4 δ 4 pq´p`η 1 P g q » 2m 1 p2πq 3 δ 3 pq´p`η 1 P g q. (2.11b) As we shall see in section 2.3, we use eqs. (2.10) and (2.11) to integrate out the virtuality of X 1 , X 2 in the radiative part of the BSF diagrams.

Potential and the running of the coupling
The interaction between X 1 and X 2 can be decomposed into irreducible representations, For eachR, the interaction is described in the non-relativistic regime by a static Coulomb potential [63] V prq "´α g {r , (2.13) where the coupling α g is related to α s according to (2.14) Here, C 2 pRq is the quadratic Casimir invariant of the representation R. The Coulomb potential (2.13) distorts the scattering-state wavefunctions and, if attractive, gives rise to bound states. The scattering-state and bound-state wavefunctions are reviewed in appendix A.
In general, the coupling α s depends on the momentum transfer Q, 15) which is different in the various interaction vertices that appear in the transitions we consider. In table 1, we list the various vertices, specify the symbols we use, and give the average Q in each case.

Amplitude for radiative transitions
Radiative transitions are represented by the diagram of fig. 1a, which can be separated into the wavefunctions of the asymptotic states and the radiative vertex. The wavefunctions resum the two-particle interactions at infinity. The long-range X 1´X2 interaction arises from the one-gluon exchange kernel, which gives rise to the static potential of eq. (2.13) in the non-relativistic regime. The low momentum transfer (" µv rel for the scattering states and " µα s for the bound states) via the exchanged gluons is responsible for the appearance of non-perturbative phenomena, the Sommerfeld effect [64,65] and the mere existence of bound states. The radiative vertex is computed perturbatively, with the leading order contributions shown in fig. 1b. We discuss them further below.
In the instantaneous approximation, the amplitude for the radiative capture into a bound state is [57] rM ν kÑtn mu s a ii 1 ,jj 1 " where the Latin indices i, i 1 , j, j 1 and a denote the colour of the initial and final state particles, as shown in fig. 1b, andφ k pqq andψ n m ppq are the scattering-state and boundstate wavefunctions in momentum space that obey the Schrödinger equation. Here, q (´q) and p (´p) are the 3-momenta of X 1 pX 2 q in the CM frame, in the scattering state and in the bound state, respectively. The scattering state wavefunction is characterised by the continuous quantum number k, which specifies the expectation value of q. In a central potential, such as eq. (2.13), the bound state wavefunction is characterised by the standard discrete principal and angular-momentum quantum numbers tn, , mu, which specify the expectation value of p. The wavefunctions in a Coulomb potential are reviewed in appendix A. 2 M ν trans pq, pq is the perturbative transition amplitude with the virtuality of the interacting particles integrated out, as follows [57, sec. 3.3] rM ν trans pq, pqs a ii 1 ,jj 1 " 1 S 0 pq; Kq S 0 pp; P q ż dq 0 2π dp 0 2π rC ν pq, p; K, P qs a ii 1 ,jj 1 .
(2.17) 2 We note that the integrand in eq. (2.16) admits corrections of higher order in q and p that arise from the relativistic normalisation of states. Here we are interested only in the leading order terms, and we shall neglect these corrections. However, these corrections become important when there is a cancellation between the lowest order contributions to Mtrans [9,57].
Here, rC ν pq, p; K, P qs a ii 1 ,jj 1 is the sum of all connected diagrams contributing to the process X 1,i pη 1 K`qq`X 2,j pη 2 K´qq Ñ X 1,i 1 pη 1 P`pq`X 2,j 1 pη 2 P´pq`g a pP g q, (2.18) with the incoming and outgoing X 1 , X 2 being off-shell and only the emitted gluon g being on-shell and amputated. We emphasise that the X 1 , X 2 incoming and outgoing legs should not be amputated in the computation of C ν ; the proper amputation is done by the prefactor in eq. (2.17). (We recall that S 0 is defined in eq. (2.9).) Note that the connected diagrams contributing to C ν may include not-fully-connected diagrams that are non-zero due to the off-shellness of the legs, such as the diagrams in which the radiation is emitted from one of the legs (cf. fig. 1b). 3 In eq. (2.18), the momenta of the particles are indicated inside the parentheses. While the 3-momenta q and p follow the probability distributions given by the wavefunctions φ k pqq andψ n m ppq that appear in eq. (2.16), q 0 and p 0 are determined by the poles of C ν , upon the integration denoted in eq. (2.17). The total 4-momenta of the scattering state, the bound state and the radiated gluon, K, P and P g respectively, essentially contain all the (discrete and continuous) quantum numbers that fully specify the system. In the non-relativistic regime, they can be expressed as where E k " k 2 {p2µq " µv 2 rel {2 is the kinetic energy of the scattering state in the CM frame, with v rel being the relative velocity of the interacting particles, and E n ă 0 is the binding energy of the bound state. Note that M n " M`E n is the mass of the bound state. For a Coulomb potential, E n "´κ 2 {p2n 2 µq, with κ " µα B s (cf. appendix A). Energy-momentum conservation, K " P`P g , implies (2.20) The leading order contributions to rC ν s a ii 1 ,jj 1 are shown in fig. 1b. We compute them next using the Feynman rules from [66].  Emission from the mediator ipC ν med q a ii 1 ,jj 1 " p´g BSF s f abc q tg ρµ rpη 1 K`q´η 1 P´pq´pη 2 K´q´η 2 P`pqs ν g νρ r´P g´p η 1 K`q´η 1 P´pqs µ`gµν rpη 2 K´q´η 2 P`pq`P g s ρ u , (2.21a) We are interested only in the spatial components of C ν , ν " 1, 2, 3, The wavefunctions in eq. (2.16) impose |q| " k " µv rel and for a Coulomb potential |p| " κ{n " µα B s {n. Noting the hierarchy of scales |q|, |p|, k, κ ! M, µ, and using eqs. (2.19) and (2.20), we find from eqs. (2.21) that to leading order in v rel and α s , in the CM frame, A few comments on C med are in order.
• The leading order contribution arises from the contraction of the momenta of the two vertices on the scalar legs.
• The factor g 2 s arises from the same vertices, where the momentum transfer is " |p´q|, with |p| " κ " µα B g and |q| " k " µv rel . This implies that g s should be evaluated at the scale Q « ?
• In eq. (2.23a), we have neglected the energy transfer along the gluon propagators (as is also done in the one-boson exchange diagrams that are resummed into the nonrelativistic potential). Upon the integration indicated in eq. (2.17), the poles of the scalar propagators, Spq; Kq Spp; P q, set q 0 " q 2 {µ " k 2 {µ and p 0 " p 2 {µ " κ 2 {µ. In contrast, the poles of the gluon propagators set q 0´p0 " |q´p| " a k 2`κ 2 . Since q 0 and p 0 indicate the off-shellness of the two scalar particles, the gluon poles -on which the off-shellness of the scalars is greater for α s , v rel ă 1 -yield subdominant contributions to the transition amplitude.
• Naively, C med may be expected to be of higher order in α s than C 1 and C 2 . However, the scaling of q 0 , p 0 , |q| and |p| with α s and v rel described above implies all three diagrams are of the same order. This will become apparent in the following.

Colour decomposition for conjugate representations
We now focus on particles transforming under conjugate representations R 1 " R and We shall not assume though that the masses of the interacting particles are equal, so that our results are more widely applicable. For the overlap integral J kÑtn mu pbq, with |b|9|P g |, the dominant contribution is independent of b [9,57]. In the following, we denote J k,tn mu " J k,tn mu pb " 0q, and the amplitude (2.25) becomes The tensor product of two conjugate representations contains always a singlet and an adjoint, and possibly other states, It is clear from eq. (2.14) that, among all irreducible representations, the singlet configuration (C 2 p1q " 0) exhibits the most attractive potential, and thus accommodates the tightest bound state. Equation (2.29) implies that at least the following capture processes are allowed by the group algebra, provided that the potential in the final state is attractive, such that the bound state exists: pX`X : q radjs Ñ BpXX : q r1s`gradjs , (2.30a) pX`X : q r1s Ñ BpXX : q radjs`gradjs , (2.30b) pX`X : q radjs Ñ BpXX : q radjs`gradjs .
Depending on the group G and the representation R, more transitions may be possible. The amplitudes for the various transitions may be computed from eq. (2.28) by projecting onto the appropriate colour representations, using the Clebsch-Gordan coefficients. Below we compute explicitly the amplitudes for the transitions (2.30a) and (2.30b) only. In the following, d R , CpRq and C 2 pRq stand for the dimension and the Casimir invariants of the representation R, and d G , CpGq and C 2 pGq are the corresponding quantities for the group G. Evidently, the wavefunctions and therefore J k,tn mu and Y k,tn mu depend on the colour representations of the scattering and bound states, R S and R B . Whenever appropriate, we shall use the notation J rR S ,R B s k,tn mu and Y rR S ,R B s k,tn mu .

Adjoint scattering states to singlet bound states
The radiative capture into colour-singlet bound states can occur only from adjoint scattering states, pX`X : q radjs Ñ BpXX : q r1s`gradjs . (2.31) It thus suffices to project only the final X´X : state onto the singlet configuration; upon summing the squared amplitude over colours, the group algebra will project the initial state onto the adjoint. The amplitude for the process (2.31) iś where in the last step we used f abc T b T c " pi{2q C 2 pGq T a . The amplitude squared, coloursummed and averaged over the colour of the initial particles is

Singlet scattering states to adjoint bound states
With the emission of a gluon, a colour-singlet scattering state can turn only into an adjoint state, pX`X : q r1s Ñ BpXX : q radjs`gradjs . (2.34) Similarly to the above, the amplitude for the process (2.34) is deduced from eq. (2.28) by projecting the initial X´X : state onto the singlet configuratioń M r1sÑradjs kÑtn mu¯a i 1 ,j 1 "

Remaining transitions
Depending on R, other transitions may be possible. In this case, appropriate projections of the amplitude (2.28) have to be computed. It is sometimes possible to obtain the amplitudesquared for a transition of interest from the total amplitude-squared by subtracting the contributions of other known transitions. For this reason, we provide here the total colouraveraged squared amplitude, Note that J k,tn mu and Y k,tn mu depend on the representations of the scattering and the bound states. Thus, if more than one transitions are possible, in which either the scattering and/or the bound states belong to different representations, then their contributions to eq. (2.37) need to be separated, and J k,tn mu , Y k,tn mu should be evaluated using the wavefunctions of the corresponding representations.

Cross-sections for capture into the ground state
The differential BSF cross-section for capture into the ground state tn mu " t100u is given by where energy-momentum conservation implies (cf. eq. (2.20)) The leading-order contributions to the amplitude are J k,t100u 9k and Y k,t100u 9k (cf. refs. [9,57] and appendix B). Thus where θ is the angle between k and P g , and |M kÑt100u | 2 is independent of θ. Thus, For convenience, in the following we shall use the parameters [cf. eqs. (A.6)] From the amplitudes of eqs. (2.33), (2.36) and (2.37), and the expressions (B.4) for the overlap integrals J and Y, we find that the colour-averaged BSF cross-sections are where f c is a numerical factor that depends on the transition, radjs Ñ r1s, r1s Ñ radjs, , rest, A few remarks are in order: • Clearly, in the radjs Ñ r1s transition, the contributions from all three diagrams of fig. 1b add up. No (partial) cancellation occurs, for any group or representation, contrary to what was found in [55,56]. We note that our result reproduces the dissociation rate via gluon absorption of the colour-singlet bound state of a particleantiparticle pair transforming in the (anti)fundamental of SU pN q, that was computed in ref. [61, eq. (19)]. 5 The potential in the singlet state is always attractive and gives rise to the tightest bound state. Thus, our results quite generally suggest that BSF can be very significant for phenomenology. 6 Moreover, the radiative transitions contribute to the self-energy of the initial state. From the optical theorem and eq. (2.42a), it follows that the forward scattering amplitude (or equivalently, the index of refraction) of the adjoint state is enhanced by the non-Abelian contribution, as is reasonable to expect.
• Since α s runs only logarithmically, to a good approximation we may set α NA s » α B s , at least in the parameter space where BSF is significant (cf. footnote 4). Then, for the radjs Ñ r1s and r1s Ñ radjs transitions, the f c factors simplify.
• We recall that the couplings α BSF s , α NA s , α B s , α B g , α S g , and thus ζ S and ζ B , depend on the colour representations of the initial and final states (cf. table 1), and are different for every transition.
• As noted in section 2.4.3, if the group algebra allows for more than one transitions in the category "rest", then their contributions to f c have to be disentangled, in order for S BSF to be computed. Note that a transition pXX : qR Ñ pXX : qR 1`g allowed by the group algebra contributes to f c even ifR 1 has a repulsive potential and cannot accommodate a bound state.
The function S BSF pζ S , ζ B q encapsulates all the velocity dependence of σ BSF v rel . The first two factors in eq. (2.42b) arise solely from the scattering-state wavefunction and coincide with the Coulomb Sommerfeld enhancement of p-wave annihilation processes S 1 pζ S q " r2πζ S {p1´e´2 πζ S qsp1`ζ 2 S q. The factors inside the square brackets in eq. (2.42b) arise from the convolution of the scattering-state and bound-state wavefunctions with the radiative vertices.
Let us now discuss the asymptotic behaviour of S BSF in various cases.
• At large velocities, |ζ S |, ζ B ! 1, BSF is very suppressed, • For an attractive interaction in the scattering state (ζ S ą 0), and at low enough v rel such that ζ S Á 1 and ζ B Á 1, is constant, S BSF exhibits the characteristic scaling S BSF 91{v rel . 7 We observe that S BSF becomes maximal for α S g {α B g " 0.5. That is, the transition probability decreases for transitions between states governed by very different potentials (α S g {α B g " 1 or ! 1).
• For a repulsive interaction in the scattering state (ζ S ă 0), and at low enough v rel such that ζ S À´1 and ζ B Á 1, At low v rel , S BSF becomes exponential suppressed. However, the exponential suppression sets at ζ B ą 1, by when BSF may already have an important effect on the DM density. It is interesting that the exponential suppression is more severe for tighter bound states (larger ζ B ), i.e. two particles that repel each other are less likely to be captured into a very deep bound state. This is consistent with the behaviour exhibited by eq. (2.43b).
3 Dark matter co-annihilating with coloured partners

Simplified model and Boltzmann equation
We assume that DM is a Majorana fermion χ of mass m χ , that co-annihilates with a complex scalar triplet under SU p3q c , denoted by X. The gauge interactions of X are specified by the Lagrangian where D µ,ij " δ ij B µ`i g s G a µ T a ij is the covariant derivative, with G a µ being the gluon fields and T a are the generators. χ and X are the lightest and next-to-lightest particles that are odd under a Z 2 symmetry which prevents χ from decaying. We also assume that the interactions between χ and X -which we shall leave unspecified -keep them in chemical equilibrium throughout the freeze-out of their annihilation processes into other species.
As long as the relative mass splitting between DM and its coannihilating partner, is small, δ ! 1, the DM density is determined by the χ´χ, χ´X, χ´X : and X´X : (co-)annhilation processes. It can be tracked by considering the sum of densities of all co-annihilating species,Ỹ where Y j " n j {s, with n j being the number density of the species j and s " p2π 2 {45q g˚S T 3 being the entropy density of the universe. Using the time parameter the evolution ofỸ is governed by the Boltzmann equation [67] dỸ dx "´c with g χ " 2 and g X " 3 being the χ and X degrees of freedom.
The effective cross-section xσ eff v rel y in eq. (3.5) includes all annihilation and co-annihilation processes weighted by the densities of the participating species. We shall assume that the dominant contribution arises from the processes that annihilate XX : , with total crosssection σ XX : , such that Both the direct annihilation and the BSF processes contribute to σ XX : , as we discuss in the following.
In this work, we shall neglect thermal effects. The thermal bath may affect the DM freeze-out in a variety of ways, including, on one hand, screening of the long-range interactions and, on the other hand, frequent (non-radiative) scattering processes that precipitate DM depletion via BSF [49]. In the context of DM coannihilation with coloured partners, the latter have been considered in Ref. [51]. The inclusion of thermal corrections for the radiative BSF processes considered here requires a comprehensive study that we leave for future work.

Colour states and the running of the coupling
The X´X : colour interaction may be decomposed as 3 b3 " 1 ' 8 . (3.7) In each irreducible representationR, the gluon exchange gives rise to the Coulomb potential of eq. (2.13) with the coupling α g given by eq. (2.14). The quadratic Casimir invariants for the SU p3q representations of interest are C 2 p1q " 0, C 2 p3q " C 2 p3q " 4{3, C 2 p8q " 3, therefore α g " α sˆ# 4{3,R " 1, 1{6,R " 8. As discussed in section 2.2, the strong coupling α s depends on the momentum transfer Q. In table 2, we list the average Q for the various vertices appearing in the annihilation and BSF processes, in this model. For the bound states, the momentum transfer depends itself on the strong coupling, Q " Qpα s q. In this case, we determine α s by solving the numerically the equation α s pQpαqq "α , (3.9) forã. We discuss further the effect of the α s running in the following.

Direct annihilation
XX : pairs annihilate dominantly into gluons (cf. fig. 2), with cross-section [68] σ XX : Ñgg v rel " 14 27 where S 0,r1s and S 0,r8s are the s-wave Sommerfeld factors of the colour-singlet and colouroctet states,  Table 2. The momentum transfer Q at which the strong coupling α s pQq is evaluated, for the various processes and the states participating in these processes, in the model of section 3. X X :¨¨F igure 2a. The XX : annihilation is influenced by the Sommerfeld effect due to gluon exchange. The black blob represents the perturbative part of the annihilation processes (hard scattering). X X : g g X X : g g X X : g g Figure 2b. Tree-level diagrams contributing to the hard process XX : Ñ gg. Besides the tchannel diagram, there is the corresponding u-channel (not shown). The s-channel yields a p-wave contribution, which we neglect. S0, [8] (2/7)S0,[1] + (5/7)S0, [8] x  =2 0 x  = 5 0 Figure 3. Top left: The Sommerfeld factor for s-wave annihilation, S 0 pζ S q, vs ζ S " α S g {v rel , for both attractive (ζ S ą 0) and repulsive (ζ S ă 0) Coulomb interaction. Top right: The running of α s in the scattering states, where the average momentum exchange is Q " pm X {2qv rel . For the colour-singlet and the colour-octet states, α S g,r1s " p4{3qα S s and α S g,r8s "´α S s {6. We show α S s , α S g,r1s and α S g,r8s in the velocity range 0.2 ă v rel ă 0.4, that is typical during the DM freeze-out. For comparison, we also show the strong coupling at the gluon emission vertices of the annihilation processes, α ann s , which corresponds to Q " m X (red diamonds). Bottom left: The parameter ζ S " α S g {v rel that determines the Sommerfeld effect. Bottom right: The thermally-averaged s-wave Sommerfeld factor,S 0 , forx " m X {T within the indicative range 20 ăx ă 50 during which the DM abundance freezes-out.
The function S 0 pζ S q is the s-wave Sommerfeld enhancement factor (cf. ref. [69] and appendix A), The annihilation XX : Ñ qq is p-wave suppressed and we neglect it for simplicity. In fig. 3, we show S 0 pζ S q, for both attractive and repulsive interactions, and depict the effect of the  α s running on the Sommerfeld factors. Because the momentum exchange in the scattering state is much smaller than on the gluon-emission vertices (cf. table 2), α S s is considerably larger than α ann s , as seen in the top right panel of fig. 3.

Bound-state formation, ionisation and decay Formation
As seen from eq. (3.8), only the colour-singlet XX : state interacts via an attractive potential and can form bound states. The only capture process via one-gluon emission is from the octet state, pX`X : q r8s Ñ BpXX : q r1s`gr8s . (3.13) Using d R " 3, Cp3q " 1{2, C 2 p3q " 4{3 and C 2 pGq " 3, and setting η 1 " η 2 " 1{2 and µ " m X {2, we find from eqs. (2.42) the colour-averaged BSF cross-section, where S BSF pζ S , ζ B q is given by eq. (2.42b), and here ζ S " α S g,r8s {v rel , The coupling α B s,r1s that determines the bound-state wavefunction, and the coupling α BSF s,r1s that corresponds to the gluon radiation vertex in the capture process, are shown in fig. 4. Due to the small momentum transfer (cf. table 2), they are considerably larger than α ann s that corresponds to the gluon vertices in the XX : Ñ gg annihilation. This enhances further the BSF cross-section with respect to the annihilation cross-section, as seen by comparing the two panels in fig. 5. In fig. 5, we compare eq. (3.14) to the cross-section for XX : Ñ gg (cf. eq. (3.10)). At large velocities, α s {v rel ! 1, the BSF cross-section scales as σ r8sÑr1s BSF v rel 9pα s {v rel q 4 and is subdominant to annihilation. At low velocities, α s {v rel " 1, it becomes exponentially suppressed due to the Coulomb repulsion in the scattering state. σ r8sÑr1s BSF v rel peaks at α s {v rel « 1, where it exceeds the annihilation cross-section by more than one order of magnitude.
The thermally-averaged BSF cross-section is where µ " m X {2 is the X´X : reduced mass. Here, f g pωq " 1{pe ω{T´1 q is the gluon occupation number, with ω being the energy of the emitted gluon, The factor 1`f g pωq accounts for the Bose enhancement due to the final-state gluon, and is necessary to ensure the detailed balance between the bound-state formation and ionisation processes at T Á ω [4], which encompasses a significant temperature range that is relevant to the DM freeze-out.

Ionisation
The ionisation and BSF cross-sections are related via the Milne relation, which we review in appendix D. From eq. (D.4), we find where g g and g B are the gluon and bound-state degrees of freedom. The ionisation rate is Using eqs. (3.16) and (3.17), we obtain the ionisation rate of the colour-singlet bound states, We compute (3.18) using g g " 8, g X " 3, g B,r1s " 1, µ " m X {2 and eq. (3.14).

Decay
The decay rate of " 0 bound states is related to the perturbative s-wave annihilation cross-section times relative velocity of the corresponding scattering states (see e.g. [57]), Γ dec " pσ s´wave ann,rRs v rel q |ψ rRs n m p0q| 2 . Note that σ s´wave ann,rRs v rel corresponds to the colour configuration of the bound state and should be averaged over the bound-state colour degrees of freedom, rather than those of an unbound XX : pair.

Effective bound-state formation cross-section
The effect of unstable bound states on the DM relic density is governed by a system of coupled Boltzmann equations for the unbound and bound particles that describe the interplay between bound-state formation, ionisation and decay processes [4]. However, it is possible to incorporate the effect of bound states in a single Boltzmann equation for the unbound particles, using an effective BSF cross-section that is weighted by the fraction of bound states that decay (rather than getting ionised), Γ dec,r1s`Γion,r1s˙. (3.23)

Relic density
The total cross-section of the processes that deplete XX : is , we obtain xσ eff v rel y which enters the Boltzmann eq. (3.5). In fig. 6, we show xσ eff v rel y as a function of the time parameter m X {T , together with the contributions it receives from direct annihilation and from BSF. At early times, the depletion of DM via BSF is impeded by the large ionisation rate of the bound states. However, BSF becomes more efficient than direct annihilation in depleting DM at m X {T Á 70 for m X " TeV, suggesting that a sizeable effect on the DM density should be expected.  Left panel: The mass splitting ∆m between DM and its coloured co-annihilating partner, that is required to obtain the observed DM density. The blue dotted band takes into account perturbative annihilation only, the purple dashed band incorporates the Sommerfeld effect on the direct annihilation, and the yellow solid band includes also the effect of bound-state formation and decay. The width of the bands arises from the 3σ uncertainty on the DM density. Right panel: The impact of the Sommerfeld effect and bound-state formation on the DM density. ∆m is fixed with respect to m χ along the yellow solid band of the left panel. We present the ratios of the relic densities predicted by perturbative annihilation only (blue dotted line) and by Sommerfeldenhanced annihilation (purple dashed line), to the relic density predicted by the full computation that includes the effect of bound states.
In fig. 7 we present the results of the relic density computation. In the left panel, we show the mass splitting ∆m " m X´mχ vs. m χ , in three different cases, (i) considering perturbative annihilation only, (ii) taking into account the Sommerfeld effect on the direct annihilation processes, and (iii) including the formation and decay of unstable bound states. In agreement with previous works [33,38], we confirm that the Sommerfeld effect has a considerable impact on the predicted ∆m and m χ . In addition, we find that BSF has a significant effect. It implies that the mass splitting can be as high as " 38 GeV, and DM can be as heavy as 3.3 TeV. For the viable m χ and ∆m values determined by the full computation, we show in the right panel of fig. 7, the depletion of DM due to the Sommerfeld enhancement of the direct annihilation and due to BSF. We find that BSF depletes DM by p40´240q%. Clearly, this far exceeds the experimental uncertainty on the DM density.

Conclusion
Long-range interactions imply that non-perturbative effects and a variety of radiative processes come into play. Here, we have considered the radiative capture of non-relativistic particles into bound states, in unbroken non-Abelian gauge theories, in the regime where the gauge coupling is perturbative. This can be important in multi-TeV WIMP DM sce-narios, in scenarios where DM co-annihilates with coloured particles, as well as in hidden sector models.
Our main results include the amplitude for the radiative formation of bound states via one-gluon emission, for arbitrary representations and masses of the interacting particles [cf. eq. (2.25)], and the complete BSF cross-sections for particles transforming in conjugate representations [cf. eqs. (2.42)], but still for arbitrary masses.
As a first application of our results, we considered a simplified model where DM coannihilates with particles transforming in the fundamental of SU p3q c , and showed that the formation and decay of particle-antiparticle bound states can affect the DM relic density very significantly. This implies larger DM mass and/or mass splitting between DM and its coannihilating particles, thereby altering the interpretation of the experimental results, and affecting the detection prospects. In particular, larger mass splittings imply the production of harder jets that can be more easily probed in collider experiments. Moreover, larger DM masses motivate indirect searches in the multi-TeV regime.
While the analytical formulae (2.42) assume a Coulomb potential, it is straightforward to generalise our results to other potentials, by computing the overlap integrals (2.26) using the wavefunctions arising from those potentials. This allows to include, for example, thermal masses for the gauge bosons, as well as the effect of multiple force mediators. The latter has been shown to be important in models where the (co-)annihilating particles possess a significant coupling to the Higgs [58,70]. We leave these extensions for future work.

Appendices A Scattering-state and bound-state wavefunctions
The non-relativistic potentials due to gluon exchange, for the scattering and bound states are where α S g may be either positive or negative, but α B g ą 0. The scattering and bound states are characterised by the momenta where v rel is the expectation value of the relative velocity in the scattering state and κ is the Bohr momentum of the bound state. The corresponding wavefunctions, φ k prq and ψ n m prq, with tn, , mu being the principal and angular-momentum quantum numbers, obey the Schrödinger equations The wavefunctions are normalised according to For convenience, we shall define and The solutions to eqs. (A.3) with the potentials (2.13), are (see e.g. [71]) where 1 F 1 is the confluent hypergeometric function of the first kind, and L a n are the generalised Laguerre polynomials of degree n. (We assume the normalisation condition ş 8 0 z a e´zL paq n pzqL paq m pzqdz " rΓ pn`a`1q{n!s δ n,m .) For the ground state, tn, , mu " t1, 0, 0u, Note that S 0 is the Sommerfeld factor for s-wave annihilation (see e.g. [69]) In section 2, we also need the Fourier transforms of the wavefunctions, defined as In deriving eq. (B.1b), we Fourier-transformedψn m ppq andφ k pqq, and used the following integral in the limit m g Ñ 0, In the following, we consider capture into the ground state only, tn mu " t100u. Following refs. [9,16,57], we compute the overlap integrals (B.1) using the identity [72] ż d 3 r e ipk´bq¨r´κr 4πr 1 F 1 riζ S , 1, ipkr´k¨rqs " where For the cross-sections of interest, we only need J k,t100u pbq evaluated at b " 0 [9,57].
We recall that ζ S and ζ B are defined in eqs. (A.6). In section 2.5, we use eqs. (B.4) to obtain analytical expressions for the BSF cross-sections.

C The non-relativistic Hamiltonian from effective field theory
Our results in section 2 differ from previous computations [55,56] in the relative sign of the Abelian and non-Abelian contributions to the radiative transition amplitude. In this appendix, we use the non-relativistic QCD (NRQCD) approach of ref. [47,73,74] to derive the effective Hamiltonian for our system, and compare it with refs. [55,56], whose computations are based on effective field theory. This offers an independent check of our computations. Before moving to NRQCD, we want to display that the discrepancy between our computations and ref. [55] originates from the expression for the transition amplitude. Indeed, starting from eq. (2.25) with η 1 " η 2 " 1{2, and using the coordinate-space expressions for the overlap integrals (B.1) where we integrate eq. (B.1a) by parts, we arrive at irM kÑtn mu s ii 1 ,jj 1 "´irM kÑtn mu s a where a is the gluon polarisation vector. Equation (C.1) may be directly compared to eqs. (41)-(43) of ref. [55]. We observe that the relative sign between the two terms is different.
This discrepancy arises from the non-relativistic Hamiltonian assumed in ref. [55]. The interactions that give rise to the radiative BSF amplitude correspond to a non-relativistic potentialṼ BSF pq, pq that can be deduced from the transition amplitude xp|iT |qy in ordinary quantum mechanics, xp|iT |qy "´iṼ BSF pq, pq , where |qy denotes a state where the two interacting particles have momentum q and´q in the CM frame. The quantum mechanical transition amplitude is related to the matrix element Mpq, pq "´M a trans pq, pq¨ a via p4m 2 X {A 0 q xp|iT |qy "´iM a trans pq, pq¨ a p2πqδpE q´Ep´ω q , where the factor 4m 2 X {A 0 accounts for the different normalization of fields in quantum field theory and quantum mechanics, with A 0 being the non-relativistic normalisation of the gauge field. E q , E p are the energies of the |qy, |py states and ω is the energy of the radiated gauge boson. From eqs. (C.2) and (C.3), we identify the non-relativistic potential in momentum space as The potential in coordinate space is 8 Using eq. (2.24) for M a trans pq, pq and the identity (B.2), we obtain .
Identifying A a pt, xq " A 0 expripωt´P g¨x qs a as the background field that induces the transition (see e.g. [75, sec. 5.7]), we rewrite the above as follows Comparing this with eq. (36) of ref. [55], we note the difference in the relative sign of the Abelian and non-Abelian contributions.
We now move on to NRQCD. At leading order, a heavy quark-antiquark system yields the same potential as a scalar particle-antiparticle system. Thus, we may compare our result of eq. (C.7) with that derived for quarkonium in NRQCD, where the "potential" gluons are integrated out, leaving four-quark operators in the effective Lagrangian in analogy to the Fermi theory. There are different formulations, e.g. pNRQCD [47,76] or vN-RQCD [73,74,77] 9 . Here, we will follow the conventions of [73,74,79]. We start from the ultrasoft NRQCD Lagrangian given in eq. (7) of ref. [73] L u Ą ÿ p ψ : where D µ " B µ`i g s A aµ T a " pD 0 ,´Dq with D 0 " B 0`i g s A a0 T a and D " ∇´ig s A a T a . The momentum p represents momenta of the soft scale. In contrast to the usual relativistic conventions, here ψ : (χ : ) annihilates (anti)particles, while ψ (χ) creates (anti)particles. Their generators are related via T a ij "´T a ji . m Q denotes the (anti)quark mass. From eq. (C.8), we can derive the Feynman rules for the fermion-antifermion system, treating the temporal and spatial component of the gauge field separately. We present them in fig. 8.
We may now compute the contributions to the radiative transition amplitudes, shown in fig. 9. In the CM frame, the momenta of the incoming/outgoing fermions are where P g is the momentum of the emitted gluon. The energies of the scattering and bound states are where k is the expectation value of q and M B " 2m Q´EB is the mass of the bound state, with E B being the binding energy (cf. appendix A). The energy of the radiated gluon, ω " |P g |, is found from the conservation of energy to be [cf. eq. (2.20)], In the following, we shall extract the leading order contributions to the radiative BSF amplitude taking into account that |q| " |k| " pm Q {2qv rel , |p| " κ " pm Q {2qα s and ω " pm Q {4qpα 2 s`v 2 rel q. From the Lagrangian of eq. (C.8) and the Feynman rules of fig. 8, it is straightforward to obtain the Abelian contributions to the transition amplitude, shown in fig. 9a, Figure 8. Feynman rules that are used to derive the non-relativistic potential that determines the radiative capture into bound states. In the gauge fields A a,µ , the first index denotes colour, while the second one if the space-time index. The first column depicts the NRQCD Feynman rules which distinguish between temporal (black dot) and spatial (black rectangle) couplings in the quark-gluon vertex. Hereby, ψ represents a fermion field. The Feynman rules for the anti-fermion field χ are obtained by simply substituting T a ij Ñ T a ij "´T a ji . The second column shows the non-relativistic 3-gluon-vertex, and the third column shows the gluon propagators. We distinguish again between temporal and spatial gauge fields. Note that according to the NRQCD Lagrangian (C.8), the NRQCD Feynman rules are defined with upper indices only, and incorporate the signs arising from the metric, g µν " diagp1,´1,´1,´1q, in 4-vector contractions done according to the relativistic conventions.
We shall now derive the potential NRQCD Lagrangian term that describes the non-Abelian contribution to the radiative transition amplitudes. We will demonstrate that (a) we recover the same sign as in [73,74,79] using our conventions eqs. (C.2) and (C.3), and that (b) this sign disagrees with the result of [55]. From the NRQCD Feynman rules of fig. 8, it is immediately evident that the first of the four diagrams shown in fig. 9b yields   9 For a comprehensive comparison, we refer to [78]. . Non-Abelian contributions to the capture of a fermion-antifermion pair into a bound state, via gluon emission. The fermion field ψ is depicted as a single solid line, the antifermion field χ as a double line. As in fig. 8, we distinguish between the temporal (circle) and spatial (rectangle) components of A aµ .
the dominant contribution (see also comment below eqs. (2.23)). The diagrams involving the spatial components of the gluon propagators are suppressed by higher orders in the momenta q and p, as shown explicitly below. Allowing for a non-zero gluon mass m g , we In order to compare eq. (C.13) with [73,79], we construct the effective action that recovers the amplitude, xp 1 , p 2 |iS eff |q 1 , q 2 y " p2πq 4 δ p4q pq 1`q2´p1´p2 q iW trans pq 1 , q 2 , p 1 , p 2 q , (C.14) such that S eff " ż d 4 q 1 d 4 q 2 d 4 p 1 d 4 p 2 p2πq 16 ψ : pp 1 qχ : pp 2 q W trans pq 1 , q 2 , p 1 , p 2 q ψpq 1 qχpq 2 q p2πq 4 δ p4q pq 1`q2´p1´p2 q . (C. 15) Exchanging the polarisation vector a for the background field A a , we arrive at the Lagrangian This Lagrangian can be compared with eqs. (6), (7), (12) and the first term of eq. (16) of ref. [73]. 10 It agrees perfectly, including the sign of the non-Abelian term. 11 Similarly to before, we may now derive the full interaction potential in position space. Starting from eq. (C.13) and following similar steps as in the beginning of this appendix, we arrive again at eq. (C.7), with the identification T a i 1 i " pT 1 q a i 1 i , T a j 1 j " pT 2 q a j 1 j "´pT 1 q a jj 1 and m Q " m X , thereby confirming the validity of the computation of section 2 and the disagreement with ref. [55]. We note that we have also compared our results to those of ref. [47] and found them in agreement. 12

D The Milne relation
For a 2-to-2 process 1`2 Ø A`B, the cross-section, averaged over the degrees of freedom of the initial state, is 10 Note that k in eq. (16) of ref. [73] is defined as k "´pq´pq. 11 The effective Lagrangian is invariant under ultrasoft gauge transformations [47,79]. 12 For comparison with eq. (1.8) in ref. [47], it is important to note the different definition of fields with respect to ref. [73]. In ref. [47], ψ (χ : ) annihilates (anti)particles, while ψ : (χ) creates (anti)particles. In that notation, one has to account for the odd Wick permutations, such that the corresponding non-Abelian amplitude M 1 NA is recovered by S eff " ż d 4 q1d 4 q2d 4 p1d 4 p2 p2πq 16 ψ : pp1qχp´p2q p´1q M 1 NA pq1, q2, p1, p2qψpq1qχ : p´q2q , (C. 17) in contrast to eq. (C.15). From M 1 NA , we may derive the non-relativistic potential, as in the beginning of this appendix, and compare it with eq. (C.7). A global sign difference is expected due to the different definition for the covariant derivative: while here we have used Dµ " Bµ`igsA a µ T a , in ref. [47] the covariant derivative is defined as Dµ " Bµ´igsA a µ T a .
where s is the first Mandelstam variable and P CM A is the momentum of A (or B) in the CM frame. Since |M 1`2ÑA`B | 2 " |M A`BÑ1`2 | 2 , the cross-sections of two 2-to-2 inverse processes are related via We now consider the radiative BSF and ionisation processes X 1`X2 Ø B n `g, where g stands generally for the massless field radiated in a BSF process with energy ω. In the non-relativistic regime, s » pm 1`m2`Ek q 2 » pm 1`m2`En `ωq 2 , (D. 3) where E k " µv 2 rel {2 is the kinetic energy of the initial state in the CM frame, and E n ă 0 is the binding energy of the bound state (cf. appendix A), with E k , |E n | ! m 1`m2 . Then, |P 1 | " |P 2 | " µv rel and |P g | " |P B | " ω " E k´En , and from eq. (D.2) we find the Milne relation σ ion " g 1 g 2 g g g Bˆµ