Dark-matter bound states from Feynman diagrams

If dark matter couples directly to a light force mediator, then it may form bound states in the early universe and in the non-relativistic environment of haloes today. In this work, we establish a field-theoretic framework for the computation of bound-state formation cross-sections, de-excitation and decay rates, in theories with long-range interactions. Using this formalism, we carry out specific computations for scalar particles interacting either via a light scalar or vector mediator. At low relative velocities of the interacting particles, the formation of bound states is enhanced by the Sommerfeld effect. For particle-antiparticle pairs, we show that bound-state formation can be faster than annihilation into radiation in the regime where the Sommerfeld effect is important. The field-theoretic formalism outlined here can be generalised to compute bound-state formation cross-sections in a variety of theories, including theories featuring non-Abelian (albeit non-confining) interactions, such as the electroweak interactions.


Introduction
Dark matter (DM) with long-range interactions, mediated by a light or massless force carrier, appears in diverse theories motivated on various grounds. Let us mention a few important examples. Dark matter with sizable self-interactions, mediated by a light particle, can explain the observed galactic structure better than collisionless DM [1][2][3][4][5][6][7][8][9][10]. Asymmetric DM [11][12][13][14] -motivated in part by the similarity of the dark and the ordinary matter abundances -resides, in most implementations, in a hidden sector that often includes long-range interactions [15][16][17][18][19][20]. Dark matter which dissipates energy via its coupling to a light force mediator may provide a dynamical explanation for some of the observed scaling relations governing haloes [21][22][23][24], and other features [25,26]. Inside haloes, the inelastic scatterings of either symmetric or asymmetric DM with long-range interactions can produce radiative signals [27][28][29][30][31][32][33][34][35], which can potentially account for anomalous excesses observed in the radiation backgrounds [27][28][29][30][31][32]. Moreover, long-range DM-nucleon scattering implies a different interpretation of the direct-detection data than the commonly assumed short-range scattering [36][37][38][39]. Notably, the long-range character of DM interactions is relevant not only for theories involving hidden sectors; even the electroweak interactions of the Standard Model manifest as long-range if DM is heavier than a few TeV [40][41][42]. Clearly, long-range interactions play a central role in the venture to identify DM. In order to extract accurate predictions for the DM phenomenology, it is then essential to fully understand their implications. Long-range interactions typically imply the existence of bound states. The formation of DM bound states in the early universe and/or in the dense environment of haloes today affects the phenomenology of DM in many important ways. In the early universe, symmetric DM may form unstable bound states, whose decay contributes to the DM annihilation rate and affects the DM relic abundance [43]. Asymmetric DM may form stable bound states [15-20, 34, 44-52], which affect all manifestations of DM today, including the DM self-scattering in haloes [19,[46][47][48], the expected indirect-detection [29][30][31][32][33][34] and direct-detection [53] signals, as well as the kinetic decoupling of DM from dark radiation [54]. Inside haloes, DM bound states -whether they are stable or unstable -may form radiatively and yield detectable signatures [28][29][30]35]. Radiative signals may also be produced in transitions between the bound-state energy levels [31][32][33], or other related inelastic processes [34]. Moreover, the formation of DM bound states may be detectable at the LHC [55]. These rich phenomenological implications strongly suggest that it is critical to accurately account for the formation of DM bound states, in order to constrain the DM properties and eventually detect DM. This work is a step toward this goal.
There are, of course, two different classes of bound states: Those arising due to non-confining interactions, such as the atomic bound states in QED, and those arising due to confining interactions, such as the hadrons in QCD. Particles charged under a confining force always combine into hadronic states, roughly once the kinetic energy in their center-of-mass frame (or the temperature of their plasma) drops below the confinement scale. On the other hand, in the case of non-confining interactions, the efficiency of bound-state formation (BSF) depends on the corresponding cross-sections and the details of the thermodynamic environment. Here we shall consider bound states due to non-confining interactions, and calculate the cross-sections for their formation in the non-relativistic regime, which is relevant for cosmology and DM indirect detection signals. 1 The formation of DM bound states in the early universe and inside the non-relativistic environment of haloes differs from BSF in colliders in some important ways. In high-energy colliders, the initial-state particles are highly relativistic. 2 However, BSF is more efficient when the relative velocity of the interacting particles is lower than the expectation value of the relative velocity of the particles inside the bound state. Equivalently, this is when the kinetic energy in the center-of-mass (CM) frame is lower than the binding energy; clearly, this lies within the non-relativistic regime. In this regime, the long-range interaction distorts the wavefunctions of the incoming particles, which cannot be approximated by momentum eigenstates (plane waves). This is the well-known Sommerfeld effect [62]. If the interaction is attractive, the Sommerfeld effect enhances the crosssection for any process the two particles may participate in, including the formation of bound states. It follows that, cosmologically, DM bound states form most efficiently after the temperature of the dark plasma drops below the binding energy. While the formation of bound states in the early universe eventually freezes out due to the cosmological expansion, bound states may again form efficiently in today's dense and non-relativistic haloes. In either case, to accurately estimate the formation of DM bound states, we must account for the Sommerfeld effect, which is a non-perturbative phenomenon, as is of course the very existence of bound states.
In this paper, we establish a field-theoretic framework for the calculation of BSF cross-sections and other level-transition rates, as well as the decay rates of unstable bound states. Then, we carry out computations for specific interactions. We organise our work as follows.
We begin, in Sec. 2, with reviewing how to determine the wavefunctions of the two-particle states and the bound states in the presence of a long-range interaction. We derive the Bethe-Salpeter equation for the wavefunctions, and reduce it to the Schrödinger equation using the instantaneous approximation in the non-relativistic regime. In Sec. 3, we determine the amplitudes for transitions between energy levels with the emission of a force mediator; this includes the radiative capture to a bound state. In the fully relativistic regime, we express the amplitudes for such processes in terms of the Bethe-Salpeter wavefunctions describing the initial and final states, and a perturbative interaction. We then employ the instantaneous and the non-relativistic approximations, to express the transition amplitudes in terms of the Schrödinger wavefunctions. In Sec. 4, we repeat the analysis for the decay of unstable bound states into radiation, as well as for the (co-)annihilation of unbound pairs of particles into radiation -two closely related processes. While we lay out our formalism in terms of scalar particles, it is straightforward to extend it to include fermionic species.
We continue by applying our formalism to specific interactions. In Sec. 5, we consider scalar particles interacting either via a scalar or a vector mediator, and calculate the cross-sections for the dominant radiative capture to a bound state. We estimate the range of validity of our computations, using the unitarity bound on the inelastic cross-section. For bound states made of particle-antiparticle pairs or pairs of self-conjugate particles, we compare BSF with annihilation, and show that BSF can be the dominant inelastic process in the regime where the Sommerfeld effect is important; we sketch this comparison in Fig. 1. In addition, we calculate the decay rates of particle-antiparticle bound states into force mediators. We cast our results in terms of a minimal parametrisation, which makes their potential implications more transparent, and summarise them in table 1. We conclude with a discussion of the phenomenological implications of DM bound-state formation in Sec. 6, and present many of the detailed calculations in the appendices.
The field-theoretic formalism developed in this work has several advantages in comparison to previous quantum mechanical calculations [29,35,51,63,64]. It can accommodate the possibility of DM coupled to non-Abelian interactions. Such interactions can convert the incoming particles into different species which may subsequently form bound states; importantly, DM coupled to the electroweak interactions of the Standard Model belongs to this category. Moreover, the field-theoretic approach allows for a systematic inclusion of higher-order corrections, both in the interaction strength and in the momentum transfer between the interacting degrees of freedom.
We emphasise that BSF is distinct from processes such as the direct annihilation into mediators or elastic scattering, in which particles coupled to a long-range interaction may participate. While all these processes are influenced by the Sommerfeld effect, the final-state particles are obviously different. Field-theoretic treatments of the annihilation processes, analogous to the formalism presented here for BSF and discrete level transitions, have been presented in Refs. [65,66].

