Dark matter bound state formation via emission of a charged scalar

The formation of stable or meta-stable bound states can dramatically affect the phenomenology of dark matter (DM). Although the capture into bound states via emission of a vector is known to be significant, the capture via scalar emission suffers from cancellations that render it important only within narrow parameter space. While this is true for neutral scalar mediators, here we show that bound-state formation via emission of a charged scalar can be extremely significant. To this end, we consider DM charged under a dark U(1) force and coupled also to a light complex scalar that is charged under the same gauge symmetry. We compute the cross-sections for bound-state formation via emission of the charged scalar, and show that they can exceed those for capture via vector emission, as well as annihilation, by orders of magnitude. This holds even for very small values of the DM coupling to the charged scalar, and remains true in the limit of global symmetry. We then compute the DM thermal freeze-out, and find that the capture into meta-stable bound states via emission of a charged scalar can cause a late period of significant DM depletion. Our results include analytical expressions in the Coulomb limit, and are readily generalisable to non-Abelian interactions. We expect them to have implications for Higgs-portal scenarios of multi-TeV WIMP DM, as well as scenarios that feature dark Higgses or (darkly-)charged inert scalars, including models of self-interacting DM.


Introduction
Light scalar bosons that mediate a long-range force between dark matter (DM) particles appear in a variety of theories, including models of self-interacting DM and other hiddensector constructions. Intriguingly, it has been recently shown that the 125 GeV Higgs boson can also mediate a long-range interaction between TeV-scale particles, that affects their annihilation rate and can bind them into bound states [1,2]. This renders the dynamics of light scalar force carriers relevant as well to the phenomenology of DM consisting of JHEP02(2020)036 Weakly interacting massive particles (WIMPs). The existence of bound levels is a generic feature of theories with light force mediators that has rich phenomenological implications. The capture of unbound particles into bound states necessitates the dissipation of energy. This may occur radiatively, typically via emission of the mediator that is responsible for the long-range force. However, for a particle-antiparticle pair or a pair of identical particles, the radiative capture via emission of a scalar boson is rather suppressed due to cancellations in the amplitude that reflect in part the angular momentum selection rules of the process (cf. section 2.2) [3,4]. These cancellations concern the contributions to the radiative part of the amplitude that arise from the trilinear DM-DM-mediator coupling alone. The couplings of the scalar potential -the self-couplings of the mediator, as well as the biquadratic couplings between DM and the mediator if DM is bosonic -also contribute to the radiative amplitude and may enhance the capture cross-sections [5]. However for natural values of the parameters, the cross-sections remain mostly small.
In this work, we point out that the situation is markedly different if the emission of the scalar boson alters the potential between the interacting particles. This may occur if the scalar is charged under either a local or a global symmetry. As we shall see, in this case, the leading-order contributions to the amplitude are proportional to the overlap of the initial-state and final-state wavefunctions, which now are not orthogonal since they are subject to different potentials. The large overlap between the incoming and outgoing states gives rise to strikingly large bound-state formation (BSF) cross-sections. This is akin to atomic transitions precipitated by "sudden perturbations", such as ionisation caused by a beta decay of the nucleus [6].
To demonstrate the phenomenological importance of the transitions we consider, we calculate the chemical decoupling of DM in the early universe taking into account the formation of particle-antiparticle bound states via charged-scalar emission, and their subsequent decay into radiation. The formation of metastable bound states in the early universe has been shown to deplete the DM abundance [7], with the effect being generally more pronounced if the bound states have sizeable binding energy. Then, they form and decay efficiently already at high temperatures, when the DM density is large [2,7,8]. Here we find that, because of the largeness of the BSF cross-sections, shallow bound states can cause a second period of rapid DM depletion at low temperatures (of the order of their binding energy), much later than the traditional freeze-out. This alters the predicted couplings of DM to other species very significantly, thereby affecting its observable signatures.
For simplicity, in the present work, we carry out our computations in an Abelian model with scalar DM that is singly-charged under a dark U(1) D gauge force and is coupled to a doubly-charged light scalar via a trilinear coupling. We shall not assume that the light scalar obtains a vacuum expectation value (VEV). Even in models where it does, our computation remains essentially valid, provided that the VEV of the scalar is not much larger than its mass. Then, if the scalar mediator is light enough to be emitted during BSF, its mass and VEV must be smaller than all other relevant scales, and the symmetry is only mildly broken. Our results are also readily applicable to non-Abelian models, whose dynamics in the unbroken phase can be reduced to the Abelian case by an appropriate decomposition of the representations of the interacting particles [9].

JHEP02(2020)036
The paper is organised as follows. In section 2, we introduce the model, compute the cross-sections for BSF via emission of a charged scalar, provide analytical results in the Coulomb limit for capture into any bound level, and discuss their features. We confront our results with partial-wave unitarity and discuss the resolution to its apparent violation where it occurs. In section 3, we consider the DM freeze-out in the presence of BSF via emission of a charged scalar, and show the effect on the DM relic density and predicted couplings. We conclude in section 4 with an outlook on the implications of our results.
2 Bound-state formation via charged scalar emission

The model
We assume that DM consists of a complex scalar field X that couples to a dark Abelian gauge force U(1) D with V being the gauge boson, as well as to a light complex scalar Φ that is doubly charged under the same force. The interaction Lagrangian is with F µν ≡ ∂ µ V ν − ∂ ν V µ and D µ j ≡ ∂ µ + iq j gV µ , where the index j denotes the particle of charge q j . The charges are q X = 1 and q Φ = 2 for X and Φ fields respectively. The quartic terms stabilize the scalar potential at large field values. It is possible that the dark sector couples also to the Standard Model (SM), via biquadratic couplings of the scalars to the Higgs, and/or kinetic mixing of V with the Hypercharge gauge boson. Here we do not attempt a detailed phenomenological study of the model, but instead focus on computing the radiative capture into bound states via emission of a charged scalar, and simply showcasing its implications. We thus do not consider any couplings to the SM and do not derive any observational constraints.
We define the parameters that will appear in the non-relativistic potential, For convenience, we also define the total mass M and the reduced mass µ of a pair of interacting particles; in our case M = 2m X and µ = m X /2. We emphasise here that we are interested in m Φ m X . This hierarchy remains stable for momentum flows Q < m X , which encompass the momentum transfer along the off-shell Φ bosons exchanged in the scattering and bound states (cf. section 2.1.3). At Q m X , X loops generate corrections to the running mass of Φ, δ(m 2 Φ ) ∝ α Φ m 2 X , that may far exceed its low-energy value. Seen from a high-energy perspective, this amounts to a near cancellation between the high-energy value of the running Φ mass and the running contribution. While this may be fine-tuned, the value of the running Φ mass at high energies does not affect our computations. Moreover, eq. (2.1) can be viewed as an effective theory valid below ∼ m X , that is potentially stabilised by additional physics at higher scales, such as supersymmetry. Figure 1. The 2PI diagrams contributing to the non-relativistic potential for XX † pairs (upper) and XX or X † X † pairs (lower). The arrows denote the flow of the U(1) D charge.

