Radiative bound-state-formation cross-sections for dark matter interacting via a Yukawa potential

We calculate the cross-sections for the radiative formation of bound states by dark matter whose interactions are described in the non-relativistic regime by a Yukawa potential. These cross-sections are important for cosmological and phenomenological studies of dark matter with long-range interactions, residing in a hidden sector, as well as for TeV-scale WIMP dark matter. We provide the leading-order contributions to the cross-sections for the dominant capture processes occurring via emission of a vector or a scalar boson. We offer a detailed inspection of their features, including their velocity dependence within and outside the Coulomb regime, and their resonance structure. For pairs of annihilating particles, we compare bound-state formation with annihilation.


Introduction
In a variety of theories, motivated on theoretical and phenomenological grounds, dark matter (DM) couples directly to light or massless force mediators, which give rise to long-range interactions. A notable example is the self-interacting DM scenario, which can currently explain the observed galactic structure better than collisionless DM. Importantly, even the Weak interactions of the Standard Model -which have long served as the prototype of short-range interactions, and as the canonical particle-physics framework for DM -exhibit a long-range behaviour if the interacting particles are heavier than a few TeV.
An important implication of long-range interactions is the existence of bound states. If bound states exist in the spectrum of the theory, then they may form efficiently in the early universe, and in the dense and non-relativistic environment of haloes today. This, in turn, may have dramatic consequences for the phenomenology and the detection signatures of DM. It is thus essential to accurately compute the formation of DM bound states, when long-range interactions are considered. In theories where the (effective) particle degrees of freedom are weakly coupled, the efficiency of bound-state formation (BSF) depends on the cross-sections of the relevant processes and on the thermodynamic environment. Here we are concerned with the former.
In Ref. [1], we established a field-theoretic framework for the computation of radiative BSF cross-sections in weakly coupled theories. We expressed the amplitudes for such processes in terms of the wavefunctions of the (initial) scattering state and the (final) bound state, and an off-shell perturbative interaction involving the radiative vertex. We then reduced the fully relativistic expressions into their non-relativistic counterparts. Finally, we focused on particles interacting via a Coulomb potential, and calculated the cross-sections for BSF with the emission of a (nearly) massless vector or scalar force mediator.
In this paper, we extend these calculations to the case of a massive force mediator giving rise to a Yukawa potential, where α parametrises the interaction strength, and m ϕ is the mediator mass. While no analytical expressions for the cross-sections of interest can be derived for m ϕ > 0, our goal is to outline how to evaluate these cross-sections, and to highlight the features that are important for DM phenomenology. The paper is organised as follows. In the next section, we offer some preliminaries that will be needed in the evaluation of the radiative BSF cross-sections. In sections 3 and 4, we compute the BSF cross-sections with emission of a vector and a scalar boson, respectively. We consider particularly capture into zero and low angular momentum bound states, and inspect in detail the velocity dependence and the resonance structure of the corresponding cross-sections. For particleantiparticle pairs, we compare BSF with direct annihilation into radiation. We conclude in section 5 with a discussion of the implications of our findings.
Many of the formulae used in sections 3 and 4 are derived in the appendices. In appendix A, we discuss the scattering-state and bound-state wavefunctions in a Yukawa potential. In appendix B, we consider the integrals that convolve these wavefunctions with the radiative vertex, and which enter in the computation of the radiative BSF cross-sections. We show how to perform an expansion in powers of the coupling, for capture into bound states of any angular momentum. For capture into zero and low angular momentum bound states, we identify the leading order contributions, which we then use in our computations of sections 3 and 4. In appendix C, we derive analytical expressions in the Coulomb limit, of the convolution integrals that enter the cross-sections for capture into zero angular momentum bound states. Throughout this work we have, in fact, applied two different methods to calculate BSF cross sections. While the main text and appendices A to C focus on a coordinate-space method, we have reproduced the calculations using a momentum-space procedure that is based on methods developed by the nuclear-physics community, for few-body problems. We outline this method in appendix D.
For easy reference, we list our results and several useful formulae in table 1, with references to the corresponding equations and figures. In table 2, we summarise the notation used throughout the paper and for succinctness, we often do not define these symbols in the text.
The radiative formation of bound states by particle-antiparticle pairs interacting via a Yukawa potential has been recently considered in Refs. [2,3], where a quantum mechanical formalism has been employed from the onset. Here, we adopt the field-theoretic formalism outlined in Ref. [1], which has a direct representation in terms of Feynman diagrams. Moreover, this formalism allows for a systematic inclusion of higher order corrections, which can be particularly important when the lowest order contributions to a specific process cancel. This, in fact, occurs in the radiative capture of a pair of identical particles or a particle-antiparticle pair, via emission of a scalar force mediator. The cancellation of the lowest-order terms implies, among else, that the BSF cross-sections may be different for bosonic and fermionic pairs of particles (c.f. section 2.4). In section 4.2, we calculate the BSF cross-sections for bosonic particle-antiparticle pairs, via emission of a scalar current. This computation is complementary to Ref. [3], which has considered fermionic DM; indeed, our results show that the two cases are different.
Besides the different formalism, the present work has broader applicability, in the following ways: (i) Our formulae are valid in the entire parameter range where bound states exist, rather than in the more limited parameter space where bound states are kinematically allowed to form with emission of the same force mediator that is responsible for their existence. Our results can therefore be easily adapted to describe BSF with emission a (lighter) vector or scalar boson that is not (primarily) responsible for the Yukawa interaction (1.1) (c.f. section 2.2). This is, in fact, relevant to DM coupled to the Weak interactions of the Standard Model (WIMPs), as has been recently discussed in Ref. [4]. (ii) We consider BSF by pairs of particles that do not necessarily belong to the same species. For capture via vector current emission, the BSF cross-sections computed for particle-antiparticle pairs can be easily adapted to describe the capture of particles belonging to different species, by an appropriate use of the reduced mass of the interacting pair. However, for capture via scalar current emission, the results are markedly different depending on whether the two interacting particles have the same or different masses and couplings to the scalar current. Lastly, we note that in our computations, we employ a minimal parametrisation that we believe facilitates phenomenological studies (c.f. section 2.1); using this parametrisation, we perform a rather thorough inspection of the features of the BSF cross-sections. Capture into = 0 states (including Coulomb limit) (4.5), (4.6) 7 Capture into = 1 states (4.7) 8 Scalar force mediator: Non-self-conjugate bosonic particle-antiparticle pairs Capture into = 0 states (4.10), (4.11) 9 -10 − Coulomb limit and comparison with annihilation (4.11), (4.16) 9 − Resonant structure 10a − Velocity dependence off-and on-resonance 10b Capture into = 1 states (4.12) 11 Scalar force mediator: Identical particles (4.15) -(4.16) 9 -11 Convolution integrals I k,n m (b), J k,n m (b), K k,n m (b) Equations Binding energy of n m bound state Kinetic energy of scattering state Phase-space suppression due to emission of a massive force mediator, during capture into the n m bound states pss n , see Eq. (2.5b) Dimensionless parameters