Scalar mediator
Vector mediator non-degenerate species Table 1. Summary of the bound-state-formation and annihilation cross-sections, the decay and deexcitation rates, computed in Sec. 5, for scalar particles interacting via a light scalar or vector mediator. The annihilation cross-sections and the rates of decay into ϕ-mediators refer to unbound and bound particle-antiparticle pairs respectively; all other formulae apply to any pair of particles. The cross-sections are normalised to σc ≡ (πα 2 /µ 2 ) × 2πζ/(1 − e −2πζ ), where α is the fine-structure constant entering the Coulomb potential, µ = m1m2/(m1 + m2) is the reduced mass of the interacting species, and ζ = α/v rel , with v rel being the relative velocity of the incoming particles. Also, η1,2 = m1,2/(m1 +m2), and g1,2, c1,2 are the couplings of the interacting particles to the force mediators, as described by the Lagrangian densities of Sec. 5. In our computations, we have neglected the mediator mass. Unitarity suggests that the range of validity of the above computations is α 0.5.
2 Bound-state and two-particle state wavefunctions We shall consider two scalar particle species χ 1 and χ 2 , interacting via a light or massless force mediator ϕ. Specific interaction Lagrangians will be introduced in Sec. 5 (c.f. Eqs. (5.13), (5.14) and (5.48)). In this section, we aim to determine the wavefunctions which describe two-particle states and bound states with the quantum numbers of χ 1 and χ 2 . Our presentation draws largely from the pedagogical discussions of Refs. [67,68]. Let us first introduce some notation. We denote the masses of χ 1 , χ 2 and ϕ by m 1 , m 2 and Vector mediator Figure 1. Comparison of the cross-sections for annihilation (dashed blue lines) and the dominant radiative capture process to a bound state (solid purple lines), for a scalar particle-antiparticle pair. The crosssections times relative velocity, σv rel , are normalised to the perturbative value for s-wave annihilation, σ0. Left: For a scalar mediator, σ0 = πα 2 /µ 2 . The leading capture process is to the excited state {210} and consists dominantly of the J = 2 partial wave. At ζ 1, the ratio of the bound-state formation and annihilation cross-sections is σ BSF /σann 0.066. Right: For a vector mediator, σ0 = πα 2 /(2µ 2 ). The leading capture process is to the ground state {100} and consists dominantly of the J = 0 and J = 2 partial waves. At ζ 1, σ BSF /σann 1.56. m ϕ respectively. We define the total and the reduced χ 1 , χ 2 masses, Obviously, µ m/4. In the following, |B Q,n stands for a χ 1 − χ 2 bound state, of total momentum Q and energy ω Q,n = Q 2 + M 2 n , where n denotes collectively all the discrete quantum numbers characterizing the bound state, and M n < m is the bound-state mass. |U Q,q stands for a χ 1 − χ 2 unbound two-particle state, with total momentum Q, expectation value of relative velocity v rel = q/µ, and energy ω Q,q m. (As is common in scattering theory, we shall often refer to the unbound states as scattering states.) Moreover, |ϕ Q stands for a ϕ particle state with momentum Q and energy ω ϕ = Q 2 + m 2 ϕ .

The Bethe-Salpeter wavefunctions
We shall now introduce the Bethe-Salpeter wavefunctions that will appear in the fully relativistic version of the transition amplitudes of sections 3 and 4. The Bethe-Salpeter wavefunctions are related to the more familiar Schrödinger wavefunctions, which we will use in evaluating the amplitudes of interest in the non-relativistic regime. We are interested in the following wavefunctions: where T is the time-ordering operator, and |Ω is the vacuum of the interacting theory. If χ 1 , χ 2 are a particle-antiparticle pair, in the above definitions we replace χ 1 → χ, χ 2 → χ † . Note that does not denote complex conjugation, for which we shall use the symbol * , as usual. In fact whereT is the anti-time-ordering operator. We define the coordinate transformation and its inverse where η 1 + η 2 = 1 for the Jacobian to be 1, and we choose specifically In the non-relativistic regime, this choice will enable us to separate the motion of the CM from the relative motion. Using the 4-momentum operatorP , we obtain Then, the wavefunction of Eqs. (2.3) becomes where it is understood that Q 0 = ω Q,n , and we defined This is the first step in separating the motion of the CM from the relative motion. For notational simplicity, we are using the same symbol for Ψ Q,n (x 1 , x 2 ) and Ψ Q,n (x). We also define the Fourier transforms We repeat the above for the amplitudes of Eqs. (2.4) -(2.6), and summarise the definitions in appendix A.

The 4-point Green's function and Dyson-Schwinger equation
Consider the 4-point Green's function and let W (x 1 , x 2 , y 1 , y 2 ) be the perturbative 4-point interaction kernel between χ 1 and χ 2 . Then, G (4) satisfies the Dyson-Schwinger equation, 18) where S 1 , S 2 are the full propagators for χ 1 , χ 2 . Symbolically, the above series can be written as Fig. 2. Due to translational invariance, W and G (4) depend only on coordinate differences. We shall take them to be x, y, X − Y , where we used the definitions of Eq. (2.9) and assumed analogous definitions for the y 1 , y 2 variables. Thus where we retained the same symbols to keep the notation simple. Equation (2.18) becomes We define the Fourier transforms of G, W, S 1 and S 2 , withS j (p) being the momentum-space propagator for χ j . From the above, we deduce the relation between the conjugate momenta of x 1 , x 2 , which we shall call here p 1 , p 2 , and the conjugate momenta of x, X, denoted above as p, Q: Analogous relations hold between the conjugate momenta of y 1 , y 2 and those of y, Y . For convenience, we also define We may now rewrite the Dyson-Schwinger Eq. (2.22) for the 4-point function, in momentum spacẽ We shall use Eq. (2.29) to derive the Bethe-Salpeter equation for the wavefunctions of Sec. 2.1.