Non-relativistic potential
The long-range potential of XX, X † X † and XX † pairs is generated by the one-bosonexchange diagrams shown in figure 1. Because the Φ-exchange diagram for XX † pairs is u-channel, the Φ-generated potential depends on the angular momentum mode of the eigenstate; we clarify this subtlety in appendix A. Combining this with well-known results for vector-mediated and scalar-mediated potentials [4,5,10], we obtain where α V and α Φ are defined in eq. (2.2). The potentials (2.3) distort the wavefunctions of pairs of unbound particles, a phenomenon known as the Sommerfeld effect [11,12]. For α V + (−1) α Φ > 0, they also give rise to XX † bound states. 1

Radiative capture processes
The capture of unbound particles into bound states can occur radiatively, via emission of a vector or scalar boson, according to the processes and We shall refer to the processes (2.4) and (2.5) as BSF V and BSF Φ respectively. The leading order Feynman diagrams are shown in figures 2 and 3. JHEP02(2020)036 While the bottom diagram appears to be naively of higher order, the momentum exchange along the two Φ propagators scales with the couplings, and renders this diagram of the same order as the two upper diagrams. We refer to appendix C for the computation. Figure 3. The capture into bound states via emission of a charged scalar (BSF Φ ), X + X → B(XX † ) + Φ. There is also the conjugate process, X † + X † → B(XX † ) + Φ † . The arrows denote the flow of the U(1) D charge. Note that in this case, the diagrams in which the incoming XX particles emit off-shell V and Φ that fuse to produce the final state Φ are of higher order, and we do not consider them here.
BSF V has been computed in [4,10] (see [8] for non-Abelian generalisations), and a number of papers have considered its effects on the DM relic density [2,7,8,[13][14][15][16] and indirect signals [13][14][15][17][18][19][20]. Here, the coupling of DM to the light scalar gives rise to an additional contribution to the BSF V amplitude at leading order, shown in figure 2. We review and adapt the computation of BSF V to the present model in appendix C.
In the rest of this section, we focus on the BSF Φ cross-sections. JHEP02(2020)036

Momentum decomposition and wavefunctions
We focus on the processes (2.5). For simplicity, in the following we neglect the mass of Φ, except in the phase-space integration. This will allow us to obtain analytical expressions for the BSF Φ cross-sections, and gain important insight. Taking fully into account the mass of Φ (and potentially also a non-zero mass for V ) requires computing the wavefunctions numerically, as done in [4] for a Yukawa potential and [1,2] for mixed Coulomb and Yukawa potentials. In the Coulomb approximation, the XX, X † X † scattering states and the XX † bound states are governed respectively by the potentials The momentum assignments for the particles participating in BSF Φ are shown in figure 4. In order to separate the motion of the center-of-momentum from the relative motion in the scattering and bound states, we decompose the momenta as follows [10] k ( ) Scattering states. In eq. (2.8a), the unprimed momenta correspond to infinite separation of XX, while the primed momenta denote the corresponding values in the XX wavepacket, which is distorted by the long-range interaction; this is the well-known Sommerfeld effect [11,12]. It is easy to see from eq. (2.8a) that in the non-relativistic regime, k = µv rel , with v rel being the relative velocity of the incoming XX pair. The nonrelativistic on-shell relations for k 0 1 and k 0 2 imply that the total energy of the scattering state is K 0 M + E k + K 2 /(2M ), with E k = k 2 /(2µ) = µv 2 rel /2. The scattering states are described by the wavefunctions φ XX k (r) in position space andφ XX k (k ) in momentum space; k is the expectation value of k . The wavefunctions φ XX k (r) obey the Schrödinger equation with the potential (2.3a) and energy eigenvalue E k .
Bound states. They are described by the wavefunctions ψ XX † n m (r) andψ XX † n m (p) in position and momentum space respectively, where n m are the standard principal and angular momentum quantum numbers in a central potential that determine the expectation value of p. The energy of the bound states is P 0 M − |E n | + P 2 /(2M ), where the binding energies can be parametrised as E n = −κ 2 B /(2n 2 µ), with κ B ≡ µα B being the Bohr momentum of the system. Note that p 0 1 and p 0 2 do not obey on-shell relations individually. The wavefunctions ψ XX † n m (r) obey the Schrödinger equation with the potential (2.3b) and energy eigenvalue E n .

JHEP02(2020)036
Hierarchy of scales. The emergence of non-perturbative phenomena -the Sommerfeld effect and the existence of bound states -is largely due to the different scales involved in the XX and XX † scattering. For the scattering states and the bound states, In our computations, we make approximations based on these hierarchies.
Energy-momentum conservation. Taking into account the above relations, the conservation of energy and momentum, K = P + P Φ , implies that Φ takes away the kinetic energy of the relative motion in the scattering state and the binding energy of the bound state [4,10], where we neglected the recoil of the bound state, as per (2.9). Equation (2.11) can be recast as |P Φ | = ω s 1/2 ps , with the phase-space suppression factor being Parametrisation. Throughout, we shall thus assume and use the well-known analytical solutions for the energy eigenstates and eigenvalues in a Coulomb potential (see e.g. ref. [21]), which we review in appendix B. We discuss the range of validity of the Coulomb approximation in section 3.4, in the context of the phenomenological application of BSF Φ on the DM freeze-out, that we present in section 3. In the Coulomb regime, we need the following two variables to parametrise the cross-sections in a minimal fashion [10], Taking the couplings (2.7) into account, ζ S and ζ B can be re-expressed in terms of (2.14) Outside the Coulomb regime, the wavefunctions can be computed as in ref. [4] (see [1,2] for results in a mixed Coulomb and Yukawa potential).

Amplitude
The amplitude for BSF Φ is [10] iM Φ k→n m where the factor 1/ √ 2µ has arisen in switching from the relativistic to the non-relativistic normalisation for the fields participating in the bound state [10]. A Φ T is the (amputated) amplitude of the radiative part of the process, (2.16) where the parentheses denote the momenta of each field; K and P are the total 4-momenta of the XX scattering state and the XX † bound state respectively. The leading order diagrams with the precise momentum assignments are shown in figure 4. Because these diagrams are not fully connected, the virtuality of the X, X † fields has to be integrated out as described in [10] (see [8] for a more recent summary). Adapting the result of [10], we find the leading order contributions to be Combining eqs. (2.15) and (2.17), Fourier transforming the wavefunctions, and taking into account that ψ n m (−r) = (−1) ψ n m (r), we obtain Amplitudes of this type can be computed by expanding in powers of P Φ · r/2 [4]. Indeed, the bound-state wavefunction is exponentially suppressed at r n/κ B . Moreover, the scattering state wavefunction oscillates at r > 1/k. Thus the integrand is significant roughly only for r 1/ κ 2 B /n 2 + k 2 . Taking eq. (2.11) into account, this implies P Φ ·r/2 κ 2 B /n 2 + k 2 /(4µ) = α 2 B /n 2 + v 2 rel /4 1. If the scattering and the bound states were JHEP02(2020)036 subject to the same potential, the zeroth order term in this expansion would vanish due to the orthogonality of the wavefunctions. 2 The essential point of our calculation is that because Φ carries away charge, the scattering and bound state wavefunctions are governed by different potentials and thus are not orthogonal. 3 Therefore, to lowest order, eq. (2.18) becomes Since the scattering state consists of a pair of identical bosons XX, the wavefunction is related to that of two distinguishable particles (DP), φ DP k (r), as follows where φ DP k, S (r) denotes the S angular mode of φ DP k (r). We now define the overlap integral of wavefunctions of distinguishable particles 4 with η1,2 ≡ m1,2/(m1 + m2). In the PΦ · r expansion, the zeroth order terms vanish due to the orthogonality of the wavefunctions, and for y1 = y2 and m1 = m2 also the first order terms cancel with each other. Then, the dominant contributions arise from the (PΦ · r) 2 terms (plus corrections of the same order that have been omitted in the above expression). Note that the second cancellation indicates the angular momentum selection rule ∆ = even. Thus, for a particle-antiparticle pair or a pair of identical particles, the capture cross-section is suppressed and becomes phenomenologically important mostly for large couplings. 3 While this is inevitable in an Abelian theory, in non-Abelian theories it is possible for a pair of particles to emit a charged boson without changing their combined representation. For example, because adj ⊗ adj ⊃ adj, two particles each transforming in the adjoint representation of a group, can begin from a combined adjoint configuration, emit an adjoint boson and end up again in an adjoint combined state. 4 In the notation of refs. [4,10], this is the overlap integral I k,n m (b), evaluated at b = 0 and up to the overall constant κ 3/2 B , here introduced to make R k,n m dimensionless.

JHEP02(2020)036
For = odd, we must keep the first order terms in the (P Φ · r) expansion of the integrand of eq. (2.18). Then the amplitude becomes non-vanishing for = odd, and is proportional to the overlap integral d 3 r (P Φ · r) φ XX k (r) [ψ XX † n m (r)] * , which imposes the selection rule | − S | = 1. As per eq. (2.20), the scattering state wavefunction contains only modes S = even. Thus, this contribution survives, and yields a cross-section that is larger than BSF via emission of a neutral scalar (cf. footnote 2) and of the equivalent order as BSF via emission of a vector boson. 5 Nevertheless, the cancellation of the zeroth order terms in (P Φ · r) for odd is a particularity of the model we are considering here, rather than a generic feature of BSF via charged scalar emission. For example, in a (coannihilation) scenario where the two incoming particles transform under different representations of the underlying symmetry (i.e. in the U(1) case, they have different charges), there in no generic cancellation between the contributing diagrams. Thus, to remain focused on the main point of our computation, we shall consider only the lowest order contributions given by eq. (2.19).
The result (2.23) should also make evident that for BSF Φ , the diagrams in which the final-state Φ is produced from the fusion of off-shell V and Φ emitted by the incoming XX pair, are subleading. Diagrams where the radiated boson is emitted from an off-shell propagator exchanged between the two interacting particles, are known to give leadingorder contributions to BSF V [22]. However, as shown here, the diagrams of figure 4 for BSF Φ yield lower order contributions than the corresponding diagrams for BSF V , where the vector-emission vertices introduce a momentum suppression in the wavefuntion overlap integral and thus in the amplitude (see appendix C for more details). Thus, with respect to the diagrams of figure 4, the Φ emission from V Φ fusion must be subleading; in fact, it turns out to be even of higher order than the corresponding diagrams in BSF V , due to the different Lorentz structure of the vertices involved.

Cross-section
The cross-section times relative velocity for the BSF Φ processes (2.5) is where the momentum of the emitted scalar P Φ is given by eq. (2.11). Collecting eqs. (2.11), (2.23), (2.24) and (B.12), we find for the capture cross-section, For fermionic DM, we find that the (PΦ · r) 0 contributions would survive for + s + 1 = even, with s =0 or 1 being the total spin. However, the XX wavefunctions contain only S + s + 1 = odd modes. Given the S = selection rule (2.22), these contributions cancel. The (PΦ · r) 1 contributions have the same fate: they would survive for + s + 1 = odd, however the selection rule now becomes | − S | = 1. Since S + s + 1 = odd, these contributions cancel as well.

JHEP02(2020)036
the following where S (ζ S ) is the Sommerfeld factor for -wave processes [23], 27) and in our model with α V and α Φ defined in eq. (2.2). 2 F 1 is the (ordinary) hypergeometric function, and s ps is the phase-space suppression factor defined in eq. (2.12). The cross-sections (2.26) is the main result of this section. They are readily generalisable to unbroken perturbative non-Abelian theories: as in ref. [8], the appropriate colour factors arising in the amplitudes (2.23) upon projection of the initial and final states onto states of definite colour must be included, α S and α B have to be chosen according to the initial and final colour representations [9], and the (anti)symmetrisation of the wavefunctions in the case of identical particles must be taken into account. Considering the couplings (2.28), we illustrate eq. (2.26) in figures 5 and 6, and compare BSF Φ with BSF V [4,10].
A few remarks are in order.
• The computed contribution (2.26) to the BSF cross-section vanishes if α B = α S , as expected from the orthogonality of the wavefunctions.
• The hypergeometric function in eq. (2.26) is a finite polynomial in its last argument (which can be also cast as 1+e i4arccot(ζ B /n) ) because its first argument is a non-positive integer, 1 + − n 0. For = n − 1, this factor reduces to 1. For an arbitrary , it tends to 1 both at large and at small velocities (ζ B , |ζ S | 1 and ζ B , |ζ S | 1). At intermediate velocities, it gives rise to cancellations, as seen for example in figure 5.
• At large velocities, v rel α B /n = (α V + α Φ )/n (i.e. ζ B /n 1), the overlap of the scattering and bound state wavefunctions is small, as seen from the term inside the square brackets in the second line of eq. (2.26). The BSF Φ cross-sections are suppressed by ζ 2( +3) B 1.