Preliminaries
We shall consider two particles X 1 and X 2 that may in general belong to different species, and interact via a vector or scalar force mediator ϕ (c.f. fig. 1a). The interaction Lagrangians will be specified in the following sections. In the non-relativistic regime, the interaction between X 1 and X 2 is described by the static Yukawa potential of eq. (1.1), which admits bound state solutions roughly if the mediator mass is less than the inverse Bohr radius, m ϕ µα ≡ κ, where µ is the X 1 − X 2 reduced mass. 1 Under this condition, we want to compute the cross-sections for the radiative formation of bound states, with emission of the force mediator ϕ, Here, U k (X 1 + X 2 ) stands for a two-particle scattering (unbound) state, characterised by the continuous vector quantum number k = µv rel , with v rel being the expectation value of the relative velocity. Because of their long-range interaction, X 1 and X 2 cannot be approximated by plane waves. This gives rise to the well-known Sommerfeld effect [5]. The U k (X 1 + X 2 ) state is instead described by a wavefunction φ k (r) that obeys the Schrödinger equation with the Yukawa potential of eq. (1.1) and a positive energy eigenvalue E k , parametrised by k. Further, B n m (X 1 X 2 ) is a bound state, whose wavefunction ψ n m (r), parametrised by the familiar principal and angular-momentum discrete quantum numbers {n, , m}, obeys the Schrödinger equation with the same potential and a negative energy eigenvalue E n . We discuss the numerical computation of ψ n m (r) and φ k (r) in appendix A. As is well-known, the Yukawa potential arises from the one-boson-exchange diagram, shown in fig. 1a (left). This is the lowest-order 2-particle-irreducible diagram contributing to the 4-point Green's function of the X 1 −X 2 pair. The (infinite) repetition of one-boson-exchange diagrams gives rise to the ladder diagrams of fig. 1a (right), whose resummation amounts to solving the Schrödinger equation. Indeed, the wavefunctions obeying the Schrödinger equation appear as multiplicative factors at (and determine the strength of) the singularities of the 4-point function: the poles corresponding to bound states, and the branch cuts corresponding to scattering states (see e.g. [1]).
It is now reasonable to wonder how do non-perturbative effects -the Sommerfeld effect and the existence of bound states -arise from the resummation of a perturbative Dyson series? Indeed, in the ladder diagrams of fig. 1a, the two vertices introduced by each boson exchange imply a suppression by one power of the coupling α. However, this suppression is cancelled by the loop momentum exchange, which also scales with α. In particular, the average momentum exchange along each virtual boson scales as |q| ∼ µα, and the off-shellness of the X 1 , X 2 propagators scales as q 0 ∝ |q| 2 ∝ α 2 ; the integration over the loop energy and momentum also yield factors of α 2 and α, respectively. It is then straightforward to see that each additional loop does not increase the order of the diagram (see e.g. [6]); instead, the ladder diagrams add up coherently.
The radiative BSF process (2.1) arises from the diagrams of fig. 1b. The initial-state ladder corresponds to a scattering state and is evaluated at center-of-momentum (CM) energy E = m 1 + m 2 + E k m 1 + m 2 , while the final-state ladder corresponds to a bound state and is evaluated at CM energy E = m 1 + m 2 + E n < m 1 + m 2 . For concreteness, we assume that the energy difference is dissipated via emission of the same particle that is responsible for the attractive interaction. However, it is straightforward to adapt our results to BSF occurring via radiation of any vector or scalar current that couples either to one or both of the interacting particles (see also section 2.2).
. Left: In the weak-coupling regime, the one-boson (either scalar or vector) exchange is the dominant contribution to the non-relativistic X1 − X2 potential. Right: The ladder diagrams -arising from the infinite repetition of the one-boson-exchange diagrams -give rise to non-perturbative effects: the Sommerfeld effect and the existence of bound states.
. The diagrams contributing in leading order, to the formation of bound states with emission of a force mediator. For this transition, the ladder to the left of the radiative vertex corresponds to the initial scattering state, while the ladder to the right corresponds to a bound state. The mediator can be either a scalar or a vector boson.

Parametrisation
We formulate our computations and present our results in terms of two dimensionless parameters, where v rel is the relative velocity of the interacting particles and µ is their reduced mass. ζ and ξ suffice to characterise the solutions of the Schrödinger equation for the Yukawa potential (1.1), for bound and scattering states (c.f. appendix A), and to parametrise the non-analytical parts of our computations that arise in the BSF cross-sections. The ζ parameter compares the average momentum transfer ∼ µv rel between two unbound particles, with the momentum transfer ∼ µα that occurs in the exchange of virtual force mediators, as explained above. In practice, ζ parametrises the velocity dependence of the cross-sections.
On the other hand, ξ parametrises the model dependence of the cross-sections. It compares the two physical scales involved: the Bohr momentum κ ≡ µα, which determines the size of the bound states and the momentum transfer along the virtual force mediators in the ladder diagrams, with the mediator mass m ϕ , which determines the range of the interaction. The interaction manifests as long-range roughly if ξ 1; this is the regime where non-perturbative phenomena, such as the Sommerfeld effect and the existence of bound states, emerge.
We believe that this minimal parametrisation, in terms of ζ and ξ, exposes the physical significance of the features of the cross-sections, and can greatly facilitate phenomenological studies.

The range of the ξ parameter
Bound states can form with emission of the force mediator that is responsible for their existence, only if the available energy from the transition to a lower energy state suffices, m ϕ < E k −E n , where E k = k 2 /(2µ) is the kinetic energy of the scattering state in the CM frame, and E n = −γ 2 n ×κ 2 /(2µ) is the binding energy of the bound state [c.f. eqs. (A.9) and (A. 13)]. In most phenomenological applications related to DM, the formation of bound states becomes important in the regime where the kinetic energy is lower than the binding energy, E k |E n | (which roughly implies v rel α/n). Then, the condition for BSF with emission of a force mediator becomes roughly m ϕ µα 2 /(2n 2 ), or, in terms of the ξ parameter, (Note that γ n (ξ) 1/n, with the equality realised in the Coulomb limit.) The condition (2.3) is significantly more stringent than the condition for the existence of bound states, γ n (ξ) > 0, which implies roughly ξ = µα/m ϕ n 2 (c.f. appendix A.2). In fact, the condition (2.3) asserts that in its regime of validity, the bound-state wavefunctions can be well approximated by their Coulomb limit 2 . This is not necessarily the case for the scattering state wavefunctions though, which depend on both ξ and ζ. As we shall see in the following, the Coulomb limit, which formally corresponds to ξ → ∞, is essentially attained for ξ ζ, (2.4) i.e. when the average momentum transfer between the interacting particles exceeds the mediator mass, µv rel m ϕ . However, in many models of interest, BSF may occur via emission of a lighter species than the force mediator that is (primarily) responsible for their existence. This species may couple only to one of the particles participating in the bound state, and would therefore not mediate a long-range interaction between them. It is also possible that a light bosonic species couples to both of the interacting particles, albeit more weakly than the heavier force mediator, whose contribution then dominates the interaction between the two particles. In either case, BSF can occur for values of ξ that do not satisfy the condition (2.3). In the following, we shall thus consider the entire range of ξ values for which bound states exist in the spectrum of a theory, even if they cannot form with emission of the force mediator responsible for the potential (1.1). To the extent possible, we will separate the phase-space suppression due to the emission of a massive particle [c.f. eq. (2.5b)], from the effect of the non-zero mediator mass on the bound-state and scattering-state wavefunctions, and through them, on the amplitude of the process. This renders it possible to adapt our results for BSF processes that occur via emission of any vector or scalar boson.