Completeness relation and decomposition of the 4-point function
To compute the Bethe-Salpeter wavefunctions of Sec. 2.1, we have to decompose the 4-point Green's function of Sec. 2.2 using the one-and two-particle completeness relation. Then, Eq. (2.29) will yield the equations which the wavefunctions Ψ Q,n and Φ Q,q satisfy. Including the one-and two-particle states with the same quantum numbers as χ 1 and χ 2 , the completeness relation is where we have assumed the standard relativistic normalization of one-particle momentum eigenstates p|k = 2E p (2π) 3 δ 3 (p − k), with E p being the energy of the state |p . To lowest (zeroth) order in the interaction strength, Next, we insert the unity operator of Eq. (2.30) in G (4) , to obtain the decomposition where G (4) n (x, y; X − Y ) and G U (x, y; X − Y ) are the contributions of the bound and the scattering states, respectively. We compute them below. We shall make use of the fact that a non-zero contribution to G (4) from a one-or two-particle state arises only when two annihilation operators act on that state to obtain the quantum numbers of the vacuum. Moreover, in order to extract the poles and the branch-cuts of G (4) n and G (4) U , we will use the integral representation of the θ-function, and

Contribution of the bound states to the 4-point function
The contribution of the nth bound state to G (4) is where in the third step we made use of the integral representation of the θ function, given in Eq. (2.34), which introduces the integration over K 0 . The Fourier transform of the above is .
(2.38) Equation (2.37) is the contribution of the nth bound state toG (4) (p, p ; Q). Evidently, the scattering amplitude has a pole at energy equal to the bound-state energy.

Contribution of two-particle scattering states to the 4-point function
Following similar steps, we find the contribution of the two-particle states to G (4) , The Fourier transform of G (4) (2.39) Clearly, the contribution of the two-particle states toG (4) (p, p ; Q) gives rise to a branch-cut in the scattering amplitude.
Summing the contributions from the bound and the scattering states, we obtain the decomposition of the momentum-space 4-point functioñ  (2.40) We shall now combine the Dyson-Schwinger Eq. (2.29) and Eq. (2.40), to obtain the Bethe-Salpeter equation for the wavefunctions Ψ Q,n and Φ Q,q .

The Bethe-Salpeter equation for bound and scattering states
We introduce the operator Then, the Dyson-Schwinger Eq. (2.29) can be cast in the form This is formally solved bỹ where C n (p; Q) and F a (p; Q) are the eigenfunctions of the discrete and the continuous spectrum of the operator A(p, q; Q), with eigenvalues c n (Q) and f a (Q) respectively, For the continuous spectrum, we identify a → q, and deduce The relations (2.47), (2.48), (2.50) and (2.51) are stipulated because c n , f a are independent of the momenta p, p ; all the p, p -dependent factors must arise from the eigenfunctions, C n and F a . The relations (2.49), (2.52) are warranted so that C n and F a are not singular in the limit Q 0 → ω Q,n and Q 0 → ω Q,q respectively; the factors [1 − ω Q,n /Q 0 ] −1 , [1 − ω Q,q /Q 0 ] −1 cannot be part of the eigenfunctions, and thus belong to the eigenvalues.
Inserting the above into the eigenvalue equations (2.44) and (2.45), and taking the limits Q 0 → ω Q,n and Q 0 → ω Q,q respectively, we obtain the Bethe-Salpeter equations for the bound and the scattering states 3Ψ Q,n (p) = S(p; Q) These are homogeneous equations and do not determine the normalisation of Ψ Q,n and Φ Q,q . Moreover, because we do not know the exact eigenvalues c n and f a , we cannot use Eq. (2.46) to obtain the normalisation of Ψ Q,n and Φ Q,q . We derive their normalisation in the next section.

Normalization of the Bethe-Salpeter wavefunctions
We derive the normalisation of the wavefunctions Ψ Q,n and Φ Q,q from the inhomogeneous Dyson-Schwinger Eq. (2.29), or equivalently from Eq. (2.42), using the method described in Ref. [68].
We define the symbolic multiplication We shall use Eq. (2.57) to obtain the normalisation of the Bethe-Salpeter wavefunctions. For later convenience, we defineÑ and their Fourier transforms, (2.61)

Bound states
Substituting the contribution to the 4-point function from the nth bound state, given in Eq. (2.38), into Eq. (2.57), and taking the limit Q 0 → ω Q,n , we obtain the normalisation condition In coordinate space, this becomes (2.63)

Two-particle states
Substituting the contribution to the 4-point function from the two-particle states, Eq. (2.39), into Eq. (2.57), we deduce the normalisation condition 4 In coordinate space, this becomes Note that in the fully relativistic case, the normalisation of the wavefunctions depends in general on the potential.

The instantaneous approximation and the Schrödinger equation
In the non-relativistic regime, it is possible to simplify the Bethe-Salpeter equations. The momentum exchange between the two unbound particles is |q| ∼ µv rel , while between two bound particles |q| ∼ µα, with α characterising the interaction strength; in either case, for α, v rel 1, the energy exchange is |q|. It is then reasonable to ignore the dependence of the kernel W (p, p ; Q) on p 0 , p 0 . This is the instantaneous approximation. 5 In fact, in the cases of interest, W (p, p ; Q) depends only on |p − p |, rather than on p and p separately, and it does not depend on Q (except perhaps for Q 2 , which, in the non-relativistic regime, we shall approximate with Q 2 m 2 ). We shall thus assume that In this approximation, we deduce from the Bethe-Salpeter Eqs.
where we choose the normalization factor such that we recover the conventional normalisation forψ Q,n andφ Q,q , as we shall see in Sec. 2.7. We calculate S 0 (p; Q) in appendix C. Multiplying both sides of Eqs. (2.68), (2.69) with S(p; Q), integrating over p 0 , and using Eq. (2.67), it follows that we obtain the exact form of Eq. (2.64). 5 For a discussion on relativistic corrections, see e.g. Ref. [69] and references within. ψ Q,n (p) andφ Q,q (p) are sometimes called the "equal-time wavefunctions". From the above, and recalling Eqs. (2.7), (2.8) we see that