JHEP02(2020)036
Unitarity: s-wave Unitarity: p-wave Figure 5. The velocity-weighted cross-section for capture into zero-angular momentum XX † bound states, n00: XX → B n00 (XX † ) + Φ for n = 1, 2, 3, and XX † → B 100 (XX † ) + V . For V emission, the capture into n > 1 states is subdominant to capture into n = 1 [4], and we do not show them here. Also shown are the s-and p-wave unitarity limits on inelastic cross-sections; for capture into n00, BSF Φ is s-wave, while BSF V is p-wave. All cross-sections have been normalised to π/µ 2 , with µ being the reduced mass of the interacting particles, and we have neglected any phase-space suppression due to the mass of Φ.
• In between, the BSF Φ cross-sections become very significant. While the velocity at which they peak depends on n, and the ratio α S /α B , it can be roughly approximated by v rel ∼ α B /n = (α V + α Φ )/n, as seen in figures 5 and 6.
• For ζ B /n = α B /(nv rel ) 1, the factor in eq. (2.26) next to S (ζ S ) yields the characteristic behaviour v 2 rel of -wave processes without Sommerfeld. Combined with S (ζ S ), we see that the velocity suppression of higher partial waves disappears, and all partial waves exhibit the velocity dependence of S 0 .
Let us also now examine two different limits of eq. (2.26).
Limit α V → 0. In the limit of vanishing gauge coupling, the symmetry that is responsible for the dynamics we are considering becomes essentially global. This limit can be also effectively attained if the gauge boson has a very large mass m V µα V such that it mediates a contact type interaction. In this case, V decouples from the low-energy effective theory, leaving a remnant unbroken global symmetry. In this regime, Since there is now no repulsion in the scattering state, this cross-section is not exponentially suppressed at very low velocities. Nevertheless, because there is also no long-range attraction in the incoming state, the cross-section scales as (2.30) Despite the very small α Φ , this can exceed the BSF V cross-section for a significant velocity range, as seen in figures 5 and 6.
Capture via scattering. While in this work we focus on radiative BSF, it is possible that the dissipation of energy necessary for the capture into bound states occurs via scattering on other particles through exchange of an off-shell mediator, if the mediator couples also to other light degrees of freedom. Although of higher order, such processes can be extremely efficient inside a relativistic thermal bath, where the density of the light particles is very high, as was recently shown in [24] and previously suggested in [25][26][27].