The cross-section
The differential cross-section times relative velocity for the 2-to-2 process where M k→n m is the transition amplitude, and |P ϕ | is the momentum of the emitted mediator.
In the CM frame, the energy dissipated during BSF is the sum of the kinetic and binding energies, which implies where we used eqs. (A.9) and (A.13). Then where pss n is the phase-space suppression factor for the emission of a force mediator, during capture into the n m bound state; in the Coulomb limit ξ → ∞ and pss n = 1. Putting everything together, we obtain In the following, we take the solid angle Ω (thus the angle θ that will appear in the transition amplitudes) to be the angle between the vectors k and P ϕ .

Initial and final state convolution integrals
The transition amplitudes M k→n m depend on the scattering-state and bound-state wavefunctions, ψ n m and φ k , and on the radiative vertex, via the convolution integrals [1] Hereψ n m andφ k are the Fourier transforms of ψ n m and φ k respectively. The momentum b is proportional to the momentum of the radiated particle, b ∝ P ϕ , as we shall see in the following.
In appendix B, we expand the integrals (2.7) in powers of |b|/κ ∝ |P ϕ |/κ ∝ α [c.f. eq. (2.5)], and identify the leading-order contributions for the transitions of interest. The range of validity of this approximation is given by the condition (B.3). We note that the formalism developed in Ref. [1] assumed that the interacting particles have zero spin. However, the lowest-order computations using the expressions of Ref. [1] for the radiative BSF cross-sections, are applicable to both fermionic and bosonic interacting species. Indeed, in the non-relativistic regime, and to lowest order in the coupling, the spin of each of the interacting particles is conserved in the capture process, and the BSF cross-sections do not depend on the spin (or the spin configuration) of the incoming particles. The computations of sections 3 and 4.1 fall in this category. However, in the radiative capture of a particle-antiparticle pair via emission of scalar force mediator, the lowest-order contributions cancel each other. As discussed in section 4.2, in this case, we are forced to consider higher-order terms, which may in general be different for fermions and bosons.

Radiative capture into a bound state
We assume X 1 and X 2 to be coupled to a gauged U(1) force, where F µν = ∂ µ ϕ ν −∂ ν ϕ µ and D µ = ∂ µ −ic j gϕ µ , with c 1 , c 2 being the charges of X 1 , X 2 . The mass m ϕ of ϕ µ may have arisen either via the Higgs or the Stückelberg mechanisms. While the details of the local U(1) breaking are not important for our purposes, we do assume in the Lagrangians (3.1) that the individual global U(1) symmetries associated with X 1 and X 2 remain unbroken. 3 In the non-relativistic regime, the ϕ µ exchange between X 1 and X 2 gives rise to the Yukawa potential of eq. (1.1), with The interaction is attractive if c 1 c 2 < 0.
Since BSF involves gauge interactions with conserved currents, as seen from eqs. (3.1), the Ward identity ensures that P µ ϕ M µ = 0. This implies M 0 = P j ϕ M j /P 0 ϕ . Then, the unpolarised amplitude is We thus need only the spatial components of the BSF amplitude, which are [1] M j k→n m = −2g 2µ where I k,n m and J j k,n m are defined in eqs. (2.7). It is evident from the momentum factor in the integrand of J j k,n m [c.f. eq. (2.7b)] and the momentum dependence of the I k,n m contribution to M k→n m , that eq. (3.4) describes a radiative process that proceeds via a derivative current interaction. Moreover, both charged particles contribute to the emitted radiation.

Capture into = 0 bound states
From eqs. (B.9), we find that the lowest-order contributions to the amplitude (3.4) arise from the J k,n00 integrals and are of zeroth order in |P ϕ |/κ. The unpolarised squared amplitude (3.3) is (Since we are now interested only in the zero-th order terms in |P ϕ |/κ, we have dropped the arguments of the J k,n00 functions.) Substituting eq. (B.9b) into eq. (3.5) yields where χ n, (x) and χ |k|, (x) are related to the bound-state and scattering-state wavefunctions as described in appendix B. From eq. (2.6) and (3.6), we find the BSF cross-section to be 4 (3.7b)