Non-relativistic approximation
Using the non-relativistic approximation described in appendix C, Eqs. (C.18) -(C. 22), and setting, in accordance to Eq. (C. 21), (2.79) These are the Schrödinger equations for the bound and the scattering states in momentum space. They are eigenvalue equations, and as such, their solutions determine E n and E q . Because in Eqs. (2.78), (2.79), all dependence on the CM momentum Q has been eliminated, we have dropped this subscript from the ψ, φ wavefunctions, but kept the same symbols in order to avoid cluttering the notation. Note that from Eq. (2.76), it follows that the mass of the bound state is It is convenient to Fourier-transform Eqs. (2.78) and (2.79) to coordinate space. We set where V (r) is the non-relativistic potential, (2.85) We quote the bound-state and scattering-state solutions of the Schrödinger equation for a Coulomb potential, in appendix F, and use them in our computations in Sec. 5.

Normalization of the Schrödinger wavefunctions
To find the normalization of ψ n and φ p , it is easiest to follow a similar procedure to that of Sec. 2.5. We first define (2.88) Integrating Eq. (2.29) with respect to p 0 , k 0 yields Following the steps of Sec. 2.5, we find the equivalent of Eq. (2.57), We may now obtain the normalisation conditions for the wavefunctions ψ n and φ q .

Bound states
Close to the pole, at Q 0 → ω Q,n , we may substitute the contribution from the nth bound state, Eq. (2.87), into (2.90). We obtain (2.92)

Two-particle states
Substituting the contribution from the two-particle states, Eq. (2.88) into (2.90), and for Q 0 → ω Q,q , we deduce the normalisation condition Note that in obtaining the normalisation conditions (2.92) and (2.93), we did not make use of the non-relativistic expansions of the factors N Q (p) and ε Q,q , given in Eqs. (3.32) and (3.34).

Radiative level transitions
In this section, we determine the amplitudes for the radiative BSF and de-excitation processes . Diagrammatic representation of the equation (3.12) for the 5-point functionG (5) . stands for the full propagator of the force mediator ϕ, which may be either a scalar or a vector boson.
in terms of the bound-state and scattering-state wavefunctions computed in Sec. 2 and a perturbative interaction which describes the emission of the force mediator. The S-matrix elements of interest are where the indices stand for the momenta and the quantum numbers of the corresponding states, as stated in the beginning of Sec. 2.

The 5-point Green's function
Since B Q,n and U Q,q are generated by the action of χ † 1 and χ † 2 on the vacuum (c.f. Eqs. (2.4), (2.6)), in order to compute the S-matrix elements of Eqs. (3.3) and (3.4) we need to consider the 5-point function We define the Fourier transform As in Eq. (2.9), we set and rewrite the above as i.e. the conjugate momenta of X, x are P, p, and the conjugate momenta of Y, y are K, k defined as The 5-point Green's function G (5) (X ϕ , x 1 , x 2 ; y 1 , y 2 ) is equal to the sum of all connected diagrams with five external points. The momentum-spaceG (5) is sketched in Fig. 3, and can be written asG (5) being the field-strength renormalisation parameter for ϕ. A (5) is defined via the relation Note that C (5) may include diagrams that are not fully connected, i.e. diagrams in which external legs are disconnected from each other, 6 but it does not, of course, include vacuum bubble diagrams.
(If only fully connected diagrams contributed to C (5) , then A (5) would simply be the sum of all connected and amputated diagrams, as conventionally defined.) For later convenience, we also define C (5) ϕ−amp as the sum of all connected diagrams with only the ϕ-leg amputated, Then, A (5) appearing in Eq. (3.12), becomes .28)). We sketch Eqs. (3.15) and (3.17) in Fig. 4.
Let us first focus on the BSF amplitude of Eq. (3.3), for which ϕ−amp is equal to C (5) with only the ϕ-leg amputated. and stand for the χ1 and χ2 full propagators, respectively.
stands for the full propagator of the force mediator ϕ, which may be either a scalar or a vector boson.
In this limit, the Lehmann-Symanzik-Zimmermann reduction formula yields Here, the ∼ sign means that the two sides have the same singularities in the limit (3.19); to compute the S-matrix element, we need to extract the residues of these singularities from both sides of Eq. (3.20).
In the above expression, the correlation functions involving the χ 1 and χ 2 fields correspond to the bound and scattering state wavefunctions (c.f. Eqs. (A.5), (A.8)). The correlation function involving the ϕ field is the ϕ field-strength renormalisation parameter (c.f. Eq. (3.14)). We Fouriertransform Eq. (3.20) with respect to x, y, to obtain The left side of the above equation isG (5) At P 0 ϕ → ω ϕ (P ϕ ) and P 0 → ω P,n , this expression has the same poles as the right side of Eq. (3.21). Identifying their residues, we obtain We still have to extract the leading singularity at K 0 → ω K,k . We multiply both sides of the above expression withÑ k (q, q ; K)Φ K,k (q ), integrate over q and q , and use the orthonormality condition (2.64), to obtain the S-matrix element for BSF Following similar steps, we obtain the S-matrix element for transition between discrete energy levels, B P,n ; ϕ Pϕ | S | B K,n .
In standard notation, we write the S-matrix elements as If non-fully connected diagrams contribute to the perturbative part of the transition amplitudes, then A (5) should be replaced by C (5) ϕ−amp using Eq. (3.18). In this case, we obtain In the case of a vector mediator ϕ µ , Z ϕ becomes the charge-renormalisation parameter and the amplitudes contain the polarisation vector µ , i.e. A (5)

Instantaneous approximation
In the instantaneous and non-relativistic approximations, we may express the transition amplitudes in terms of the Schrödinger wavefunctions defined in Eqs. (2.68), (2.69), as follows where we took Z ϕ (P ϕ ) 1 to lowest order, and set It is sufficient for our purposes, and consistent with our approximation (see footnote 14), to expand the normalisation factors up to first order in p 2 , q 2 , as follows The p 2 , q 2 terms in Eq. (3.33) introduce corrections of order α 2 and v 2 rel (see appendix F), where α parametrises the strength of the interaction (for a Coulomb potential, it is the fine-structure constant) and gives the expectation value of the relative velocity inside the bound state. Similar corrections arise also in M trans (see appendix E). We shall retain such corrections only where the dominant term in the respective expansion vanishes, as in the case of degenerate particles interacting via scalar boson exchange (see Sec. 5.1). Moreover, from Eqs. (2.31) and (2.77), we find that to zeroth order in the relative velocity, ε K,k µ . (3.34) Because ε K,k factors out of the integrals, as seen in Eq. (3.29), we neglect k 2 corrections, which always produce subdominant terms in v 2 rel . In addition, in our computations, we consistently ignore corrections that involve at least one power of the total momentum of any of the initial or final states (denoted typically with capital letters). In the CM frame, these momenta are of order ∼ O(α 2 +v 2 rel ), which renders their scalar products with any other momenta, of higher order in α and v rel than the p 2 , q 2 corrections.