JHEP02(2020)036
Reference [24] found that the rate of capture via scattering factorises into the radiative cross-section (albeit without any phase-space suppression due to the mass of the emitted scalar), times a part that includes the kinematics and dynamics of the bath particles. The cross-sections (2.26) can then be recast to calculate BSF via off-shell exchange of a charged scalar. A corollary of this and the largeness of the radiative cross-sections (2.26) is that the corresponding bath scattering processes must also be very significant in the early universe in the presence of light relativistic particles coupled to the charged scalar. In the case of an unbroken or mildly broken gauge symmetry, the gauge bosons and charged scalars already provide the relativistic bath necessary for such scattering processes to occur, In contrast to the radiative capture, BSF via bath scattering is not kinematically blocked if the mediator mass is larger than the energy available to be dissipated [cf. eq. (2.11)]. This implies that bound-state effects are not only enhanced, but also relevant to a broader parameter space.

Partial-wave unitarity
The unitarity of the S matrix implies an upper bound on the partial wave inelastic crosssection [28] σ where J denotes the partial wave of the scattering state wavefunction that participates in the process. For BSF Φ , this is the same as that of the bound state formed, thus we consider the ratio where we have neglected the phase-space suppression factor s 1/2 ps . As v rel varies, ζ S and ζ B scan a range of values, but the ratio r ≡ ζ S /ζ B of course remains constant (neglecting the possible running of the coupling, which may render α S mildly dependent on v rel ). Unitarity must be respected for all v rel . We thus define (2.34) Then eq. (2.32) implies that unitarity is respected provided that α Φ is sufficiently small, The coupling α Φ must be sufficiently small, α Φ < max(f n ;r ), such that unitarity is respected for all velocities. This is possible because f n ;r is bounded from above. Figure 8. The maximum value of α Φ vs. r ≡ α S /α B , for which the BSF Φ cross-sections (2.26) remain below the unitarity limit for all velocities. In the model considered in this work In figure 7, we present f n ;r (ζ B ) vs. ζ B for n = 1, = 0 and various values of r. In the present model, However, in any model where transitions of the type considered here occur, the BSF amplitudes will be proportional to the overlap integrals (2.21), and the cross-sections will be similar to eq. (2.26), up to a possible numerical factor. Thus, to get a broader insight into the implications of unitarity, in figure 7 we consider a wider range of r values. As seen, f n ;r is bounded from above; this remains true for all n, . It is then indeed possible to find a maximum value for α Φ , below which our calculation is consistent with unitarity, but above which it evidently fails. We determine this numerically and present it in figure 8. Notably, for α S /α B < 0, our computation fails JHEP02(2020)036 already at rather small values of α Φ . This is a consequence of the very large overlap between the initial and final states. The high peak of the BSF Φ cross-sections at v rel ∼ α B /n, explained in section 2.3, results in a rather stringent upper bound on α Φ .
What is the underlying reason for this apparent violation of unitarity, and how can unitarity be restored in the computation of the BSF Φ cross-sections? At such low values of α Φ it is unlikely that higher order corrections to the perturbative part of the amplitude A Φ T [cf. eq. (2.17)] may have any significant effect on the cross-section. Moreover, it has been pointed out that the breakdown of unitarity in perturbative calculations at low energies suggests that the two-particle interactions at infinity must be resummed [13, section 5]. In our computation, this has been done at leading order, by the resummation of the one-boson exchange diagrams of figure 1 that give rise to the potentials (2.3). However, by the optical theorem, all the inelastic processes to which the two interacting particles may participate also contribute to the self-energy of this two-particle state. Such contributions are typically neglected because they are of higher order than the one-boson exchange diagrams, and give rise to shorter-range (or contact) potentials that may have only limited impact on the large-distance behaviour of the wavefunctions. Still, the fact that the BSF Φ cross-sections can become so large suggests that their contribution to the two-particle self-energy may be significant, thus it must be resummed. The effect of this resummation will likely be significant mostly for incoming momenta around the peak of the BSF Φ cross-sections, k peak . While the corrected BSF Φ cross-sections should be consistent with unitarity, we expect them to remain very significant for k ∼ k peak , and essentially unaffected for k k peak and k k peak . Therefore, we still expect significant phenomenological implications. We leave this computation for future work.

Freeze-out of thermal-relic dark matter
To showcase the phenomenological applications of the above, we consider the effect of BSF Φ on the density of thermal-relic DM. Below, we list the pertinent cross-sections and rates, and present the Boltzmann equations that govern the evolution of the unbound and bound DM particle densities. We then describe how freeze-out is modified due to BSF Φ , and compute the couplings that reproduce the observed DM density. For simplicity, we assume that the DM particles and the radiation to which they couple are at the same temperature as the SM plasma, and use the standard time parameter The generalisation to different dark sector and SM temperatures is straightforward, see e.g. [14].