Coulomb limit
At m ϕ → 0, the above becomes where, using eq. (C.5b), we find Expanding the sum, S {n00} BSF,C takes the form where n (ζ) is a rational function of ζ 2 (and specifically, a polynomial of degree n − 2, for n 2), with lim ζ→0 n (ζ) = 1. We give the explicit expressions of n (ζ) for 1 n 5 in table 3. For capture into the ground state in particular, We illustrate in detail the features of S {n00} BSF (ζ, ξ) and S {n00} BSF,C (ζ) in figs. 2 to 5, and we discuss them in section 3.5. 4 Here and in the following, we typically choose to factorise the various BSF cross-sections (i.e. separate out the S BSF factors) in a way that facilitates the comparison between BSF and annihilation processes for particleantiparticle pairs. However, we note that the S BSF factors do not carry any physical significance on their own. This is in contrast to the Sann factors appearing in the various annihilation cross-sections [c.f. sections 3.4 and 4.2.3], which represent the enhancement of the corresponding processes due to the non-perturbative Sommerfeld effect, and are often referred to as "Sommerfeld enhancement factors". At ζ, ξ 1, the non-perturbative effects switch off, and the annihilation cross-sections reduce to their perturbative values (Sann 1). In the same regime, S BSF → 0, as the very existence of bound states is a non-perturbative effect. The BSF processes do not have a perturbative limit, and the S BSF factors are not enhancement factors of otherwise perturbative processes, by a non-perturbative effect. n n (ζ) Table 3.
The functions n(ζ ) appearing in the Coulomb limit of the cross-sections for the radiative formation of zero angular momentum states of principal quantum number n, Let us now recount the origin of the various factors in eqs. (3.7a) and (3.8c). In eq. (3.7a), the factor ( emanates from the vertices of the radiated gauge boson in the Feynman diagrams of fig. 1b, and asserts that lighter particles radiate more easily; it becomes 1 for c 1 = −c 2 . The phase-space suppression factor pss 1/2 n,0 is due to the limited energy available for the radiation of a massive force mediator. The factor (3 − pss n,0 ) accounts for the contribution of the three polarizations of the emitted massive gauge boson; it reduces to 2 for a massless mediator.
The function S BSF (ζ, ξ) measures the overlap of the initial and final states. Its Coulomb limit (3.8c) is illustrative. The first factor, S 0 (ζ) = 2πζ/(1−e −2πζ ), is an overall multiplicative constant in the scattering-state wavefunction [c.f. eqs. (C.1)], and is responsible for the characteristic scaling of long-range inelastic processes, σ BSF v rel ∝ 1/v rel at ζ 1. The other factors in eq. (3.8c) depend on the details (space/momentum dependence) of the scattering-state and bound-state wavefunctions, as well as the radiative vertex. At ζ < 1, S BSF (ζ, ξ) is very small and BSF is typically inefficient; however, at ζ 1, S BSF,C (ζ) 3.13 × 2πζ. This point will become important in section 3.4, where we compare BSF and annihilation for particle-antiparticle pairs.
The BSF cross-sections we calculate in the following sections have similar structure to the one described above.

Capture into = 1 bound states
Similarly to before, the zeroth order terms in |P ϕ |/κ, of the J k,n1m integrals yield the dominant contribution to the amplitude (3.4) for capture into = 1 bound states, with the unpolarised squared amplitude (3.3) being For convenience, we define Then, from eqs. (B.10) we find and To calculate the total cross-section for capture to any = 1 state, we sum over all values of m.
Note thatk ·P ϕ = cos θ k cos θ Pϕ + sin θ k sin θ Pϕ cos(φ k − φ Pϕ ), and the term proportional to Im(A 2 A * 0 ) in eq. (3.12) gives a vanishing contribution when integrated over dΩ Pϕ . Using the above and eq. (2.6), we find the total cross-section for capture into an = 1 state (for fixed n), with A 0 and A 2 defined in eqs. (3.10).

Annihilation of particle-antiparticle pairs
Bound states of particle-antiparticle pairs are unstable and decay into radiation. This effectively provides an extra annihilation channel. It is then instructive to compare BSF with the direct annihilation into radiation. The (spin-averaged) annihilation cross-section times relative velocity of a particle-antiparticle pair of scalars and fermions respectively, is where µ = m X /2 = M/4 and Note that the perturbative value of the annihilation cross-section (recovered in the limit S (0) ann → 1) is different for bosons and fermions. In contrast, the leading order BSF cross-section does not depend on the spin of the particles, as already discussed in section 2.4.
We compare annihilation and BSF in figs. 2 and 3, and offer our comments in the next section.

Discussion
We now discuss the main features of the BSF cross-sections computed in this section.

Coulomb limit
For a massless force mediator, and at large enough ζ, the BSF cross-sections scale as σv rel ∝ ζ. This behaviour is realised at ζ n for capture into an {n m} bound state, as evident in fig. 2; for capture into = 0 bound states, it can also be confirmed analytically using eq. (3.8c). Importantly, the 1/v rel scaling is anticipated due to the upper bound on the inelastic cross-sections imposed by unitarity [1,7], which we discuss below. The radiative capture into the ground state dominates over capture into excited levels, n > 1. On the other hand, for n > 1, the formation of = 0 states (summed over − m ) dominates over capture into the = 0 state, in part due to the larger multiplicity of the former.
The above points imply that for a given value of ζ, the total BSF cross-section is dominated by capture into the n < ζ levels, with the > 0 states yielding a significant contribution. Indeed, the capture into excited states gives rise to a logarithmic enhancement of the total cross-section for radiative BSF [2,8], described by Kramer's formula,  .15) is not always relevant for the phenomenology of DM, since the observable implications of bound states do not depend only on the total rate at which bound states form, but also on their features. For example, the formation of unstable (particle-antiparticle) bound states in the early universe, and their subsequent decay into radiation, can reduce the DM relic density [7]; however, the efficiency with which DM is depleted depends on the balance between BSF, ionisation and decay, which in turn depends sensitively on the quantum numbers of the bound states that form. Similarly, the cosmological formation of stable bound states by asymmetric DM [9], and their survival until today, depends typically on a rather complex interplay between formation, ionisation, excitation and de-excitation processes [10]. Moreover, the detectability of the radiation emitted inside halos today during the formation of stable bound states of asymmetric DM [11][12][13], obviously depends on the energy release in the specific transition that takes place. In all these cases, the phenomenological importance of BSF cannot be assessed based solely on eq. (3.15), even in the Coulomb regime. On the other hand, the high-energy signals arising from the decay of unstable bound states of symmetric DM inside halos may reflect the logarithmic enhancement of eq. (3.15), provided that the excited states get de-excited or decay into radiation promptly enough in astrophysical timescales [2]. Note though that away from the Coulomb regime, this enhancement is curtailed due the phase-space suppression [c.f. eq. (2.5b)] that becomes more severe for capture into excited states [2].
The Coulomb limit is attained at ξ ζ, as can be observed in figs. 3, 4 and 6 and figs. 5, and was already noted in section 2.1. The physical interpretation of this condition is that the momentum transfer between the two incoming particles exceeds the mediator mass, µv rel m ϕ .

Resonance structure
The BSF cross-sections exhibit a rich resonance structure, as seen in figs. 3a, 4 and 6. The resonances are features of the scattering-state wavefunction, which determines the strength of the branch-cut singularity of the X 1 − X 2 4-point function (see e.g. [1]). The resonances arise at the points of the parameter space where a pole of the 4-point function lives at zero energy, E n → 0, and thus overlaps with the branch-cut, which lives at E k 0 (see e.g. Ref. [14, section 7.7]). The locations of resonances therefore denote the thresholds for the existence of genuine bound states with E n < 0. These thresholds signify that if the mediator is massive, the potential has to be sufficiently strong for bound states to exist. In contrast, for a massless mediator, bound states exist independently of the strength of the coupling, as is the case with QED. For every bound-state energy level, this implies a minimum value for ξ ≡ µα/m ϕ .
Unlike the Coulomb potential, the Yukawa potential does not imply the conservation of the Laplace-Runge-Lenz vector. The energy eigenvalues of the discrete spectrum depend on both the n and quantum numbers (c.f. appendix A.2). As a result, the resonances that appear in the scattering-state wavefunction depend on the mode.
Angular momentum conservation implies that the modes of the scattering-state wavefunction participating in a capture process depend on the angular momentum of the bound state formed, the spin of the emitted particle, and the orbital angular momentum of the final state (c.
ann , the Sommerfeld enhancement factor of p-wave annihilation processes, As seen in eq. (3.16), S ann depends only on the = 1 mode of the scattering-state wavefunction [see eq.  fig. 5a attests that the resonances in the BSF crosssections emanate from the scattering state. However, the capture to the {200} bound state exhibits also anti-resonances (c.f. fig. 4), whose origin is clearly not the scattering state alone. Instead, the anti-resonances arise from the convolution of the scattering-state wavefunction with the {200} bound-state wavefunction that contains a node.

Velocity dependence away from the Coulomb regime
At low velocities, v rel < m ϕ /µ, the inelastic cross-sections depart from the Coulombic scaling σ inel v rel ∝ 1/v rel . Their velocity dependence in this regime, is determined by the modes of the scattering-state wavefunction that participate in the corresponding processes.
For ξ values away from resonances, the contribution of an mode of the scattering-state wavefunction to σ inel v rel scales as v 2 rel at v rel < m ϕ /µ. That is, as the velocity decreases, a contribution from the = 0 mode saturates to its Coulomb value at v rel ≈ m ϕ /µ, while contributions from > 0 modes drop below their Coulomb value at v rel ≈ m ϕ /µ. The radiative capture into {n00} bound states is dominated by the = 1 mode of the scattering-state wavefunction, and scales as v 2 rel at v rel < m ϕ /µ, as can been seen in figs. 3b and 4 (bottom left panel).
For ξ values near or on resonances, σ inel v rel grows faster than 1/v rel at v rel < m ϕ /µ, and raises above its Coulomb value for a given velocity. This growth is sustained for a range of velocities that depends on how close ξ is to a resonance value. Then, at sufficiently low velocities, the resonant growth stops, and the v 2 rel scaling ensues. 5 This behaviour can be observed in figs. 3c and 4 (bottom right panel).
If more than one modes of the scattering-state wavefunction participate in a process (at the same order in the coupling), then more complex patterns arise. This is the case with the radiative capture to = 1 bound states, which receives contributions from the = 0 and = 2 modes of the scattering-state wavefunction (c.f. fig. 6).
Because of the different location of the resonances, as well as the different velocity dependence on-and off-resonance that the various BSF cross-sections exhibit, the relative strength of these processes at low velocities can be very different than in the Coulomb regime. In fact, at sufficiently low velocities, the capture to = 1 bound states, if kinematically allowed, should always dominate, since these are the only mono-photon transitions to which the = 0 mode of the scattering-state wavefunction participates.

Near-threshold behaviour
As mentioned above, ξ values near the thresholds for the existence of bound states imply resonances in the scattering-state wavefunction. These resonances then appear in the cross-sections of the processes in which this scattering state participates. However, the radiative formation of a bound state for ξ values close to the threshold for the existence of the same bound state, is suppressed [c.f. figs. 3d and 4 (top right panel)]. In this limit, the bound state wavefunction becomes very spatially extended, and approaches zero.
Note that this suppression is independent of the phase-space suppression due to the emission of a massive mediator. In fact, radiative BSF near threshold is kinematically allowed only if it occurs with emission of a nearly massless particle that therefore cannot be the force mediator itself.
Comparison of bound-state formation and annihilation, for particle-antiparticle pairs At ζ 1, the BSF processes are rather suppressed. However, at ζ 1, i.e. in the regime where the Sommerfeld effect is important, BSF can be comparable to and even more significant than annihilation [7].
In the Coulomb regime, we may easily compare the radiative capture to the ground state with the direct annihilation into force mediators, using eqs. (3.8) and (3.14), (3.17) At ζ 1, this becomes σ {100} BSF,C /σ f ann,C = 2σ {100} BSF,C /σ s ann,C 3.13. In fact, the radiative formation of n = 2, = 1 bound states is also faster than the annihilation of a fermion-antifermion pair, for ζ 4. Figure 2 compares the dominant BSF processes with annihilation, in the Coulomb limit.
Since BSF is the dominant inelastic process for particle-antiparticle pairs, whose bound states are unstable and decay into radiation, the BSF via emission of a vector boson can significantly affect the relic density 6 [7] and enhance the indirect detection signals [2,15,16] of symmetric DM. Away from the Coulomb regime, the comparison of annihilation and BSF becomes more complex. Figures 3 illustrate the main features, which are related to the discussion offered above, and which we now summarise: • The formation of zero-angular-momentum bound states and the direct annihilation into radiation exhibit resonances at different locations (ξ values), as clearly seen in fig. 3a. This is due to the different modes of the scattering-state wavefunction that contribute to each process: = 1 for the former and = 0 for the latter. In fact, the locations of the = 1 resonances exhibit a mild velocity dependence, in contrast to the = 0 resonances.
• The modes of the scattering-state wavefunction also determine the velocity dependence away from the Coulomb limit. At low velocities, σ ann v rel saturates to a constant value, while σ {100} BSF v rel scales as v 2 rel . This scaling, together with the relative strength of the two processes in the Coulomb regime that we discussed above, imply that the radiative capture to the ground state dominates over annihilation within a finite range of velocities, as can be seen in figs. 3b and 3c. This range is roughly 1 ζ ξ (or equivalently m ϕ µv rel µα) for non-resonant ξ values, but it may extend to much lower velocities (by orders of magnitude) for ξ values near an = 1 resonance.
• Near the threshold for the existence of the ground state, i.e. for ξ ≈ 1, the radiative BSF is always suppressed with respect to its Coulomb value, while annihilation is on resonance. This behaviour can be seen in fig. 3d. (We repeat that in this regime, BSF is kinematically allowed to occur only via emission a nearly massless particle, which cannot therefore be the force mediator itself.)

Unitarity
Unitarity implies an upper bound on the partial-wave inelastic cross-sections, which in the nonrelativistic regime is [17], where J is the partial wave.
In the Coulomb limit, inelastic cross-sections have the same velocity scaling as eq. (3.18); setting σ inel σ uni then implies an upper bound on α that does not depend on any other physical parameter, and is roughly α 0.5 [1,7]. Around this upper bound on α, higher order corrections (of perturbative or non-perturbative origin) need to be considered. Notably, α ∼ 0.5 is well below the naive perturbativity limit, α ∼ 4π.
However, the resonances that appear away from the Coulomb limit, imply that the leadingorder computations presented here, may violate the unitarity bound even at much lower values of α that now also depend on ζ and ξ. Evidently, the peaks of the resonances can be unphysical. 6 Reference [2] argued that BSF cannot affect the DM relic density due to the rapid ionization of the bound states. However, ionisation was fully taken into account in a proper analysis in Ref. [7], which employed a set of coupled Boltzmann equations that incorporate bound-state formation, ionisation and decay processes. This analysis showed that, in a dark QED scenario, BSF reduces the DM density by a factor greater than 2 if the DM mass is m DM 15 TeV (and up to factor of 4 at the unitarity limit, m DM ≈ 140 TeV). As pointed out in [7], BSF depletes efficiently the DM density only after the ionisation rate drops below the decay rate of the bound states. The detailed timeline shows that this occurs around or before freeze-out for m DM 20 TeV. Note that Sommerfeld-enhanced processes -either annihilations or BSF -remain important even after the DM freeze-out (conventionally defined as the time of departure of the DM density from its equilibrium value). This explains the significant effect of BSF on the relic density even for m DM 20 TeV.
The resonances can be regulated by taking into account the short-range elastic scattering of the interacting particles (see e.g. [18,19]). For particle-antiparticle pairs, the short-range inelastic scattering -in particular the annihilation processes -may also contribute to taming the unphysical behaviour [20].  Figure 2. Velocity dependence of the cross-sections for the radiative formation of bound states, by a pair of particles charged under an unbroken dark U(1) force. The solid lines correspond to capture into zero angular momentum bound states, with n denoting the principal quantum number. The dashed line is the total cross-section for capture into any n = 2, = 1 state (i.e. summed over all possible projections of the bound-state angular momentum on the z axis). We also show the cross-sections for the annihilation of a particle-antiparticle pair of fermions and scalars (red dashed lines). All cross-section have been normalised to σ f 0 ≡ πα 2 /(4µ 2 ), where µ is the reduced mass of the interacting particles. σ f 0 is the spin-averaged perturbative annihilation cross-section times relative velocity of a Dirac fermion-antifermion pair; for a complex scalar, this quantity is σ s 0 = 2σ f 0 . The bound-state formation cross-sections do not depend on the spin of the incoming particles, which is conserved in the non-relativistic regime (at leading order), as particles get captured in a bound state. For particle-antiparticle pairs, at ζ = α/v rel 1, the radiative capture into the ground state is the dominant inelastic process, . Resonances due to bound states at threshold (zero binding energy) appear at discrete values of ξ = µα/mϕ, with a magnitude that increases at low energies (large ζ = α/v rel ). They arise from the wavefunction of the (initial) two-particle scattering state, and depend on the modes that contribute to the process of interest.
Top left: The factor S {100} BSF (ζ, ξ), appearing in the cross-section for capture into the ground state, with emission of a vector force mediator. The = 1 mode of the scattering state wavefunction yields the dominant contribution to this process, and results in sharp resonances whose precise location has a mild ζ ann, for non-resonant values of ξ = µα/mϕ. When the momentum transfer between the interacting particles drops below the mediator mass, µv rel mϕ (equivalently, ζ ξ), the velocity dependence of the cross-sections departs from the Coulombic behaviour.
ann, for values of ξ = µα/mϕ near the n = 2 resonances. As in the non-resonant case, the velocity dependence of the cross-sections departs from the Coulombic behaviour when the momentum transfer between the interacting particles drops below the mediator mass, µv rel mϕ (ζ ξ). In contrast to non-resonant ξ values though, for ξ near a resonance and at ζ ξ, the interaction cross-sections become larger than their Coulomb values at the same velocity, within a finite range of velocities.
Top left: The cross-section for radiative capture into the ground state with emission of a vector mediator rises above its Coulomb limit, and then drops approximately as ζ −2 ∝ v 2 rel at ζ ξ (see also fig. 5b).
Top right: For s-wave annihilation processes, the Sommerfeld enhancement rises monotonically with ζ; at ζ ξ, it asymptotes to a value that can be much larger than its value at ζ ≈ ξ.  Vector mediator: Formation of zero angular momentum bound states vs. p-wave annihilation    Figure 5b.
Comparison of the velocity dependence of radiative capture processes with emission of a vector force mediator, to p-wave annihilation processes. The factor ζ −2 S ann] tends to a ξ-dependent constant, which itself saturates to a fixed value (its Coulomb limit) at large ξ. For the ground state and first excited state, Right: At large ζ, the contribution from the = 2 mode decreases with decreasing velocity as v 4 rel , while the contribution from the = 0 mode saturates to a constant value and dominates.
The superposition of two modes -which have different velocity dependence and resonances -gives rise to the various features of the total cross-section.

Scalar force mediator
We now consider the interaction Lagrangians In eqs. (4.1a), (4.1b), and (4.1c), X 1 and X 2 are real scalar fields, complex scalar fields, and Dirac fermions, respectively. ϕ is a real scalar boson, and g 1 , g 2 are dimensionless couplings. The interaction between X 1 , X 2 via ϕ exchange is described by the Yukawa potential of eq. (1.1), with α = α sc or α = α f depending on whether the interacting particles are scalars or fermions respectively, where [1] α sc = g 1 g 2 16π and The interaction is attractive if g 1 g 2 > 0. In the following, the parameters ζ ≡ α/v rel and ξ = αµ/m ϕ are always defined using the appropriate α.

Pairs of non-degenerate particles
The BSF amplitude is [1] M k→n m −M 2µ [g 1 I k,n m (η 2 P ϕ ) + g 2 I k,n m (−η 1 P ϕ )] . The dominant contributions to the amplitudes M k→n00 and M k→n1m that we will consider below, are of order |P ϕ |/κ.

Capture into = 1 bound states
Similarly to above, we combine eqs. (2.6), (4.3) and (B.5a), and find the BSF cross-section to be We inspect the features of S {21} BSF (ζ, ξ) in fig. 8.  Right: At large ζ, the contribution from the = 2 mode decreases with decreasing velocity as v 4 rel , while the contribution from the = 0 mode saturates to a constant value and dominates.

Bosonic particle-antiparticle pairs (non-identical)
We now consider particle-antiparticle pairs of a non-self-conjugate species X. In this case, g 1 = g 2 = g, η 1 = η 2 = 1/2 and µ = M/4. As seen in eqs. (4.5) and (4.7), the lowest-order terms cancel, and we are forced to consider the next order contributions. For this reason, the computations that follow, which are based on the formalism of Ref. [1] that was developed for bosonic species, will be valid for bosonic particle-antiparticle pairs only. A computation of BSF cross-sections for fermion-antifermion pairs with emission of a scalar force mediator, can be found in Ref. [3].
In the radiative BSF amplitude, we shall now include higher-order contributions from (i) the relativistic normalisation of states, and (ii) the off-shellness of the incoming and outgoing fields in the perturbative part of the amplitude that includes the radiative vertex (i.e. the part of the diagrams of fig. 1b that remains when the incoming and outgoing ladders are amputated). 7 Then, the BSF amplitude for a particle-antiparticle pair of non-self-conjugate bosons is [1] M k→n m −M 2µ [g 1 I k,n m (η 2 P ϕ ) + g 2 I k,n m (−η 1 P ϕ )+

Capture into = 0 bound states
Keeping the terms of order (|P ϕ |/κ) 2 from the I k,n00 integral [c.f. eq. (B.9a)] and the zeroth-order terms from the K k,n00 integral [c.f. eq. (B.9c)], the amplitude (4.8) becomes Note that here, M k→n00 receives its dominant contribution from the = 0 and = 2 components of the scattering state wavefunction. This is in contrast to the formation of zero angular momentum bound states via vector emission [c.f. eq. (3.6)], or via scalar emission but by non-degenerate particles [c.f. eq. (4.4)], where the = 1 mode of the scattering-state wavefunction dominates. Substituting eq. (4.9) into eq. (2.6), we find the corresponding cross-section, where S {n00} BSF includes terms suppressed by different powers of pss n,0 , (4.10e) Note that the dependence of S {n00} BSF (ζ, ξ; α) on α (independently of ζ and ξ) arises from the phasespace suppression factor pss n,0 in eq. (4.10b).

Capture into = 1 bound states
Similarly to the previous section, we keep the terms of order (|P ϕ |/κ) 2 from the I k,n1m integrals and the zeroth-order terms from the K k,n1m integrals [c.f. eqs. (B.5)]. Then, from eqs. (2.6) and (4.8), we find the total cross-section for capture to any = 1 state, for a fixed n, to be 8 n,1 , (4.12a) 8 We note that the computation of σ where S {n1} BSF includes terms suppressed by different powers of pss n,1 , We showcase the resonant features of S {21} BSF (ζ, ξ; α) in fig. 11.

Annihilation
We now consider the annihilation of a bosonic particle-antiparticle pair into two scalar force mediators, χ * + χ → ϕ + ϕ. The dominant contribution to the annihilation cross-section arises from the = 0 component of the scattering state wavefunction, where S (0) ann and its Coulomb limit are given in eqs. (3.14c) and (3.14d); we repeat them here for convenience (4.13c)

Identical particles
If the interacting particles are identical, then their total wavefunction is either symmetric or antisymmetric in their interchange, depending on whether the particles are bosons or fermions, respectively. For a pair of fermions, the spatial wavefunction depends on their spin state. A pair of spin-1/2 particles may be either in the antisymmetric spin-singlet state, or in the symmetric spin-triplet state. Their spatial wavefunction should then be symmetric or antisymmetric, respectively. Thus, for a pair of identical particles, the scattering-state spatial wavefunctions are Bosons, Fermions with total spin 0:

14a)
Fermions with total spin 1: The wavefunction (4.14a) implies that the contribution of the even-modes participating in a process is doubled, while the contribution of the odd-modes vanishes. The opposite is true for the wavefunction (4.14b). For a pair of bosonic identical particles (IP), the BSF and annihilation cross-sections are related to those for distinguishable particles (DP), computed in section 4.2, as follows σ ann for IP = 2 × [σ ann for DP], c.f. eq. (4.13a) . (4.15c) Note that the vanishing result in eq. (4.15b) holds to working order in α; contributions of higher order in α, which we have not computed here, will yield a non-zero cross-section.

Discussion
The general aspects discussed in section 3.5 in the context of BSF via emission of a vector boson, are pertinent also for BSF with emission of a scalar boson. Here, we point out some features that are specific to the latter. For BSF via emission of a scalar boson, the dominant transition modes are different for particleantiparticle or identical-particle pairs, than for pairs of particles with different masses and couplings to the emitted scalar boson. The capture to the ground state is dominated by the monopole and quadrupole modes in the first case [ = 0 and = 2 modes of the scattering-state wavefunction, respectively, c.f. eq. (4.10)], and by the dipole mode in the second case [ = 1 mode, c.f. eq. (4.5)]. The monopole and quadrupole modes contribute also to the latter case, but at higher order in the coupling than the dipole transition.
For bosonic particle-antiparticle pairs and pairs of annihilating identical bosons, BSF via scalar emission is significantly slower than annihilation into two scalar bosons. In the Coulomb regime, using eqs. (4.11), (4.13) and (4.15), we find, for both self-conjugate and non-self-conjugate species, (4.16) At ζ 1, this ratio becomes ∼ 0.16α 2 . Away from the Coulomb regime, the relative significance of BSF with respect to annihilation is further diminished. Indeed, as we have seen, the contribution of an mode of the scattering-state wavefunction to an inelastic process scales as σ inel v rel ∝ v 2 rel at low velocities. For BSF, the = 2 contribution will thus diminish, leaving ultimately only the = 0 mode at sufficiently low v rel . Equation (4.16) and the above discussion imply that the formation and decay of unstable bound states via emission of a scalar mediator cannot deplete significantly the DM density in the early universe, or enhance the indirect detection signals today, in contrast to the case of BSF via vector emission [2,7].
Nevertheless, BSF may be important for the capture of asymmetric DM into stable bound states. Moreover, because of the different velocity scaling of the various modes, even for pairs of particles with different masses/couplings to the radiated scalar boson, the = 0 mode may dominate the formation of zero-angular momentum bound states at low enough velocities. This is despite the contribution of the = 0 mode being suppressed by a higher order in the coupling with respect to that of the = 1 mode. Scalar mediator, non-self-conjugate bosonic particle-antiparticle pairs Coulomb limit contains also a contribution from the = 2 mode; away from the Coulomb regime, this contribution decreases with decreasing velocity, and eventually becomes subdominant, as can been seen in the left panel. Scalar mediator, non-self-conjugate bosonic particle-antiparticle pairs Capture into n = 2, = 1 states: Resonances

Conclusion
We have computed the cross-sections for the radiative formation of bound states by particles whose interaction is described in the non-relativistic regime by a Yukawa potential. We considered capture processes via emission of either a vector or a scalar boson, and inspected in detail the features of the cross-sections in the entire parametric regime where bound states exist. Bound-state effects can be important both for hidden-sector scenarios in which DM couples directly to light force mediators, as well as for TeV-scale WIMP models where non-perturbative effects due to the long-range nature of the interactions have already been shown to be significant [21][22][23][24][25][26][27].
The formation of DM bound states has multifaceted implications. The formation of unstable bound states in the early universe, and their subsequent decay, can deplete the density of symmetric or self-conjugate thermal-relic DM, and therefore affect the predictions for its mass and couplings [7,[28][29][30]. The same chain of processes taking place in the dense environment of haloes today [2,15,31], or in the interior of stars where DM may be captured [16], enhances the expected rate of the indirect detection signals, and, importantly, gives rise to correlated spectral features. Indeed, besides the high-energy radiation produced in the decay of the bound states, their formation is accompanied by the emission of low-energy radiation that dissipates the binding energy, and which may be detectable. In fact, the low-energy radiation emitted during capture, or in various level transitions between bound-state energy levels, is a potential source of indirect detection signals from asymmetric DM with long-range self-interactions, which forms stable bound states [11-13, 32, 33]. The cosmological formation of stable bound states by asymmetric DM typically screens or curtails the DM self-interactions, and has to be properly accounted for in any consistent phenomenologically study [34]. This is particularly important in the context of the self-interacting DM scenario [35][36][37][38], as well as in scenarios that feature a dissipative hidden sector [39][40][41][42][43][44][45][46]. Finally, DM bound states can give rise to distinct signatures in direct detection experiments [47,48] Appendices A Wavefunctions

A.1 The Schrödinger equation
The Schrödinger equations for the bound and scattering states are with the wavefunctions normalised as follows For a central potential V (r) = V (r) potential, we perform the standard separation of variables We set x = κr and γ ≡ −2µE/κ .
Then, for the Yukawa potential of eq. (1.1), the radial Schrödinger equations read where we temporarily dropped the indices in the wavefunctions and energy eigenvalues for generality. At x → 0, and for > 0, the second term of eq. (A.6) is dominated by the centrifugal contribution. In this region, the two independent solutions of eq. (A.6) scale as x +1 (regular) and x − (irregular). Here, we are interested in regular solutions, for E = E n < 0 and E = E k > 0, representing bound and scattering states respectively. This implies the boundary condition The condition (A.7) will be valid also for = 0.   Table 4. The fit parameters ξc and ρ, in the fitting formula γ n (ξ) (1/n)(1 − n 2 ξc/ξ) ρ , for the energy eigenvalues E n = −γ 2 n (ξ) × κ 2 /(2µ). In the Coulomb limit, lim ξ→∞ γ n (ξ) = 1/n.

A.2 Bound states
For E < 0, we seek solutions of eq. (A.6) that vanish at infinity, The boundary conditions (A.7) and (A.8), and the normalisation condition (A.4a) completely specify the discrete spectrum of wavefunctions and energy eigenvalues. For a Yukawa potential, the discrete energy eigenvalues depend on the principal quantum number n, as well as on , The lifting of the well-known -degeneracy of the energy eigenvalues of the Coulomb limit, is due to the non-conservation of the Laplace-Runge-Lenz vector by the Yukawa potential. We determine γ n (ξ) numerically, and present it in fig. 12, for n = 1, 2, 3. We find that it can be well fit by the formula The best fit parameters ξ c and ρ are given in table 4, for 1 n 3. 9 Equation (A.10) reproduces the Coulomb limit, lim ξ→∞ γ n (ξ) = 1/n. Away from the Coulomb limit, the existence of bound states implies an n-and -dependent lower bound on ξ.

Coulomb limit
In the limit ξ → ∞, eq. (A.6) admits analytic solutions with where L a n are the generalised Laguerre polynomials of degree n. (We assume the normalisation condition

A.3 Two-particle scattering state
For the continuous spectrum with E = E k > 0, we set 13) or equivalently γ 2 = −1/ζ 2 . The wavefunctions are specified by the boundary condition (A.7), and the asymptotic behaviour at x → ∞. At large x, the wavefunction χ |k|, behaves as (see e.g. Ref. [14, chapter 7]) where the phase shifts δ depend on ζ and ξ. This implies that

Coulomb limit
At ξ → ∞, the analytical solutions of eq. (A.6) with E = E k > 0, are The sum (A.3b) over the modes can be expressed in closed form, The cross-sections for radiative BSF depend on the following integrals [1] where b is proportional to the momentum of the radiated particle, b ∝ P ϕ . Using the Schrödinger eq. (A.1a), and eq. (A.9), K k,n m (b) takes also the form The first term in eq. (B.1d) yields a higher order correction to the contribution from I k,n m (b) in a given process, and we shall typically ignore it. The second term above is the leading relativistic correction.

B.2 Useful identities for angular integration
The plane waves e −ib·r , which appear in the integrals (B.1) and essentially stand for the wavefunction of the radiated boson, can be expanded in terms of Legendre polynomials using the identity where j is the spherical Bessel function. Derivatives of j that will arise in the following, can be re-expressed in terms of Bessel functions using the identity The Legendre polynomials can be expanded in spherical harmonics as follows Using the expansion (B.2a) and the identity (B.2c), will give rise to angular integrals that involve three spherical harmonics, and which can be expressed in terms of the Wigner-3j symbol, (Ω) = (2 1 + 1)(2 2 + 1)(2 3 + 1) 4π (B.2d) We always assume that the spherical harmonics are normalised according to

B.3 Expansion in the momentum of the emitted radiation
The exponential decay of the bound-state wavefunction ψ n m (r) at large r ensures that the integrands in eqs. (B.1) are significant only for κr γ n (ξ) few (c.f. section A.2). In this range, the argument of the Bessel function in eq. (B.2a) is where here η = η 1 or η 2 , and we used eq. (2.5a) for |P ϕ |. From this we deduce that br The condition (B.3) covers the range of interest. Provided that it is satisfied, we may expand the Bessel function of eq. (B.2a), and keep only leading-order terms. For this purpose, we shall use the expansion where for our purposes, z = br.
In computing the integrals (B.1), we express the wavefunctions as in eqs. (A.3), and carry out the radial integration in the variable x = κr. The expansion of the Bessel function over z = (b/κ)x amounts thus to an expansion over b/κ. Since b ∝ |P ϕ |, with |P ϕ | given in eq. (2.5a), this is ultimately an expansion in α.

B.4 Capture into bound states of arbitrary angular momentum
In order to evaluate the integrals (B.1), we first perform the angular integration using the identities (B.2), and then expand in powers of b/κ, using eq. (B.4). We find In the above, { I , m I } and { R , m R } are the orbital angular momentum quantum numbers of the incoming state and the radiated particle, respectively. In the ∇ x operator of eq. (B.5b), the radial coordinate should be understood to be x ≡ κr. We also note that (−1) m Y * m (Ω) = Y ,−m (Ω). The expansions (B.5) can be used to evaluate the amplitudes for the capture processes of interest, by keeping the leading order terms, as appropriate.

B.5 Capture into = 0 bound states
We may evaluate the integrals needed for capture into zero angular momentum bound states, directly from eqs. (B.5). Instead, here we shall perform the angular integration independently, and then expand in powers of the radiated momentum.

B.5.1 Angular integration
We will need the following angular integrals × k dP dy j (bx/κ) bx/κ +b P (y) − y dP (y) dy j (bx/κ) bx/κ − P (y) j +1 (bx/κ) For I k,n00 , the zero-th order contribution vanishes, due to the orthogonality of the ψ n m and φ k (or particularly, the χ n and χ |k|, ) wavefunctions. In our computations in sections 3 and 4, we shall need terms up to O[(b/κ) 2 ]. For J k,n00 and K k,n00 , the zero-th order terms in b/κ suffice. Using eqs. (C.3), we expand the integrals (C.4), keeping terms up to order (b/κ) 2 for I k,n00 , and only zero-th order terms for J k,n00 and K k,n00 . We obtain the following 11 11 Note that the contribution from the b-independent term of eq. (C.3c) to the I C k→n00 integral vanishes (as expected) when the summation over s is performed, i.e.

D Bound-state formation in momentum space
In this appendix, we describe a momentum-space procedure to compute BSF cross-sections, that is based on methods developed originally for few-body problems, by the nuclear-physics community.
We have adjusted codes that were written for low-energy proton-proton collisions [50] and deuteron formation via neutron capture on a proton target [51], in order to calculate the BSF cross sections discussed in the rest of this paper. Below, we outline the procedure for the formation of scalar DM bound states via a vector mediator. More details can be found in Ref. [50,52,53]. All BSF cross sections obtained in this work have been checked by both the coordinate-and momentum-space routine. We point out that the notation and some conventions in this appendix are somewhat disjoint from those in the main text, and some symbols are used here for different purposes. where E is the center-of-mass energy, p = |p| and p = |p | are the relative momenta of the incoming and outgoing DM particles in the center-of-mass frame, and T l l s s j denotes the T -matrix (scattering matrix) element corresponding to conserved total angular momentum j for states with initial and final orbital angular momentum (spin) l (s) and l (s ). V l l s s j (p , p) denotes a partialwave-decomposition of the DM potential. For scalar DM we remove the spin indices and use j = l = l , such that T l (p , p, E) = V l (p , p) + ∞ 0 dp V l (p , p ) p 2 E − p 2 /(2µ) + i T l (p , p, E) . (D.1) The potential for a vector mediator in momentum space is given by (D. 2) The partial-wave-decomposed potential that appears in eq. (D.1) is defined as V l (p , p) = 1 (2π) 3 p l||V (p, p )||p l where P l (x) denotes the Legendre polynomials.
To numerically solve the LS equation, we need to deal with the i in the numerator of eq. (D.1). We write where P denotes the principal value integral, and we introduced E ≡ q 2 0 /(2µ). The LS equation can then be written a T l (p , p, E) = V l (p , p) + (2µ)P pmax 0 dp V l (p , p ) p 2 q 2 0 − p 2 T l (p , p, E) −i π(2µ)q 0 2 V l (p , q 0 )T l (q 0 , p, E) , (D.5) where we introduced p max which corresponds to the maximum momentum of the momentum grid that is applied in the actual numerical solution. The main problem is the divergence at p = q 0 which we therefore subtract and add to get T l (p , p, E) = V l (p , p) + (2µ) pmax 0 dp V l (p , p ) p 2 q 2 0 − p 2 T l (p , p, E) −V l (p , q 0 ) q 2 0 q 2 0 − p 2 T l (q 0 , p, E) +(2µ)q 2 0 V l (p , q 0 )T l (q 0 , p, E)P pmax 0 dp 1 q 2 0 − p 2 −i π(2µ)q 0 2 l s V l (p , q 0 )T l (q 0 , p, E) , (D. 6) such that the first integral is no longer singular. The second integral can be done analytically and does not depend on the form of the potential The LS equation can now be discretized on a momentum grid and written as a complex eigenvalue equation which we solve using the LAPACK library [54].
To obtain the momentum-space bound-state wave function we solve the homogeneous part of the LS equation where E b nlm is the (negative) binding energy and the real wave function is normalized as dp p 2 Ψ nlm (p) 2 = 1 . (D.9) Equation (D.8) can be immediately discretized and written as an eigenvalue equation which we again solve with a LAPACK routine [54]. The binding energy is varied until we find a consistent solution.
Although for simplicity we discussed only scalar DM and vector mediator, the above routines can be easily extended to solve coupled-channel LS and bound-state equations.

D.2 The BSF cross sections
The goal is to obtain the amplitude for the BSF process X 1 + X 2 → B nlm (X 1 X 2 ) + ϕ. We write this amplitude as A nlm,λ (p in , P ϕ ) = Ψ nlm | O λ (P ϕ ) | p in (+) , (D. 10) where λ denotes the polarization of the vector mediator and p in (p in = |p in | µv rel ) and |P ϕ | = P ϕ are, respectively, the incoming relative momentum of the DM pair in the c.o.m. frame and the outgoing mediator momentum. O λ describes the current. The (+) superscript on the incoming state implies that this is the fully scattered state obtained from applying the T -matrix to a free state.
We insert a complete set of states 1 = in terms of Clebsch-Gordan coefficients (j 1 j 2 j; m 1 m 2 m 1 + m 2 ), we introduced a nine-J symbol, andl = 2l + 1. The function g nlm f G (p, P ϕ ) denotes a single numerical angular integral 15) where k(x ) = p 2 + P 2 ϕ /4 − pP ϕ x . We are now in the position to evaluate the scattering amplitude in eq. (D.10). Although the obtained expressions are valid in any coordinate frame, it is convenient to specify P ϕ = P ϕẑ . In this frame, the differential BSF cross section is given by |A nlm,λ (p in , P ϕ = P ϕẑ )| 2 + m 2 ϕ P 2 ϕ + m 2 ϕ |A nlm,λ=0 (p in , P ϕ = P ϕẑ )| 2 , (D.16) such that only the transverse polarizations contribute in the Coulomb limit.