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

Kalliopi Petraki, Marieke Postma, Jordy de Vries

► To cite this version:

Kalliopi Petraki, Marieke Postma, Jordy de Vries. Radiative bound-state-formation cross-sections for dark matter interacting via a Yukawa potential. Journal of High Energy Physics, 2017, 4, pp.077. 10.1007/JHEP04(2017)077. hal-01522604

HAL Id: hal-01522604

https://hal.sorbonne-universite.fr/hal-01522604v1

Submitted on 15 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RECEIVED: November 14, 2016

REVISED: January 31, 2017

ACCEPTED: March 27, 2017

PUBLISHED: April 13, 2017

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

Kalliopi Petraki,ª,ᵇ Marieke Postmaᵇ and Jordy de Vriesᵇ

ªLPTHE, CNRS, UMR 7589,
4 Place Jussieu, F-75252, Paris, France

ᵇNikhef,
Science Park 105, 1098 XG Amsterdam, The Netherlands

E-mail: kpetraki@lpthe.jussieu.fr, mpostma@nikhef.nl,
jordy.de.vries@nikhef.nl

**ABSTRACT:** 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.

**KEYWORDS:** Beyond Standard Model, Cosmology of Theories beyond the SM, Nonperturbative Effects

ARXIV ePRINT: 1611.01394
# Contents

1 Introduction 1

2 Preliminaries 6

2.1 Parametrisation 7

2.2 The range of the $\xi$ parameter 8

2.3 The cross-section 9

2.4 Initial and final state convolution integrals 9

3 Vector force mediator 10

3.1 Radiative capture into a bound state 10

3.2 Capture into $\ell = 0$ bound states 11

3.3 Capture into $\ell = 1$ bound states 13

3.4 Annihilation of particle-antiparticle pairs 15

3.5 Discussion 15

4 Scalar force mediator 27

4.1 Pairs of non-degenerate particles 28

4.1.1 Capture into $\ell = 0$ bound states 28

4.1.2 Capture into $\ell = 1$ bound states 29

4.2 Bosonic particle-antiparticle pairs (non-identical) 31

4.2.1 Capture into $\ell = 0$ bound states 32

4.2.2 Capture into $\ell = 1$ bound states 33

4.2.3 Annihilation 34

4.3 Identical particles 34

4.4 Discussion 35

5 Conclusion 38

A Wavefunctions 41

A.1 The Schrödinger equation 41

A.2 Bound states 42

A.3 Two-particle scattering state 43

B Convolution integrals 44

B.1 Definition 44

B.2 Useful identities for angular integration 44

B.3 Expansion in the momentum of the emitted radiation 45

B.4 Capture into bound states of arbitrary angular momentum 45

B.5 Capture into $\ell = 0$ bound states 46

B.5.1 Angular integration 46

B.5.2 Leading-order contributions 47
<table><tbody><tr><td>B.6</td><td>Capture into <i>ℓ</i> = 1 bound states</td><td>48</td></tr><tr><td>C</td><td>Coulomb limit for capture into <i>ℓ</i> = 0 bound states</td><td>48</td></tr><tr><td></td><td>C.1 Wavefunctions</td><td>48</td></tr><tr><td></td><td>C.2 Identities</td><td>49</td></tr><tr><td></td><td>C.3 Convolution integrals</td><td>49</td></tr><tr><td></td><td>C.4 Convolution integrals: expansion in the momentum of the emitted radiation</td><td>50</td></tr><tr><td>D</td><td>Bound-state formation in momentum space</td><td>50</td></tr><tr><td></td><td>D.1 Solution of the Lippmann-Schwinger and bound-state equation</td><td>51</td></tr><tr><td></td><td>D.2 The BSF cross sections</td><td>52</td></tr></tbody></table>

# 1 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,