Annihilation
The tree-level annihilation channels for XX, X † X † and XX † pairs are shown in figure 9.
The annihilation processes are affected by the Sommerfeld effect as depicted in figure 10. Figure 9. The tree level diagrams contributing to the annihilation of XX † , XX and X † X † pairs. Note that for the XX † → V V and XX → ΦV annihilation processes, there are both t-and uchannel diagrams. In contrast, there is no u-channel diagram for XX † → ΦΦ † . The arrows denote the flow of the U(1) D charge. Figure 10. The long-range interaction affects the rate of the annihilation processes, and necessitates the resummation of the 2PI interactions at infinity. The 2PI kernels are shown in figure 1.

JHEP02(2020)036
The black blob stands for the tree-level annihilation diagrams of figure 9.
We consider only s-wave contributions, at leading order in each coupling and zeroth order in v rel . The full velocity-weighted cross-sections are where we recall that ζ V ≡ α V /v rel and ζ Φ = α Φ /v rel , and S 0 (ζ S ) ≡ 2πζ S /(1 − e −2πζ S ) is the s-wave Sommerfeld factor [cf. eqs. (2.14) and (2.27)]. Note that the cross-section (3.2a) for XX † → V V is twice as large as the spin-averaged cross-section for the annihilation of a fermionic particle-antiparticle pair into two Abelian vector bosons [4,7]. For the XX † → ΦΦ † annihilation, the s-channel diagram (annihilation via off-shell V ) is p-wave and we have neglected it in eq. (3.2b). For simplicity, in the following we shall also ignore the λ XΦ contribution. This coupling does not affect BSF Φ , and ignoring it will allow us to compare more easily the strength of the processes that arise from the essential couplings of the model, α Φ and α V . Thus, the total velocity-weighted annihilation cross-section we will consider is 3)

JHEP02(2020)036
with its thermal average being 3.1.2 Bound-state formation, ionisation and decay Formation. As already discussed, XX † bound states can form via emission of a V or a Φ boson, according to the processes (2.4) and (2.5), with the Feynman diagrams shown in figures 2 and 3. For simplicity, we shall consider the capture into the ground state only, n = 1, = m = 0, for both BSF V and BSF Φ . The larger binding energy and decay rate of the ground state render the ionisation processes unimportant earlier on, and imply that the capture into the ground state has a higher efficiency in depleting DM than the other BSF processes. Moreover, for BSF V , the capture into the ground state is the dominant contribution [4, figure 2]. On the other hand, for BSF Φ , the rate of capture into excited states may exceed that of capture into the ground state in some velocity range, as seen in figures 5 and 6. While we do expect that the capture into excited states plays an important role, here we only aim at showcasing the effect of BSF via emission of a charged scalar. We leave more detailed phenomenological studies for future work.
The effect of BSF V on the DM relic density was shown in [7], in a setup where V was the sole mediator; the corresponding cross-sections have been computed in [4,8,10]. In appendix C, we review the computation and adapt it to the present model.
Taking the above into account, the BSF cross-sections we will consider are where S 0 (ζ S ) ≡ 2πζ S /(1 − e −2πζ S ) is the s-wave Sommerfeld factor [cf. eq. (2.27)], and we neglect any phase-space suppression due to the mass of the Φ (cf. section 3.4). The total BSF cross-section is then The thermally averaged BSF cross-section is where the last factor accounts for the Bose enhancement due to the low-energy boson (V or Φ) emitted in the capture process; including the Bose enhancement is necessary in order to ensure detailed balance at large temperatures [7].

JHEP02(2020)036
Ionization. The ionisation rate of the bound states can be found either by using the Milne relation between the capture and ionization cross-sections (see [8, appendix D] for the proof), or more directly, by invoking detailed balance, where s is the entropy density of the universe, and the equilibrium yields of the unbound particles and the bound states, Y eq X ≡ n eq X /s and Y eq B ≡ n eq B /s, are given in section 3.2 below. Using these densities, we obtain where E B = E 10 is the binding energy of the ground state.
Decay into radiation. The dominant decays of the ground state are where (σ ann v rel ) pert 0 is the perturbative s-wave velocity-weighted annihilation cross-section (to zeroth order in v rel ), which is contained in eq. (3.3). Then,

Boltzmann equations and effective depletion rate
Let Y X = n X /s and Y B = n B /s be the yields of the unbound X particles and the bound states respectively. The Boltzmann equations that govern the evolution of the densities are [7] where σ ann v rel , σ BSF v rel and Γ dec JHEP02(2020)036 with g * and g * S being the energy and entropy relativistic degrees of freedom. We will take g * = g * S = g SM * + 4 to account for the SM plus the V and Φ degrees of freedom, during the DM freeze-out. We recall that the entropy density of the universe is s = (2π 2 /45)g * S T 3 . In the non-relativistic regime, the equilibrium yields Y eq X and Y eq The relic density of the X, X † particles is where s 0 2840 cm −3 and ρ c 3.67 × 10 −47 GeV 4 are the entropy and critical energy density of the universe today [29]. We require that Ω X = Ω DM , where the observed DM density is Ω DM 0.264 [29].
We note that in the minimal setup considered here, the Φ particles are stable, being the lightest degrees of freedom charged under U(1) D . As such, they contribute to the DM density. However, due to their lightness, even a small α V suffices to ensure that they annihilate into vector bosons efficiently, via the t/u-channel process ΦΦ † → V V , thereby reaching a cosmologically negligible relic abundance. Moreover, the radiation density due to V evades the CMB and BBN constraints provided that the dark sector temperature is somewhat lower than that of the SM plasma at the corresponding times. If Φ acquires a VEV and/or the dark sector couples to the SM, then there are more possibilities for the cosmological abundance of Φ and V , as well as constraints. A complete phenomenological study is beyond the scope of this work.
Effective depletion cross-section. Instead of the system of coupled eqs. (3.13), the effect of bound states on the relic density can be described by a single Boltzmann equation for the DM particles and an effective annihilation cross-section that includes BSF weighted by the fraction of bound states that decay rather than being ionised. We define first the effective BSF cross-section (3.17) The effective DM depletion cross-section is We may compute the X relic density by solving the Boltzmann equation [30] dY X dx = − π 45 m Pl m X g 1/2 * ,eff