On-shell approximation
Let us now consider the case when the perturbative part of the transition amplitude C (5) consists only of fully connected diagrams. 7 Then, A (5) (P ϕ , η 1 P + p, η 2 P − p; η 1 K + q, η 2 K − q) is the perturbative amplitude for the 2 → 3 transition (with no on-shell conditions imposed). Equation (3.31) becomes M trans (q; p) dp 0 2π dq 0 2π S(p; P ) S(q; K) S 0 (p; P ) S 0 (q; K) A (5) (P ϕ , η 1 P + p, η 2 P − p; η 1 K + q, η 2 K − q). (3.35) Provided that A (5) has no singularities in p 0 and q 0 , 8 the integrations over p 0 , q 0 force the evaluation of A (5) on the poles of S(p; P ) and S(q; K) that are located in either the lower or upper p 0 and q 0 complex planes (depending on the choice of integration contours). As described in appendix C, each integration picks out two poles: One physical pole, which corresponds to setting one of the particles on-shell, and one unphysical pole, where the energy of the other particle is negative (c.f. Eqs. (C.11), (C.12)). In the non-relativistic regime, the contribution from the physical pole dominates. For concreteness, let us take these poles to be in the lower p 0 , q 0 complex planes (as in Eq. (C.12)), Fixing p 0 and q 0 to the pole values means that the energies of the χ 1 , χ 2 particles in the bound state and in the two-particle states are specified as functions of the 3-momenta p, q, P, K and the quantum numbers n and k (note that P 0 = ω P,n and K 0 = ω K,k ), as follows where we used Eqs. (2.76), (2.77) and (C.18), (C. 19). Evidently, in both the bound state and the two-particle state, the χ 2 degree of freedom is off-shell, by E n − p 2 /2µ and E k − q 2 /2µ respectively. However, p 2 /(2µ) ∼ −E n and q 2 /(2µ) ∼ E k ; provided that E n , E k µ, 9 we may ignore this small deviation from the on-shell condition and evaluate A (5) on-shell. Then, from Eq. (3.31), we obtain This is the approximation presented in Sec. 5.3 of Ref. [70]. Note that, for consistency, when using Eq. (3.38) inside Eqs. (3.29) and (3.30), the normalisation factor of Eq. (3.33) should be approximated to zeroth order in p 2 , q 2 . Indeed, corrections of the order p 2 , q 2 arise not only due to the normalisation factor, but also due to the off-shellness of the amplitude A (5) . When important, such corrections should be included self-consistently, by making use both of the full expansion of Eq. (3.33), and of the off-shell momenta of Eqs. (3.36), (3.37) instead of the on-shell conditions. We will not make use of Eq. (3.38) in our computations in Sec. 5.

Bound-state formation cross-sections
In the CM frame (K = 0), the differential cross-section for radiative BSF, U K=0,k → B P,n + ϕ −P , is where we used M n = m + E n (c.f. Eq. (2.80)). In addition,

Partial-wave decomposition and unitarity
It will be useful to decompose the amplitude M k→n in partial waves where P J are the Legendre polynomials, and with the partial-wave cross-section given by (3.46) Unitarity implies an upper limit on the partial-wave inelastic cross-sections. In the nonrelativistic regime, for the Jth partial wave, this is [71] For a given inelastic process and associated cross-section, this bound yields an estimate for the value of the coupling at which the probability for inelastic scattering saturates. In Sec. 5, we use the unitarity bound to deduce the range of validity of our calculation.

Bound-state de-excitation rates
The differential rate for the radiative de-excitation of a bound state, B K=0,n → B P,n + ϕ −P , is where |M n →n | 2 is found from Eq. (3.30) and |P| is given by Eq. (3.40) with the replacement E k → E n . Setting M n = m + E n m, we obtain

Non-perturbative amplitude
If an unbound χ 1 , χ 2 pair can (co-)annihilate into a number of light particles f 1 , · · · f N , then the χ 1 − χ 2 bound states are unstable against decay into the same final states (provided that this is allowed by angular momentum conservation). For example, unbound and bound particleantiparticle pairs can annihilate and decay, respectively, into force mediators. (Co-)annihilation and bound-state decay have the same diagrammatic representation, shown in Fig. 5; the difference A ann Figure 5. Diagrammatic representation of equation (4.5). represent the full propagators of annihilation/decay products fj.
in evaluating these two processes is which initial state is singled out from the G (4) function. In this section, we express the (co-)annihilation and the bound-state decay amplitudes, in terms of the initial state wavefunction and the perturbative interaction that gives rise to these processes. 10 To calculate the S-matrix elements (4.1) and (4.2), we need to consider the Green's function and its Fourier transform Let A ann (p 1 , · · · p N ; k 1 , k 2 ) be the sum of all connected and amputated diagrams contributing the (co-)annihilation of a χ 1 , χ 2 pair with momenta k 1 , k 2 , into f 1 , · · · , f N particles with momenta p 1 , · · · , p N . 11 Then, as sketched in Fig. 5, whereS fj (p j ) is the propagator of the f j particle with momentum p j . As always, energy-momentum conservation implies that where M pert ann is the perturbative annihilation amplitude, with no on-shell conditions imposed on the incoming and outgoing degrees of freedom. We follow a similar procedure as in Sec. 3.2, and determine the S-matrix elements of interest to be dec , (4.8) 10 As already mentioned, for annihilation processes, similar analyses have been carried out in previous works, e.g. [40-42, 65, 66, 72-74]. 11 Note that, in contrast to level transitions processes, all diagrams contributing to the perturbative part of the (co-)annihilation processes are fully connected. The amputation of fully connected diagrams is well defined, and it is thus sensible to express the full amplitudes for the (co-)annihilation and decay processes of interest, in terms of the sum of amputated diagrams.

with
where K 0 = ω K,k and K 0 = ω K,n , respectively. In the above, Z fj (p) ≡ | Ω|f j (0)|f j,p | 2 , with f j,p being a f j particle with momentum p. In the second line Eqs. (4.9) and (4.10), we have used the instantaneous approximation for the wavefunctions, and set Z fj (p) 1.

On-shell approximation
Following the same arguments as in Sec. 3.3, we may evaluate the perturbative amplitude M pert ann on-shell. This enables us to express the annihilation and decay amplitudes as follows whereM pert ann is the on-shell perturbative (co-)annihilation amplitude. Note that, as discussed below Eq. (3.38), the integrands in Eqs. (4.11), (4.12) admit q 2 and higher order corrections from the normalisation factor of Eq. (3.32) and from the off-shellness of the perturbative amplitude M pert ann .