$$V_Y(\mathbf{r}) = -\frac{\alpha e^{-m_{\phi}r}}{r}, \qquad (1.1)$$
where $\alpha$ parametrises the interaction strength, and $m_\phi$ is the mediator mass. While no analytical expressions for the cross-sections of interest can be derived for $m_\phi > 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 particle-antiparticle 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 (cf. 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 fol-
lowing 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 kinemati-
cally 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 interac-
tion (1.1) (cf. 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 (cf. section 2.1); using this parametrisation, we
perform a rather thorough inspection of the features of the BSF cross-sections.
<table><thead><tr><th>Cross-sections for the radiative formation of bound states {nℓm}, by particles interacting via a Yukawa potential</th><th>Equations</th><th>figures</th></tr></thead><tbody><tr><td><strong>Vector force mediator</strong></td><td></td><td></td></tr><tr><td>■ Capture into ℓ = 0 bound states</td><td>(3.7), (3.8)<br>(3.8),</td><td>2–5</td></tr><tr><td>• Coulomb limit and comparison with annihilation</td><td>(3.14),<br>(3.17)</td><td>2</td></tr><tr><td>• Capture into the ground state {100}, and comparison with annihilation</td><td>(3.7),<br>(3.8d),<br>(3.14)</td><td>3</td></tr><tr><td>– Resonant structure</td><td></td><td>3a</td></tr><tr><td>– Velocity dependence off-resonance</td><td></td><td>3b</td></tr><tr><td>– Velocity dependence on-resonance</td><td></td><td>3c</td></tr><tr><td>– Velocity dependence near threshold</td><td></td><td>3d</td></tr><tr><td>• Capture into the first excited state n = 2, ℓ = 0</td><td>(3.7), (3.8)</td><td>4</td></tr><tr><td>• Comparison: the resonant structure of bound-state formation and p-wave annihilation</td><td></td><td>5</td></tr><tr><td>■ Capture into ℓ = 1 bound states</td><td>(3.13)</td><td>2, 6</td></tr><tr><td><strong>Scalar force mediator: non-degenerate particles</strong></td><td></td><td></td></tr><tr><td>■ Capture into ℓ = 0 states (including Coulomb limit)</td><td>(4.5), (4.6)</td><td>7</td></tr><tr><td>■ Capture into ℓ = 1 states</td><td>(4.7)</td><td>8</td></tr><tr><td><strong>Scalar force mediator: non-self-conjugate bosonic particle-antiparticle pairs</strong></td><td></td><td></td></tr><tr><td>■ Capture into ℓ = 0 states</td><td>(4.10), (4.11)</td><td>9–10</td></tr><tr><td>– Coulomb limit and comparison with annihilation</td><td>(4.11), (4.16)</td><td>9</td></tr><tr><td>– Resonant structure</td><td></td><td>10a</td></tr><tr><td>– Velocity dependence off- and on-resonance</td><td></td><td>10b</td></tr><tr><td>■ Capture into ℓ = 1 states</td><td>(4.12)</td><td>11</td></tr><tr><td><strong>Scalar force mediator: identical particles</strong></td><td>(4.15)–(4.16)</td><td>9–11</td></tr></tbody></table>

<table><caption>**Table 1.** References to the main results and formulae.</caption><thead><tr><th>Convolution integrals I<sub>k,nℓm</sub>(b), J<sub>k,nℓm</sub>(b), K<sub>k,nℓm</sub>(b)</th><th>Equations</th></tr></thead><tbody><tr><td>Definition</td><td>(2.7)</td></tr><tr><td>Expansion in powers of radiated momentum |b|: validity of approximation</td><td>(B.3)</td></tr><tr><td>Capture into {nℓm} bound states, expansion in |b|</td><td>(B.5)</td></tr><tr><td>Capture into ℓ = 0 bound states</td><td></td></tr><tr><td>– Expansion in |b|, leading-order terms</td><td>(B.9)</td></tr><tr><td>– Coulomb limit (no expansion)</td><td>(C.4)</td></tr><tr><td>– Coulomb limit: expansion in |b|, leading-order terms</td><td>(C.5)</td></tr><tr><td>Capture into ℓ = 1 bound states: expansion in |b|, leading-order terms</td><td>(B.10)</td></tr></tbody></table>
<table><thead><tr><th>Particles and masses</th><th>Symbols</th></tr></thead><tbody><tr><td>Interacting particles</td><td>X<sub>1</sub>, X<sub>2</sub></td></tr><tr><td>Masses of interacting particles</td><td>m<sub>1</sub>, m<sub>2</sub></td></tr><tr><td>Total mass of interacting particles</td><td>M ≡ m<sub>1</sub> + m<sub>2</sub></td></tr><tr><td>Reduced mass of interacting particles</td><td>μ ≡ <sup>m<sub>1</sub>m<sub>2</sub></sup>&frasl;<sub>m<sub>1</sub> + m<sub>2</sub></sub></td></tr><tr><td>Force mediator (scalar or vector)</td><td>φ</td></tr><tr><td>Mass of force mediator</td><td>m<sub>φ</sub></td></tr></tbody></table>

<table><thead><tr><th>Description</th><th>Symbol</th></tr></thead><tbody><tr><td>Dark fine structure constant</td><td>α</td></tr><tr><td>Expectation value of relative velocity of interacting particles in the scattering state</td><td><b>v</b><sub>rel</sub>,<br>v<sub>rel</sub> = |<b>v</b><sub>rel</sub>|</td></tr><tr><td>Bohr momentum</td><td>κ ≡ μα</td></tr><tr><td>Momentum of reduced system of interacting particles in the scattering state</td><td><b>k</b> ≡ μ<b>v</b><sub>rel</sub>,<br>k ≡ |<b>k</b>|</td></tr><tr><td>Wavefunction of nlm bound state</td><td>ψ<sub>nlm</sub>(<b>r</b>) = κ<sup>3/2</sup> [χ<sub>nℓ</sub>(κr) / κr] Y<sub>ℓm</sub>(Ω<sub><b>r</b></sub>)</td></tr><tr><td>Wavefunction of scattering state</td><td>φ<sub><b>k</b></sub>(<b>r</b>) = ∑<sub>ℓ=0</sub><sup>∞</sup> (2ℓ + 1) [χ<sub>|<b>k</b>|,ℓ</sub>(κr) / κr] P<sub>ℓ</sub>(hat(<b>k</b> · ˆr))</td></tr><tr><td>Binding energy of nlm bound state</td><td>ε<sub>nℓ</sub> = -κ<sup>2</sup> / (2μ) γ<sub>nℓ</sub><sup>2</sup>(ξ)</td></tr><tr><td>Kinetic energy of scattering state</td><td>ε<sub><b>k</b></sub> = <sup>k<sup>2</sup></sup>&frasl;<sub>2μ</sub></td></tr><tr><td>Phase-space suppression due to emission of a massive force mediator, during capture into the nlm bound states</td><td>pssl<sub>nℓ</sub>, see eq. (2.5b)</td></tr></tbody></table>

<table><caption>**Table 2.** Notation.</caption><thead><tr><th>Dimensionless parameters</th></tr></thead><tbody><tr><td>η<sub>1,2</sub> ≡ <sup>m<sub>1,2</sub></sup>&frasl;<sub>m<sub>1</sub> + m<sub>2</sub></sub></td></tr><tr><td>ζ ≡ κ/k = α/v<sub>rel</sub></td></tr><tr><td>ξ ≡ κ/m<sub>φ</sub> = μα/m<sub>φ</sub></td></tr></tbody></table>
## 2 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 $\phi$ (cf. figure 1a). The interaction La-grangians 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_\phi \lesssim \mu\alpha \equiv \kappa$, where $\mu$ is the $X_1-X_2$ reduced mass.<fnref n="1" /> Under this condition, we want to compute the cross-sections for the radiative formation of bound states, with emission of the force mediator $\phi$,

$$ \mathcal{U}_{\mathbf{k}}(X_1 + X_2) \rightarrow \mathcal{B}_{n\ell m}(X_1 X_2) + \phi. \qquad (2.1) $$

Here, $\mathcal{U}_{\mathbf{k}}(X_1 + X_2)$ stands for a two-particle scattering (unbound) state, characterised by the continuous vector quantum number $\mathbf{k} = \mu\mathbf{v}_{\text{rel}}$, with $\mathbf{v}_{\text{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 $\mathcal{U}_{\mathbf{k}}(X_1 + X_2)$ state is instead described by a wavefunction $\phi_{\mathbf{k}}(\mathbf{r})$ that obeys the Schrödinger equation with the Yukawa potential of eq. (1.1) and a positive energy eigenvalue $\mathcal{E}_{\mathbf{k}}$, parametrised by $\mathbf{k}$. Further, $\mathcal{B}_{n\ell m}(X_1 X_2)$ is a bound state, whose wavefunction $\psi_{n\ell m}(\mathbf{r})$, parametrised by the familiar principal and angular-momentum discrete quantum numbers $\{n, \ell, m\}$, obeys the Schrödinger equation with the same potential and a negative energy eigenvalue $\mathcal{E}_{nl}$. We discuss the numerical computation of $\psi_{n\ell m}(\mathbf{r})$ and $\phi_{\mathbf{k}}(\mathbf{r})$ in appendix A.

As is well-known, the Yukawa potential arises from the one-boson-exchange diagram, shown in figure 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 figure 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 figure 1a, the two vertices introduced by each boson exchange imply a suppression by one power of the coupling $\alpha$. However, this suppression is cancelled by the loop momentum exchange, which also scales with $\alpha$. In particular, the average momentum exchange along each virtual boson scales as $|\mathbf{q}| \sim \mu\alpha$, and the off-shellness of the $X_1, X_2$ propagators scales as $q^0 \propto |\mathbf{q}|^2 \propto \alpha^2$; the integration over the loop energy and momentum also yield factors of $\alpha^2$ and $\alpha$, 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 figure 1b. The initial-state ladder corresponds to a scattering state and is evaluated at center-of-momentum

<footnote n="1">¹More precise conditions for the existence of bound states are discussed in appendix A.2.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1a.** Left: in the weak-coupling regime, the one-boson (either scalar or vector) exchange is the dominant contribution to the non-relativistic $X_1 - X_2$ 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.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 1b.** 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.</figcaption></figure>

(CM) energy $E = m_1 + m_2 + \mathcal{E}_k \ge 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 + \mathcal{E}_{nl} < 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).

## 2.1 Parametrisation

We formulate our computations and present our results in terms of two dimensionless parameters,

$$ \zeta \equiv \frac{\mu\alpha}{\mu v_{\text{rel}}} = \frac{\alpha}{v_{\text{rel}}} \quad \text{and} \quad \xi \equiv \frac{\mu\alpha}{m_{\varphi}}, \qquad (2.2) $$

where $v_{\text{rel}}$ is the relative velocity of the interacting particles and $\mu$ is their reduced mass. $\zeta$ and $\xi$ suffice to characterise the solutions of the Schrödinger equation for the Yukawa potential (1.1), for bound and scattering states (cf. appendix A), and to parametrise the non-analytical parts of our computations that arise in the BSF cross-sections.

The $\zeta$ parameter compares the average momentum transfer $\sim \mu v_{\text{rel}}$ between two un-bound particles, with the momentum transfer $\sim \mu\alpha$ that occurs in the exchange of virtual force mediators, as explained above. In practice, $\zeta$ parametrises the velocity dependence of the cross-sections.
On the other hand, $\xi$ parametrises the model dependence of the cross-sections. It compares the two physical scales involved: the Bohr momentum $\kappa \equiv \mu\alpha$, 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_\phi$, which determines the range of the interaction. The interaction manifests as long-range roughly if $\xi \gtrsim 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 $\zeta$ and $\xi$, exposes the physical significance of the features of the cross-sections, and can greatly facilitate phenomenological studies.

## 2.2 The range of the $\xi$ 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_\phi < \mathcal{E}_k - \mathcal{E}_{n\ell}$, where $\mathcal{E}_k = k^2/(2\mu)$ is the kinetic energy of the scattering state in the CM frame, and $\mathcal{E}_{n\ell} = -\gamma_{n\ell}^2 \times \kappa^2/(2\mu)$ is the binding energy of the bound state [cf. 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, $\mathcal{E}_k \lesssim |\mathcal{E}_{n\ell}|$ (which roughly implies $v_{rel} \lesssim \alpha/n$). Then, the condition for BSF with emission of a force mediator becomes roughly $m_\phi \lesssim \mu\alpha^2/(2n^2)$, or, in terms of the $\xi$ parameter,

$$ \xi > \frac{2}{\alpha \gamma_{n\ell}^2(\xi)} \gg n^2. \qquad (2.3) $$

(Note that $\gamma_{n\ell}(\xi) \le 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, $\gamma_{n\ell}(\xi) > 0$, which implies roughly $\xi = \mu\alpha/m_\phi \gtrsim n^2$ (cf. 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.<fnref n="2" /> This is not necessarily the case for the scattering state wavefunctions though, which depend on both $\xi$ and $\zeta$. As we shall see in the following, the Coulomb limit, which formally corresponds to $\xi \to \infty$, is essentially attained for

$$ \xi \gtrsim \zeta, \qquad (2.4) $$

i.e. when the average momentum transfer between the interacting particles exceeds the mediator mass, $\mu v_{rel} \gtrsim m_\phi$.

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 $\xi$ that do not satisfy the condition (2.3). In the

<footnote n="2">²This approximation was employed in refs. [2, 3].</footnote>
following, we shall thus consider the entire range of $\xi$ 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 [cf. 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.

## 2.3 The cross-section

The differential cross-section times relative velocity for the 2-to-2 process (2.1) is

$$v_{\text{rel}} \frac{d\sigma_{\text{BSF}}^{\{n\ell m\}}}{d\Omega} = \frac{|\mathbf{P}_{\varphi}|}{64\pi^2 M^2 \mu} |\mathcal{M}_{\mathbf{k} \to n\ell m}|^2,$$

where $\mathcal{M}_{\mathbf{k} \to n\ell m}$ is the transition amplitude, and $|\mathbf{P}_{\varphi}|$ 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

$$\sqrt{|\mathbf{P}_{\varphi}|^2 + m_{\varphi}^2} = \mathcal{E}_{\mathbf{k}} - \mathcal{E}_{n\ell} = \frac{\mathbf{k}^2}{2\mu} + \frac{\kappa^2}{2\mu} \gamma_{n\ell}^2(\xi),$$

where we used eqs. (A.9) and (A.13). Then

$$\frac{|\mathbf{P}_{\varphi}|}{\kappa} = \frac{\alpha}{2} \left[ \frac{1 + \zeta^2 \gamma_{n\ell}^2(\xi)}{\zeta^2} \right] \times \text{pss}_{n\ell}^{1/2}, \qquad (2.5a)$$

$$\text{pss}_{n\ell} = 1 - \frac{4\zeta^4}{\alpha^2\xi^2[1 + \zeta^2\gamma_{n\ell}^2(\xi)]^2}, \qquad (2.5b)$$

where $\text{pss}_{n\ell}$ is the phase-space suppression factor for the emission of a force mediator, during capture into the $n\ell m$ bound state; in the Coulomb limit $\xi \to \infty$ and $\text{pss}_{n\ell} = 1$. Putting everything together, we obtain

$$v_{\text{rel}} \frac{d\sigma_{\text{BSF}}^{\{n\ell m\}}}{d\Omega} = \frac{\alpha^2}{2^7 \pi^2 M^2} \text{pss}_{n\ell}^{1/2} \left( \frac{1 + \zeta^2 \gamma_{n\ell}^2(\xi)}{\zeta^2} \right) |\mathcal{M}_{\mathbf{k} \to n\ell m}|^2. \qquad (2.6)$$

The general expression for the amplitude $\mathcal{M}_{\mathbf{k} \to n\ell m}$ in terms of the initial and final state wavefunctions and the radiative vertex can be found in ref. [1, section 3]. In the following, we shall reproduce only the expressions that are relevant for the interactions we are considering.

## 2.4 Initial and final state convolution integrals

The transition amplitudes $\mathcal{M}_{\mathbf{k} \to n\ell m}$ depend on the scattering-state and bound-state wavefunctions, $\psi_{n\ell m}$ and $\phi_{\mathbf{k}}$, and on the radiative vertex, via the convolution integrals [1]
$$ \mathcal{I}_{\mathbf{k},nlm}(\mathbf{b}) \equiv \int \frac{d^3p}{(2\pi)^3} \tilde{\psi}_{nlm}^*(\mathbf{p}) \tilde{\phi}_{\mathbf{k}}(\mathbf{p}+\mathbf{b}) = \int d^3r \, \psi_{nlm}^*(\mathbf{r}) \phi_{\mathbf{k}}(\mathbf{r}) e^{-i\mathbf{b}\cdot\mathbf{r}}, \quad (2.7a) $$

$$ \mathcal{J}_{\mathbf{k},nlm}(\mathbf{b}) \equiv \int \frac{d^3p}{(2\pi)^3} \mathbf{p} \, \tilde{\psi}_{nlm}^*(\mathbf{p}) \tilde{\phi}_{\mathbf{k}}(\mathbf{p}+\mathbf{b}) = i \int d^3r \, [\nabla \psi_{nlm}^*(\mathbf{r})] \phi_{\mathbf{k}}(\mathbf{r}) e^{-i\mathbf{b}\cdot\mathbf{r}}, \quad (2.7b) $$

$$ \mathcal{K}_{\mathbf{k},nlm}(\mathbf{b}) \equiv \int \frac{d^3p}{(2\pi)^3} \mathbf{p}^2 \, \tilde{\psi}_{nlm}^*(\mathbf{p}) \tilde{\phi}_{\mathbf{k}}(\mathbf{p}+\mathbf{b}) = - \int d^3r \, [\nabla^2 \psi_{nlm}^*(\mathbf{r})] \phi_{\mathbf{k}}(\mathbf{r}) e^{-i\mathbf{b}\cdot\mathbf{r}}. \quad (2.7c) $$

Here $\tilde{\psi}_{nlm}$ and $\tilde{\phi}_{\mathbf{k}}$ are the Fourier transforms of $\psi_{nlm}$ and $\phi_{\mathbf{k}}$ respectively. The momentum $\mathbf{b}$ is proportional to the momentum of the radiated particle, $\mathbf{b} \propto \mathbf{P}_{\varphi}$, as we shall see in the following. In appendix B, we expand the integrals (2.7) in powers of $|\mathbf{b}|/\kappa \propto |\mathbf{P}_{\varphi}|/\kappa \propto \alpha$ [cf. 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.

# 3 Vector force mediator

## 3.1 Radiative capture into a bound state

We assume $X_1$ and $X_2$ to be coupled to a gauged U(1) force,

$$ \mathcal{L} = (D_{\mu}X_{1})^{\dagger}(D^{\mu}X_{1}) + (D_{\mu}X_{2})^{\dagger}(D^{\mu}X_{2}) - \frac{1}{4}F_{\mu\nu}F^{\mu\nu} \\ - m_{1}^{2}|X_{1}|^{2} - m_{2}^{2}|X_{2}|^{2} + \frac{1}{2}m_{\varphi}^{2}\varphi_{\mu}\varphi^{\mu}, \quad (3.1a) $$

$$ \mathcal{L} = \bar{X}_1 i \not D X_1 + \bar{X}_2 i \not D X_2 - m_1 \bar{X}_1 X_1 - m_2 \bar{X}_2 X_2 - \frac{1}{4} F_{\mu\nu} F^{\mu\nu} + \frac{1}{2} m_{\varphi}^2 \varphi_{\mu} \varphi^{\mu}, \quad (3.1b) $$

where $F^{\mu\nu} = \partial^{\mu}\varphi^{\nu} - \partial^{\nu}\varphi^{\mu}$ and $D^{\mu} = \partial^{\mu} - ic_jg\varphi^{\mu}$, with $c_1, c_2$ being the charges of $X_1, X_2$. The mass $m_{\varphi}$ of $\varphi^{\mu}$ 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.<fnref n="3" /> In the non-relativistic regime, the $\varphi^{\mu}$ exchange between $X_1$

<footnote n="3">³In the Higgs mechanism, this can be ensured by an appropriate choice of the charge of the scalar field breaking the gauged U(1) symmetry. In the Stückelberg mechanism, the massive gauge boson couples to a conserved current.</footnote>
and $X_2$ gives rise to the Yukawa potential of eq. (1.1), with

$$ \alpha = -\frac{c_1 c_2 g^2}{4\pi}. \qquad (3.2) $$

The interaction is attractive if $c_1c_2 < 0$.

Since BSF involves gauge interactions with conserved currents, as seen from eq. (3.1), the Ward identity ensures that $P_{\phi}^{\mu}M_{\mu} = 0$. This implies $M^0 = P_{\phi}^{j}M^j/P_{\phi}^0$. Then, the unpolarised amplitude is

$$ \begin{aligned} \sum_{\epsilon} |\mathcal{M}_{\mathbf{k} \to n\ell m}|^2 &= - \left( g_{\mu\nu} - \frac{P_{\phi,\mu} P_{\phi,\nu}}{m_{\phi}^2} \right) \mathcal{M}_{\mathbf{k} \to n\ell m}^{\mu} \mathcal{M}_{\mathbf{k} \to n\ell m}^{\nu*} \\ &= \mathcal{M}_{\mathbf{k} \to n\ell m}^{j} \mathcal{M}_{\mathbf{k} \to n\ell m}^{j*} - \frac{|P_{\phi}^{j} \mathcal{M}_{\mathbf{k} \to n\ell m}^{j}|^2}{P_{\phi}^2 + m_{\phi}^2}. \end{aligned} \qquad (3.3) $$

We thus need only the spatial components of the BSF amplitude, which are [1]

$$ \mathcal{M}_{\mathbf{k} \to n\ell m}^{j} = -2g\sqrt{2\mu} \left\{ \frac{c_1}{\eta_1} \mathcal{J}_{\mathbf{k},n\ell m}^{j}(\eta_2 \mathbf{P}_{\phi}) - \frac{c_2}{\eta_2} \mathcal{J}_{\mathbf{k},n\ell m}^{j}(-\eta_1 \mathbf{P}_{\phi}) \right. \\ \left. + \left[ c_1 \left( K^j - \frac{\eta_1 - \eta_2}{2\eta_1} P_{\phi}^j \right) \mathcal{I}_{\mathbf{k},n\ell m}(\eta_2 \mathbf{P}_{\phi}) + c_2 \left( K^j + \frac{\eta_1 - \eta_2}{2\eta_2} P_{\phi}^j \right) \mathcal{I}_{\mathbf{k},n\ell m}(-\eta_1 \mathbf{P}_{\phi}) \right] \right\}, \qquad (3.4) $$

where $\mathcal{I}_{\mathbf{k},n\ell m}$ and $\mathcal{J}_{\mathbf{k},n\ell m}^j$ are defined in eq. (2.7). It is evident from the momentum factor in the integrand of $\mathcal{J}_{\mathbf{k},n\ell m}^j$ [cf. eq. (2.7b)] and the momentum dependence of the $\mathcal{I}_{\mathbf{k},n\ell m}$ contribution to $\mathcal{M}_{\mathbf{k} \to n\ell 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.

## 3.2 Capture into $\ell = 0$ bound states

From eq. (B.9), we find that the lowest-order contributions to the amplitude (3.4) arise from the $\mathcal{J}_{\mathbf{k},n00}$ integrals and are of zeroth order in $|\mathbf{P}_{\phi}|/\kappa$. The unpolarised squared amplitude (3.3) is

$$ \sum_{\epsilon} |\mathcal{M}_{\mathbf{k} \to n00}|^2 \approx \frac{2^5 \pi \mu \alpha}{\eta_1^2 \eta_2^2} \left[ \frac{(\eta_2 c_1 - \eta_1 c_2)^2}{-c_1 c_2} \right] \left( |\mathcal{J}_{\mathbf{k},n00}|^2 - pss_{n,0} |\hat{\mathbf{P}}_{\phi} \cdot \mathcal{J}_{\mathbf{k},n00}|^2 \right). \qquad (3.5) $$

(Since we are now interested only in the zero-th order terms in $|\mathbf{P}_{\phi}|/\kappa$, we have dropped the arguments of the $\mathcal{J}_{\mathbf{k},n00}$ functions.) Substituting eq. (B.9b) into eq. (3.5) yields

$$ \begin{aligned} \sum_{\epsilon} |\mathcal{M}_{\mathbf{k} \to n00}|^2 &\approx \frac{2^7 \pi^2 M^2}{\mu^2} \left[ \frac{(\eta_2 c_1 - \eta_1 c_2)^2}{-c_1 c_2} \right] (1 - pss_{n,0} \cos^2 \theta) \times \\ &\quad \left| \int_0^{\infty} dx \left[ \frac{d\chi_{n,0}^*(x)}{dx} - \frac{\chi_{n,0}^*(x)}{x} \right] \chi_{|\mathbf{k}|,1}(x) \right|^2, \end{aligned} \qquad (3.6) $$

where $\chi_{n,\ell}(x)$ and $\chi_{|\mathbf{k}|,\ell}(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<fnref n="4" />

$$ \sigma_{\text{BSF}}^{\{n00\}} v_{\text{rel}} \approx \frac{\pi \alpha^2}{4\mu^2} \left[ \frac{(c_1\eta_2 - c_2\eta_1)^2}{-c_1c_2} \right] pss_{n,0}^{1/2} \left( \frac{3-pss_{n,0}}{2} \right) \times S_{\text{BSF}}^{\{n00\}}(\zeta, \xi), \quad (3.7a) $$

where

$$ S_{\text{BSF}}^{\{n00\}}(\zeta, \xi) = \frac{2^5}{3} \left( \frac{1 + \zeta^2 \gamma_{n,0}^2(\xi)}{\zeta^2} \right) \left| \int_0^\infty dx \left[ \frac{d\chi_{n,0}^*(x)}{dx} - \frac{\chi_{n,0}^*(x)}{x} \right] \chi_{|k|,1}(x) \right|^2 . \quad (3.7b) $$

**Coulomb limit.** At $m_\phi \to 0$, the above becomes

$$ \lim_{\xi \to \infty} \sigma_{\text{BSF}}^{\{n00\}} v_{\text{rel}} \approx \frac{\pi \alpha^2}{4\mu^2} \left[ \frac{(c_1\eta_2 - c_2\eta_1)^2}{-c_1c_2} \right] S_{\text{BSF,C}}^{\{n00\}}(\zeta), \quad (3.8a) $$

where, using eq. (C.5b), we find

$$ S_{\text{BSF,C}}^{\{n00\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \frac{2^9}{3n^7} \zeta^4 (1 + \zeta^2) (1 + \zeta^2/n^2) \times \\ \times \left[ \sum_{s=0}^{n-1} \frac{n!(2n-s)}{(n-s-1)! (s+2)! s!} \frac{2^s \zeta_n^s}{d\zeta_n^s} \left[ \frac{e^{-2\zeta\arccot\zeta_n}}{(1 + \zeta_n^2)^2} \right] \right]_{\zeta_n = \zeta/n}^2 . \quad (3.8b) $$

Expanding the sum, $S_{\text{BSF,C}}^{\{n00\}}$ takes the form

$$ S_{\text{BSF,C}}^{\{n00\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \frac{2^9}{3n^3} \frac{\zeta^4 (1 + \zeta^2)}{(1 + \zeta^2/n^2)^{2n-1}} [\varrho_n(\zeta)]^2 e^{-4\zeta \arccot(\zeta/n)}, \quad (3.8c) $$

where $\varrho_n(\zeta)$ is a rational function of $\zeta^2$ (and specifically, a polynomial of degree $n-2$, for $n \ge 2$), with $\lim_{\zeta \to 0} \varrho_n(\zeta) = 1$. We give the explicit expressions of $\varrho_n(\zeta)$ for $1 \le n \le 5$ in table 3. For capture into the ground state in particular,

$$ S_{\text{BSF,C}}^{\{100\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \frac{2^9}{3} \frac{\zeta^4}{(1 + \zeta^2)^2} e^{-4\zeta \arccot \zeta}. \quad (3.8d) $$

We illustrate in detail the features of $S_{\text{BSF}}^{\{n00\}}(\zeta, \xi)$ and $S_{\text{BSF,C}}^{\{n00\}}(\zeta)$ in figures 2 to 5, and we discuss them in section 3.5.

Let us now recount the origin of the various factors in eqs. (3.7a) and (3.8c). In eq. (3.7a), the factor $(c_1\eta_2 - c_2\eta_1)^2/(-c_1c_2) = (c_1/\eta_1 - c_2/\eta_2)^2 \times [\eta_1^2\eta_2^2/(-c_1c_2)]$ emanates from the vertices of the radiated gauge boson in the Feynman diagrams of figure 1b, and

<footnote n="4">⁴Here and in the following, we typically choose to factorise the various BSF cross-sections (i.e. separate out the $S_{\text{BSF}}$ factors) in a way that facilitates the comparison between BSF and annihilation processes for particle-antiparticle pairs. However, we note that the $S_{\text{BSF}}$ factors do not carry any physical significance on their own. This is in contrast to the $S_{\text{ann}}$ factors appearing in the various annihilation cross-sections [cf. 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 $\zeta, \xi \ll 1$, the non-perturbative effects switch off, and the annihilation cross-sections reduce to their perturbative values ($S_{\text{ann}} \simeq 1$). In the same regime, $S_{\text{BSF}} \to 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_{\text{BSF}}$ factors are not enhancement factors of otherwise perturbative processes, by a non-perturbative effect.</footnote>
<table><caption>**Table 3.** The functions $\rho_n(\zeta)$ appearing in the Coulomb limit of the cross-sections for the radiative formation of zero angular momentum states of principal quantum number $n$, (i) with emission of a vector force mediator [cf. eq. (3.8c)], and (ii) with emission of a scalar mediator, by particles of different species [cf. eq. (4.6c)].</caption><thead><tr><th>n</th><th>ρ<sub>n</sub>(ζ)</th></tr></thead><tbody><tr><td>1</td><td><sup>1</sup>&frasl;<sub>1 + ζ<sup>2</sup></sub></td></tr><tr><td>2</td><td>1</td></tr><tr><td>3</td><td>1 + <sup>7 ζ<sup>2</sup></sup>&frasl;<sub>3<sup>3</sup></sub></td></tr><tr><td>4</td><td>1 + <sup>3 ζ<sup>2</sup></sup>&frasl;<sub>2<sup>3</sup></sub> + <sup>23 ζ<sup>4</sup></sup>&frasl;<sub>2<sup>8</sup>3</sub></td></tr><tr><td>5</td><td>1 + <sup>11 ζ<sup>2</sup></sup>&frasl;<sub>5<sup>2</sup></sub> + <sup>509 ζ<sup>4</sup></sup>&frasl;<sub>3 · 5<sup>5</sup></sub> + <sup>7 · 13 ζ<sup>6</sup></sup>&frasl;<sub>3 · 5<sup>6</sup></sub></td></tr></tbody></table>

asserts that lighter particles radiate more easily; it becomes 1 for $c_1 = -c_2$. The phase-space suppression factor $pss_{n,0}^{1/2}$ 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_{\text{BSF}}(\zeta, \xi)$ measures the overlap of the initial and final states. Its Coulomb limit (3.8c) is illustrative. The first factor, $S_0(\zeta) = 2\pi\zeta/(1 - e^{-2\pi\zeta})$, is an overall multiplicative constant in the scattering-state wavefunction [cf. eq. (C.1)], and is responsible for the characteristic scaling of long-range inelastic processes, $\sigma_{\text{BSF}} v_{\text{rel}} \propto 1/v_{\text{rel}}$ at $\zeta \gg 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 $\zeta < 1$, $S_{\text{BSF}}(\zeta, \xi)$ is very small and BSF is typically inefficient; however, at $\zeta \gtrsim 1$, $S_{\text{BSF,C}}(\zeta) \simeq 3.13 \times 2\pi\zeta$. 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.

### 3.3 Capture into $\ell = 1$ bound states

Similarly to before, the zeroth order terms in $|\mathbf{P}_{\varphi}|/\kappa$, of the $\mathcal{J}_{\mathbf{k},n1m}$ integrals yield the dominant contribution to the amplitude (3.4) for capture into $\ell = 1$ bound states, with the unpolarised squared amplitude (3.3) being

$$ \sum_{\epsilon} |\mathcal{M}_{\mathbf{k} \to n1m}|^2 \approx \frac{2^5 \pi \mu\alpha}{\eta_1^2 \eta_2^2} \left[ \frac{(\eta_2 c_1 - \eta_1 c_2)^2}{-c_1 c_2} \right] \left( |\mathcal{J}_{\mathbf{k},n1m}|^2 - pss_{n,0} |\hat{\mathbf{P}}_{\varphi} \cdot \mathcal{J}_{\mathbf{k},n1m}|^2 \right). \quad (3.9) $$
For convenience, we define

$$A_0(\zeta, \xi) \equiv \int_0^\infty dx \left[ \chi'_{n,1}(x) + \frac{\chi_{n,1}(x)}{x} \right]^* \chi_{|k|,0}(x), \quad (3.10a)$$

$$A_2(\zeta, \xi) \equiv \int_0^\infty dx \left[ \chi'_{n,1}(x) - \frac{2\chi_{n,1}(x)}{x} \right]^* \chi_{|k|,2}(x). \quad (3.10b)$$

Then, from eq. (B.10) we find

$$|\mathcal{J}_{\mathbf{k},n10}|^2 = \frac{4\pi}{\kappa} \left[ \frac{1}{3}(|A_0|^2 + |A_2|^2) + |A_2|^2 \cos^2 \theta_{\mathbf{k}} + (A_0 A_2^* + A_0^* A_2) \left( \cos^2 \theta_{\mathbf{k}} - \frac{1}{3} \right) \right], \quad (3.11a)$$

$$|\mathcal{J}_{\mathbf{k},n1\pm 1}|^2 = \frac{2\pi}{\kappa} \left[ \frac{2}{3}(|A_0|^2 + |A_2|^2) + |A_2|^2 \sin^2 \theta_{\mathbf{k}} + (A_0 A_2^* + A_0^* A_2) \left( \sin^2 \theta_{\mathbf{k}} - \frac{2}{3} \right) \right], \quad (3.11b)$$

and

$$|\mathcal{J}_{\mathbf{k},n10} \cdot \hat{\mathbf{P}}_{\varphi}|^2 = \frac{12\pi}{\kappa} \left[ \frac{|A_0 - A_2|^2}{9} \cos^2 \theta_{\mathbf{P}_{\varphi}} + (\hat{\mathbf{k}} \cdot \hat{\mathbf{P}}_{\varphi})^2 |A_2|^2 \cos^2 \theta_{\mathbf{k}} \right. \\ \left. + (\hat{\mathbf{k}} \cdot \hat{\mathbf{P}}_{\varphi}) \cos \theta_{\mathbf{k}} \cos \theta_{\mathbf{P}_{\varphi}} \frac{A_2(A_0^* - A_2^*) + \text{c.c.}}{3} \right], \quad (3.11c)$$

$$|\mathcal{J}_{\mathbf{k},n1\pm 1} \cdot \hat{\mathbf{P}}_{\varphi}|^2 = \frac{6\pi}{\kappa} \left[ \frac{|A_0 - A_2|^2}{9} \sin^2 \theta_{\mathbf{P}_{\varphi}} + (\hat{\mathbf{k}} \cdot \hat{\mathbf{P}}_{\varphi})^2 |A_2|^2 \sin^2 \theta_{\mathbf{k}} + \right. \\ \left. + (\hat{\mathbf{k}} \cdot \hat{\mathbf{P}}_{\varphi}) \sin \theta_{\mathbf{k}} \sin \theta_{\mathbf{P}_{\varphi}} \frac{A_2(A_0^* - A_2^*) \exp[i(\phi_{\mathbf{k}} - \phi_{\mathbf{P}_{\varphi}})] + \text{c.c.}}{3} \right]. \quad (3.11d)$$

To calculate the total cross-section for capture to any $\ell = 1$ state, we sum over all values of $m$.

$$\sum_{m=-1}^{1} \sum_{\epsilon} |\mathcal{M}_{\mathbf{k} \to n1m}|^2 \approx \frac{2^7 \pi^2}{\eta_1^2 \eta_2^2} \left[ \frac{(\eta_2 c_1 - \eta_1 c_2)^2}{-c_1 c_2} \right] \times \\ \times \left( |A_0|^2 + 2|A_2|^2 - pss_{n,1} \left[ \frac{|A_0 - A_2|^2}{3} + (|A_2|^2 + A_2 A_0^* + A_2^* A_0) (\hat{\mathbf{k}} \cdot \hat{\mathbf{P}}_{\varphi})^2 \right. \right. \\ \left. \left. - 2 \operatorname{Im}(A_2 A_0^*) \sin \theta_{\mathbf{k}} \sin \theta_{\mathbf{P}_{\varphi}} \sin(\phi_{\mathbf{k}} - \phi_{\mathbf{P}_{\varphi}}) (\hat{\mathbf{k}} \cdot \hat{\mathbf{P}}_{\varphi}) \right] \right). \quad (3.12)$$

Note that $\hat{\mathbf{k}} \cdot \hat{\mathbf{P}}_{\varphi} = \cos \theta_{\mathbf{k}} \cos \theta_{\mathbf{P}_{\varphi}} + \sin \theta_{\mathbf{k}} \sin \theta_{\mathbf{P}_{\varphi}} \cos(\phi_{\mathbf{k}} - \phi_{\mathbf{P}_{\varphi}})$, and the term proportional to $\operatorname{Im}(A_2 A_0^*)$ in eq. (3.12) gives a vanishing contribution when integrated over $d\Omega_{\mathbf{P}_{\varphi}}$. Using the above and eq. (2.6), we find the total cross-section for capture into an $\ell = 1$ state (for fixed $n$),

$$\sum_{m=-1}^{1} \sigma_{\text{BSF}}^{\{n1m\}} v_{\text{rel}} = \frac{\pi \alpha^2}{4\mu^2} \left[ \frac{(c_1\eta_2 - c_2\eta_1)^2}{-c_1c_2} \right] pss_{n,1}^{1/2} \left( \frac{3 - pss_{n,1}}{2} \right) \times S_{\text{BSF}}^{\{n1\}}(\zeta, \xi), \quad (3.13a)$$

where

$$S_{\text{BSF}}^{\{n1\}}(\zeta, \xi) = \frac{2^5}{3} \left( \frac{1 + \zeta^2 \gamma_{n,1}^2(\xi)}{\zeta^2} \right) [|A_0(\zeta, \xi)|^2 + 2|A_2(\zeta, \xi)|^2], \quad (3.13b)$$

with $A_0$ and $A_2$ defined in eq. (3.10).
## 3.4 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

$$ \sigma_{\text{ann}}^{s} v_{\text{rel}} = \frac{\pi \alpha^2}{2 \mu^2} \times S_{\text{ann}}^{(0)}(\zeta, \xi), \qquad (3.14\text{a}) $$

$$ \sigma_{\text{ann}}^{f} v_{\text{rel}} = \frac{\pi \alpha^2}{4 \mu^2} \times S_{\text{ann}}^{(0)}(\zeta, \xi), \qquad (3.14\text{b}) $$

where $\mu = m_X/2 = M/4$ and

$$ S_{\text{ann}}^{(0)}(\zeta, \xi) \equiv |\phi_{\mathbf{k}}(\mathbf{r} = 0)|^2 = \lim_{\mathbf{x} \to 0} \left[ \frac{\chi_{|\mathbf{k}|, \ell=0}(\mathbf{x})}{\mathbf{x}} \right]. \qquad (3.14\text{c}) $$

The Coulomb limit of eq. (3.14c) is [cf. eq. (C.1)]

$$ S_{\text{ann,C}}^{(0)}(\zeta) = \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \equiv S_0(\zeta). \qquad (3.14\text{d}) $$

Note that the perturbative value of the annihilation cross-section (recovered in the limit $S_{\text{ann}}^{(0)} \to 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 figures 2 and 3, and offer our comments in the next section.

## 3.5 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 $\zeta$, the BSF cross-sections scale as $\sigma v_{\text{rel}} \propto \zeta$. This behaviour is realised at $\zeta \gtrsim n$ for capture into an $\{nlm\}$ bound state, as evident in figure 2; for capture into $\ell = 0$ bound states, it can also be confirmed analytically using eq. (3.8c). Importantly, the $1/v_{\text{rel}}$ scaling is anticipated due to the upper bound on the inelastic cross-sections imposed by unitarity [1, 7, 8], 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 $\ell \neq 0$ states (summed over $-\ell \le m \le \ell$) dominates over capture into the $\ell = 0$ state, in part due to the larger multiplicity of the former.

The above points imply that for a given value of $\zeta$, the total BSF cross-section is dominated by capture into the $n < \zeta$ levels, with the $\ell > 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, 9], described by Kramer's formula,

$$ \sigma_{\text{BSF}}^{\text{tot}} v_{\text{rel}} \approx \frac{\pi \alpha^2}{4 \mu^2} \times \frac{2^7}{3\sqrt{3}} \zeta [\ln \zeta + 0.16 + \mathcal{O}(\zeta^{-1})], \qquad (3.15) $$
which is valid for $\zeta \gtrsim 2$. The enhancement with respect to capture into the ground state [cf. eqs. (3.8a) and (3.8d)], is $\sigma_{\text{BSF}}^{\text{tot}}/\sigma_{\text{BSF}}^{\{100\}} \simeq 1 + 1.25 \ln(\zeta/1.89)$.

However, the enhanced BSF rate implied by eq. (3.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 [10], and their survival until today, depends typically on a rather complex interplay between formation, ionisation, excitation and de-excitation processes [11]. Moreover, the detectability of the radiation emitted inside halos today during the formation of stable bound states of asymmetric DM [12–14], 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 [cf. eq. (2.5b)] that becomes more severe for capture into excited states [2].

The Coulomb limit is attained at $\xi \gtrsim \zeta$, as can be observed in figures 3 to 6, 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, $\mu v_{\text{rel}} \gtrsim m_{\phi}$.

**Resonance structure.** The BSF cross-sections exhibit a rich resonance structure, as seen in figures 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, $\mathcal{E}_{n\ell} \to 0$, and thus overlaps with the branch-cut, which lives at $\mathcal{E}_k \ge 0$ (see e.g. ref. [15, section 7.7]). The locations of resonances therefore denote the thresholds for the existence of genuine bound states with $\mathcal{E}_{n\ell} < 0$. These thresholds signify that if the mediator is massive, the potential has to be sufficiently strong for bound states to exist. For every bound-state energy level, this implies a minimum value for $\xi \equiv \mu\alpha/m_\phi$. In contrast, for a massless mediator, bound states exist independently of the strength of the coupling, as is the case with QED.

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 $\ell$ quantum numbers (cf. appendix A.2). As a result, the resonances that appear in the scattering-state wavefunction depend on the $\ell$ mode.

Angular momentum conservation implies that the $\ell$ 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 (cf. appendix B.4). The formation of $\ell = 0$ bound states with emission of a vector boson is dominated by the $\ell = 1$ mode of the scattering-state wavefunction [cf. eq. (3.7b)]. On the other hand, the formation of $\ell = 1$ bound states with vector emission is dominated by the $\ell = 0$ and $\ell = 2$ modes of the scattering-state wavefunction [cf. eqs. (3.13b) and (3.10)]. This explains the resonant patterns observed in figures 3a, 4 and 6.

In figure 5a, we present the ratio of the factors $S_{\text{BSF}}^{\{100\}}$ and $S_{\text{BSF}}^{\{200\}}$ to $S_{\text{ann}}^{(1)}$, the Sommerfeld enhancement factor of p-wave annihilation processes,

$$S_{\text{ann}}^{(1)} = \left| \frac{3\zeta}{2\kappa} \frac{d}{dr} \int_{-1}^{1} d(\cos\theta_{\mathbf{r}}) P_{1}(\cos\theta_{\mathbf{r}}) \phi_{\mathbf{k}}(\mathbf{r}) \right|_{\mathbf{r}=0}^{2} = 9\zeta^{2} \left| \frac{\chi_{|\mathbf{k}|,1}(x)}{x^{2}} \right|_{x \to 0}^{2} . \quad (3.16)$$

As seen in eq. (3.16), $S_{\text{ann}}^{(1)}$ depends only on the $\ell = 1$ mode of the scattering-state wavefunction; in contrast $S_{\text{BSF}}^{\{100\}}$ and $S_{\text{BSF}}^{\{200\}}$ depend on the overlap of the wavefunctions of the scattering state, the bound state and the emitted vector boson, as seen in eq. (2.7). The smoothness of the curves in figure 5a attests that the resonances in the BSF cross-sections emanate from the scattering state. However, the capture to the $\{200\}$ bound state exhibits also anti-resonances (cf. figure 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_{\text{rel}} < m_{\varphi}/\mu$, the inelastic cross-sections depart from the Coulombic scaling $\sigma_{\text{inel}}v_{\text{rel}} \propto 1/v_{\text{rel}}$. Their velocity dependence in this regime, is determined by the $\ell$ modes of the scattering-state wavefunction that participate in the corresponding processes.

For $\xi$ values away from resonances, the contribution of an $\ell$ mode of the scattering-state wavefunction to $\sigma_{\text{inel}}v_{\text{rel}}$ scales as $v_{\text{rel}}^{2\ell}$ at $v_{\text{rel}} < m_{\varphi}/\mu$. That is, as the velocity decreases, a contribution from the $\ell = 0$ mode saturates to its Coulomb value at $v_{\text{rel}} \approx m_{\varphi}/\mu$, while contributions from $\ell > 0$ modes drop below their Coulomb value at $v_{\text{rel}} \approx m_{\varphi}/\mu$. The radiative capture into $\{n00\}$ bound states is dominated by the $\ell = 1$ mode of the scattering-state wavefunction, and scales as $v_{\text{rel}}^{2}$ at $v_{\text{rel}} < m_{\varphi}/\mu$, as can been seen in figures 3b and 4 (bottom left panel).

For $\xi$ values near or on resonances, $\sigma_{\text{inel}}v_{\text{rel}}$ grows faster than $1/v_{\text{rel}}$ at $v_{\text{rel}} < m_{\varphi}/\mu$, and raises above its Coulomb value for a given velocity. This growth is sustained for a range of velocities that depends on how close $\xi$ is to a resonance value. Then, at sufficiently low velocities, the resonant growth stops, and the $v_{\text{rel}}^{2\ell}$ scaling ensues.<fnref n="5" /> This behaviour can be observed in figures 3c and 4 (bottom right panel).

If more than one $\ell$ 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 $\ell = 1$ bound states, which receives contributions from the $\ell = 0$ and $\ell = 2$ modes of the scattering-state wavefunction (cf. figure 6).

<footnote n="5">⁵For $\ell = 0$ and $\xi$ values exactly on resonance, $\sigma_{\text{inel}}v_{\text{rel}}$ grows as $1/v_{\text{rel}}^2$ at $v_{\text{rel}} \lesssim m_{\varphi}/\mu$, indefinitely. This is unphysical behaviour that needs to be regulated (see comments on unitarity below).</footnote>
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 $\ell = 1$ bound states, if kinematically allowed, should always dominate, since these are the only mono-photon transitions to which the $\ell = 0$ mode of the scattering-state wavefunction participates.

**Near-threshold behaviour.** As mentioned above, $\xi$ 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 $\xi$ values close to the threshold for the existence of the same bound state, is suppressed [cf. figures 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 $\zeta \ll 1$, the BSF processes are rather suppressed. However, at $\zeta \gtrsim 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),

$$ \frac{\sigma_{\text{BSF,C}}^{\{100\}}}{\sigma_{\text{ann,C}}^{f}} = \frac{2\sigma_{\text{BSF,C}}^{\{100\}}}{\sigma_{\text{ann,C}}^{s}} = \frac{2^9 \zeta^4}{3 (1 + \zeta^2)^2} e^{-4\zeta\arccot\zeta}. \qquad (3.17) $$

At $\zeta \gtrsim 1$, this becomes $\sigma_{\text{BSF,C}}^{\{100\}}/\sigma_{\text{ann,C}}^{f} = 2\sigma_{\text{BSF,C}}^{\{100\}}/\sigma_{\text{ann,C}}^{s} \simeq 3.13$. In fact, the radiative formation of $n = 2, \ell = 1$ bound states is also faster than the annihilation of a fermion-antifermion pair, for $\zeta \gtrsim 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<fnref n="6" /> [7] and enhance the indirect detection signals [2, 16–18] of symmetric DM.

<footnote n="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_{\text{DM}} \gtrsim 15$ TeV (and up to a factor of 4 for $m_{\text{DM}} \gtrsim 100$ 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_{\text{DM}} \gtrsim 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_{\text{DM}} \lesssim 20$ TeV.</footnote>
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 ($\xi$ values), as clearly seen in figure 3a. This is due to the different $\ell$ modes of the scattering-state wavefunction that contribute to each process: $\ell = 1$ for the former and $\ell = 0$ for the latter. In fact, the locations of the $\ell = 1$ resonances exhibit a mild velocity dependence, in contrast to the $\ell = 0$ resonances.

* The $\ell$ modes of the scattering-state wavefunction also determine the velocity dependence away from the Coulomb limit. At low velocities, $\sigma_{\text{ann}}v_{\text{rel}}$ saturates to a constant value, while $\sigma_{\text{BSF}}^{100}v_{\text{rel}}$ scales as $v_{\text{rel}}^2$. 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 figures 3b and 3c. This range is roughly $1 \lesssim \zeta \lesssim \xi$ (or equivalently $m_{\phi} \lesssim \mu v_{\text{rel}} \lesssim \mu\alpha$) for non-resonant $\xi$ values, but it may extend to much lower velocities (by orders of magnitude) for $\xi$ values near an $\ell = 1$ resonance.

* Near the threshold for the existence of the ground state, i.e. for $\xi \approx 1$, the radiative BSF is always suppressed with respect to its Coulomb value, while annihilation is on resonance. This behaviour can be seen in figure 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.)

**Partial-wave unitarity.** Partial-wave unitarity implies an upper bound on the inelastic cross-sections, which in the non-relativistic regime is [19],

$$ (\sigma_{\text{uni}})_J v_{\text{rel}} = \frac{(2J+1)\pi}{\mu^2 v_{\text{rel}}}, \qquad (3.18) $$

where $J$ is the partial wave.

In the Coulomb limit, inelastic cross-sections have the same velocity scaling as eq. (3.18); setting $\sigma_{\text{inel}} \le \sigma_{\text{uni}}$ then implies an upper bound on $\alpha$ that does not depend on any other physical parameter, and is roughly $\alpha \lesssim 0.85$ [7, 8]. Around this upper bound on $\alpha$, higher order corrections (of perturbative or non-perturbative origin) need to be considered. Notably, $\alpha \sim 0.85$ is well below the naive perturbativity limit, $\alpha \sim 4\pi$.

However, the resonances that appear away from the Coulomb limit, imply that the leading-order computations presented here, may violate the unitarity bound even at much lower values of $\alpha$ that now also depend on $\zeta$ and $\xi$. Evidently, the peaks of the resonances can be unphysical.

The resonances can be regulated by taking into account the short-range elastic scattering of the interacting particles (see e.g. [20, 21]). For particle-antiparticle pairs, the short-range inelastic scattering — in particular the annihilation processes — may also contribute to taming the unphysical behaviour [22].
<figure><img src="image_4.png" /><figcaption>**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 $\sigma_0^f \equiv \pi\alpha^2/(4\mu^2)$, where $\mu$ is the reduced mass of the interacting particles. $\sigma_0^f$ is the spin-averaged perturbative annihilation cross-section times relative velocity of a Dirac fermion-antifermion pair; for a complex scalar, this quantity is $\sigma_0^s = 2\sigma_0^f$. 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 $\zeta = \alpha/v_{rel} \gg 1$, the radiative capture into the ground state is the dominant inelastic process, $\sigma_{\text{BSF}}^{\{100\}}/\sigma_{\text{ann}}^f \simeq 2\sigma_{\text{BSF}}^{\{100\}}/\sigma_{\text{ann}}^s \simeq 3.13$ [cf. eqs. (3.8), (3.14) and table 3].</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Vector mediator: resonances</figcaption></figure>

**Figure 3a.** Resonances due to bound states at threshold (zero binding energy) appear at discrete values of $\xi = \mu\alpha/m_\phi$, with a magnitude that increases at low energies (large $\zeta = \alpha/v_{\text{rel}}$). They arise from the wavefunction of the (initial) two-particle scattering state, and depend on the $\ell$ modes that contribute to the process of interest. Top left: the factor $S_{\text{BSF}}^{\{100\}}(\zeta, \xi)$, appearing in the cross-section for capture into the ground state, with emission of a vector force mediator. The $\ell = 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 $\zeta$ (velocity) dependence [cf. eqs. (3.7), (3.8)]. Top right: the Sommerfeld enhancement factor $S_{\text{ann}}^{(0)}$ of the annihilation of a particle-antiparticle pair into two dark photons. The $\ell = 0$ mode of the scattering state wavefunction yields the dominant contribution to this process. The resonances are less sharp than those of the $\ell = 1$ mode, and their location does not depend on $\zeta$ [cf. eqs. (3.14)]. Bottom: comparison of $S_{\text{BSF}}^{\{100\}}$ (solid lines), $S_{\text{ann}}^{(0)}$ (dashed lines), $2S_{\text{ann}}^{(0)}$ (dotted lines), for values of $\zeta$ shown in the top panels. For a particle-antiparticle pair, this is a comparison between the strength of the radiative capture into the ground state (ignoring the phase-space suppression), and the spin-averaged annihilation of a pair of fermions, or the annihilation of a pair of scalars, respectively. The location and the magnitude of the $S_{\text{BSF}}^{\{100\}}$ and $S_{\text{ann}}^{(0)}$ resonances is clearly seen to be different.
<figure><img src="image_6.png" /><figcaption>Vector mediator: $\xi$ values away from $\ell = 0$ and $\ell = 1$ resonances</figcaption></figure>

**Figure 3b.** Velocity dependence of $S_{\text{BSF}}^{\{100\}}$ and $S_{\text{ann}}^{(0)}$, for non-resonant values of $\xi = \mu\alpha/m_\varphi$. When the momentum transfer between the interacting particles drops below the mediator mass, $\mu v_{\text{rel}} \lesssim m_\varphi$ (equivalently, $\zeta \gtrsim \xi$), the velocity dependence of the cross-sections departs from the Coulombic behaviour. **Top left:** the cross-section for radiative capture into the ground state with emission of a vector mediator drops approximately as $\zeta^{-2} \propto v_{\text{rel}}^2$ (see also figure 5b). **Top right:** for s-wave annihilation processes, the Sommerfeld enhancement saturates to a constant value. **Bottom:** comparison of the velocity dependence of $S_{\text{BSF}}^{\{100\}}$ (solid lines), $S_{\text{ann}}^{(0)}$ (dashed lines), $2S_{\text{ann}}^{(0)}$ (dotted lines), for various values of $\xi$. For a non-zero mediator mass, and for a particle-antiparticle pair, the bound-state formation cross-section is larger than the annihilation cross-section within a finite range of velocities, roughly $m_\varphi \lesssim \mu v_{\text{rel}} \lesssim \mu\alpha$, or equivalently $1 \lesssim \zeta \lesssim \xi$. (This comparison ignores the phase-space suppression of the bound-state formation process.)
Vector mediator

<figure><img src="image_7.png" /><figcaption>$\xi$ values near the $n = 2, \ell = 0$ resonance</figcaption></figure>

<figure><img src="image_8.png" /></figure>

**Figure 3c.** Velocity dependence of $S_{\text{BSF}}^{\{100\}}$ and $S_{\text{ann}}^{(0)}$, for values of $\xi = \mu\alpha/m_{\varphi}$ 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, $\mu v_{\text{rel}} \lesssim m_{\varphi}$ ($\zeta \gtrsim \xi$). In contrast to non-resonant $\xi$ values though, for $\xi$ near a resonance and at $\zeta \gtrsim \xi$, 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 $\zeta^{-2} \propto v_{\text{rel}}^2$ at $\zeta \gg \xi$ (see also figure 5b). **Top right:** for s-wave annihilation processes, the Sommerfeld enhancement rises monotonically with $\zeta$; at $\zeta \gg \xi$, it asymptotes to a value that can be much larger than its value at $\zeta \approx \xi$. **Bottom:** comparison of the velocity dependence of $S_{\text{BSF}}^{\{100\}}$ (solid lines), $S_{\text{ann}}^{(0)}$ (dashed lines), $2S_{\text{ann}}^{(0)}$ (dotted lines), for various values of $\xi$ near resonances arising in the $\ell = 1$ (left) and $\ell = 0$ (right) modes of the scattering state wavefunction.
<figure><img src="image_9.png" /><figcaption>Vector mediator: $\xi$ values near the $n = 1, \ell = 0$ threshold/resonance</figcaption></figure>

**Figure 3d.** *Left:* for $\xi$ values close to the threshold for the existence of a bound state, the cross-section for the radiative formation of this bound state is suppressed with respect to its Coulomb value, for all velocities. This suppression arises from the bound-state wavefunction, which becomes very extended for $\xi$ values near threshold, and is independent of the phase-space suppression due to the emission of a massive force mediator. *Right:* the same values of $\xi$ yield a resonance in the annihilation processes.
<figure><img src="image_10.png" /><figcaption>Vector mediator: capture into first excited state {200}</figcaption></figure>

**Figure 4.** Resonance structure and velocity dependence of the $S_{\text{BSF}}^{\{200\}}$ factor of the cross-section for the radiative formation of the first excited state {200}, with emission of a vector force mediator [cf. eqs. (3.7)]. In the bottom right panel, we show the velocity dependence near the $n = 3, \ell = 1$ resonance. Besides the resonances, which we have already seen in the cross-section for capture to the ground state, the formation of excited states may feature anti-resonances.
<figure><img src="image_11.png" /><figcaption>Vector mediator:
Formation of zero angular momentum bound states vs. p-wave annihilation</figcaption></figure>

**Figure 5a.** The ratios $S_{\text{BSF}}^{\{100\}}/S_{\text{ann}}^{(1)}$ (left) and $S_{\text{BSF}}^{\{200\}}/S_{\text{ann}}^{(1)}$ (right), where $S_{\text{ann}}^{(1)}(\zeta, \xi)$ is the Sommerfeld enhancement factor of p-wave annihilation processes, which depends on the $\ell = 1$ component of the scattering state wavefunction only [cf. eq. (3.16)]. The absence of any resonances in these ratios implies that the resonances in $S_{\text{BSF}}^{\{100\}}$ and $S_{\text{BSF}}^{\{200\}}$ arise solely from the scattering-state wavefunction. On the other hand, the anti-resonances arise from the interference of the scattering-state and bound-state wavefunctions.

<figure><img src="image_12.png" /></figure>

**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 $\zeta^{-2}S_{\text{ann}}^{(1)}$ captures the entire velocity dependence of the latter. At large $\zeta$, the ratio $S_{\text{BSF}}^{\{n00\}}/[\zeta^{-2}S_{\text{ann}}^{(1)}]$ tends to a $\xi$-dependent constant, which itself saturates to a fixed value (its Coulomb limit) at large $\xi$. For the ground state and first excited state, $S_{\text{BSF}}^{\{100\}}/[\zeta^{-2}S_{\text{ann}}^{(1)}] \simeq 2^9/(3e^4) \simeq 3.13$ and $S_{\text{BSF}}^{\{200\}}/[\zeta^{-2}S_{\text{ann}}^{(1)}] \simeq 2^{12}/(3e^8) \simeq 0.46$ [cf. eq. (3.8c) and table 3].
Vector mediator: capture into states with non-zero angular momentum
$n = 2, \ell = 1$

<figure><img src="image_13.png" /><figcaption>**Figure 6.** Resonance structure and velocity dependence of the $S_{\text{BSF}}^{\{21\}}$ factor of the total cross-section for the radiative capture into any $n = 2, \ell = 1$ state, with emission of a vector force mediator [cf. eqs. (3.13)]. Left: $S_{\text{BSF}}^{\{21\}}$ receives contributions from the $\ell = 0$ and the $\ell = 2$ modes of the scattering state wavefunction. The $\ell = 0$ mode dominates for small $\xi$, while the $\ell = 2$ mode dominates for larger $\xi$ values (including the Coulomb regime) and gives rise to sharper resonances. Right: at large $\zeta$, the contribution from the $\ell = 2$ mode decreases with decreasing velocity as $v_{\text{rel}}^4$, while the contribution from the $\ell = 0$ mode saturates to a constant value and dominates. The superposition of two $\ell$ modes — which have different velocity dependence and resonances — gives rise to the various features of the total cross-section.</figcaption></figure>

## 4 Scalar force mediator

We now consider the interaction Lagrangians

$$ \mathcal{L} = \frac{1}{2} \partial_{\mu} X_{1} \partial^{\mu} X_{1} + \frac{1}{2} \partial_{\mu} X_{2} \partial^{\mu} X_{2} + \frac{1}{2} \partial_{\mu} \varphi \partial^{\mu} \varphi - \frac{1}{2} m_{1}^{2} X_{1}^{2} - \frac{1}{2} m_{2}^{2} X_{2}^{2} - \frac{1}{2} m_{\varphi}^{2} \varphi^{2} \\ - \frac{1}{2} g_{1} m_{1} \varphi X_{1}^{2} - \frac{1}{2} g_{2} m_{2} \varphi X_{2}^{2}, \qquad (4.1a) $$

$$ \mathcal{L} = \partial_{\mu} X_{1}^{\dagger} \partial^{\mu} X_{1} + \partial_{\mu} X_{2}^{\dagger} \partial^{\mu} X_{2} + \frac{1}{2} \partial_{\mu} \varphi \partial^{\mu} \varphi - m_{1}^{2} |X_{1}|^{2} - m_{2}^{2} |X_{2}|^{2} - \frac{1}{2} m_{\varphi}^{2} \varphi^{2} \\ - g_{1} m_{1} \varphi |X_{1}|^{2} - g_{2} m_{2} \varphi |X_{2}|^{2}, \qquad (4.1b) $$

$$ \mathcal{L} = \bar{X}_1 i \not D X_1 + \bar{X}_2 i \not D X_2 + \frac{1}{2} \partial_{\mu} \varphi \partial^{\mu} \varphi - m_1 \bar{X}_1 X_1 - m_2 \bar{X}_2 X_2 - \frac{1}{2} m_{\varphi}^2 \varphi^2 \\ - g_1 \varphi \bar{X}_1 X_1 - g_2 \varphi \bar{X}_2 X_2. \qquad (4.1c) $$

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. $\varphi$ is a real scalar boson, and $g_1, g_2$ are dimensionless couplings. The interaction between $X_1, X_2$ via $\varphi$ exchange is described by the Yukawa potential of eq. (1.1), with $\alpha = \alpha_{sc}$ or $\alpha = \alpha_f$ depending on whether the interacting particles are scalars
or fermions respectively, where [1]

$$ \alpha_{sc} = \frac{g_1 g_2}{16\pi} \quad \text{and} \quad \alpha_f = \frac{g_1 g_2}{4\pi}. \qquad (4.2) $$

The interaction is attractive if $g_1g_2 > 0$. In the following, the parameters $\zeta \equiv \alpha/v_{\text{rel}}$ and $\xi = \alpha\mu/m_{\varphi}$ are always defined using the appropriate $\alpha$.

## 4.1 Pairs of non-degenerate particles

The BSF amplitude is [1]

$$ \mathcal{M}_{\mathbf{k} \to n\ell m} \approx -M\sqrt{2\mu} [g_1 \mathcal{I}_{\mathbf{k},n\ell m}(\eta_2 \mathbf{P}_{\varphi}) + g_2 \mathcal{I}_{\mathbf{k},n\ell m}(-\eta_1 \mathbf{P}_{\varphi})]. \qquad (4.3) $$

The dominant contributions to the amplitudes $\mathcal{M}_{\mathbf{k} \to n00}$ and $\mathcal{M}_{\mathbf{k} \to n1m}$ that we will consider below, are of order $|\mathbf{P}_{\varphi}|/\kappa$.

### 4.1.1 Capture into $\ell = 0$ bound states

Using eqs. (4.3) and (B.9a), we find

$$ |\mathcal{M}_{\mathbf{k} \to n00}|^2 \approx \frac{2^5 \pi^2 M^2}{\mu^2} \left[ \frac{(g_1 \eta_2 - g_2 \eta_1)^2}{16\pi\alpha} \right] pss_{n,0} \times \\ \times \left( \frac{1 + \zeta^2 \gamma_{n,0}^2}{\zeta^2} \right)^2 \left| \int_0^\infty dx \, x \, \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,1}(x) \right|^2 \cos^2 \theta. \qquad (4.4) $$

In the limit $g_1 = g_2$, $\eta_1 \gg \eta_2$, the factor in the square brackets is $(g_1\eta_2 - g_2\eta_1)^2/(16\pi\alpha) = 1$ or $1/4$, for a bosonic and fermionic pair respectively. From eqs. (2.6) and (4.4), we find the corresponding cross-section to be

$$ \sigma_{\text{BSF}}^{\{n00\}} v_{\text{rel}} \approx \frac{\pi \alpha^2}{\mu^2} \left[ \frac{(g_1 \eta_2 - g_2 \eta_1)^2}{16\pi\alpha} \right] pss_{n,0}^{3/2} \times S_{\text{BSF}}^{\{n00\}}(\zeta, \xi), \qquad (4.5a) $$

where

$$ S_{\text{BSF}}^{\{n00\}}(\zeta, \xi) = \frac{1}{3} \left( \frac{1 + \zeta^2 \gamma_{n,0}^2(\xi)}{\zeta^2} \right)^3 \left| \int_0^\infty dx \, x \, \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,1}(x) \right|^2. \qquad (4.5b) $$

**Coulomb limit.** The cross-section (4.5) becomes

$$ \lim_{\xi \to \infty} \sigma_{\text{BSF}}^{\{n00\}} v_{\text{rel}} \approx \frac{\pi \alpha^2}{\mu^2} \left[ \frac{(g_1 \eta_2 - g_2 \eta_1)^2}{16\pi\alpha} \right] S_{\text{BSF,C}}^{\{n00\}}(\zeta, \xi), \qquad (4.6a) $$

where from eq. (C.5a), and keeping terms of order $|\mathbf{P}_{\varphi}|/\kappa$, we find

$$ S_{\text{BSF,C}}^{\{n00\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \frac{2^6}{3n^5} \zeta^2 (1 + \zeta^2) (1 + \zeta^2/n^2)^3 \times \\ \times \left| \sum_{s=0}^{n-1} \frac{n! 2^s \zeta_n^s}{(n-s-1)!(s+1)! s!} \frac{d^s}{d\zeta_n^s} \left[ \frac{(\zeta - 2\zeta_n) e^{-2\zeta} \operatorname{arccot} \zeta_n}{(1 + \zeta_n^2)^3} \right] \right|^2_{\zeta_n = \zeta/n} . \qquad (4.6b) $$
Expanding the sum, eq. (4.6b) becomes

$$S_{\text{BSF,C}}^{\{n00\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \frac{2^6}{3n^3} \frac{\zeta^4(1 + \zeta^2)}{(1 + \zeta^2/n^2)^{2n-1}} [\varrho_n(\zeta)]^2 e^{-4\zeta \arccot(\zeta/n)}, \quad (4.6c)$$

where $\varrho_n(\zeta)$ is a rational function of $\zeta^2$ with $\lim_{\zeta \to 0} \varrho_n(\zeta) = 1$. The explicit expressions of $\varrho_n(\zeta)$ for $1 \le n \le 5$ are given in table 3. For $n=1$,

$$S_{\text{BSF,C}}^{\{100\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \times \frac{2^6}{3} \left( \frac{\zeta^2}{1 + \zeta^2} \right)^2 e^{-4\zeta \arccot \zeta}, \quad (4.6d)$$

We inspect the features of $S_{\text{BSF}}^{\{100\}}(\zeta, \xi)$ and $S_{\text{BSF,C}}^{\{100\}}(\zeta)$ in figure 7.

### 4.1.2 Capture into $\ell = 1$ bound states

Similarly to above, we combine eqs. (2.6), (4.3) and (B.5a), and find the BSF cross-section to be

$$\sigma_{\text{BSF}}^{\{n1\}} v_{\text{rel}} \equiv \sum_{m=-1}^{1} \sigma_{\text{BSF}}^{\{n1m\}} v_{\text{rel}} \approx \frac{\pi \alpha^2}{\mu^2} \left[ \frac{(g_1 \eta_2 - g_2 \eta_1)^2}{16 \pi \alpha} \right] p s_{n,1}^{3/2} \times S_{\text{BSF}}^{\{n1\}}(\zeta, \xi), \quad (4.7a)$$

where

$$S_{\text{BSF}}^{\{n1\}}(\zeta, \xi) = \frac{1}{3} \left( \frac{1 + \zeta^2 \gamma_{n,1}^2(\xi)}{\zeta^2} \right)^3 \times \left( \left| \int_0^\infty dx \, x \, \chi_{n,1}^*(x) \chi_{|k|,0}(x) \right|^2 + 2 \left| \int_0^\infty dx \, x \, \chi_{n,1}^*(x) \chi_{|k|,2}(x) \right|^2 \right). \quad (4.7b)$$

We inspect the features of $S_{\text{BSF}}^{\{21\}}(\zeta, \xi)$ in figure 8.
<figure><img src="image_14.png" /><figcaption>Scalar mediator, non-degenerate particles</figcaption></figure>

**Figure 7.** Resonance structure (upper two panels) and velocity dependence (lower four panels) of the $S_{\text{BSF}}$ factors determining the cross-sections for the formation of bound states with emission of a scalar force mediator, by two particles with different masses and/or couplings [cf. eqs. (4.5) and (4.6)].
Scalar mediator, non-degenerate particles: capture into $n = 2, \ell = 1$ bound states

<figure><img src="image_15.png" /><figcaption>**Figure 8.** Resonance structure and velocity dependence of the $S_{\text{BSF}}^{\{21\}}$ factor of the total cross-section for the radiative capture into any $n = 2, \ell = 1$ state of two particles with different masses and/or couplings, with emission of a scalar force mediator [cf. eqs. (4.7)]. Left: $S_{\text{BSF}}^{\{21\}}$ receives contributions from the $\ell = 0$ and the $\ell = 2$ modes of the scattering-state wavefunction. The $\ell = 0$ mode dominates for small $\xi$, while the $\ell = 2$ mode dominates for larger $\xi$ values (including the Coulomb regime) and gives rise to sharper resonances. Right: at large $\zeta$, the contribution from the $\ell = 2$ mode decreases with decreasing velocity as $v_{\text{rel}}^4$, while the contribution from the $\ell = 0$ mode saturates to a constant value and dominates.</figcaption></figure>

## 4.2 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$, $\eta_1 = \eta_2 = 1/2$ and $\mu = 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 figure 1b that remains when the incoming and outgoing ladders are amputated).<fnref n="7" /> Then, the BSF amplitude for a particle-antiparticle pair of non-

<footnote n="7">⁷We believe that the difference between the bosonic case considered here, and the fermionic case considered in ref. [3], is the corrections due to (ii), which are not present in the latter case.</footnote>
self-conjugate bosons is [1]

$$ \mathcal{M}_{\mathbf{k} \to n\ell m} \approx -M\sqrt{2\mu} \left[ g_1 \mathcal{I}_{\mathbf{k},n\ell m}(\eta_2 \mathbf{P}_{\varphi}) + g_2 \mathcal{I}_{\mathbf{k},n\ell m}(-\eta_1 \mathbf{P}_{\varphi}) + \frac{g_1 \mathcal{K}_{\mathbf{k},n\ell m}(\eta_2 \mathbf{P}_{\varphi}) + g_2 \mathcal{K}_{\mathbf{k},n\ell m}(-\eta_1 \mathbf{P}_{\varphi})}{2M\mu} \right]. \quad (4.8) $$

### 4.2.1 Capture into $\ell = 0$ bound states

Keeping the terms of order $(|\mathbf{P}_{\varphi}|/\kappa)^2$ from the $\mathcal{I}_{\mathbf{k},n00}$ integral [cf. eq. (B.9a)] and the zeroth-order terms from the $\mathcal{K}_{\mathbf{k},n00}$ integral [cf. eq. (B.9c)], the amplitude (4.8) becomes

$$ |\mathcal{M}_{\mathbf{k} \to n00}|^2 \approx 2^7 \pi^2 \alpha^2 \times \left| - \int_0^\infty dx \left( \frac{2e^{-x/\xi}}{x} \right) \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,0}(x) \right. \\ \left. + \frac{pss_{n,0}}{12} \left( \frac{1 + \zeta^2 \gamma_{n,0}^2(\xi)}{\zeta^2} \right)^2 \int_0^\infty dx \, x^2 \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,0}(x) \right. \\ \left. + P_2(\hat{\mathbf{k}} \cdot \hat{\mathbf{P}}_{\varphi}) \frac{pss_{n,0}}{6} \left( \frac{1 + \zeta^2 \gamma_{n,0}^2(\xi)}{\zeta^2} \right)^2 \int_0^\infty dx \, x^2 \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,2}(x) \right|^2. \quad (4.9) $$

Note that here, $\mathcal{M}_{\mathbf{k} \to n00}$ receives its dominant contribution from the $\ell = 0$ and $\ell = 2$ components of the scattering state wavefunction. This is in contrast to the formation of zero angular momentum bound states via vector emission [cf. eq. (3.6)], or via scalar emission but by non-degenerate particles [cf. eq. (4.4)], where the $\ell = 1$ mode of the scattering-state wavefunction dominates.

Substituting eq. (4.9) into eq. (2.6), we find the corresponding cross-section,

$$ \sigma_{\text{BSF}}^{\{n00\}} v_{\text{rel}} \approx \left(\frac{\pi\alpha^2}{\mu^2}\right) \alpha^2 S_{\text{BSF}}^{\{n00\}}(\zeta, \xi; \alpha) pss_{n,0}^{1/2}, \quad (4.10a) $$

where $S_{\text{BSF}}^{\{n00\}}$ includes terms suppressed by different powers of $pss_{n,0}$,

$$ S_{\text{BSF}}^{\{n00\}}(\zeta, \xi; \alpha) = \Sigma_{0}^{\{n00\}}(\zeta, \xi) + pss_{n,0} \times \Sigma_{1}^{\{n00\}}(\zeta, \xi) + pss_{n,0}^2 \times \Sigma_{2}^{\{n00\}}(\zeta, \xi), \quad (4.10b) $$

with

$$ \Sigma_{0}^{\{n00\}}(\zeta, \xi) = \frac{1}{4} \left( \frac{1 + \zeta^2 \gamma_{n,0}^2(\xi)}{\zeta^2} \right) \times \left| \int_{0}^{\infty} dx \left( \frac{2e^{-x/\xi}}{x} \right) \chi_{n,0}^{*}(x) \chi_{|k|,0}(x) \right|^2, \quad (4.10c) $$

$$ \Sigma_{1}^{\{n00\}}(\zeta, \xi) = -\frac{1}{2^4 \cdot 3} \left( \frac{1 + \zeta^2 \gamma_{n,0}^2(\xi)}{\zeta^2} \right)^3 \times \left( \left[ \int_{0}^{\infty} dx \left( \frac{2e^{-x/\xi}}{x} \right) \chi_{n,0}^{*}(x) \chi_{|k|,0}(x) \right] \right. \\ \left. \times \left[ \int_{0}^{\infty} dx \, x^2 \chi_{n,0}(x) \chi_{|k|,0}^{*}(x) \right] + \text{c.c.} \right), \quad (4.10d) $$
$$ \Sigma_2^{\{n00\}}(\zeta, \xi) = \frac{1}{2^4 3^2} \left( \frac{1 + \zeta^2 \gamma_{n,0}^2(\xi)}{\zeta^2} \right)^5 \times \\ \times \left( \frac{1}{4} \left| \int_0^\infty dx \, x^2 \chi_{n,0}^*(x) \chi_{|k|,0}(x) \right|^2 \right. \\ \left. + \frac{1}{5} \left| \int_0^\infty dx \, x^2 \chi_{n,0}^*(x) \chi_{|k|,2}(x) \right|^2 \right). \qquad (4.10e) $$

Note that the dependence of $S_{\text{BSF}}^{\{n00\}}(\zeta, \xi; \alpha)$ on $\alpha$ (independently of $\zeta$ and $\xi$) arises from the phase-space suppression factor $pss_{n,0}$ in eq. (4.10b).

**Coulomb limit.** The BSF cross-section (4.10) becomes

$$ \lim_{\xi \to \infty} \sigma_{\text{BSF}}^{\{n00\}} v_{\text{rel}} \approx \left( \frac{\pi \alpha^2}{\mu^2} \right) \alpha^2 S_{\text{BSF,C}}^{\{n00\}}(\zeta), \qquad (4.11a) $$

where the analytic expression for $S_{\text{BSF,C}}^{\{n00\}}$ can be found using eqs. (2.6), (4.8), (C.5a) and (C.5c). Since the general expression is rather lengthy, here we give the explicit form only for $n = 1, 2, 3$,

$$ S_{\text{BSF,C}}^{\{100\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \times \frac{2^6 \zeta^2 (3 + 2\zeta^2)}{15 (1 + \zeta^2)^2} e^{-4\zeta\arccot\zeta}, \qquad (4.11b) $$

$$ S_{\text{BSF,C}}^{\{200\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \times \frac{2^5 \zeta^2 (192 + 272\zeta^2 + 100\zeta^4 + 15\zeta^6)}{15 (4 + \zeta^2)^4} e^{-4\zeta\arccot(\zeta/2)}, \qquad (4.11c) $$

$$ S_{\text{BSF,C}}^{\{300\}}(\zeta) = \left( \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}} \right) \times \frac{2^6 \zeta^2 (3^3 + 7\zeta^2)(3^7 + 2^2 3^6 \zeta^2 + 915\zeta^4 + 122\zeta^6)}{3^3 5 (9 + \zeta^2)^5} e^{-4\zeta\arccot(\zeta/3)}. \qquad (4.11d) $$

At $\zeta \gg 1$, we find $S_{\text{BSF,C}}^{\{100\}}(\zeta) \simeq 0.16 \times 2\pi\zeta$.

We illustrate the features of $S_{\text{BSF}}^{\{100\}}(\zeta, \xi; \alpha)$ and $S_{\text{BSF,C}}^{\{100\}}(\zeta)$ in figures 9 and 10.

### 4.2.2 Capture into $\ell = 1$ bound states

Similarly to the previous section, we keep the terms of order $(|\mathbf{P}_\phi|/\kappa)^2$ from the $\mathcal{I}_{\mathbf{k},n1m}$ integrals and the zeroth-order terms from the $\mathcal{K}_{\mathbf{k},n1m}$ integrals [cf. eq. (B.5)]. Then, from eqs. (2.6) and (4.8), we find the total cross-section for capture to any $\ell = 1$ state, for a fixed $n$, to be<fnref n="8" />

$$ \sigma_{\text{BSF}}^{\{n1\}} v_{\text{rel}} \equiv \sum_{m=-1}^{1} \sigma_{\text{BSF}}^{\{n1m\}} v_{\text{rel}} \approx \left( \frac{\pi \alpha^2}{\mu^2} \right) \alpha^2 S_{\text{BSF}}^{\{n1\}}(\zeta, \xi; \alpha) pss_{n,1}^{1/2}, \qquad (4.12a) $$

where $S_{\text{BSF}}^{\{n1\}}$ includes terms suppressed by different powers of $pss_{n,1}$,

$$ S_{\text{BSF}}^{\{n1\}}(\zeta, \xi; \alpha) = \Sigma_0^{\{n1\}}(\zeta, \xi) + pss_{n,1} \times \Sigma_1^{\{n1\}}(\zeta, \xi) + pss_{n,1}^2 \times \Sigma_2^{\{n1\}}(\zeta, \xi), \qquad (4.12b) $$

<footnote n="8">⁸We note that the computation of $\sigma_{\text{BSF}}^{\{210\}}$ with emission of a (nearly massless) scalar mediator in ref. [1] is incorrect, due to an error in eq. (F.29).</footnote>
with

$$ \Sigma_0^{\{n1\}}(\zeta, \xi) = \frac{3}{4} \left( \frac{1 + \zeta^2 \gamma_{n,1}^2(\xi)}{\zeta^2} \right) \times \left| \int_0^\infty dx \left( \frac{2e^{-x/\xi}}{x} \right) \chi_{n,1}^*(x) \chi_{|k|,1}(x) \right|^2, \quad (4.12c) $$

$$ \Sigma_1^{\{n1\}}(\zeta, \xi) = -\frac{1}{2^4} \left( \frac{1 + \zeta^2 \gamma_{n,1}^2(\xi)}{\zeta^2} \right)^3 \times \left( \left[ \int_0^\infty dx \left( \frac{2e^{-x/\xi}}{x} \right) \chi_{n,1}^*(x) \chi_{|k|,1}(x) \right] \times \left[ \int_0^\infty dx \, x^2 \chi_{n,1}(x) \chi_{|k|,1}^*(x) \right] + \text{c.c.} \right), \quad (4.12d) $$

$$ \Sigma_2^{\{n1\}}(\zeta, \xi) = \frac{1}{2^4 \cdot 5^2} \left( \frac{1 + \zeta^2 \gamma_{n,1}^2(\xi)}{\zeta^2} \right)^5 \times \left( 11 \left| \int_0^\infty dx \, x^2 \chi_{n,1}^*(x) \chi_{|k|,1}(x) \right|^2 + 4 \left| \int_0^\infty dx \, x^2 \chi_{n,1}^*(x) \chi_{|k|,3}(x) \right|^2 \right). \quad (4.12e) $$

We showcase the resonant features of $S_{\text{BSF}}^{\{21\}}(\zeta, \xi; \alpha)$ in figure 11.

### 4.2.3 Annihilation

We now consider the annihilation of a bosonic particle-antiparticle pair into two scalar force mediators, $\chi^* + \chi \to \varphi + \varphi$. The dominant contribution to the annihilation cross-section arises from the $\ell = 0$ component of the scattering state wavefunction,

$$ \sigma_{\text{ann}}^{\text{sc}} v_{\text{rel}} = \frac{\pi \alpha^2}{\mu^2} \times S_{\text{ann}}^{(0)}(\zeta, \xi), \quad (4.13a) $$

where $S_{\text{ann}}^{(0)}$ and its Coulomb limit are given in eqs. (3.14c) and (3.14d); we repeat them here for convenience

$$ S_{\text{ann}}^{(0)}(\zeta, \xi) = |\phi_{\mathbf{k}}(\mathbf{r} = 0)|^2 = \lim_{x \to 0} \left[ \frac{\chi_{|\mathbf{k}|, \ell=0}(x)}{x} \right], \quad (4.13b) $$

$$ S_{\text{ann},C}^{(0)}(\zeta) = \lim_{\xi \to \infty} S_{\text{ann}}^{(0)}(\zeta, \xi) = \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}}. \quad (4.13c) $$

### 4.3 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

$$ \text{Bosons, Fermions with total spin 0:} \qquad \frac{1}{\sqrt{2}} [\phi_{\mathbf{k}}(\mathbf{r}) + \phi_{-\mathbf{k}}(\mathbf{r})], \qquad (4.14a) $$

$$ \text{Fermions with total spin 1:} \qquad \frac{1}{\sqrt{2}} [\phi_{\mathbf{k}}(\mathbf{r}) - \phi_{-\mathbf{k}}(\mathbf{r})]. \qquad (4.14b) $$

The wavefunction (4.14a) implies that the contribution of the even-$\ell$ modes participating in a process is doubled, while the contribution of the odd-$\ell$ 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

$$ \sigma_{\text{BSF}}^{\{n00\}} \text{ for IP} = 2 \times [\sigma_{\text{BSF}}^{\{n00\}} \text{ for DP}], \quad \text{cf. eq. (4.10a)}, \qquad (4.15a) $$

$$ \sigma_{\text{BSF}}^{\{n1\}} \text{ for IP} = 0, \qquad (4.15b) $$

$$ \sigma_{\text{ann}} \text{ for IP} = 2 \times [\sigma_{\text{ann}} \text{ for DP}], \quad \text{cf. eq. (4.13a)}. \qquad (4.15c) $$

Note that the vanishing result in eq. (4.15b) holds to working order in $\alpha$; contributions of higher order in $\alpha$, which we have not computed here, will yield a non-zero cross-section.

## 4.4 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 particle-antiparticle 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 [$\ell = 0$ and $\ell = 2$ modes of the scattering-state wavefunction, respectively, cf. eq. (4.10)], and by the dipole mode in the second case [$\ell = 1$ mode, cf. 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,

$$ \frac{\sigma_{\text{BSF,C}}^{\{100\}}}{\sigma_{\text{ann,C}}} = \alpha^2 \frac{2^6 \zeta^2 (3 + 2\zeta^2)}{15(1 + \zeta^2)^2} e^{-4\zeta\arccot\zeta}. \qquad (4.16) $$

At $\zeta \gtrsim 1$, this ratio becomes $\sim 0.16\alpha^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 $\ell$ mode of the scattering-state wavefunction to an inelastic process scales as $\sigma_{\text{inel}}v_{\text{rel}} \propto v_{\text{rel}}^{2\ell}$ at low velocities. For BSF, the $\ell = 2$ contribution will
Scalar mediator, non-self-conjugate bosonic particle-antiparticle pairs
Coulomb limit

<figure><img src="image_16.png" /><figcaption>**Figure 9.** Velocity dependence of the $S_{\text{BSF}}$ factors for the radiative bound-state formation by a bosonic particle-antiparticle pair coupled to a nearly massless scalar. $n$ and $\ell$ denote the principal and the angular-momentum quantum number of the bound states formed. For capture into $\ell = 1$ states, we have summed over all $m$ (the angular-momentum projections on one axis). Annihilation (not shown here) is the dominant inelastic process; in the Coulomb regime, $\sigma_{\text{BSF}}^{\{100\}}/\sigma_{\text{ann}}^{\text{sc}} = \alpha^2 S_{\text{BSF}}^{\{100\}}(\zeta)/S_{\text{ann}}^{(0)}(\zeta) \simeq 0.16\alpha^2$, at $\zeta = \alpha/v_{\text{rel}} \gg 1$ [cf. eqs. (4.11), (4.13), (4.16)].</figcaption></figure>

thus diminish, leaving ultimately only the $\ell = 0$ mode at sufficiently low $v_{\text{rel}}$. Eq. (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 $\ell$ modes, even for pairs of particles with different masses/couplings to the radiated scalar boson, the $\ell = 0$ mode may dominate the formation of zero-angular momentum bound states at low enough velocities. This is despite the contribution of the $\ell = 0$ mode being suppressed by a higher order in the coupling with respect to that of the $\ell = 1$ mode.
<figure><img src="image_17.png" /><figcaption>Scalar mediator, non-self-conjugate bosonic particle-antiparticle pairs
Capture into the ground state: resonances</figcaption></figure>

**Figure 10a.** The factors $\Sigma_0^{\{100\}}(\zeta, \xi)$, $\Sigma_1^{\{100\}}(\zeta, \xi)$ and $\Sigma_2^{\{100\}}(\zeta, \xi)$ that contribute to the cross-section for the radiative capture of a particle-antiparticle pair, into the ground state, with emission of a scalar force mediator [cf. eqs. (4.10)]. The bottom panel shows their sum, which corresponds to $S_{\text{BSF}}^{\{100\}}$ if the phase-space suppression is negligible. $\Sigma_0^{\{100\}}$ and $\Sigma_1^{\{100\}}$ involve only the $\ell=0$ mode of the scattering-state wavefunction, while $\Sigma_2^{\{100\}}$ includes contributions both from the $\ell=0$ and $\ell=2$ components. The $\ell$ modes determine the resonances.

<figure><img src="image_18.png" /><figcaption>Capture into the ground state: velocity dependence</figcaption></figure>

**Figure 10b.** Velocity dependence of the various contributions to the cross-section for the radiative capture of a particle-antiparticle pair into the ground state, with emission of a scalar force mediator: $\Sigma_0^{\{100\}}(\zeta, \xi)$ (dashed lines), $\Sigma_1^{\{100\}}(\zeta, \xi)$ (dot-dashed lines), $\Sigma_2^{\{100\}}(\zeta, \xi)$ (dot-dot-dashed lines), and their sum (solid lines) [cf. eq. (4.10)]. Because $\Sigma_0^{\{100\}}$, $\Sigma_1^{\{100\}}$ and $\Sigma_2^{\{100\}}$ involve the $\ell=0$ mode of the scattering-state wavefunction, they saturate to a constant value at low velocities. $\Sigma_2^{\{100\}}$ contains also a contribution from the $\ell=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, \ell = 1$ states: resonances

<figure><img src="image_19.png" /><figcaption>**Figure 11.** The factors $\Sigma_0^{\{21\}}(\zeta, \xi)$, $\Sigma_1^{\{21\}}(\zeta, \xi)$ and $\Sigma_2^{\{21\}}(\zeta, \xi)$ that contribute to the cross-section for the radiative capture of a particle-antiparticle pair into any $n = 2, \ell = 1$ state, with emission of a scalar force mediator [cf. eqs (4.12)]. The bottom right panel shows their sum, which corresponds to $S_{\text{BSF}}^{\{21\}}$, provided that the phase-space suppression is negligible. $\Sigma_0^{\{21\}}$ and $\Sigma_1^{\{21\}}$ involve only the $\ell = 1$ mode of the scattering state wavefunction, while $\Sigma_2^{\{21\}}$ includes contributions both from the $\ell = 1$ and $\ell = 3$ components. The latter is responsible for the sharper resonances at large $\xi$.</figcaption></figure>

# 5 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 [23–29].

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, 30–32]. The same chain of processes taking place in the dense environment of haloes today [1, 15, 18, 33], or in the interior of stars where DM may be captured [17], enhances the expected rate of the indirect detection signals and results in stronger constraints [18]. Moreover, it gives rise to correlated spectral features; 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.
We showcase some of the above in figure 12, for a minimal model of fermionic DM coupled to a light but massive dark photon that mixes kinetically with hypercharge [34, 35]. Models of this kind are frequently invoked in the literature [36–41], for example in the context of self-interacting DM [42–44], as well as a plausible explanation of various astrophysical anomalies [16, 45–48]. A thorough investigation of its phenomenology — including, for first time in the literature, a self-consistent treatment of bound-state effects both in the DM relic density determination and the indirect detection signals — has been recently carried out in ref. [18], where it was demonstrated that the formation and decay of bound states strengthen the constraints derived from γ-ray observations of the Milky Way and its Dwarf Spheroidal galaxies.

Beyond symmetric or self-conjugate DM, bound-state effects can be even more signifi-
cant in asymmetric DM models. Asymmetric DM with long-range self-interactions can form
stable bound states. Inside haloes today, the low-energy radiation emitted during capture
into a bound state, or in various level transitions between bound-state energy levels, can
give rise to signals observable by indirect searches [12–14, 49, 50]. The cosmological forma-
tion of stable bound states typically screens or curtails the DM self-interactions, and has to
be properly accounted for in any consistent phenomenological study [51]. This is particu-
larly important in the context of the self-interacting DM scenario [42, 44, 52, 53], as well as
in scenarios that feature a dissipative hidden sector [54–61]. Moreover, stable bound states
of asymmetric DM can give rise to distinct signatures in direct detection experiments [62–
64]. Finally, DM bound states may result in detectable collider signals [65–68].

Acknowledgments

We thank Franz Herzog, Bira van Kolck, and Andreas Nogga for useful discussions. This work was supported by the Netherlands Foundation for Fundamental Research of Matter (FOM) and the Netherlands Organisation for Scientific Research (NWO). J.d.V. and K.P. acknowledge support by NWO in the form of, respectively, the VENI and VIDI grants. K.P. was supported by the European Research Council (ERC) under the EU Seventh Framework Programme (FP7/2007-2013)/ERC Starting Grant (agreement n. 278234 — ‘NewDark’ project), and by the ANR ACHN 2015 Grant (‘TheIntricateDark’ project).
<figure><img src="image_20.png" /><figcaption>**Figure 12.** (Plots reproduced from ref. [18].) We consider dark matter consisting of Dirac Fermions $X, \bar{X}$, of mass $M_{DM}$, that couple to a dark photon of mass $m_{V_D}$. We assume that the DM relic density arose from thermal freeze-out in the hidden sector, and determine the dark fine-structure constant $\alpha_D$ as a function of $M_{DM}$ and $m_{V_D}$, taking into account both the direct DM annihilation into radiation and the formation and decay of unstable bound states, wherever applicable, according to ref. [7]. Left: below the blue dashed line, the Bohr momentum of the $X-\bar{X}$ pair is larger than the dark photon mass ($\xi \gtrsim 1$); bound states exist (cf. appendix A.2) and the Sommerfeld effect is significant. However, bound states can form with emission of a dark photon only in the parameter space below the blue solid line, where the mediator mass is less than the binding energy. Below the red dotted lines, the average momentum transfer between $X$ and $\bar{X}$ is larger than the mediator mass, $(M_{DM}/2)v_{\text{rel}} \gtrsim m_{V_D}$ [or $\xi \gtrsim \zeta$, cf. eq. (2.4)], for typical velocities during the DM chemical decoupling in the early universe, in the Milky Way, and in the Dwarf Spheroidal galaxies, as indicated in the plot. (We take $v_{\text{rel}} = \sqrt{2}v_0$.) In this regime, the interaction cross-sections are well approximated by their Coulomb limit. In the Coulomb limit, bound-state formation is faster than annihilation, $\sigma_{\text{BSF}}^{100} \simeq 3.13 \sigma_{\text{ann}}$ [cf. figure 2 and eqs. (3.8d), (3.14d)]. Right: the cross-sections times relative velocity, for annihilation (solid lines) and radiative capture to the ground state (dotted lines), averaged over the DM velocity distribution in the Milky Way ($v_0 = 220$ km/s, $v_{\text{esc}} = 533$ km/s), for two different values of the dark photon mass, $m_{V_D} = 10$ MeV (green smooth lines) and $m_{V_D} = 10$ GeV (orange lines with resonances). For $m_{V_D} = 10$ MeV, the Coulomb limit pertains to $M_{DM} \gtrsim 10$ GeV, as can be seen in the left panel. Bound-state formation is faster than annihilation wherever it is kinematically allowed, $M_{DM} \gtrsim 400$ GeV. For $M_{DM} \lesssim 1$ TeV, the Sommerfeld effect is not significant during freeze-out, which implies roughly $\alpha_D \propto M_{DM}$. However, in the Milky Way, where the average velocity is lower, the Sommerfeld effect enhances the cross-sections proportionally to $\alpha_D/v_{\text{rel}}$, for $M_{DM} \gtrsim 10$ GeV. Then, in the mass range $M_{DM} \sim 10$ GeV–1 TeV, the annihilation cross-section scales roughly as $\langle \sigma v_{\text{rel}} \rangle \propto M_{DM}$. For $M_{DM} \gtrsim 1$ TeV, the Sommerfeld enhancement is operative both during freeze-out and in the Milky Way, which causes the annihilation and BSF cross-sections to vary only slowly with $M_{DM}$. For $m_{V_D} = 10$ GeV, the Coulomb limit pertains only to $M_{DM} \gtrsim 30$ TeV. For lower values of $M_{DM}$, both the annihilation and the BSF cross-sections exhibit resonances, which are however located at different values of the DM mass. Bound-state formation may still dominate over annihilation in the mass range where it is kinematically possible. At the BSF threshold ($M_{DM} \simeq 5.5$ TeV), a disruption appears in the annihilation cross-section, and the resonances are shifted to lower values of the DM mass; this is due to the effect of BSF on the DM relic density, which reduces the estimated coupling $\alpha_D$.</figcaption></figure>
# A Wavefunctions

## A.1 The Schrödinger equation

The Schrödinger equations for the bound and scattering states are

$$ \left[-\frac{\nabla^2}{2\mu} + V(\mathbf{r})\right] \psi_{n\ell m}(\mathbf{r}) = \mathcal{E}_{n\ell} \psi_{n\ell m}(\mathbf{r}), \qquad (A.1a) $$

$$ \left[-\frac{\nabla^2}{2\mu} + V(\mathbf{r})\right] \phi_{\mathbf{k}}(\mathbf{r}) = \mathcal{E}_{\mathbf{k}} \phi_{\mathbf{k}}(\mathbf{r}), \qquad (A.1b) $$

with the wavefunctions normalised as follows

$$ \int d^3\mathbf{r} \, \psi_{n\ell m}^*(\mathbf{r}) \psi_{n\ell m}(\mathbf{r}) = 1, \qquad (A.2a) $$

$$ \int d^3\mathbf{r} \, \phi_{\mathbf{k}}^*(\mathbf{r}) \phi_{\mathbf{k}'}(\mathbf{r}) = (2\pi)^3 \delta^3(\mathbf{k} - \mathbf{k}'). \qquad (A.2b) $$

For a central potential $V(\mathbf{r}) = V(r)$ potential, we perform the standard separation of variables

$$ \psi_{n\ell m}(\mathbf{r}) = \kappa^{3/2} \left[ \frac{\chi_{n\ell}(\kappa r)}{\kappa r} \right] Y_{\ell m}(\Omega_{\mathbf{r}}), \qquad (A.3a) $$

$$ \phi_{\mathbf{k}}(\mathbf{r}) = \sum_{\ell=0}^{\infty} (2\ell + 1) \left[ \frac{\chi_{|\mathbf{k}|,\ell}(\kappa r)}{\kappa r} \right] P_{\ell}(\hat{\mathbf{k}} \cdot \hat{\mathbf{r}}), \qquad (A.3b) $$

with

$$ \int_0^\infty dx |\chi_{n\ell}(x)|^2 = 1, \qquad (A.4a) $$

$$ \int_0^\infty dx \chi_{n\ell}^*(x) \chi_{|k|,\ell}(x) = 0. \qquad (A.4b) $$

We set $x = \kappa r$ and

$$ \gamma \equiv \sqrt{-2\mu\mathcal{E}}/\kappa. \qquad (A.5) $$

Then, for the Yukawa potential of eq. (1.1), the radial Schrödinger equations read

$$ \chi''(x) + \left[ -\frac{\ell(\ell+1)}{x^2} - \gamma^2 + \frac{2e^{-x/\xi}}{x} \right] \chi(x) = 0, \qquad (A.6) $$

where we temporarily dropped the indices in the wavefunctions and energy eigenvalues for generality.

At $x \to 0$, and for $\ell > 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^{\ell+1}$ (regular) and $x^{-\ell}$ (irregular). Here, we are interested in regular solutions, for $\mathcal{E} = \mathcal{E}_n < 0$ and $\mathcal{E} = \mathcal{E}_{\mathbf{k}} > 0$, representing bound and scattering states respectively. This implies the boundary condition

$$ \lim_{x \to 0} \chi'(x) = (\ell + 1) \lim_{x \to 0} [\chi(x)/x]. \qquad (A.7) $$

The condition (A.7) will be valid also for $\ell = 0$.
<figure><img src="image_21.png" /><figcaption>Figure 13. $\gamma_{nl}(\xi) = \sqrt{-2\mu\mathcal{E}_{nl}/\kappa}$, for $n = 1$ (blue), $n = 2$ (purple), $n = 3$ (yellow), and $\ell = 0$ (solid), $\ell = 1$ (dashed) and $\ell = 2$ (dotted), with $0 \le \ell \le n-1$.</figcaption></figure>

## A.2 Bound states

For $\mathcal{E} < 0$, we seek solutions of eq. (A.6) that vanish at infinity,

$$ \lim_{x \to \infty} \chi(x) = 0. \qquad (A.8) $$

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 $\ell$,

$$ \mathcal{E} = \mathcal{E}_{nl}(\xi) \equiv -\gamma_{nl}^2(\xi) \times \frac{\kappa^2}{2\mu}. \qquad (A.9) $$

The lifting of the well-known $\ell$-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 $\gamma_{nl}(\xi)$ numerically, and present it in figure 13, for $n = 1, 2, 3$. We find that it can be well fit by the formula

$$ \gamma_{nl}(\xi) \approx \frac{1}{n} \left( 1 - \frac{n^2 \xi_c}{\xi} \right)^{\rho}. \qquad (A.10) $$

The best fit parameters $\xi_c$ and $\rho$ are given in table 4, for $1 \le n \le 3$.<fnref n="9" /> Equation (A.10) reproduces the Coulomb limit, $\lim_{\xi \to \infty} \gamma_{nl}(\xi) = 1/n$. Away from the Coulomb limit, the existence of bound states implies an $n$- and $\ell$-dependent lower bound on $\xi$.

In figure 14, we show the wavefunctions $\chi_{1,0}(x)$ and $\chi_{2,0}(x)$, for various values of $\xi$.

**Coulomb limit.** In the limit $\xi \to \infty$, eq. (A.6) admits analytic solutions with

$$ \gamma_{nl}^C = \lim_{\xi \to \infty} \gamma_{nl}(\xi) = 1/n, \qquad (A.11) $$

$$ \chi_{nl}^C(x) = \lim_{\xi \to \infty} \chi_{nl}(x) = \frac{1}{n} \left[ \frac{(n-\ell-1)!}{(n+\ell)!} \right]^{1/2} e^{-x/n} (2x/n)^{\ell+1} L_{n-\ell-1}^{2\ell+1}(2x/n) \qquad (A.12) $$

<footnote n="9">⁹Note that we are using a higher precision numerical fit for the computation of cross-sections of bound-state related processes.</footnote>
<table><caption>**Table 4.** The fit parameters $\xi_c$ and $\rho$, in the fitting formula $\gamma_{nl}(\xi) \approx (1/n)(1 - n^2\xi_c/\xi)^\rho$, for the energy eigenvalues $\mathcal{E}_{nl} = -\gamma_{nl}^2(\xi) \times \kappa^2/(2\mu)$. In the Coulomb limit, $\lim_{\xi \to \infty} \gamma_{nl}(\xi) = 1/n$.</caption><thead><tr><th>n</th><th>1</th><th colspan="2">2</th><th colspan="3">3</th></tr><tr><th>ℓ</th><th>0</th><th>0</th><th>1</th><th>0</th><th>1</th><th>2</th></tr></thead><tbody><tr><th>ξ<sub>c</sub></th><td>0.8399</td><td>0.8059</td><td>1.1195</td><td>0.79678</td><td>0.96883</td><td>1.1991</td></tr><tr><th>ρ</th><td>1.1129</td><td>1.1597</td><td>0.81847</td><td>1.1746</td><td>0.94518</td><td>0.7638</td></tr></tbody></table>

<figure><img src="image_22.png" /><figcaption>**Figure 14.** Left: the bound-state wavefunction $\chi_{n=1,\ell=0}(x)$, for various values of $\xi \equiv \mu\alpha/m_\phi$: analytical Coulomb limit $\xi \to \infty$ (solid), $\xi = 3$ (dashed), $\xi = 1.5$ (dotdashed), and $\xi = 1$ (dotted). Right: the bound-state wavefunction $\chi_{n=2,\ell=0}(x)$, in the Coulomb limit $\xi \to \infty$ (solid), and for $\xi = 6$ (dashed), $\xi = 5$ (dotdashed), and $\xi = 4$ (dotted).</figcaption></figure>

where $L_n^a$ are the generalised Laguerre polynomials of degree $n$. (We assume the normalization condition $\int_0^\infty x^a e^{-x} L_n^{(a)} L_m^{(a)} dx = [\Gamma(n+a+1)/n!] \delta_{n,m}$.)

### A.3 Two-particle scattering state

For the continuous spectrum with $\mathcal{E} = \mathcal{E}_k > 0$, we set

$$ \mathcal{E}_k = \frac{\mathbf{k}^2}{2\mu} = \frac{1}{2}\mu v_{\text{rel}}^2, \qquad (A.13) $$

or equivalently $\gamma^2 = -1/\zeta^2$. The wavefunctions are specified by the boundary condition (A.7), and the asymptotic behaviour at $x \to \infty$. At large $x$, the wavefunction $\chi_{|k|,\ell}$ behaves as (see e.g. ref. [15, chapter 7])

$$ \chi_{|k|,\ell}(x) \xrightarrow{x \to \infty} \frac{\zeta}{2i} \left[ e^{i(x/\zeta + \delta_\ell)} - e^{-i(x/\zeta - \ell\pi)} \right], \qquad (A.14a) $$

where the phase shifts $\delta_\ell$ depend on $\zeta$ and $\xi$. This implies that

$$ \left| \chi_{|k|,\ell}(x) \right|^2 + \left| \chi_{|k|,\ell}(x - \pi\zeta/2) \right|^2 = \zeta^2. \qquad (A.14b) $$

**Coulomb limit.** At $\xi \to \infty$, the analytical solutions of eq. (A.6) with $\mathcal{E} = \mathcal{E}_k > 0$, are

$$ \chi_{|k|,\ell}^{C}(x) = e^{\pi\zeta/2} \frac{\Gamma(1+\ell-i\zeta)}{(2\ell+1)!} x(2ix/\zeta)^{\ell} e^{-ix/\zeta} {}_1F_1(1+\ell+i\zeta; 2\ell+2; 2ix/\zeta). \qquad (A.15a) $$
The sum (A.3b) over the $\ell$ modes can be expressed in closed form,

$$ \phi_{\mathbf{k}}^{C}(\mathbf{r}) = \lim_{\xi \to \infty} \phi_{\mathbf{k}}(\mathbf{r}) = e^{\pi\zeta/2} \Gamma(1 - i\zeta) {}_1F_1[i\zeta; 1; i(kr - \mathbf{k} \cdot \mathbf{r})] e^{i\mathbf{k} \cdot \mathbf{r}}. \quad (A.15b) $$

# B Convolution integrals

## B.1 Definition

The cross-sections for radiative BSF depend on the following integrals [1]

$$ \mathcal{I}_{\mathbf{k},n\ell m}(\mathbf{b}) = \int \frac{d^3 p}{(2\pi)^3} \tilde{\psi}_{n\ell m}^*(\mathbf{p}) \tilde{\phi}_{\mathbf{k}}(\mathbf{p} + \mathbf{b}) = \int d^3 r \psi_{n\ell m}^*(\mathbf{r}) \phi_{\mathbf{k}}(\mathbf{r}) e^{-i\mathbf{b} \cdot \mathbf{r}}, \quad (B.1a) $$

$$ \mathcal{J}_{\mathbf{k},n\ell m}(\mathbf{b}) = \int \frac{d^3 p}{(2\pi)^3} \mathbf{p} \tilde{\psi}_{n\ell m}^*(\mathbf{p}) \tilde{\phi}_{\mathbf{k}}(\mathbf{p} + \mathbf{b}) = i \int d^3 r [\nabla \psi_{n\ell m}^*(\mathbf{r})] \phi_{\mathbf{k}}(\mathbf{r}) e^{-i\mathbf{b} \cdot \mathbf{r}}, \quad (B.1b) $$

$$ \mathcal{K}_{\mathbf{k},n\ell m}(\mathbf{b}) = \int \frac{d^3 p}{(2\pi)^3} \mathbf{p}^2 \tilde{\psi}_{n\ell m}^*(\mathbf{p}) \tilde{\phi}_{\mathbf{k}}(\mathbf{p} + \mathbf{b}) = - \int d^3 r [\nabla^2 \psi_{n\ell m}^*(\mathbf{r})] \phi_{\mathbf{k}}(\mathbf{r}) e^{-i\mathbf{b} \cdot \mathbf{r}}, \quad (B.1c) $$

where $\mathbf{b}$ is proportional to the momentum of the radiated particle, $\mathbf{b} \propto \mathbf{P}_{\varphi}$. Using the Schrödinger eq. (A.1a), and eq. (A.9), $\mathcal{K}_{\mathbf{k},n\ell m}(\mathbf{b})$ takes also the form

$$ \mathcal{K}_{\mathbf{k},n\ell m}(\mathbf{b}) = -\kappa^2 \gamma_{nl}^2(\xi) \mathcal{I}_{\mathbf{k},n\ell m}(\mathbf{b}) - 2\mu \int d^3 r V(\mathbf{r}) \psi_{n\ell m}^*(\mathbf{r}) \phi_{\mathbf{k}}(\mathbf{r}) e^{-i\mathbf{b} \cdot \mathbf{r}}. \quad (B.1d) $$

The first term in eq. (B.1d) yields a higher order correction to the contribution from $\mathcal{I}_{\mathbf{k},n\ell m}(\mathbf{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^{-i\mathbf{b} \cdot \mathbf{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

$$ e^{-i\mathbf{b} \cdot \mathbf{r}} = \sum_{\ell=0}^{\infty} (2\ell + 1) (-i)^{\ell} j_{\ell}(br) P_{\ell}(\hat{\mathbf{b}} \cdot \hat{\mathbf{r}}), \quad (B.2a) $$

where $j_{\ell}$ is the spherical Bessel function. Derivatives of $j_{\ell}$ that will arise in the following, can be re-expressed in terms of Bessel functions using the identity

$$ \frac{d j_{\ell}(z)}{d z} = \frac{\ell j_{\ell}(z)}{z} - j_{\ell+1}(z). \quad (B.2b) $$

The Legendre polynomials can be expanded in spherical harmonics as follows

$$ P_{\ell}(\hat{\mathbf{x}} \cdot \hat{\mathbf{y}}) = \frac{4\pi}{2\ell + 1} \sum_{m=-\ell}^{\ell} Y_{\ell m}^{*}(\Omega_{\mathbf{x}}) Y_{\ell m}(\Omega_{\mathbf{y}}), \quad (B.2c) $$

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,

$$ \int d\Omega Y_{\ell_1 m_1}(\Omega) Y_{\ell_2 m_2}(\Omega) Y_{\ell_3 m_3}(\Omega) = \sqrt{\frac{(2\ell_1+1)(2\ell_2+1)(2\ell_3+1)}{4\pi}} \begin{pmatrix} \ell_1 & \ell_2 & \ell_3 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_1 & \ell_2 & \ell_3 \\ m_1 & m_2 & m_3 \end{pmatrix}. \qquad (B.2d) $$

We always assume that the spherical harmonics are normalised according to

$$ \int d\Omega Y_{\ell m}^*(\Omega) Y_{\ell' m'}(\Omega) = \delta_{\ell\ell'} \delta_{mm'}. \qquad (B.2e) $$

### B.3 Expansion in the momentum of the emitted radiation

The exponential decay of the bound-state wavefunction $\psi_{n\ell m}(\mathbf{r})$ at large $r$ ensures that the integrands in eq. (B.1) are significant only for $\kappa r \gamma_{n\ell}(\xi) \lesssim$ few (cf. section A.2). In this range, the argument of the Bessel function in eq. (B.2a) is

$$ br = (\eta|\mathbf{P}_{\varphi}|) r \lesssim \text{few} \times \frac{\eta|\mathbf{P}_{\varphi}|}{\kappa\gamma_{n\ell}(\xi)} = \text{few} \times \frac{\eta\alpha}{2\gamma_{n\ell}(\xi)} \frac{1 + \zeta^2\gamma_{n\ell}^2(\xi)}{\zeta^2} pss_{n\ell}^{1/2}, $$

where here $\eta = \eta_1$ or $\eta_2$, and we used eq. (2.5a) for $|\mathbf{P}_{\varphi}|$. From this we deduce that $br \ll 1$, if

$$ v_{\text{rel}}^2 \ll \alpha \gamma_{n\ell}(\xi). \qquad (B.3) $$

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

$$ j_{\ell}(z) = \sum_{s=0}^{\infty} \frac{(-1)^s z^{\ell+2s}}{2^s s! (2s + 2\ell + 1)!!}, \qquad (B.4) $$

where for our purposes, $z = br$.

In computing the integrals (B.1), we express the wavefunctions as in eq. (A.3), and carry out the radial integration in the variable $x = \kappa r$. The expansion of the Bessel function over $z = (b/\kappa)x$ amounts thus to an expansion over $b/\kappa$. Since $b \propto |\mathbf{P}_{\varphi}|$, with $|\mathbf{P}_{\varphi}|$ given in eq. (2.5a), this is ultimately an expansion in $\alpha$.

### 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/\kappa$, using eq. (B.4). We find

$$ \begin{aligned} \mathcal{I}_{\mathbf{k},n\ell m}(\mathbf{b}) = & \left(\frac{4\pi}{\kappa}\right)^{3/2} \sum_{s=0}^{\infty} \sum_{\ell_R=0}^{\infty} \left(\frac{b}{\kappa}\right)^{\ell_R+2s} \frac{(-1)^{\ell_R+s} i^{\ell_R}}{2^s s! (2s + 2\ell_R + 1)!!} \\ & \times \sum_{m_R=-\ell_R}^{\ell_R} \sum_{\ell_I=0}^{\infty} \sum_{m_I=-\ell_I}^{\ell_I} (-1)^{m_I} Y_{\ell_I m_I}^*(\Omega_{\mathbf{k}}) Y_{\ell_R m_R}(\Omega_{\mathbf{b}}) \\ & \times \sqrt{(2\ell+1)(2\ell_R+1)(2\ell_I+1)} \begin{pmatrix} \ell & \ell_R & \ell_I \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell & \ell_R & \ell_I \\ -m & -m_R & m_I \end{pmatrix} \end{aligned} \qquad (B.5a) $$
$$ \times \int_0^\infty dx \, x^{\ell_R+2s} \chi_{n\ell}^*(x) \chi_{|\mathbf{k}|, \ell_I}(x), $$

$$ \mathcal{J}_{\mathbf{k},nlm}(\mathbf{b}) = \frac{(4\pi)^2 i}{\kappa^{1/2}} \sum_{s=0}^{\infty} \sum_{\ell_R=0}^{\infty} \left(\frac{b}{\kappa}\right)^{\ell_R+2s} \frac{(-1)^{\ell_R+s} i^{\ell_R}}{2^s s! (2s+2\ell_R+1)!!} \quad (\text{B.5b}) $$

$$ \times \sum_{m_R=-\ell_R}^{\ell_R} \sum_{\ell_I=0}^{\infty} \sum_{m_I=-\ell_I}^{\ell_I} Y_{\ell_I m_I}^*(\Omega_{\mathbf{k}}) Y_{\ell_R m_R}(\Omega_{\mathbf{b}}) $$

$$ \times \int d\Omega \, Y_{\ell_I m_I}(\Omega) Y_{\ell_R m_R}^*(\Omega) \int_0^\infty dx \, x^{1+\ell_R+2s} \nabla_{\mathbf{x}} \left[ \frac{\chi_{n\ell}^*(x)}{x} Y_{\ell m}^*(\Omega) \right] \chi_{|\mathbf{k}|, \ell_I}(x), $$

$$ \mathcal{K}_{\mathbf{k},nlm}(\mathbf{b}) = \sqrt{(4\pi)^3 \kappa} \sum_{s=0}^{\infty} \sum_{\ell_R=0}^{\infty} \left(\frac{b}{\kappa}\right)^{\ell_R+2s} \frac{(-1)^{\ell_R+s} i^{\ell_R}}{2^s s! (2s+2\ell_R+1)!!} \quad (\text{B.5c}) $$

$$ \times \sum_{m_R=-\ell_R}^{\ell_R} \sum_{\ell_I=0}^{\infty} \sum_{m_I=-\ell_I}^{\ell_I} (-1)^{m_I} Y_{\ell_I m_I}^*(\Omega_{\mathbf{k}}) Y_{\ell_R m_R}(\Omega_{\mathbf{b}}) $$

$$ \times \sqrt{(2\ell+1)(2\ell_R+1)(2\ell_I+1)} \begin{pmatrix} \ell & \ell_R & \ell_I \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell & \ell_R & \ell_I \\ -m & -m_R & m_I \end{pmatrix} $$

$$ \times \int_0^\infty dx \left[ -\gamma_{n\ell}^2(\xi) + \frac{2e^{-x/\xi}}{x} \right] x^{\ell_R+2s} \chi_{n\ell}^*(x) \chi_{|\mathbf{k}|, \ell_I}(x). $$

In the above, {$\ell_I, m_I$} and {$\ell_R, m_R$} are the orbital angular momentum quantum numbers of the incoming state and the radiated particle, respectively. In the $\nabla_{\mathbf{x}}$ operator of eq. (B.5b), the radial coordinate should be understood to be $x \equiv \kappa r$. We also note that $(-1)^m Y_{\ell m}^*(\Omega) = Y_{\ell, -m}(\Omega)$.

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 $\ell = 0$ bound states

We may evaluate the integrals needed for capture into zero angular momentum bound states, directly from eq. (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

$$ \int d\Omega_{\mathbf{r}} P_{\ell}(\hat{\mathbf{k}} \cdot \hat{\mathbf{r}}) e^{-i\mathbf{b} \cdot \mathbf{r}} \equiv \varpi, \quad (\text{B.6a}) $$

$$ \int d\Omega_{\mathbf{r}} \hat{\mathbf{r}} P_{\ell}(\hat{\mathbf{k}} \cdot \hat{\mathbf{r}}) e^{-i\mathbf{b} \cdot \mathbf{r}} = -\frac{\nabla_{\mathbf{b}}\varpi}{ir}. \quad (\text{B.6b}) $$

Using the identities (B.2), we find

$$ \varpi = \int d\Omega_{\mathbf{r}} P_{\ell}(\hat{\mathbf{k}} \cdot \hat{\mathbf{r}}) e^{-i\mathbf{b} \cdot \mathbf{r}} = 4\pi(-i)^{\ell} j_{\ell}(br) P_{\ell}(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}}), \quad (\text{B.7a}) $$
and

$$ \int d\Omega_r \hat{r} P_\ell(\hat{k} \cdot \hat{r}) e^{-i\mathbf{b} \cdot \mathbf{r}} = -4\pi(-i)^{1+\ell} \times \left\{ \hat{\mathbf{k}} \frac{j_\ell(br)}{br} \frac{dP_\ell(y)}{dy} + \hat{\mathbf{b}} \left[ \frac{j_\ell(br)}{br} \left( \ell P_\ell(y) - y \frac{dP_\ell(y)}{dy} \right) - j_{\ell+1}(br) P_\ell(y) \right] \right\}_{y=\hat{\mathbf{k}} \cdot \hat{\mathbf{b}}} . \quad (B.7b) $$

Combining the wavefunction decompositions (A.3) and the angular integrals (B.7), the integrals (B.1) for capture into {$n00$} bound states, become

$$ \mathcal{I}_{\mathbf{k},n00}(\mathbf{b}) = \sqrt{\frac{4\pi}{\kappa^3}} \sum_{\ell=0}^{\infty} P_\ell(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}}) (2\ell + 1)(-i)^\ell \int_0^\infty dx \, \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,\ell}(x) j_\ell(bx/\kappa), \quad (B.8a) $$

$$ \mathcal{J}_{\mathbf{k},n00}(\mathbf{b}) = -\sqrt{\frac{4\pi}{\kappa}} \sum_{\ell=0}^{\infty} (2\ell + 1)(-i)^\ell \int_0^\infty dx \left[ \frac{d\chi_{n,0}^*(x)}{dx} - \frac{\chi_{n,0}^*(x)}{x} \right] \chi_{|\mathbf{k}|,\ell}(x) \times \\ \times \left\{ \hat{\mathbf{k}} \frac{dP_\ell}{dy} \frac{j_\ell(bx/\kappa)}{bx/\kappa} + \hat{\mathbf{b}} \left[ \left( \ell P_\ell(y) - y \frac{dP_\ell(y)}{dy} \right) \frac{j_\ell(bx/\kappa)}{bx/\kappa} - P_\ell(y) j_{\ell+1}(bx/\kappa) \right] \right\}_{y=\hat{\mathbf{b}} \cdot \hat{\mathbf{k}}} \quad (B.8b) $$

and

$$ \mathcal{K}_{\mathbf{k},n00}(\mathbf{b}) = -\sqrt{4\pi\kappa} \sum_{\ell=0}^{\infty} P_\ell(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}}) (2\ell+1)(-i)^\ell \int_0^\infty dx \frac{d^2\chi_{n,0}^*(x)}{dx^2} \chi_{|\mathbf{k}|,\ell}(x) j_\ell(bx/\kappa) . \quad (B.8c) $$

## B.5.2 Leading-order contributions

Using eq. (B.4), we now expand the integrals (B.8) in powers of $b/\kappa$, and keep the leading-order terms. For $\mathcal{I}_{\mathbf{k},n00}$, the zero-th order contribution vanishes, due to the orthogonality of the $\psi_{n\ell m}$ and $\phi_{\mathbf{k}}$ (or particularly, the $\chi_{n\ell}$ and $\chi_{|\mathbf{k}|,\ell}$) wavefunctions. In our computations in sections 3 and 4, we shall need terms up to $\mathcal{O}[(b/\kappa)^2]$. For $\mathcal{J}_{\mathbf{k},n00}$ and $\mathcal{K}_{\mathbf{k},n00}$, the zero-th order terms in $b/\kappa$ suffice.

$$ \mathcal{I}_{\mathbf{k},n00}(\mathbf{b}) \approx -\sqrt{\frac{4\pi}{\kappa^3}} \left\{ \left(\frac{b}{\kappa}\right) i P_1(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}}) \int_0^\infty dx \, x \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,1}(x) \right. \\ \left. + \left(\frac{b}{\kappa}\right)^2 \left[ \frac{P_0(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}})}{6} \int_0^\infty dx \, x^2 \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,0}(x) \right. \right. \\ \left. \left. + \frac{P_2(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}})}{3} \int_0^\infty dx \, x^2 \chi_{n,0}^*(x) \chi_{|\mathbf{k}|,2}(x) \right] \right\}, \quad (B.9a) $$

$$ \mathcal{J}_{\mathbf{k},n00}(\mathbf{b}) \approx -\hat{\mathbf{k}} \sqrt{\frac{4\pi}{\kappa}} \int_0^\infty dx \left[ \frac{d\chi_{n,0}^*(x)}{dx} - \frac{\chi_{n,0}^*(x)}{x} \right] \chi_{|\mathbf{k}|,1}(x), \quad (B.9b) $$

$$ \begin{aligned} \mathcal{K}_{\mathbf{k},n00}(\mathbf{b}) &\approx -\sqrt{4\pi\kappa} \int_0^\infty dx \frac{d^2\chi_{n,0}^*(x)}{dx^2} \chi_{\mathbf{k},0}(x) \\ &= \sqrt{4\pi\kappa} \int_0^\infty dx \left(\frac{2e^{-x/\xi}}{x}\right) \chi_{n,0}^*(x) \chi_{\mathbf{k},0}(x). \end{aligned} \quad (B.9c) $$
## B.6 Capture into $\ell = 1$ bound states

It is straightforward to obtain the leading order terms of the $\mathcal{I}_{k,n1m}$ and $\mathcal{K}_{k,n1m}$ integrals, using eq. (B.5). Here, we give explicitly only the leading order contributions to the $\mathcal{J}_{k,n1m}$ integrals. To zero-th order in $b$,

$$ \mathcal{J}_{k,n10}(\mathbf{b}) \approx i \sqrt{\frac{12\pi}{\kappa}} \left\{ \left( \hat{\mathbf{k}} \cos \theta_{\mathbf{k}} - \frac{\hat{\mathbf{e}}_z}{3} \right) \int_0^\infty dx \left[ \chi_{n,1}'(x) - \frac{2\chi_{n,1}(x)}{x} \right]^* \chi_{|\mathbf{k}|,2}(x) \right. \\ \left. + \frac{\hat{\mathbf{e}}_z}{3} \int_0^\infty dx \left[ \chi_{n,1}'(x) + \frac{\chi_{n,1}(x)}{x} \right]^* \chi_{|\mathbf{k}|,0}(x) \right\}, \quad (B.10a) $$

$$ \mathcal{J}_{k,n11}(\mathbf{b}) \approx -i \sqrt{\frac{6\pi}{\kappa}} \left\{ \left( \hat{\mathbf{k}} \sin \theta_{\mathbf{k}} e^{i\phi_{\mathbf{k}}} - \frac{\hat{\mathbf{e}}_x + i\hat{\mathbf{e}}_y}{3} \right) \int_0^\infty dx \left[ \chi_{n,1}'(x) - \frac{2\chi_{n,1}(x)}{x} \right]^* \chi_{|\mathbf{k}|,2}(x) \right. \\ \left. + \frac{\hat{\mathbf{e}}_x + i\hat{\mathbf{e}}_y}{3} \int_0^\infty dx \left[ \chi_{n,1}'(x) + \frac{\chi_{n,1}(x)}{x} \right]^* \chi_{|\mathbf{k}|,0}(x) \right\}, \quad (B.10b) $$

and

$$ \mathcal{J}_{k,n1-1}(\mathbf{b}) = -\mathcal{J}_{k,n11}(\mathbf{b}), \quad (B.10c) $$

where $\hat{\mathbf{k}} = \hat{\mathbf{e}}_x \sin \theta_{\mathbf{k}} \cos \phi_{\mathbf{k}} + \hat{\mathbf{e}}_y \sin \theta_{\mathbf{k}} \sin \phi_{\mathbf{k}} + \hat{\mathbf{e}}_z \cos \theta_{\mathbf{k}}$.

# C Coulomb limit for capture into $\ell = 0$ bound states

## C.1 Wavefunctions

For our analytical computations in the Coulomb limit, we shall use the closed form of the scattering state wavefunction [cf. eq. (A.15b)],

$$ \phi_{\mathbf{k}}^{C}(\mathbf{r}) = \sqrt{S_0(\zeta)} {}_1F_1[i\zeta; 1; i(kr - \mathbf{k} \cdot \mathbf{r})] e^{i\mathbf{k} \cdot \mathbf{r}}, \quad (C.1a) $$

where<fnref n="10" />

$$ S_0(\zeta) \equiv \frac{2\pi\zeta}{1 - e^{-2\pi\zeta}}. \quad (C.1b) $$

We will consider only capture to zero angular momentum bound states. The $\ell = 0$ bound-state wavefunctions and their derivatives, are [cf. eq. (A.12)]

$$ \psi_{n00}^{C}(\mathbf{r}) = \sqrt{\frac{\kappa^3}{\pi n^5}} e^{-\kappa r/n} \sum_{s=0}^{n-1} \frac{n! (-2\kappa r/n)^s}{(n-s-1)!(s+1)! s!}, \quad (C.1c) $$

$$ \nabla \psi_{n00}^{C}(\mathbf{r}) = -\hat{\mathbf{r}} \sqrt{\frac{\kappa^5}{\pi n^7}} e^{-\kappa r/n} \sum_{s=0}^{n-1} \frac{n! (2n-s) (-2\kappa r/n)^s}{(n-s-1)!(s+2)! s!}, \quad (C.1d) $$

$$ \nabla^2 \psi_{n00}^{C}(\mathbf{r}) = \sqrt{\frac{\kappa^7}{\pi n^9}} e^{-\kappa r/n} \sum_{s=0}^{n-1} \frac{n! (2n-s) (-2\kappa r/n)^s}{(n-s-1)!(s+2)! s!} \left(1 - \frac{2+s}{\kappa r/n}\right), \quad (C.1e) $$

where we expanded the Laguerre polynomials for later convenience.

<footnote n="10">¹⁰The factor $\sqrt{S_0(\zeta)}$ in eq. (A.15b) often appears in the literature as $e^{\pi\zeta/2}\Gamma(1 - i\zeta)$.</footnote>
## C.2 Identities

We shall use the identity [69]

$$ \int d^3 r \frac{e^{-\bar{\kappa} r}}{4\pi r} {}_1F_1[i\zeta; 1; i(kr - \mathbf{k} \cdot \mathbf{r})] e^{i(\mathbf{k}-\mathbf{b}) \cdot \mathbf{r}} = \frac{[\mathbf{b}^2 + (\bar{\kappa} - ik)^2]^{-i\zeta}}{[(\mathbf{k}-\mathbf{b})^2 + \bar{\kappa}^2]^{1-i\zeta}} \equiv f_{\mathbf{k},\mathbf{b}(\bar{\kappa})}, \quad (C.2) $$

where ${}_1F_1$ is the confluent hypergeometric functions of the first kind. From the eq. (C.2), we find

$$ f_{\mathbf{k},\mathbf{b}=0}(\bar{\kappa}) = \frac{1}{k^2} \frac{e^{-2\zeta \arccot(\bar{\kappa}/k)}}{1 + (\bar{\kappa}/k)^2}, \quad (C.3a) $$

$$ [\nabla_{\mathbf{b}} f_{\mathbf{k},\mathbf{b}(\bar{\kappa})}]_{\mathbf{b}=0} = \hat{\mathbf{k}} \frac{2(1 - i\zeta)}{k^3} \frac{e^{-2\zeta \arccot(\bar{\kappa}/k)}}{[1 + (\bar{\kappa}/k)^2]^2}, \quad (C.3b) $$

and, keeping up to $b^2$ terms,

$$ \frac{d}{d\bar{\kappa}} f_{\mathbf{k},\mathbf{b}(\bar{\kappa})} \approx \frac{2}{k^3} \frac{e^{-2\zeta \arccot(\bar{\kappa}/k)}}{[1 + (\bar{\kappa}/k)^2]^2} \left[ \zeta - (\bar{\kappa}/k) + \frac{2b \cos(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}})}{k} \frac{(1 - i\zeta)(\zeta - 2\bar{\kappa}/k)}{1 + (\bar{\kappa}/k)^2} + \frac{2b^2 \cos^2(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}})}{k^2} \frac{(\zeta - 3\bar{\kappa}/k)(1 - i\zeta)(2 - i\zeta) + i(1 - i\bar{\kappa}/k)[(\zeta - \bar{\kappa}/k)^2 - i(\bar{\kappa}/k)(1 - i\zeta)]}{[1 + (\bar{\kappa}/k)^2]^2} \right]. \quad (C.3c) $$

## C.3 Convolution integrals

Starting from the definitions (B.1), and using the Coulomb wavefunctions (C.1) and the identity (C.2), we find

$$ \mathcal{I}_{\mathbf{k} \to n00}^{C}(\mathbf{b}) = -\sqrt{\frac{16\pi\kappa^3}{n^5}} S_0(\zeta) \left[ \sum_{s=0}^{n-1} \frac{n! (2\bar{\kappa})^s}{(n-s-1)! (s+1)! s!} \frac{d^{s+1}}{d\bar{\kappa}^{s+1}} f_{\mathbf{k},\mathbf{b}(\bar{\kappa})} \right]_{\bar{\kappa}=\kappa/n}, \quad (C.4a) $$

$$ \mathcal{J}_{\mathbf{k} \to n00}^{C}(\mathbf{b}) = -\sqrt{\frac{16\pi\kappa^5}{n^7}} S_0(\zeta) \nabla_{\mathbf{b}} \left[ \sum_{s=0}^{n-1} \frac{n! (2n-s)!}{(n-s-1)! (s+2)! s!} \frac{d^s}{d\bar{\kappa}^s} f_{\mathbf{k},\mathbf{b}(\bar{\kappa})} \right]_{\bar{\kappa}=\kappa/n}, \quad (C.4b) $$

and

$$ \mathcal{K}_{\mathbf{k} \to n00}^{C}(\mathbf{b}) = \sqrt{\frac{16\pi\kappa^7}{n^9}} S_0(\zeta) \left[ \sum_{s=0}^{n-1} \frac{n! (2n-s)!}{(n-s-1)! (s+2)! s!} \times \left( \frac{d^{s+1}}{d\bar{\kappa}^{s+1}} f_{\mathbf{k},\mathbf{b}(\bar{\kappa})} + \frac{n(2+s)}{\kappa} \frac{d^s}{d\bar{\kappa}^s} f_{\mathbf{k},\mathbf{b}(\bar{\kappa})} \right) \right]_{\bar{\kappa}=\kappa/n}. \quad (C.4c) $$
### C.4 Convolution integrals: expansion in the momentum of the emitted radiation

Using eq. (C.3), we expand the integrals (C.4), keeping terms up to order $(b/\kappa)^2$ for $\mathcal{I}_{k,n00}$, and only zero-th order terms for $\mathcal{J}_{k,n00}$ and $\mathcal{K}_{k,n00}$. We obtain the following$^{11}$

$$ \mathcal{I}_{k \to n00}^{C}(\mathbf{b}) \approx -\sqrt{\frac{n^3 \pi}{\kappa^3}} S_0(\zeta) \frac{b}{\kappa} \sum_{s=0}^{n-1} \frac{n! (2\zeta_n)^{4+s}}{(n-s-1)!(s+1)! s!} \times \\ \times \frac{d^s}{d\zeta_n^s} \left\{ \frac{e^{-2\zeta \arccot \zeta_n}}{(1+\zeta_n^2)^3} \left[ \cos(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}}) (1-i\zeta)(\zeta-2\zeta_n) \right. \right. \\ \left. \left. + \frac{b\zeta \cos^2(\hat{\mathbf{k}} \cdot \hat{\mathbf{b}})(\zeta-3\zeta_n)(1-i\zeta)(2-i\zeta)+i(1-i\zeta)[(\zeta-\zeta_n)^2-i\zeta_n(1-i\zeta)]}{1+\zeta_n^2} \right] \right\}_{\zeta_n=\zeta/n}, \quad (C.5a) $$

$$ \mathcal{J}_{k \to n00}^{C}(\mathbf{b}) \approx -\hat{\mathbf{k}} \sqrt{\frac{\pi}{\kappa n}} S_0(\zeta) (1-i\zeta) \sum_{s=0}^{n-1} \frac{n! (2n-s)}{(n-s-1)! (s+2)! s!} \frac{(2\zeta_n)^{3+s}}{(s+2)! s!} \\ \times \frac{d^s}{d\zeta_n^s} \left[ \frac{e^{-2\zeta \arccot \zeta_n}}{(1+\zeta_n^2)^2} \right]_{\zeta_n=\zeta/n}, \quad (C.5b) $$

$$ \mathcal{K}_{k \to n00}^{C}(\mathbf{b}) \approx \sqrt{\frac{\pi \kappa}{n^3}} S_0(\zeta) \sum_{s=0}^{n-1} \frac{n! (2n-s)}{(n-s-1)! (s+2)! s!} \times \\ \times \frac{d^s}{d\zeta_n^s} \left[ \frac{e^{-2\zeta \arccot \zeta_n}}{1+\zeta_n^2} \left( \frac{\zeta-\zeta_n}{1+\zeta_n^2} + \frac{n(2+s)}{2\zeta} \right) \right]_{\zeta_n=\zeta/n}. \quad (C.5c) $$

We use eq. (C.5) in our computations of the Coulomb limit of BSF cross-sections in eqs. (3.8), (4.6) and (4.11b).

### 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 [70] and deuteron formation via neutron capture on a proton target [71], 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 refs. [70, 72, 73]. 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.

$$ \sum_{s=0}^{n-1} \frac{n! 2^s \zeta_n^s}{(n-s-1)!(s+1)!s!} \frac{d^s}{d\zeta_n^s} \left[ \frac{e^{-2\zeta \arccot \zeta_n}}{(1+\zeta_n^2)^2} (\zeta - \zeta_n) \right]_{\zeta_n = \zeta/n} = 0. $$

<footnote>$^{11}$Note that the contribution from the b-independent term of eq. (C.3c) to the $\mathcal{I}_{k \to n00}^C$ integral vanishes (as expected) when the summation over s is performed, i.e.</footnote>
## D.1 Solution of the Lippmann-Schwinger and bound-state equation

The starting point is the non-relativistic Lippmann-Schwinger (LS) equation which in its general form is written as

$$ T_j^{l'l's's}(p', p, E) = V_j^{l'l's's}(p', p) + \sum_{l''s''} \int_0^\infty dp'' V_j^{l'l''s's''}(p', p'') \left( \frac{p''^2}{E - p''^2/(2\mu) + i\epsilon} \right) T_j^{l''l's's}(p'', p, E), $$

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_j^{l'l's's}$ 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_j^{l'l's's}(p', p)$ denotes a partial-wave-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) + \int_0^\infty dp'' V^l(p', p'') \left( \frac{p''^2}{E - p''^2/(2\mu) + i\epsilon} \right) T^l(p'', p, E) . \quad (D.1) $$

The potential for a vector mediator in momentum space is given by

$$ V(\mathbf{p}, \mathbf{p}') = - \frac{4\pi\alpha}{(\mathbf{p} - \mathbf{p}')^2 + m_{\varphi}^2}. \quad (D.2) $$

The partial-wave-decomposed potential that appears in equation (D.1) is defined as

$$ \begin{aligned} V^l(p', p) &= \frac{1}{(2\pi)^3} \langle p' l || V(\mathbf{p}, \mathbf{p}') || p l \rangle \\ &= -\frac{4\pi\alpha}{(2\pi)^2} \int_{-1}^{+1} dx P_l(x) \frac{1}{p^2 + p'^2 - 2pp'x + m_{\varphi}^2}, \end{aligned} \quad (D.3) $$

where $P_l(x)$ denotes the Legendre polynomials.

To numerically solve the LS equation, we need to deal with the $i\epsilon$ in the numerator of equation (D.1). We write

$$ \frac{1}{E - p''^2/(2\mu) + i\epsilon} = \frac{2\mu}{q_0^2 - p''^2 + i\epsilon} = \frac{2\mu}{q_0 + p''} \left( \frac{\mathcal{P}}{q_0 - p''} - i\pi\delta(q_0 - p'') \right), \quad (D.4) $$

where $\mathcal{P}$ denotes the principal value integral, and we introduced $E \equiv q_0^2/(2\mu)$. The LS equation can then be written a

$$ \begin{aligned} T^l(p', p, E) = V^l(p', p) &+ (2\mu)\mathcal{P} \int_0^{p_{\text{max}}} dp'' V^l(p', p'') \left( \frac{p''^2}{q_0^2 - p''^2} \right) T^l(p'', p, E) \\ &- i \frac{\pi(2\mu)q_0}{2} V^l(p', q_0) T^l(q_0, p, E), \end{aligned} \quad (D.5) $$

where we introduced $p_{\text{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\mu) \int_0^{p_{\max}} dp'' \left[ \begin{aligned} & V^l(p', p'') \left( \frac{p''^2}{q_0^2 - p''^2} \right) T^l(p'', p, E) \\ & - V^l(p', q_0) \left( \frac{q_0^2}{q_0^2 - p''^2} \right) T^l(q_0, p, E) \end{aligned} \right] \\ + (2\mu) q_0^2 V^l(p', q_0) T^l(q_0, p, E) \mathcal{P} \int_0^{p_{\max}} dp'' \frac{1}{q_0^2 - p''^2} \\ - i \frac{\pi(2\mu)q_0}{2} \sum_{l'' s''} V^l(p', q_0) T^l(q_0, p, E), \qquad (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

$$ \mathcal{P} \int_0^{p_{\max}} dp'' \frac{1}{q_0^2 - p''^2} = \frac{1}{2q_0} \ln \left( \frac{p_{\max} + q_0}{p_{\max} - q_0} \right). \qquad (D.7) $$

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 [74].

To obtain the momentum-space bound-state wave function we solve the homogeneous part of the LS equation

$$ \Psi_{nlm}(p) = \frac{1}{E_{nlm}^b - \frac{p^2}{2\mu}} \int dp' p'^2 V^l(p, p') \Psi_{nlm}(p'), \qquad (D.8) $$

where $E_{nlm}^b$ is the (negative) binding energy and the real wave function is normalized as

$$ \int dp p^2 \Psi_{nlm}(p)^2 = 1. \qquad (D.9) $$

Equation (D.8) can be immediately discretized and written as an eigenvalue equation which we again solve with a LAPACK routine [74]. 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 \rightarrow B_{nlm}(X_1 X_2) + \varphi$. We write this amplitude as

$$ A^{nlm, \lambda}(\mathbf{p}_{\text{in}}, \mathbf{P}_{\varphi}) = \langle \Psi_{nlm} | O^{\lambda}(\mathbf{P}_{\varphi}) | \mathbf{p}_{\text{in}} \rangle^{(+)}, \qquad (D.10) $$

where $\lambda$ denotes the polarization of the vector mediator and $\mathbf{p}_{\text{in}}$ ($p_{\text{in}} = |\mathbf{p}_{\text{in}}| \simeq \mu v_{rel}$) and $|\mathbf{P}_{\varphi}| = P_{\varphi}$ are, respectively, the incoming relative momentum of the DM pair in the c.o.m. frame and the outgoing mediator momentum. $O^{\lambda}$ 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 = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \int dp p^2 |plm\rangle\langle plm|$ and use

$$|\mathbf{P}_{\text{in}}\rangle = \sum_{l'=0}^{\infty} \sum_{m'=-l'}^{l'} Y_{l'm'}^{*}(\hat{p}_{\text{in}}) |p_{\text{in}} l'm'\rangle, \qquad (\text{D.11})$$

to write

$$A^{nlm,\lambda}(\mathbf{P}_{\text{in}}, \mathbf{P}_{\varphi}) = A_{\text{free}}^{nlm,\lambda}(\mathbf{P}_{\text{in}}, \mathbf{P}_{\varphi}) + A_{\text{scat}}^{nlm,\lambda}(\mathbf{P}_{\text{in}}, \mathbf{P}_{\varphi}), \qquad (\text{D.12})$$

where

$$A_{\text{free}}^{nlm,\lambda}(\mathbf{P}_{\text{in}}, \mathbf{P}_{\varphi}) = \sum_{l'=0}^{\infty} \sum_{m'=-l'}^{l'} \langle \Psi_{nlm} | O^{\lambda}(\mathbf{P}_{\varphi}) | p^{\text{in}} l'm' \rangle Y_{l'm'}^{*}(\hat{p}^{\text{in}}),$$

$$A_{\text{scat}}^{nlm,\lambda}(\mathbf{P}_{\text{in}}, \mathbf{P}_{\varphi}) = \sum_{l'=0}^{\infty} \sum_{m'=-l'}^{l'} \int dp p^2 \langle \Psi_{nlm} | O^{\lambda}(\mathbf{P}_{\varphi}) | p l'm' \rangle Y_{l'm'}^{*}(\hat{p}_{\text{in}}) \frac{2\mu}{p_{\text{in}}^2 - p^2 + i\epsilon} T^{l'}(p, p_{\text{in}}),$$

where $T$ denotes the $T$-matrix obtained above. The integral appearing in $A_{\text{scat}}^{nlm,\lambda}$ is numerically solved in the same way as was done for the $T$-matrix by adding and subtracting the divergence at $p = p_{\text{in}}$.

The next step requires the calculation of $O_{nlm}^{\lambda l'm'}(p, \mathbf{P}_{\varphi}) \equiv \langle \Psi_{nlm} | O^{\lambda}(\mathbf{P}_{\varphi}) | p l'm' \rangle$. For simplicity we investigate the current for two DM scalars with equal mass and opposite charge $c_1 = -c_2 = c$ such that $O^{\lambda} = cg/(4\mu)(P + P')^{\lambda}\tau_3$ where $P$ and $P'$ are, respectively, the in- and outgoing momentum of the scalar interacting with the outgoing mediator. We have introduced an 'isospin' operator, $\tau_3 = \text{diag}(1, -1)$, that indicates that the two DM scalars carry opposite charge. Of course, different charge configurations can be considered as well. At the same time the bound- and scattering state have been assigned<fnref n="12" /> a total isospin $(t$ and $t')$ and third component of total isospin $(m_t$ and $m'_t)$ which are useful bookkeeping devices. For scalars with opposite charge we have $m_t = m'_t=0$.

Using momentum conservation and the fact that the currents under consideration only couple to one of the DM particles (two-body current appear at higher order and can be included along the same lines) we obtain

$$O_{nlm}^{\lambda l'm'}(p, \mathbf{P}_{\varphi}) = \int d\Omega(\hat{p}) \Psi_{nlm}(k) Y_{lm}^{*}(\hat{k}) Y_{l'm'}(\hat{p}) * \frac{cg}{4\mu}(2k^{\lambda}) * 2 \langle tm_t | \tau_3 | t' m'_t \rangle \quad (\text{D.13})$$

where $\mathbf{k} = \mathbf{p} - \mathbf{P}_{\varphi}/2$ ($k = |\mathbf{k}|$ and $\hat{k} = \mathbf{k}/k$). As each scalar carries isospin 1/2 the combination of 2 scalars gives total isospin $t = 0, 1$ and $m_t = 0, \pm 1$. For the case at hand, we have $m_t = 0$ such that the total isospin can be both $t = 0$ or $t = 1$. The Pauli principle requires a symmetric wave function for two scalars such that $t = 0$ ($t = 1$) implies odd (even) orbital angular momentum. Since $\langle t m_t = 0 | \tau_3 | t' m'_t = 0 \rangle = (1 - \delta^{tt'})$, we see that the total isospin flip requires $|l - l'|$ to be odd.<fnref n="13" />

<footnote n="12">¹²We have not written these isospin factors in the potential in appendix D.1 for simplicity. We could have replaced $-4\pi\alpha \to (4\pi\alpha)\tau_3^{(1)}\tau_3^{(2)}$ where $\tau_3^{(i)}$ indicates the isospin of DM particle (i). As $\langle t m_t = 0 | \tau_3^{(1)} \tau_3^{(2)} | t' m'_t = 0 \rangle = -\delta^{tt'}$, we obtain the same scattering equations for $t = t' = 0$ and $t = t' = 1$ and we therefore dropped the isospin indices.</footnote>

<footnote n="13">¹³For identical scalars, however, we have $m_t = \pm 1$ so that only $t = 1$ is allowed. This implies that $|l - l'|$ must be even, as discussed in section 4.3.</footnote>
We can either solve equation (D.13) numerically or, as we do here, perform the angular integral analytically. Using standard angular momentum techniques we can write

$$O_{nlm}^{\lambda'm'}(p, \mathbf{P}_{\varphi}) = \frac{cg}{\mu}(1 - \delta^{tt'})\sqrt{\pi}(-1)^{m+m'}\sqrt{\hat{l}}\sum_{G=0}^{\infty}\sum_{m_G=-G}^{G} (l1G; 000)(l1G; -m \lambda m_G) \\ \times \sum_{\lambda_1=0}^{G} \sum_{\lambda_2=0}^{G} \delta^{\lambda_1+\lambda_2; G} \sqrt{\hat{\lambda}_1 \hat{\lambda}_2} \sqrt{\frac{\hat{G}!}{\hat{\lambda}_1! \hat{\lambda}_2!}} (p)^{\lambda_1} (-P_{\varphi}/2)^{\lambda_2} \sum_{f=0}^{\infty} (-1)^f \hat{f}^{3/2} g_{fG}^{nlm}(p, P_{\varphi}) \\ \times \sum_{g=0}^{\infty} \sum_{m_g=-g}^{m_g} \begin{pmatrix} f & f & 0 \\ \lambda_1 & \lambda_2 & G \\ l' & g & G \end{pmatrix} (f\lambda_1 l'; 000) (f\lambda_2 g; 000) \\ \times (l'gG; -m' m_g m_G) Y_{gm_g}(\hat{P}_{\varphi}), \qquad (D.14)$$

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, and $\hat{l} = 2l+1$. The function $g_{fG}^{nlm}(p, P_{\varphi})$ denotes a single numerical angular integral

$$g_{fG}^{nlm}(p, P_{\varphi}) = \int_{-1}^{1} dx' \frac{|k(x')|}{|k(x')|} \Psi_{nlm}(k(x')) P_f(x'), \qquad (D.15)$$

where $k(x') = p^2 + P_{\varphi}^2/4 - pP_{\varphi}x'$.

We are now in the position to evaluate the scattering amplitude in equation (D.10). Although the obtained expressions are valid in any coordinate frame, it is convenient to specify $\mathbf{P}_{\varphi} = P_{\varphi}\hat{z}$. In this frame, the differential BSF cross section is given by

$$\frac{d\sigma_{\text{BSF}}^{nlm}}{d\Omega} = \frac{\pi\mu P_{\varphi}}{p_{\text{in}}} \left( \sum_{\lambda=\pm 1} |A^{nlm,\lambda}(\mathbf{p}_{\text{in}}, \mathbf{P}_{\varphi} = P_{\varphi}\hat{z})|^2 + \frac{m_{\varphi}^2}{P_{\varphi}^2 + m_{\varphi}^2} |A^{nlm,\lambda=0}(\mathbf{p}_{\text{in}}, \mathbf{P}_{\varphi} = P_{\varphi}\hat{z})|^2 \right), \qquad (D.16)$$

such that only the transverse polarizations contribute in the Coulomb limit.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] K. Petraki, M. Postma and M. Wiechers, Dark-matter bound states from Feynman diagrams, JHEP 06 (2015) 128 [arXiv:1505.00109] [INSPIRE].

[2] H. An, M.B. Wise and Y. Zhang, Effects of bound states on dark matter annihilation, Phys. Rev. D 93 (2016) 115020 [arXiv:1604.01776] [INSPIRE].

[3] H. An, M.B. Wise and Y. Zhang, Strong CMB constraint on P-wave annihilating dark matter, arXiv:1606.02305 [INSPIRE].

[4] P. Asadi, M. Baumgart, P.J. Fitzpatrick, E. Krupczak and T.R. Slatyer, Capture and decay of electroweak WIMPonium, JCAP 02 (2017) 005 [arXiv:1610.07617] [INSPIRE].
[5] A. Sommerfeld, Über die Beugung und Bremsung der Elektronen (in German), *Ann. Phys.* **403** (1931) 257.

[6] P. Hoyer, *Bound states — from QED to QCD*, arXiv:1402.5005 [INSPIRE].

[7] B. von Harling and K. Petraki, *Bound-state formation for thermal relic dark matter and unitarity*, JCAP **12** (2014) 033 [arXiv:1407.7874] [INSPIRE].

[8] I. Baldes and K. Petraki, *Asymmetric thermal-relic dark matter: Sommerfeld-enhanced freeze-out, annihilation signals and unitarity bounds*, arXiv:1703.00478 [INSPIRE].

[9] K.M. Belotsky, E.A. Esipova and A.A. Kirillov, *On the classical description of the recombination of dark matter particles with a Coulomb-like interaction*, Phys. Lett. **B** **761** (2016) 81 [arXiv:1506.03094] [INSPIRE].

[10] K. Petraki and R.R. Volkas, *Review of asymmetric dark matter*, Int. J. Mod. Phys. **A** **28** (2013) 1330028 [arXiv:1305.4939] [INSPIRE].

[11] F.-Y. Cyr-Racine and K. Sigurdson, *Cosmology of atomic dark matter*, Phys. Rev. **D** **87** (2013) 103515 [arXiv:1209.5752] [INSPIRE].

[12] L. Pearce, K. Petraki and A. Kusenko, *Signals from dark atom formation in halos*, Phys. Rev. **D** **91** (2015) 083532 [arXiv:1502.01755] [INSPIRE].

[13] L. Pearce and A. Kusenko, *Indirect detection of self-interacting asymmetric dark matter*, Phys. Rev. **D** **87** (2013) 123531 [arXiv:1303.7294] [INSPIRE].

[14] J.M. Cline, Y. Farzan, Z. Liu, G.D. Moore and W. Xue, 3.5 keV X-rays as the "21 cm line" of dark atoms and a link to light sterile neutrinos, Phys. Rev. **D** **89** (2014) 121302 [arXiv:1404.3729] [INSPIRE].

[15] J.J. Sakurai, *Modern quantum mechanics*, (1985).

[16] M. Pospelov and A. Ritz, *Astrophysical signatures of secluded dark matter*, Phys. Lett. **B** **671** (2009) 391 [arXiv:0810.1502] [INSPIRE].

[17] C. Kouvaris, K. Langæble and N.G. Nielsen, *The spectrum of darkonium in the sun*, JCAP **10** (2016) 012 [arXiv:1607.00374] [INSPIRE].

[18] M. Cirelli, P. Panci, K. Petraki, F. Sala and M. Taoso, *Dark matter’s secret liaisons: phenomenology of a dark U(1) sector with bound states*, arXiv:1612.07295 [INSPIRE].

[19] K. Griest and M. Kamionkowski, *Unitarity limits on the mass and radius of dark matter particles*, Phys. Rev. Lett. **64** (1990) 615 [INSPIRE].

[20] X. Kong and F. Ravndal, *Proton proton scattering lengths from effective field theory*, Phys. Lett. **B** **450** (1999) 320 [Erratum ibid. **B** **458** (1999) 565] [nucl-th/9811076] [INSPIRE].

[21] X. Kong and F. Ravndal, *Proton proton fusion in leading order of effective field theory*, Nucl. Phys. **A** **656** (1999) 421 [nucl-th/9902064] [INSPIRE].

[22] K. Blum, R. Sato and T.R. Slatyer, *Self-consistent calculation of the Sommerfeld enhancement*, JCAP **06** (2016) 021 [arXiv:1603.01383] [INSPIRE].

[23] M. Cirelli, N. Fornengo and A. Strumia, *Minimal dark matter*, Nucl. Phys. **B** **753** (2006) 178 [hep-ph/0512090] [INSPIRE].

[24] M. Cirelli, A. Strumia and M. Tamburini, *Cosmology and astrophysics of minimal dark matter*, Nucl. Phys. **B** **787** (2007) 152 [arXiv:0706.4071] [INSPIRE].
[25] M. Cirelli and A. Strumia, *Minimal dark matter: model and results*, New J. Phys. 11 (2009) 105005 [arXiv:0903.3381] [INSPIRE].

[26] A. Hryczuk, R. Iengo and P. Ullio, *Relic densities including Sommerfeld enhancements in the MSSM*, JHEP 03 (2011) 069 [arXiv:1010.2172] [INSPIRE].

[27] A. Hryczuk, I. Cholis, R. Iengo, M. Tavakoli and P. Ullio, *Indirect detection analysis: wino dark matter case study*, JCAP 07 (2014) 031 [arXiv:1401.6212] [INSPIRE].

[28] M. Beneke, C. Hellmann and P. Ruiz-Femenia, *Heavy neutralino relic abundance with Sommerfeld enhancements — a study of pMSSM scenarios*, JHEP 03 (2015) 162 [arXiv:1411.6930] [INSPIRE].

[29] M. Cirelli, F. Sala and M. Taoso, *Wino-like minimal dark matter and future colliders*, JHEP 10 (2014) 033 [Erratum ibid. 01 (2015) 041] [arXiv:1407.7058] [INSPIRE].

[30] J. Ellis, F. Luo and K.A. Olive, *Gluino coannihilation revisited*, JHEP 09 (2015) 127 [arXiv:1503.07142] [INSPIRE].

[31] S. Kim and M. Laine, *Rapid thermal co-annihilation through bound states in QCD*, JHEP 07 (2016) 143 [arXiv:1602.08105] [INSPIRE].

[32] S. Kim and M. Laine, *On thermal corrections to near-threshold annihilation*, JCAP 01 (2017) 013 [arXiv:1609.00474] [INSPIRE].

[33] J.D. March-Russell and S.M. West, *WIMPonium and boost factors for indirect dark matter detection*, Phys. Lett. B 676 (2009) 133 [arXiv:0812.0559] [INSPIRE].

[34] B. Holdom, *Two U(1)'s and $\epsilon$ charge shifts*, Phys. Lett. B 166 (1986) 196 [INSPIRE].

[35] R. Foot and X.-G. He, *Comment on ZZ' mixing in extended gauge theories*, Phys. Lett. B 267 (1991) 509 [INSPIRE].

[36] B. Körs and P. Nath, *A Stückelberg extension of the Standard Model*, Phys. Lett. B 586 (2004) 366 [hep-ph/0402047] [INSPIRE].

[37] D. Feldman, B. Körs and P. Nath, *Extra-weakly interacting dark matter*, Phys. Rev. D 75 (2007) 023503 [hep-ph/0610133] [INSPIRE].

[38] M. Pospelov, A. Ritz and M.B. Voloshin, *Secluded WIMP dark matter*, Phys. Lett. B 662 (2008) 53 [arXiv:0711.4866] [INSPIRE].

[39] P. Fayet, *U-boson production in $e^+e^-$ annihilations, $\psi$ and $\Upsilon$ decays and light dark matter*, Phys. Rev. D 75 (2007) 115017 [hep-ph/0702176] [INSPIRE].

[40] M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, *Naturally light hidden photons in LARGE volume string compactifications*, JHEP 11 (2009) 027 [arXiv:0909.0515] [INSPIRE].

[41] P. Fayet, *The light U-boson as the mediator of a new force, coupled to a combination of Q, B, L and dark matter*, Eur. Phys. J. C 77 (2017) 53 [arXiv:1611.05357] [INSPIRE].

[42] D.N. Spergel and P.J. Steinhardt, *Observational evidence for selfinteracting cold dark matter*, Phys. Rev. Lett. 84 (2000) 3760 [astro-ph/9909386] [INSPIRE].

[43] J.L. Feng, M. Kaplinghat, H. Tu and H.-B. Yu, *Hidden charged dark matter*, JCAP 07 (2009) 004 [arXiv:0905.3039] [INSPIRE].

[44] A. Loeb and N. Weiner, *Cores in dwarf galaxies from dark matter with a Yukawa potential*, Phys. Rev. Lett. 106 (2011) 171302 [arXiv:1011.6374] [INSPIRE].
[45] N. Arkani-Hamed, D.P. Finkbeiner, T.R. Slatyer and N. Weiner, *A theory of dark matter*, Phys. Rev. D **79** (2009) 015014 [arXiv:0810.0713] [INSPIRE].

[46] I. Cholis, D.P. Finkbeiner, L. Goodenough and N. Weiner, *The PAMELA positron excess from annihilations into a light boson*, JCAP **12** (2009) 007 [arXiv:0810.5344] [INSPIRE].

[47] M. Abdullah, A. DiFranzo, A. Rajaraman, T.M.P. Tait, P. Tanedo and A.M. Wijangco, *Hidden on-shell mediators for the galactic center γ-ray excess*, Phys. Rev. D **90** (2014) 035004 [arXiv:1404.6528] [INSPIRE].

[48] A. Berlin, P. Gratia, D. Hooper and S.D. McDermott, *Hidden sector dark matter models for the galactic center gamma-ray excess*, Phys. Rev. D **90** (2014) 015032 [arXiv:1405.5204] [INSPIRE].

[49] K.K. Boddy, J.L. Feng, M. Kaplinghat, Y. Shadmi and T.M.P. Tait, *Strongly interacting dark matter: self-interactions and keV lines*, Phys. Rev. D **90** (2014) 095016 [arXiv:1408.6532] [INSPIRE].

[50] W. Detmold, M. McCullough and A. Pochinsky, *Dark nuclei I: cosmology and indirect detection*, Phys. Rev. D **90** (2014) 115013 [arXiv:1406.2276] [INSPIRE].

[51] K. Petraki, L. Pearce and A. Kusenko, *Self-interacting asymmetric dark matter coupled to a light massive dark photon*, JCAP **07** (2014) 039 [arXiv:1403.1077] [INSPIRE].

[52] A. Kusenko and P.J. Steinhardt, *Q ball candidates for selfinteracting dark matter*, Phys. Rev. Lett. **87** (2001) 141301 [astro-ph/0106008] [INSPIRE].

[53] J.L. Feng, H. Tu and H.-B. Yu, *Thermal relics in hidden sectors*, JCAP **10** (2008) 043 [arXiv:0808.2318] [INSPIRE].

[54] R. Foot and Z.K. Silagadze, *Thin disk of co-rotating dwarfs: a fingerprint of dissipative (mirror) dark matter?*, Phys. Dark Univ. **2** (2013) 163 [arXiv:1306.1305] [INSPIRE].

[55] J. Fan, A. Katz, L. Randall and M. Reece, *Dark-disk universe*, Phys. Rev. Lett. **110** (2013) 211302 [arXiv:1303.3271] [INSPIRE].

[56] R. Foot, *Tully-Fisher relation, galactic rotation curves and dissipative mirror dark matter*, JCAP **12** (2014) 047 [arXiv:1307.1755] [INSPIRE].

[57] R. Foot, *A dark matter scaling relation from mirror dark matter*, Phys. Dark Univ. **5-6** (2014) 236 [arXiv:1303.1727] [INSPIRE].

[58] R. Foot and S. Vagnozzi, *Dissipative hidden sector dark matter*, Phys. Rev. D **91** (2015) 023512 [arXiv:1409.7174] [INSPIRE].

[59] R. Foot and S. Vagnozzi, *Diurnal modulation signal from dissipative hidden sector dark matter*, Phys. Lett. B **748** (2015) 61 [arXiv:1412.0762] [INSPIRE].

[60] R. Foot, *Dissipative dark matter and the rotation curves of dwarf galaxies*, JCAP **07** (2016) 011 [arXiv:1506.01451] [INSPIRE].

[61] K.K. Boddy, M. Kaplinghat, A. Kwa and A.H.G. Peter, *Hidden sector hydrogen as dark matter: small-scale structure formation predictions and the importance of hyperfine interactions*, Phys. Rev. D **94** (2016) 123017 [arXiv:1609.03592] [INSPIRE].

[62] R. Laha and E. Braaten, *Direct detection of dark matter in universal bound states*, Phys. Rev. D **89** (2014) 103510 [arXiv:1311.6386] [INSPIRE].

[63] R. Laha, *Directional detection of dark matter in universal bound states*, Phys. Rev. D **92** (2015) 083509 [arXiv:1505.02772] [INSPIRE].
[64] A. Butcher, R. Kirk, J. Monroe and S.M. West, *Can tonne-scale direct detection experiments discover nuclear dark matter?*, arXiv:1610.01840 [INSPIRE].

[65] W. Shepherd, T.M.P. Tait and G. Zaharijas, *Bound states of weakly interacting dark matter*, Phys. Rev. D 79 (2009) 055022 [arXiv:0901.2125] [INSPIRE].

[66] H. An, B. Echenard, M. Pospelov and Y. Zhang, *Probing the dark sector with dark matter bound states*, Phys. Rev. Lett. 116 (2016) 151801 [arXiv:1510.05020] [INSPIRE].

[67] X.-J. Bi, Z. Kang, P. Ko, J. Li and T. Li, *ADMonium: asymmetric dark matter bound state*, Phys. Rev. D 95 (2017) 043540 [arXiv:1602.08816] [INSPIRE].

[68] F. Nozzoli, *A balance for dark matter bound states*, Astropart. Phys. 91 (2017) 22 [arXiv:1608.00405] [INSPIRE].

[69] A.I. Akhiezer and N.P. Merenkov, *The theory of lepton bound-state production*, J. Phys. B 29 (1996) 2135.

[70] J. de Vries, U.-G. Meißner, E. Epelbaum and N. Kaiser, *Parity violation in proton-proton scattering from chiral effective field theory*, Eur. Phys. J. A 49 (2013) 149 [arXiv:1309.4711] [INSPIRE].

[71] J. de Vries, N. Li, U.-G. Meißner, A. Nogga, E. Epelbaum and N. Kaiser, *Parity violation in neutron capture on the proton: determining the weak pion-nucleon coupling*, Phys. Lett. B 747 (2015) 299 [arXiv:1501.01832] [INSPIRE].

[72] W. Glöckle, *The quantum mechanical few-body problem*, Springer, Germany, (1983).

[73] E. Epelbaum, W. Glöckle and U.-G. Meißner, *The two-nucleon system at next-to-next-to-next-to-leading order*, Nucl. Phys. A 747 (2005) 362 [nucl-th/0405048] [INSPIRE].

[74] E. Anderson et al., LAPACK users' guide, third ed., Society for Industrial and Applied Mathematics, Philadelphia PA U.S.A., (1999).