JHEP02(2020)036
In figures 11 and 12, we show the DM annihilation and BSF cross-sections, and their thermal averages. The BSF Φ cross-sections can exceed both the BSF V and annihilation cross-sections by orders of magnitude, even for very small values of the couplings. However, the effect on the DM density depends on the interplay among the bound-state formation, ionisation and decay processes. To anticipate the result, it is useful to discern between two phases during the DM chemical decoupling.
(a) While the temperature is large enough to ensure that Γ ion B Γ dec B , the system is in a state of ionization equilibrium [31], where the effective DM depletion rate due to bound states is essentially independent of the BSF cross-section. Indeed, combining eqs. (3.9), (3.11) and (3.17) under the aforementioned condition, we obtain Note that eq. (3.20) does not rely on the specific couplings or interactions of the model considered here, but is rather general. It is easy to check that (3.20) is small in comparison to the annihilation cross-section, unless or until the temperature approaches or drops below the binding energy. In a U(1) model where DM couples only to the gauge boson, the ionization equilibrium ends at a temperature somewhat higher than the binding energy, thus (3.20) remains mostly small and most of the BSF effect on the DM density arises after that the end of ionisation equilibrium [7]. However, in the present model, the largeness of the BSF Φ cross-section sustains ionization equilibrium down to temperatures below the binding energy (cf. figure 12), thereby rendering the DM depletion significant during this phase. Clearly, while eq. (3.20) is independent of σ BSF v rel , the duration of ionisation equilibrium depends on it.
(b) The ionisation equilibrium ends when the ionisation rate drops below the decay rate, Γ ion B Γ dec B . Then, the DM depletion rate approaches rapidly the actual BSF rate, and is therefore sensitive to the exact BSF cross-section, σ BSF v rel eff σ BSF v rel .
Taking into account the above considerations, and in order to gain insight on whether BSF Φ may affect the DM density, in figures 13 and 14 we present the effective DM depletion cross-section for the following four cases: AnnP: perturbative annihilation only (diagrams shown in figure 9).  α V = 10 -3 α Φ = 10 -2 Figure 11. The BSF and annihilation velocity-weighted cross-sections that we consider in the computation of the DM freeze-out, vs. v rel . All cross-sections have been normalised to π/µ 2 , with µ = m X /2 being the reduced mass of the interacting particles. Also shown, the unitarity limits on s-wave and p-wave inelastic processes, which have to be respected by BSF Φ and BSF V , respectively, for capture into the ground state. α V = 10 -3 α Φ = 10 -2 Figure 12. The thermally averaged BSF and annihilation velocity-weighted cross-sections that we consider in the computation of the DM freeze-out, vs. x ≡ m X /T . All cross-sections have been normalised to π/µ 2 , with µ = m X /2 being the reduced mass of the interacting particles. We also mark two important mileposts: (i) the time when the temperature equals the binding energy T B = m X (α V + α Φ ) 2 /4, at and below which the equilibrium occupation number of the bound states becomes very significant, and (ii) the end of ionisation equilibrium, below which the DM depletion rate saturates to the BSF rate and thus becomes sensitive to the BSF cross-section.

AnnS
JHEP02(2020)036   Figure 14. The value of the thermally averaged velocity-weighted effective cross-section at its peak as a function of α V , for different values of α Φ . Note that the time at which the peak occurs, x = x peak , depends on α V and α Φ , and is chosen accordingly. We observe that, by varying α V , σv rel peak eff rises at α V α Φ , and becomes most significant at the limit of global symmetry, α V → 0.
In figure 13 we show the evolution of σv rel eff as the temperature drops. We choose a small value for the DM coupling to the charged scalar, α Φ = 10 −3 , to be well within the range that is consistent with unitarity (cf. section 2.4). We observe that the BSF V and BSF Φ contributions to the effective cross-section begin to rise at T ∼ |E B |, as implied by eq. (3.20). Later on, the ionisation equilibrium ends, and σ BSF v rel eff saturates to σ BSF v rel ; because this occurs at T < |E B |, when v 2 rel < α B = α V + α Φ < α V , the effective BSF cross-

JHEP02(2020)036
H 〈σv rel 〉 eff n X α V = 10 -5 〈σv rel 〉 eff n X α V = 10 -3 〈σv rel 〉 eff n X α V = 10 -2 〈σv rel 〉 eff n X α V = 10 -1. section σ BSF v rel eff decreases beyond this point, due to the repulsive potential in the XX and X † X † scattering states. The largeness of σ BSF v rel and consequently of σ BSF v rel eff around its peak, suggest that the DM depletion processes may recouple, even at a very low temperature, as we shall see in section 3.3.
Comparing the two plots of figure 13, corresponding to α V = 10 −2 and α V = 10 −4 , we also note that the rise and the peak of BSF Φ become more pronounced for smaller α V values, which allow for the effective BSF Φ cross-sections to grow larger before the suppression due to the repulsion in the scattering state settles in. To investigate this further, in figure 14 we show the dependence of σv rel eff evaluated at its peak, on α V . We see that σv rel peak eff rises at α V α Φ , and becomes most significant at α V → 0, i.e. in the limit where the symmetry becomes global.

Solutions of the Boltzmann equations
We now solve the Boltzmann eq. (3.19), discuss the qualitative features of the solutions and present numerical results. We focus mostly on small values of α Φ , roughly α Φ 10 −2 , in order to be consistent with unitarity (cf. section 2.4), and mark any parameter space where the BSF Φ cross-section violates it.

Freeze-out and recoupling of DM depletion processes
The Boltzmann eq. (3.19) describes the balance between the DM depletion processes and the expansion of the universe. Motivated by the sharp increase of σv rel eff at low temperatures seen in figure 13, we compare the X depletion rate Γ X ≡ n X σv rel eff with the expansion rate of the universe H = 4π 3 g * /45 T 2 /m Pl in the left plot of figure 15. Then in figure 16, we show the evolution of the DM density for different sets of parameters, in the four cases defined in section 3.2. The DM chemical decoupling is marked by two important events.  Figure 16. The evolution of the DM yield Y X ≡ n X /s vs. x ≡ m X /T for the four cases defined in section 3.2. We have fixed m X = 10 3 GeV.

JHEP02(2020)036
First freeze-out. As is standard, for x x FO1 ≈ 30, the DM depletion and creation processes are in equilibrium, i.e. Γ X H. Beyond this point, the X, X † densities depart from their equilibrium values, their exponential drop is stalled, and they begin to freezeout. However, because the annihilation and BSF cross-sections increase with decreasing temperature, the depletion of DM continues to be important until somewhat later, and may lead to the reduction of the DM density by a factor of a few. The Sommerfeld enhancement of the annihilation processes is important for v rel 10(α V + α Φ ), which upon thermal averaging implies x 10 −2 (α V + α Φ ) −2 . Thus if α V + α Φ 10 −2 , the Sommerfeld enhancement becomes significant already at x ∼ 10 2 , i.e. soon after freeze-out, while the DM density is still quite large. For BSF, a somewhat larger coupling, α V + α Φ 0.03, is required. In this range of couplings, the DM chemical decoupling is prolonged beyond freeze-out, as clearly seen in the bottom right plot of figure 16.