Two-body (co-)annihilation cross-sections and bound-state decay rates
Let us now focus on the case of decays and (co-)annihilations into two final-state particles. In the CM frame (K = 0), the momenta of the final particles are |p 1 | = |p 2 | = |p|, with |p| ω K=0,k /2 = (m + E k )/2 m/2 in the case of (co-)annihilation, and |p| M n /2 = (m + E n )/2 m/2 in the case of decay. We ignore the masses of the final-state particles for simplicity. The (co-)annihilation and decay amplitudes can be expanded in partial waves as follows where M ann, ≡ dΩ p P (cos θ p )M ann (Ω p ),  Here and in the following, the indices in the angle variables specify the vector to which this angle refers; a double index denotes the angle between the two vectors. The (co-)annihilation cross-section times relative velocity and the decay rate are  In general, the expansion coefficientsã J may depend on q; in the non-relativistic regime, they can be expanded asã (q) a + F (q 2 , A · q) , (4.20) where A stands for the polarisation vectors of possible final-state vector bosons, and F is a polynomial function of the scalar products q 2 and A · q that vanishes at q = 0. Note that |p| is determined by energy conservation. a and F may depend on scalar products such as p 2 , A ·p and A · B . In the following, we consider only the a contribution toã ; any corrections arising from the q-dependent terms of Eq. (4.20) may be included only in conjunction with similar corrections arising from the normalisation factor of Eq. (3.32) and from the off-shellness of the perturbative amplitude M pert ann . We may now insert Eq. (4.19) into Eqs. (4.11) and (4.12), and use the formula , (4.21) and similarly for ψ n . We prove Eq. (4.21) in appendix D. Keeping only the q-independent term from the expansion of Eq. (4.20), we find M ann, a |p| (mµ) . (4.28) In the limit where the interaction in the two-particle state can be neglected, φ k (r) = e ik·r and S ,ann = (|k|/µ) 2 = v 2 rel . Similar analyses to the above for the non-perturbative annihilation cross-section, have been performed in Refs. [65,66], where also the Sommerfeld enhancement factors of Eq. (4.27) have been computed for a Yukawa potential.

Bound-state formation, de-excitation and decay rates for specific interactions
We now focus on specific interactions and apply the formalism of the previous sections to calculate the BSF cross-sections, and the rates for de-excitation or decay of bound states, where relevant. We consider the interaction of two scalar particles (i) via a light scalar boson (Sec. 5.1), and (ii) via an Abelian gauge vector boson (Sec. 5.2).
In the instantaneous approximation, these interactions are described in general by a Yukawa potential. 12 (Of course, an unbroken gauge symmetry gives rise to a Coulomb potential.) A Yukawa potential admits bound state solutions if m ϕ < αµ, where α is the fine-structure constant of the interaction. On the other hand, the radiative formation of bound states via emission of a force mediator is kinematically possible if m ϕ < (α 2 + v 2 rel )µ/2; for v rel < α -which is when the Sommerfeld effect renders bound-state formation efficient -this is a much stronger condition. Provided that this condition holds, the distortion of the wavefunctions due to the non-zero mediator mass, from their Coulomb limit, is expected to be negligible. For simplicity, we shall thus perform our computations in the Coulomb limit.
As is well known, in the presence of an attractive Coulomb potential there is a discrete spectrum and a continuous spectrum of energy eigenstates. The continuous spectrum corresponds to the two-particle states, and is characterised by a continuous quantum number that stands for the expectation value of the momentum of the reduced system, k = µv rel , with v rel being the expectation value of the relative velocity. The discrete spectrum corresponds to the bound states, and is characterised by the integer-valued quantum numbers {n m}. In the discrete spectrum, the expectation value of the momentum of the reduced system is κ/n, with κ ≡ µα being the Bohr momentum. As we shall see, the parameter that essentially determines the efficiency of BSF is the ratio of the momentum expectation values of the bound and the scattering states, The energies of the states of the discrete and the continuous spectra are where P, K are the momenta of the CM of the bound and the two-particle states respectively, and The wavefunctions are given in appendix F.1.

Useful integrals
For the calculation of the amplitudes M k→n and M n →n , we will find it useful to define the integrals We evaluate Ξ 1 and Ξ 2 in appendix E. We will also need the following integrals involving the initial and final state wavefunctions and We evaluate I, J and K in appendix F, for the initial and final states of interest. We shall use the integrals (5.5) -(5.12) in sections 5.1 and 5.2.

Scalar mediator
We consider the interaction Lagrangians In Eq. (5.13), χ 1 and χ 2 are real scalar fields, while in Eq. (5.14) they are complex. φ is a real scalar boson, with mass m ϕ m 1 , m 2 , and g 1 , g 2 are dimensionless couplings. 13 To lowest order, the interaction between χ 1 and χ 2 is mediated by one-boson exchange, as shown in Fig. 6. ThenW In the instantaneous approximation 14 (5.16) From Eq. (2.85), we find the non-relativistic potential, That is, The interaction is attractive if g 1 g 2 > 0, i.e. it is always attractive between particles of the same species or particles and antiparticles, but it can be either attractive or repulsive between particles of different species. As already mentioned, for our computations, we shall consider the limit m ϕ → 0.

Bound-state formation amplitudes
The lowest order contribution to the perturbative part of the radiative BSF amplitude arises from the diagrams shown in Fig. 7. In our approximation, the entire BSF amplitude corresponds to the ladder diagrams of Fig. 8. Recalling Eqs. (3.16) and (3.17), the diagrams of Fig. 7 evaluate to 14 Note that the leading-order correction to the approximation of Eq.  Figure 7. The lowest order contribution to level transition amplitudes, including bound-state formation.
where we used Eq. (F.22). We shall drop the α 2 correction in the coefficient of the I k,{100} integrals. The mediator momentum is · · · · · · + · · · · · · where θ is the angle between k and P ϕ , and R(ζ) is given in Eq. (F.35). (We emphasise that the above expression is not a consistent expansion in α, but rather only in v 2 rel ). We discern the following cases: • For a particle-antiparticle pair, or for identical particles, g 1 = g 2 = g, η 1 = η 2 = 1/2 and µ = m/4; the first term in Eq. (5.22) vanishes, and we obtain where in the second line, we decomposed the amplitude in partial waves.
• For non-degenerate particles, the first term in Eq. (5.22) dominates, and (The factor inside the square brackets becomes equal to 1 in the limit g 1 = g 2 , η 1 η 2 .)

Capture in excited state with non-zero angular momentum
Since for a pair of degenerate particles, the cross-section for radiative capture to the ground state is either v 2 rel or α 2 suppressed (as seen by comparing Eqs. (5.23) and (5.24)), we shall now calculate the amplitude for capture in the {210} state. In this case, |P ϕ | = E k − E {210} = (4 + ζ 2 )k 2 /(8µ). From Eq. (5.21), and keeping only the leading order terms of Eq. (F.42) for the I k,{210} integrals, we find The J = 0 and J = 2 contributions to the above yield the squared amplitudes (5.26) (We remind that M J is defined in Eq. (3.44), and note that the factor inside the square brackets becomes equal to 1 in the limit g 1 = g 2 , independently of η 1 , η 2 .)