JHEP02(2020)036
Recoupling of DM depletion and second freeze-out. If α V + α Φ 10 −2 , the ionisation processes impede the DM depletion via BSF until quite late, when the DM density is rather low. However, the largeness of the BSF Φ cross-section may compensate for the smallness of the DM density, and result in the recoupling of the DM depletion processes around the time when σv rel eff peaks, at T |E B |. 7 We may estimate if and when this occurs as follows. The X, X † yield after the first freeze-out is estimated by the standard result, m X Y FO1 X σ ann v rel FO1 ∼ 8 × 10 −19 GeV −1 (see e.g. [34]), where we assumed that the direct annihilation dominates the DM depletion rate at that time; this is indeed true for the range of α V , α Φ where the recoupling may occur. The DM depletion recouples if Y FO1 X σv rel eff H/s, which implies σv rel eff / σ ann v rel stands for the thermally averaged Sommerfeld factor of the annihilation processes around the time of the first freeze-out. For large α V + α Φ , this condition yields a time close to the first freeze-out, and corresponds to the case when the DM chemical decoupling is simply delayed due to the BSF processes, as discussed above (cf. bottom right plot in figure 16). However, for smaller couplings, we obtain the following estimate for the time of recoupling and approximately the peak of σv rel eff , where we kept the leading order logarithmic correction in x peak . Note that x peak is independent of the DM mass m X . In the right plot of figure 15, we compare the semi-analytical prediction (3.21) with values of x peak determined numerically, and we find them in very good agreement. For α V + α Φ 10 −2 , (3.21) occurs much after the first chemical decoupling. In the top and the bottom left plots of figure 16, this manifests as a second plateau of the DM yield at large x. Clearly, the recoupling of the depletion processes at low temperatures results in very significant decrease of the DM abundance. This impels the re-determination of the couplings that give rise to the observed DM density.

Mass-coupling relation
We solve the Boltzmann eq. (3.19) numerically and determine the relation between α V , α Φ and m X that reproduces the observed DM density. In figures 17 to 19, we present our results.
From the previous discussion, we expect that the effect of BSF Φ is more pronounced at small α V , and in particular for α V α Φ . For this reason, in figure 17, we focus on the limit of global symmetry, α V = 0, and determine the relation between m X and α Φ , while in figures 18 and 19 we consider also the dependence on α V . 7 Recoupling of the DM depletion at late times may occur also due to the strong velocity dependence of Sommerfeld-enhanced cross-sections at (or close to) parametric resonance points [32,33]. The recoupling observed here is not due to resonant features, and applies to broader parameter space. In all cases, we see that taking BSF Φ into account changes the predicted couplings or mass very significantly, even by an order of magnitude. The effect on the DM density is illustrated in the right plot of figure 17 and the bottom plots of figure 18. We pick the combination of parameters that reproduce the observed DM abundance when considering AnnS+BSF V +BSF Φ , and then calculate the final density attained if only AnnP, AnnS, or AnnS+BSF V are taken into account. We observe the BSF Φ can deplete the DM density by more than two orders of magnitude.
Because of the interplay of the couplings α V and α Φ in σv rel eff , the DM relic density does not always vary monotonically with the parameters. In particular, for fixed m X and α Φ within some range, we observe in figures 18 and 19 that there are two values of α V that reproduce the observed DM density: a value in the range α V > α Φ where BSF Φ has little effect, and a value in the range α V < α Φ where BSF Φ has significant impact.

Validity of the Coulomb approximation
Throughout this paper, we have neglected the mass of the charged scalar Φ (as well as the possibility of a non-vanishing V mass). The calculation of section 2 can be generalised to include non-zero masses for Φ and V by evaluating numerically the wavefunctions and the overlap integral, as in refs. [1,2,4]. Here we examine the validity of the Coulomb approximation in the computation of the DM relic density.
In the present model, m Φ affects BSF Φ via the bound-state wavefunction and the phasespace suppression due to Φ emission. Note that m Φ does not affect the XX and X † X † scattering states that participate in BSF Φ . The conditions for the Coulomb approximation to be valid are as follows.  Figure 18. Top: the relation between the DM mass m X and the coupling to the gauge boson α V that reproduces the observed DM density, when considering different contributions to the DM depletion. In the gray-shaded region, our computation of the BSF Φ cross-section violates unitarity within a range of velocities. Bottom: the X, X † relic density when taking into account only some of the DM depletion processes, normalised to the observed DM density. For each value of α V we have chosen m X such that the observed DM density is reproduced by the AnnS+BSF V +BSF Φ calculation (red line in plots above).
(i) For a pure Yukawa potential, the ground state is Coulombic if the mediator mass is much smaller than the Bohr momentum. In the absence of V , this would imply m Φ µα Φ [4]. The presence of the V -mediated attractive Coulomb potential relaxes this condition [2]. Indicatively we note that, neglecting the V -generated potential, the binding energy is larger than 90% of its Coulomb value if m Φ < µα Φ /10 [4, figure 13]. On the other hand, for α V = α Φ , this occurs if m Φ < µα Φ /2 [2, figure 6]. For simplicity, JHEP02(2020)036  Figure 19. The combination of the couplings α V , α Φ that reproduces the observed DM abundance, for fixed values of the DM mass, when considering different contributions to the DM depletion.
we shall thus assume the following condition for the Coulomb approximation In the thermal bath, µv 2 rel /2 = 3T /2. During the first freeze-out the temperature is large, T |E B |, and E k dominates the energy available to be dissipated. However, the recoupling of the DM depletion processes occurs at T ∼ |E B | (if at all), when E k ∼ |E B |. Therefore, we require (3.23)

JHEP02(2020)036
The condition (3.23) is stronger than (3.22). Particularly for the small values of α Φ and α V we have considered here, it ensures that the bound states are very nearly Coulombic.
For BSF V , there is no kinematic blocking, provided that V is massless. However, a non-vanishing m Φ affects the XX † scattering state, as well as the XX † bound state. For the latter, the condition for the Coulomb approximation is (3.22). We briefly discuss the scattering state. In the case of a pure attractive Yukawa potential, the Coulomb limit is obtained if the mediator mass is lower than the average momentum transfer, m Φ µv rel [4]. For a pure repulsive Yukawa potential the condition is somewhat stronger. On the other hand, this condition is relaxed by the superposition of the V -mediated attractive Coulomb force [1, figure 2]. Since BSF V can be important only at early times, when the average kinetic energy is still fairly large, and provided that α V + α Φ is sufficiently large, the Coulomb approximation is typically justified (see e.g. discussion in ref. [20]). Regardless of the validity of the approximation, BSF V is not the focus of this paper, thus we do not elaborate on this issue further.
Finally, we note that if the charged scalar obtains a VEV, v Φ , the symmetry-breaking phase transition is expected to occur at temperature T PT v Φ . Then, if v Φ < |E B |, the DM chemical decoupling -including both the first freeze-out and the recoupling epoch -takes place essentially in the unbroken phase, and the computation of this section is applicable. Assuming that (3.23) However, the DM chemical decoupling may occur in the unbroken phase even for v Φ m Φ , if both v Φ and m Φ are much lower than |E B |.