Bound-state formation cross-sections and partial-wave unitarity
Combining Eq. (3.46) for the partial-wave cross-section and the amplitudes (5.23) -(5.26), we obtain • For g 1 = g 2 = g and η 1 = η 2 = 1/2 (i.e. µ = m/4), • For non-degenerate particles, (5.30) • For capture to the {210} state, It is interesting to note that in this low-velocity regime, the velocity dependence of all partial waves is the same. 15 This is in fact expected by unitarity, since ζ 1 is both the large coupling and the low-velocity limit. Indeed, the unitarity bounds on the partial-wave inelastic cross-sections, shown in Eq. (3.47), all have the same velocity dependence. They are realised when the factors to the right of the × symbols in Eqs. (5.34) -(5.38) become ≈ 1. The validity of our calculation is thus limited to at most α α uni , with the strongest bound, α uni ≈ 1, obtained from σ {100} BSF,J=0 .
· · · Figure 9. Annihilation of an unbound particle-antiparticle pair or decay of a bound particle-antiparticle pair into two force mediators. The mediator may be either a scalar or a vector boson.

De-excitation rate
The radiative capture to the {210} state is the dominant BSF process for particle-antiparticle pairs and pairs of self-conjugate particles. Moreover, for non-degenerate particles, it is slower but comparable to the capture to the ground state. Here, we shall thus compute the de-excitation rate of the {210} state. The radiative de-excitation of a bound state arises from the same diagrams as the radiative capture to a bound state, albeit for different initial and final states. In our approximation, these are the ladder diagrams shown in Fig. 8. Inserting Eq. (5.19) into Eq. (3.30), we find M n →n −m g 1 I n ,n (η 2 P ϕ ) + g 2 I n ,n (−η 1 P ϕ ) + g 1 K n ,n (η 2 P ϕ ) + g 2 K n ,n (−η 1 P ϕ ) 2mµ .
We may now compare the BSF and annihilation cross-sections,  16 We compare σ ann and σ {210} BSF in Fig. 1. Since annihilation is the dominant inelastic process, it lowers the value of α at which the unitarity bound appears to be realised, to α uni ≈ 0.54.

Particle-antiparticle bound-state decay rates
From (4.25), (4.26) and (4.28), we find that the decay of the {100} particle-antiparticle bound state into two mediators, is dominated by the s-wave contribution, For the {210} bound state, the s-wave decay mode vanishes, since ψ 210 (0) = 0. However, a nonvanishing contribution arises from the p-wave mode. From Eq. (4.26), we find σ 1 = f s |a 1 | 2 /(2 7 3πmµ) = πα 2 /(12µ 2 ). From Eq. (4.28) and the wavefunction (F.9), we obtain S 1,dec = κ 5 /(32πµ 2 ). Then The decay rates into three mediators are expected to be suppressed by one additional power of α with respect to the above. Recalling Eq. (5.41), this suggests that for the excited state {210}, the transition to the ground state is the dominant decay mode.

Vector mediator
We now consider two scalar particles χ 1 , χ 2 coupled to a gauged U (1) force. The interaction Lagrangian is where χ 1 , χ 2 are complex scalar bosons, F µν = ∂ µ ϕ ν − ∂ ν ϕ µ and D µ = ∂ µ − ic j gϕ µ , with c 1 , c 2 being the charges of χ 1 , χ 2 . 13 The one-boson exchange diagram gives In the non-relativistic regime, we shall approximate the above with In the instantaneous approximation, and setting Q 2 m 2 and That is, The interaction is attractive if c 1 c 2 < 0. 16 Note that for a Dirac fermion-antifermion pair, the annihilation into scalars is dominantly p-wave. The spinaveraged annihilation cross-section times relative velocity is σannv rel σ 1 S 1,ann , with σ 1 3πα 2 /(2µ 2 ) and S 1,ann = v 2 rel (1 + ζ 2 ) 2πζ/(1 − e −2πζ ). On the other hand, the BSF cross-sections do not depend on the spins of the interacting particles, since in the non-relativistic regime, the spin is conserved separately from the orbital angular momentum. Thus, for Dirac fermions, at ζ 1, σ {100} BSF /σann 0.1 and σ {210} BSF /σann 0.044/α 2 . This means that for α 0.2, BSF is faster than annihilation in the regime where the Sommerfeld effect is important.

Bound-state formation amplitude, cross-section and partial-wave unitarity
The perturbative part of the level-transition amplitudes is C (5) ϕ−amp = µ C µ ϕ−amp , with the lowest order contribution depicted in Fig. 7, From this, we obtain Then, M trans = µ M µ trans , where for j = 1, 2, 3, We remind that Ξ 1 , Ξ 2 are defined in Eqs. (5.5), (5.6). Using their non-relativistic approximations (E.23), (E.24), we find the amplitude of Eq. (3.29) to be M k→n = µ M µ k→n , where for j = 1, 2, 3, We may rewrite the above in terms of the I, J integrals, defined in Eqs. (5.7) and (5.8), as follows Because the vector boson ϕ µ is transverse, the µ = 0 component and the component parallel to P ϕ do not contribute to the amplitude M k→n = µ M µ k→n . Dropping those components, we obtain M j k→n →M j k→n . In the rest frame, K = 0, and for capture in the ground state {100}, using Eq. (F.40) for the J k,{100} integrals, and keeping only the leading term, we find where the sin θ factor arises from the projection of k on the plane vertical to P ϕ , and R(ζ) is defined in Eq. (F. 35). (Recall that c 1 c 2 < 0 for an attractive potential.) Note that the partial wave decomposition of the θ-dependent factor of Eq. (5.55) is sin θ = 1 − cos 2 θ = π 4 P 0 (cos θ) − 5π 2 5 P 2 (cos θ) − 3 2 π 2 8 P 4 (cos θ) + . . . . (5.56) The sum over the vector-boson polarisations is Note that for c 1 = −c 2 = 1, which includes the case of a particle-antiparticle pair, the factor in the square brackets in the above expression becomes [(η 2 c 1 − η 1 c 2 ) 2 /(−c 1 c 2 )] = 1.

Particle-antiparticle bound-state decay rate
From Eq. (4.25), we find the unpolarised decay rate of a particle-antiparticle bound state into two vector mediators to be (5.64)

Discussion
The formation of bound states affects the phenomenology of dark matter in a variety of ways.
Computing the rates for bound-state formation and other related processes is essential in calculating the cosmology of DM and accurately estimating the expected DM signals and detection prospects.
In the non-relativistic regime, the formation of bound states is enhanced by the Sommerfeld effect. The Sommerfeld effect has already been incorporated in computations of the DM annihilation rate, in the context of various theories, and has been shown to have important phenomenological implications. Besides enhancing the total DM annihilation rate, it may also modify -depending on the nature of the DM interactions -the relative strength of the various annihilation channels, thus changing the spectrum of the annihilation products [73]. Our results demonstrate that, for particleantiparticle pairs or pairs of self-conjugate particles, the radiative formation of bound states can be faster than annihilation, in the entire regime where the Sommerfeld effect is important. This suggests that bound-state formation and decay may affect the annihilation signals of symmetric thermal-relic dark matter, as well as its relic abundance, well beyond the experimental uncertainty in the DM density [43]. Bound-state dynamics should then be incorporated in any relevant analyses.
The importance of this point is underscored by present experimental results, which strongly constrain sub-TeV DM with electroweak interactions, and thus motivate investigations in the multi-TeV regime. As is well known, for symmetric (or self-conjugate) thermal-relic DM heavier than a few TeV, including WIMP DM, the Sommerfeld effect is important, both in the determination of the relic abundance and in the estimation of the expected indirect-detection signals. For the indirect detection of hidden-sector DM, the Sommerfeld effect -and therefore the formation of bound states -can be important even for lower DM masses.
Asymmetric dark matter can couple even more strongly to light force mediators than symmetric DM; indeed, in the presence of a particle-antiparticle asymmetry, the very efficient annihilation which such a coupling would imply, cannot destroy the DM relic abundance. It follows that, for a much larger range of masses, asymmetric DM may efficiently form stable bound states in the early universe. This has important implications for its phenomenology. On one hand, the formation of bound states typically curtails the DM self-interactions and hastens the kinetic decoupling of DM from dark radiation in the early universe; consequently, it regulates the potential effect of the DM dynamics on the galactic structure. On the other hand, DM may participate in a variety of radiative processes inside haloes, such as excitations and de-excitations of bound states, or outright formation of bound states. In addition, the scattering of DM on nucleons may involve a variety of interactions, including both elastic and inelastic processes. This interplay between cosmology and the fundamental interactions of the DM constituents, determines all manifestations of DM today, and can be calculated only with precise knowledge of the rates governing bound-state-related processes.
In this work, we established a field-theoretic framework for computing rates for processes involving bound states. This framework can be employed in future investigations of related effects, in a variety of theories. In particular, the computation of bound-state formation rates in theories which involve non-Abelian interactions -including the electroweak theory of the Standard Model -necessitates adopting a field-theoretic formalism. Moreover, this framework allows for systematic expansions in the interaction strength and in the momentum exchange between the interacting degrees of freedom; these higher-order corrections are important when the leading-order terms cancel, as was explicitly shown in our computations.
The significance of long-range interactions -and therefore, the importance of comprehending their implications -is affirmed by unitarity. Unitarity sets an upper bound on the partial-wave inelastic cross-section, (σ inel,J ) max , shown in Eq. (3.47). This, in turn, yields an upper bound on the mass of thermal-relic DM [43,71]. Notably, the velocity dependence of (σ inel,J ) max suggests that the unitarity bound can be realised only if the underlying interactions are long-ranged [43]. However, in the presence of long-range interactions, the formation of bound states may be the dominant inelastic process, as shown in the present work. The realisation of the unitarity bound, and its phenomenological implications, are thus largely determined by the dynamics of bound states, which should be fully incorporated in any related study. For example, the DM self-destruction in the early universe via bound-state formation and decay involves an interplay between capture, disassociation and decay processes that is absent in the case of direct annihilation into radiation [43].
Moreover, our results show that in the large-coupling (or low-velocity) regime, the dominant inelastic channel often belongs to a higher partial wave than usually assumed. Since the unitarity bound on higher partial waves is more relaxed, this implies that thermal-relic DM may be significantly heavier than previously estimated.
Lastly, using the computed bound-state-formation cross-sections, we may estimate the interaction strength for which the unitarity bound is seemingly realised. Our leading-order computations show that this is at α ∼ 0.5, i.e. well below what is often considered to be the perturbativity limit, α ∼ π or 4π.

Non-relativistic approximation
In the non-relativistic regime, P, p P 0 , m 1 , m 2 . Then Note that in the last expression, the cancellation of the mixed terms proportional to p · P, can be traced to Eq. (2.11). This reflects the fact that in the non-relativistic regime, the relative motion can be separated from the motion of the CM. For convenience, we set With these approximations, Eq. (C.14) becomes D Partial-wave analysis for (co-)annihilation and decay processes We will prove the following two relations d 3 q (2π) 3φ k (q) |q| P (cos θ q,p ) = We shall use the addition theorem of spherical harmonics, wherex,ŷ are unit vectors, and Y m are the spherical harmonics. From Eq. (D.3) and the orthonormality of Y m , it follows that We will also need the expansion where j is the spherical Bessel function, which satisfies We begin with the right side of Eq. (D.1). Using Eqs. (D.3) -(D.6) and the Fourier transform ofφ k (q), we find d dr dΩ r P (cos θ p,r ) φ k (r) This proves Eq. (D.1). Acting on the left side of Eq. (D.1) with dΩ p P (cos θ p ), we find

Non-relativistic approximation
In the following, we consider the CM frame, K = 0, as in Sec. 5. We shall evaluate Ξ 1 (q, p; K, P ) and Ξ 2 (q, p; K, P ) at next-to-leading order in the momenta q, p, applying the non-relativistic approximations of Eqs. (C.18) -(C.20), and setting, according to Eq. (C.21), where E n and E k are given by Eqs. (F.7) and (F.13). The next-to-leading order corrections in q, p become important when the leading term in a v 2 rel expansion cancels, as is the case for the interaction of two degenerate scalar particles via a scalar mediator (c.f. Sec. 5.1). Note though that we drop subleading terms in the couplings; such corrections do not change the structure of the wavefunction convolution integrals (c.f. Eqs. (5.7), (5.8)) that enter into the transition amplitudes of Sec. 5, and thus do not avert the cancellation of the leading term in the v 2 rel expansion. The same holds for P 2 corrections. In addition, as seen from Eq. (3.40), P 2 corrections are of order v 4 rel , α 4 and α 2 v 2 rel , while p 2 corrections are only of order v 2 rel . Similarly, P · p corrections are suppressed with respect to p 2 . As mentioned, we are interested, in particular, in evaluating Ξ 1 (q, p; K, P ) at q − p − η 2 P ϕ = 0 and Ξ 2 (q, p; K, P ) at q − p + η 1 P ϕ = 0. Then, according to the above, we shall keep only the p 2 corrections.

F Schrödinger wavefunctions and convolution integrals
In the following, we shall consider the attractive Coulomb potential V (r) = − α r . (F.1)