Conclusion
The existence of bound states is a generic feature of theories with light force mediators. The formation of stable or metastable bound states has severe implications for the phenomenology of DM today. Scalar force mediators have been invoked in a variety of theories, including self-interacting DM and Higgs portal models residing in the multi-TeV regime.
Here, we computed the cross-sections for the radiative capture of non-relativistic particles into bound states via emission of a scalar that is charged under either a local or global symmetry. The emission of a charged scalar alters the Hamiltonian between the interacting particles, and precipitates extremely rapid transitions. We have provided analytical formulae in the Coulomb approximation for the capture into any bound level [cf. eq. (2.26)]. While we carried out our calculations in the context of a minimal U(1) model, our results are readily generalisable to more complex models, including perturbative non-Abelian theories, and can thus be relevant to the phenomenology of various scenarios, e.g. [1,2,[35][36][37]. Importantly, our results can be recast to compute BSF via scattering on a bath of relativistic particles, through exchange of a charged scalar, according to ref. [24]. This can be particularly important for the chemical decoupling of multi-TeV WIMP DM coupled to the 125 GeV Higgs [37]; in this regime the Higgs can indeed act as a light mediator [1,2], even if its on-shell emission in capture processes is not kinematically allowed.

JHEP02(2020)036
The phenomenological implications of the processes we computed can be striking. Here, we demonstrated that the formation of particle-antiparticle bound states via emission of a charged scalar and their subsequent decay can deplete the DM density by as much as two orders of magnitude. While for simplicity we considered only capture into the ground state, the computed cross-sections strongly suggest that the capture into excited states during the DM chemical decoupling should also be significant, thereby producing an even more important effect. The depletion of DM via these processes in the early universe alters the predicted relation between the DM mass and couplings rather dramatically. This in turn implies very different predictions for the DM signals in collider, direct and indirect detection experiments.
For indirect detection, the modification in the predicted relation between the DM mass and couplings implies that the signals arising from the direct annihilation processes and BSF via vector emission are very suppressed with respect to what expected when neglecting BSF via charged scalar emission during freeze-out. This essentially invalidates any existing constraints. On the other hand, BSF via charged scalar emission occurring during CMB or inside halos today may itself produce very significant radiative signals that result in strong constraints. For direct detection experiments, the implications are again varied. The larger predicted DM mass can bring a model previously thought to reside in the sub-GeV regime, within the threshold of current detectors. On the other hand, it can relax existing constraints for models already within the experimental sensitivity. Finally, the large BSF Φ cross-sections may imply late kinetic decoupling of DM from radiation in the early universe, as well as strong DM self-interactions inside halos today; both features can potentially affect the galactic structure very significantly.
However, as discussed earlier, the magnitude of the computed cross-sections eventuates in the apparent violation of unitarity already at rather low values of the relevant couplings. Reliable phenomenological studies therefore necessitate first the appropriate treatment of this issue.
Note added. While we were finalising this manuscript, ref. [38] appeared on the arXiv, which considers a similar setup: fermionic DM coupled to a dark U(1) gauge force and to a doubly charged complex scalar, which is assumed to obtain a VEV. Reference [38] points out the non-cancellation of the leading order term in the overlap integral governing the transitions between eigenstates of different potentials that occur via emission of the Goldstone mode, which, in that setup, has been absorbed by the gauge boson. In the present work, we do not assume a VEV for the scalar mediator, and our computation demonstrates the effect also in the limit of the underlying symmetry being global. Moreover, we provide analytical formulas in the Coulomb limit for the cross-sections of BSF via charged scalar emission, which are readily generalisable to non-Abelian theories. We note that the fermionic and scalar models exhibit different symmetry properties that determine whether the contributions under consideration survive (cf. section 2.2). Both ref. [38] and the present work find that transitions proportional to the overlap of the incoming and outgoing wavefunctions without any momentum suppression, result in the recoupling of the DM depletion processes at late times.

JHEP02(2020)036
X X Figure 20. Momentum decomposition of the 2PI diagrams of the XX † interaction.

JHEP02(2020)036
We observe that the u-channel contribution depends on the angular momentum mode of the eigenstate. Equation (A.4) can now be rewritten in the familiar order We consider scattering and bound states in two different Coulomb potentials The expectation value of the momentum of each particle in the CM frame in the scattering state, and the Bohr momenta for the scattering and bound states are For convenience, we define the parameters as well as the space variables The energy eigenvalues of the scattering and bound states are and the corresponding wavefunctions are 8

JHEP02(2020)036
where S 0 (ζ S ) ≡ 2πζ S /(1 − e −2πζ S ), and in eq. (B.6b), we have expressed the bound state wavefunction in terms of the confluent hypergeometric function 1 F 1 rather than the Laguerre polynomials. Note that the wavefunctions (B.6) assume distinguishable interacting particles. We include the necessary (anti)symmetrization factors for identical particles in sections 2 and 3, where we discuss the processes of interest. 9

B.2 Integral
We are interested in computing the overlap integral for the wavefunctions of appendix B.1. Note that the prefactor in eq. (B.7) has been chosen such that R k,n m is dimensionless. Substituting eqs. (B.6) into (B.7) and performing the angular integration, picks out the S = mode of the scattering state. Then, setting t ≡ 2x B /n = 2ζ B x S /n, we obtain Note that the hypergeometric function in eq. (B.11) is a finite polynomial in its last argument because its first argument of is a non-positive integer, 1 + − n 0. The last factor in eq. (B.11) is an unimportant overall phase. Combining eqs. (B.8) and (B.11), we find m=− dΩ k |R k,n m | 2 = 2 4( +2) π 2 n 4 2 + 1 (B.12) R k,n m vanishes if ζ B = ζ S , as expected from the orthogonality of the wavefunctions.

C Bound-state formation via vector emission
The BSF V amplitude is [10] iM V k→n m where the leading order contributions to the perturbative transition amplitude A V T (k , p) are shown in figure 21. Note that the third diagram does not appear in more minimal U(1) theories where the interacting particles do not couple to a doubly charged scalar. This diagram is akin to the one that appears in non-Abelian theories, where the final-state gluon is radiated from a gluon exchanged between the interacting particles, via the trilinear gluon vertex. Such diagrams seem naively to be of higher order than those with emission directly from one of the interacting particles. However, the momenta exchanged along the propagators scale with the couplings and render these diagrams of the same order as those with emission from the legs [22]. See also ref. [5] for analogous contributions arising from couplings in the scalar potential.

JHEP02(2020)036
Note that eq. (C.9) holds also for fermionic DM. This has been already established for the two diagrams on the left of figure 21. In the third diagram, the contraction of two spinor pairs due to the two fermion-antifermion-scalar vertices gives rise to an extra factor 2 2 . Combined with the fact that y 2 = 4πα Φ for fermions [5, appendix A], reproduces eq. (C.9). Despite the different potentials in the initial and final states -which in fact occur commonly in BSF V in non-Abelian theories where the emitted gauge vector boson carries away non-Abelian charge [8] -BSF V is suppressed by α 2 B with respect to the BSF Φ processes computed in section 2, due to the momentum dependence of the vector emission vertices. (For the third diagram of figure 21, and its interference with the other two diagrams, the suppression is actually of order α 2 Φ and α Φ α B respectively, at the level of the cross-section.) 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.