Dark matter bound states via emission of scalar mediators

If dark matter (DM) couples to a force carrier that is much lighter than itself, then it may form bound states in the early universe and inside haloes. While bound-state formation via vector emission is known to be efficient and have a variety of phenomenological implications, the capture via scalar emission typically requires larger couplings and is relevant to more limited parameter space, due to cancellations in the radiative amplitude. However, this result takes into account only the trilinear DM-DM-mediator coupling. Theories with scalar mediators include also a scalar potential, whose couplings may participate in the radiative transitions. We compute the contributions of these couplings to the radiative capture, and determine the parameter space in which they are important.


Introduction
In a variety of theories, dark matter (DM) is hypothesised to couple to scalar mediators. Such a coupling may determine the DM density via thermal freeze-out, and/or connect the dark sector to the Standard Model (SM) particles. It is a particularly compelling possibility in view of the discovery of the Higgs boson, which may itself be the mediator, or provide a portal to DM via an extended scalar sector. If DM possesses a sizeable coupling to the SM Higgs, then it is constrained by current experiments to be significantly heavier than the Higgs. On the other hand, if DM couples to a non-SM scalar mediator that is much lighter than itself, then this coupling may result in significant DM self-scattering inside haloes that can affect the galactic structure and bring theoretical predictions in better agreement with observations [1]. Even outside the self-interacting DM regime, the coupling of DM to a lighter non-SM scalar mediator is a generic possibility within dark sector models.
In scenarios where DM couples directly to a light force carrier, non-perturbative effects associated with the long-range nature of the interaction impact the DM phenomenology.
It is well known that the Sommerfeld effect [2,3] can influence the DM annihilation and self-scattering rates. This can be the case either if the mediator is a non-SM scalar or the SM Higgs [4]. More recently, it has been realised that the cosmological and astrophysical formation of DM bound states is a generic implication of theories with light mediators. The formation and subsequent decay of unstable bound states can deplete the DM density [5], and contribute to the DM indirect detection signals [6][7][8][9][10][11][12][13][14][15]. The formation of stable bound states may quell the DM self-scattering inside haloes [16], and give rise to novel radiative [17][18][19][20] and direct detection signatures [21,22]. Attractive interactions may result in the formation of large bound states [23][24][25] or non-topological solitons [26][27][28][29].
While the importance of bound-state formation (BSF) for DM is now well established for weakly coupled theories with vector mediators [5-13, 16-18, 30-39], and of course in confining DM theories [40][41][42][43][44], in models with scalar mediators bound-state effects appear to be either less severe or relevant to a limited parameter space. In such models, the radiative capture into bound states suffers from two cancellations: (i) The lowest order s-wave contribution vanishes due to the orthogonality of the incoming and outgoing wavefunctions, thus abdicating the leading order to the p-wave. (ii) The leading order p-wave terms cancel for particle-antiparticle or identical-particle pairs, yielding to s-and d-wave contributions that are suppressed by higher orders in the coupling with respect to the capture via vector emission and to the annihilation cross-section [9,24,31]. Even so, BSF has some important implications. Since it contains an s-wave component, it results in CMB constraints for fermionic DM, whose direct annihilation is p-wave and thus unconstrained by indirect probes [15]. Moreover, asymmetric DM coupled to a light scalar may form stable multiparticle bound states, provided that the coupling is large enough to overcome the two-particle capture bottleneck [24,25,45] (or BSF is facilitated by a parametric resonance [46]).
The above results take into account only the trilinear DM-DM-mediator coupling. This coupling is responsible for the long-range interaction between the DM particles that determines the scattering and bound state wavefunctions. The same coupling contributes to the radiative part of the capture process. However, there may be also other couplings that contribute to this piece. Since the radiative vertex enters the overlap integral between the initial and final states, its strength and momentum dependence are critical in determining the efficacy of BSF. It is thus important to consider all relevant contributions, particularly if they are generic within a theory. Theories with scalar mediators include a scalar potential whose couplings can contribute to the radiative part of the capture process. In this work, we investigate the impact of the various couplings in the scalar potential to the radiative BSF in theories with scalar mediators.
The paper is organised as follows. In section 2, we introduce the interaction Lagrangians, review the computation of radiative transitions, and summarise the past results on BSF with scalar emission via the trilinear DM-DM-mediator coupling. In section 3, we compute the contributions to radiative BSF from other scalar couplings. We consider both scalar and fermionic DM, and compute BSF via one and two scalar emission (BSF1 and BSF2 respectively). We discuss the features of the resulting cross-sections and compare them with the past results. We conclude in section 4, with a discussion of their phenomenological implications. Various technical computations are included in the appendices. 2 Bound-state formation via emission of scalar mediators

Preliminaries
We consider two particles X 1 and X 2 with masses m 1 and m 2 respectively, that interact via a light scalar mediator ϕ of mass m ϕ . We shall allow X 1 and X 2 to be real or complex scalars, or Dirac fermions. The relevant interactions are described by the following Lagrangians, and Note that, since we are interested in the application of our results to DM, we shall assume that the interacting particles X 1 , X 2 carry a Z 2 symmetry. For later convenience, we define the total and the reduced mass of the two interacting particles M " m 1`m2 , µ " m 1 m 2 m 1`m2 , (2.4) and the dimensionless factors In the non-relativistic regime, the interaction between X 1 and X 2 is described to leading order by a static Yukawa potential that arises from the resummation of the oneboson-exchange diagrams, with α " α sc or α " α f , depending on whether the interacting particles are scalars or fermions, where α sc " g 1 g 2 16π and α f " (2.7) We derive the Yukawa potential and α sc , α f in appendix A. The long-range interaction between X 1 and X 2 described by the potential (2.6) distorts the wavefunction of the scattering (unbound) states -a phenomenon known as the Sommerfeld effect [2,3] -and gives rise to bound states. Bound states exist if the mediator is sufficiently light. For the ground state to exist, µα{m ϕ ą 0.84 , (2.8) while stronger conditions apply for excited states [9]. The condition (2.8) also roughly marks the regime where the Sommerfeld effect is significant. The capture into bound states necessitates the dissipation of the binding energy and the kinetic energy of the relative motion of the X 1 X 2 pair, which may occur radiatively. In section 2.2, we review the computation of radiative BSF amplitudes. The capture via emission of one scalar mediator, X 1`X2 Ñ BpX 1 X 2 q`ϕ, has been previously considered in Refs. [9,15,24,31], where only the trilinear ϕX 2 j couplings were taken into account. We review the main results in section 2.3. As we shall see, for a particle-antiparticle pair or a pair of identical particles, the dipole contribution -which is the leading order term for X 1 , X 2 with different masses and couplings -vanishes identically.

Radiative transition amplitude
We consider the radiative transitions X 1 pk 1 q`X 2 pk 2 q Ñ X 1 pp 1 q`X 2 pp 2 q`radiation , (2.9) where the parentheses denote the 4-momenta of the incoming and outgoing X 1 and X 2 fields. We will be interested in particular in the case where the incoming X 1 , X 2 particles form a scattering state, while the outgoing X 1 , X 2 are captured into a bound state. In order to separate the motion of the CM from the relative motion, we make the following transformation in the momenta [31,47] k 1 " η 1 K`q, k 2 " η 2 K´q , (2.10) where η 1,2 are defined in eq. (2.5).
In the presence of a long-range interaction, the relative motion of X 1 , X 2 is not well approximated by a plane wave. It is described more generally by wavefunctions, which in momentum space we shall denote asφ k pqq andψ n m ppq for the scattering and the bound states respectively. The wavefunctions obey the Schrödinger equation with the potential (2.6). The continuous spectrum is characterised by the momentum k " µv rel , which is the expectation value of q and parametrises the energy of the relative motion in the scattering states, k " k 2 {p2µq " µv 2 rel {2. The bound states are characterised by the standard discrete principal and angular momentum quantum numbers tn mu, which determine the expectation value of p and the binding energy n . As is well known, for a Coulomb potential, n " κ 2 {p2n 2 µq " µα 2 {p2n 2 q, where κ " µα is the Bohr momentum; a non-negligible mediator mass suppresses n and introduces a dependence on . We review the wavefunctions in appendix B (see Ref. [9] for a more detailed discussion). For the purpose of evaluating the leading order contributions to the transition amplitude, we shall keep in mind that the wavefunctions impose |q| " |k| " µv rel and |p| " κ " µα.
In the non-relativistic regime, the total 4-momenta of the scattering and the bound states are We will work in the CM frame, K " 0. Then, taking into account that k , n ! M (or equivalently α, v rel ! 1), the total energy available to be dissipated is (2.14) Evidently, the bound state acquires momentum |P| " ω. The full amplitude for the radiative capture into a bound state depends on the overlap of the initial (scattering) and final (bound) state wavefunctions and the radiative vertex. The diagrammatic representation for transitions with emission of one or two scalars is shown in fig. 1. In the instantaneous and non-relativistic approximations, the amplitude is [31] M kÑtn mu » 1 ? 2µ 15) where M T pq, pq is the radiative amplitude for the (off-shell) transition (2.9), under the transformations of eqs. (2.10) and (2.11). The fully connected diagrams contributing to M T can be evaluated at leading order by setting the incoming and outgoing X 1 , X 2 onshell. For any non-fully-connected diagrams contributing to M T , the virtuality of X 1 and X 2 has to be integrated out. This can be done starting from the off-shell amplitude as described in [31] (see also [34, section 2.3] for a brief summary), or by adopting an effective field theory approach [48][49][50][51][52][53][54][55][56][57][58][59]. 1 The dominant contributions to the capture with emission of a vector or scalar mediator via the trilinear coupling arise from non-fully-connected diagrams [9,31]. (However, in non-Abelian theories, a leading order contribution to capture via gluon emission arises also from a fully connected diagram [8].) In the computations of this paper in section 3, we will consider only fully connected diagrams.
In eq. (2.15), the factor inside the square brackets includes the leading order corrections in p 2 , q 2 arising from the relativistic normalisation of states [31]. Upon the convolution with the wavefunctions and integration over p and q, these terms amount to corrections in α 2 and v 2 rel . They become important when the leading order contribution from M T alone cancels, as is the case for the capture of particle-antiparticle or identical-particle pairs with emission of a scalar mediator via the trilinear coupling.

Capture with scalar emission via the trilinear coupling
The leading contributions to capture with emission of a scalar via the trilinear coupling only are shown in fig. 2. Starting from eq. (2.15), and neglecting the correction arising Figure 1. The amplitude for the radiative capture into bound states via emission of one or two scalars consists of the initial and final state wavefunctions, and the perturbative radiative amplitude M T that includes the radiative vertices. from the normalisation of states, these diagrams yield 2 [9, 31] where P ϕ is the momentum of the emitted boson, with P 2 ϕ`m 2 ϕ " ω 2 and ω given in eq. (2.14). The wavefunctions imply that the integrand is significant in the region r À 1{ maxpµα{n, µv rel q ! 1{ω, therefore we may evaluate eq. (2.16) by expanding in P ϕ¨r " maxpα{n, v rel q (cf. ref. [9, appendix B]). Clearly, the zeroth order terms vanish due to the orthogonality of the wavefunctions.
For a pair of particles with different masses and/or couplings to the scalar mediator, the leading order contributions arise from the P ϕ¨r terms in the expansion of eq. (2.16). The cross-section for capture into the ground state tn mu " t100u is [9,31] where S TC BSF1 depends on the dimensionless parameters α{v rel and µα{m ϕ . In the Coulomb approximation, which is valid at m ϕ À µv rel , S TC BSF1 depends only on the parameter ζ " α{v rel as follows [9,31] S TC BSF1 pζq "ˆ2 πζ 1´e´2 πζ˙ˆ2 For the more general case that includes the effect of the mediator mass, we refer to [9]. We note that eqs. (2.17) and (2.18) correspond to dipole emission ( S " 1 mode of the scattering state wavefunction). In eq. (2.17), the factor in the square brackets reduces to 1 for a pair of scalar particles with g 1 " g 2 and η 1 " η 2 . In eq. (2.18), the first factor inside the brackets is responsible for the characteristic σv rel 9 1{v rel scaling of the Sommerfeld-enhanced processes at low velocities (ζ Á 1), while the remaining factors tend to a constant. Evidently, for a particle-antiparticle pair or a pair of identical particles (g 1 " g 2 " g and η 1 " η 2 " 1{2), the two terms proportional to P ϕ in the expansion of eq. (2.16) cancel. This cancellation persists even for capture into bound states of non-zero angular momentum ( ą 0), and implies that the next order contributions should be considered.
For a scalar particle-antiparticle pair the next order terms, which include also the leading order correction from the normalisation of states shown in eq. (2.15), yield [9] where in the Coulomb regime S TC BSF1,XX˚d epends on ζ " α sc {v rel as follows (2.20) (For a pair of identical scalars, an extra factor of 2 arises from the symmetrization of the scattering state wavefunction.) The case of fermionic DM was considered in ref. [24]. In the Coulomb regime and for ζ Á 1, the cross-section was found to be Note that both eq. (2.19) and eq. (2.21) are suppressed by α 2 with respect to eq. (2.17).

The contribution of scalar couplings to the radiative capture
In this section, we investigate how the scalar couplings λ 1ϕ , λ 2ϕ , λ ϕ and ρ ϕ in eqs. (2.1) to (2.3) -which do not affect the long-range potential between X 1 and X 2 -contribute to the radiative capture of X 1 , X 2 pairs into bound states. In order to exhibit the leading order contribution from all couplings, we consider BSF via one and two scalar emission, It is important to note that in most diagrams we will consider, the trilinear couplings g j also participate in the radiative part of the process (cf. figs. 3 to 6). In fact, it may naively seem that these diagrams are of the same or higher order in α than those giving rise to the cross-sections of section 2.3 (even ignoring the additional suppression introduced by the new couplings). However, the momentum transfer along the mediators exchanged in these diagrams and the off-shellness of the interacting particles scale also with α, thereby reducing the order of dependence of the diagrams on α. 3

Capture via one scalar emission
For simplicity, in this section we will not include the contributions to BSF1 that arise from the trilinear DM-DM-mediator coupling alone, which were reviewed in section 2.3. For clarity, we first carry out our computation assuming that X 1 and X 2 are distinguishable scalars, and then discuss how the result is modified for fermionic and/or identical species.

Amplitude
The leading order diagrams in λ 1ϕ , λ 2ϕ and ρ ϕ that contribute to the radiative part of the transition amplitude due to the scalar couplings in eqs. (2.1) to (2.3) are shown in fig. 3. We find Then according to the discussion in section 2.2 on the scaling of the momenta, Inserting this in eq. (2.15), and neglecting the corrections arising from the relativistic normalisation of states, we find 4 (3.7) The prefactors in eqs. (3.6) and (3.7) have been chosen such that V k,tn mu and R k,tn mu are dimensionless. We derive analytical expressions for these integrals in appendix C, using an appropriate approximation. They can be expressed in terms of two dimensionless parameters As stated in section 2.1, for the Yukawa potential (2.6) the ground state exists if ξ ą 1 [9]. 5 Moreover, BSF1 is kinematically possible only if where we used the Hulthen approximation for the binding energy (cf. appendix B). In the regime where BSF is important, v rel ă α (see e.g. section 2.3), this condition is much stronger than the requirement ξ ą 1, and reduces roughly to m ϕ À µα 2 {2 or ξ Á ξ min » 2.4{α. For capture into the ground state, we find [cf. eqs. (C.10) and (C.11)] In the Coulomb limit, eq. (3.12) reduces to S C 0 " 2πζ{p1´e´2 πζ q. It is important to note that the contributions in eqs. (3.10) and (3.11) arise from the S " 0 mode of the scattering state wavefunction.

Cross-section for capture into the ground state
The cross-section for capture into the ground state is The momentum of the emitted scalar is found from the conservation of energy where we adopted the Coulomb value for the binding energy since ξ " 1. Using the amplitude (3.5) with the overlap integrals (3.10) and (3.11), we obtain where s BSF1 ps is the phase-space suppression factor for BSF1, We shall now simplify and adapt the above expression -which has been derived for a pair of distinguishable scalars -to various cases.
• For a pair of distinguishable scalars, taking into account that α " g 1 g 2 {p16πq, (3.17) • For a particle-antiparticle pair, g 1 " g 2 and λ 1ϕ " λ 2ϕ " λ Xϕ ; the above reduces to • For a pair of identical scalars, the total wavefunction has to be symmetric in the interchange of the two particles. Thus the scattering state wavefunction is This implies that the contribution of the even-S modes participating in a process is doubled with respect to the case of distinguishable scalars, while the contribution of the odd-S modes vanishes. Similarly, there are only " even bound states of two identical bosons. Since the contributions included in eq. (3.18) are s-wave, the cross-section for a pair of identical scalars is twice as large as that given by eq. (3.18).
For fermions, there are no renormalisable λ jϕ couplings and only the second term in eq. (3.5) contributes. The amplitude is related to that for scalars as follows. The Dirac propagators S D j can be expressed in terms of the scalar propagators S j ppq and the spinors u,ū where r denotes the spin and j " 1, 2 refers to the particle species. In order to compute the fermionic BSF diagrams, we insert a factor ř r i u r i j pp i qū r i j pp i q for each propagator and contractū r i`1 j pp i`1 qu r i j pp i q across each vertex i. Since all the fermion-fermion-scalar vertices in the BSF diagrams are either soft or ultrasoft, we can use the identitȳ as also in appendix A. Then: -The spin Kronecker deltas in eq. (3.21), upon summation over the internal spin indices, ensure that the spin of each particle is conserved across the entire diagram, including both the ladders and the vertices in the radiative parts of the diagrams.
-The factors m j combine with the fermionic vertex values p´ig j q to yield the scalar vertices p´ig j m j q.
-The factors of 2 along the ladder parts of the BSF2 diagrams are responsible for the difference between α sc and α f , as shown in appendix A. The two factors of 2 arising from the two vertices in the radiative parts of the diagrams imply that for fermionic pairs, M T is larger than for scalar pairs by a factor 4.
Based on the above and taking into account that α " g 1 g 2 {p4πq for fermions, the ρ ϕ contribution to the amplitude (3.5) turns out to be the same for scalars and fermions, when expressed in terms of α. Upon squaring the amplitude and averaging over the spin of the initial particles, we find the following.
• For a pair of distinguishable fermions, including a particle-antiparticle pair, • For a pair of identical fermions, the total wavefunction has to be antisymmetric in the interchange of the two particles. This implies that the spatial wavefunction depends on their total spin. A pair of identical spin-1/2 particles may be either in the antisymmetric spin-0 state, or in the symmetric spin-1 state. Their spatial wavefunction should then be symmetric or antisymmetric, respectively, fermions with total spin 0: As for identical bosons, the wavefunction (3.23) implies that the contribution of the even-S modes participating in a process is doubled with respect to the case of distinguishable particles, while the contribution of the odd-S modes vanishes. The opposite holds for the wavefunction (3.24). Similarly, there are only " even bound states of two identical total-spin-0 fermions, and only " odd bound states of two identical total-spin-one fermions. Since eq. (3.22) is an s-wave process, the spinaveraged cross-section for capture of two identical fermions is half of that given in eq. (3.22), with the entire contribution arising from the spin-0 state.

Discussion
We shall now discuss the features of the contributions from the scalar couplings to the BSF1 cross-sections, and compare them with the past results quoted in section 2.3.

Velocity scaling
In the Coulomb limit m ϕ Ñ 0, and at low velocities or equivalently for large coupling α ą v rel , the cross-sections of section 3.1.2 exhibit the characteristic velocity scaling of the Sommerfeld enhanced processes in the Coulomb regime, σ SC BSF1 v rel 9 1{v rel . 6 The Coulomb limit is a satisfactory approximation as long as the average momentum transfer is larger than the mediator mass µv rel Á m ϕ , or equivalently ξ Á ζ (cf. e.g. ref. [9]).
At v rel ă m ϕ {µ ă α, the cross-sections exhibits parametric resonances at discrete values of µα{m ϕ that correspond to the thresholds for the existence of bound states; in the Hulthen approximation (cf. appendix B) they occur at ξ " n 2 with n being a positive integer. For non-resonant parametric points σ SC BSF1 v rel saturates to its value at v rel « m ϕ {µ as the velocity decreases. Close to resonant points, it grows as σ BSF1 v rel 9 1{v 2 rel , to eventually saturate at a constant value at a lower velocity that depends on the proximity to the resonant point. The resonant features and velocity dependence of the cross-section are encapsulated in the s-wave Sommerfeld enhancement factor S 0 pζ, ξq of eq. (3.12).
The leading order contribution to BSF1 due to the trilinear coupling alone for particleantiparticle or identical particle pairs [cf. eqs. (2.19) and (2.21)] contains an s-wave component that exhibits the same velocity scaling for v rel ă α. However, the lower order contribution (2.17) to BSF1 for particles with different masses and/or couplings is p-wave and dwindles at low velocities.

Scaling with the couplings
We now compare the scaling of the cross-sections of sections 2.3 and 3.1.2 with the various couplings. We always refer to the regime v rel ă α where BSF is significant, and we set aside the common s-wave Sommerfeld factor S 0 and the other ζ-dependent factors that result in the same velocity dependence at v rel ă α. We focus on particle-antiparticle or identical-particle pairs. 6 This scaling also appears in the partial-wave unitarity limit on the inelastic cross-sections, which in the non-relativistic regime reads σ pJq inel v rel ď σ pJq uni v rel " πp2J`1q{pµ 2 v rel q with J being the partial wave [61]. This observation implies that the unitarity limit may be approached or realised only by Sommerfeld enhanced processes [5,12].
• For scalar DM, the contribution of the DM-mediator quartic coupling to BSF1 scales as σ SC BSF1 v rel 9 α 2 λ 2 Xϕ [cf. eq. (3. 18)]. Comparing this with the contribution of the trilinear coupling alone [cf. eq. (2.19)], we find that it dominates for λ Xϕ Á 18α . (3.25) • The contribution to BSF1 from the trilinear mediator self-coupling -which arises both for bosonic and fermionic DM -scales as σ SC BSF1 v rel 9 p1{αqpρ ϕ {µq 2 . This is clearly of lower order in α than the cross-sections of section 2.3, however it is suppressed by the square of the ratio of the mediator scale to the DM scale. Assuming that ρ ϕ " m ϕ , the kinematic threshold m ϕ À µα 2 {2 implies at best σ SC BSF1 v rel 9 α 3 . This is still of lower order in α than eqs. (2.19) and (2.21). It is of course important to note that ρ ϕ may differ significantly from m ϕ and/or the binding energy. 7 We find that the ρ ϕ contribution dominates over those of eqs. (2.19) and (2.21) if which encompasses significant parameter space.

Capture via two scalar emission
The scalar couplings appearing in eqs. (2.1) to (2.3) raise the possibility that the radiative part of the capture into bound states may carry a lower suppression in α if two mediators are emitted. On the other hand, the emission of two mediators, which share the available energy ω » µpα 2`v2 rel q{2, implies that the BSF2 cross-section picks up additional suppression in powers of α (for α ą v rel ) due the phase space of the second mediator. To determine which diagrams may contribute significantly, we first work out the phase-space integration. 8

Phase-space integration
The cross-section for radiative capture via emission of two scalar mediators with momenta P a and P b is (3.27) 7 We note that the third diagram of fig. 3, from where the ρϕ contribution to BSF1 arises, resembles the diagram that appears in the radiative capture via one gluon emission in non-Abelian theories due to the trilinear gluon coupling [8]. That coupling is momentum dependent, with the relevant momentum scale in the capture process being κ. We may recover the scaling of the non-Abelian BSF cross-section on α [34] from eq. (3.22) by mapping ρϕ Ñ g κ 9 α 3{2 . 8 In a previous version of this paper, it was stated that the contributions to BSF2 from the trilinear DM-DM-mediator couplings gj alone do not suffer from the cancellations of the gj contributions to BSF1 shown in fig. 2. However, the resummation of the one-scalar-exchange diagrams between the two emission vertices in the corresponding BSF2 diagrams reveals that the same cancellations are present. Thus, these contributions to BSF2 are not of interest, and we do not consider them here. We thank Joan Soto for pointing this out.
where the factor p1{2q is due to the two identical particles in the final state. We work in the CM frame, K " 0, and use the three-momentum delta function to integrate over the momentum of the bound state P "´pP a`Pb q. The energy delta function yields the For convenience, we define the dimensioneless parameters Then, combining energy and momentum conservation, we obtain the phase space condition for the momenta of the emitted scalars, where θ ab is the angle between P a and P b . Because δ ! 1, the phase-space encompassed by eq. (3.31) extends on the x a´xb plane essentially along the line with 0 ď x a ď ? 1´2d, and a small width along the x b direction, due to cos θ ab ranging in r´1, 1s. Note that for the capture process to be kinematically allowed, 2d ă 1. The x b width can be estimated by differentiating eq. (3.31) with respect to x b and cos θ ab . We find (3.33) For the diagrams of interest, the amplitude is independent of P a,b (cf. section 3.2.2). Thus, putting everything together, eq. (3.27) yields where s BSF2 ps is the phase-space suppression factor due to the kinematic threshold for BSF2, with s BSF2 ps " 1 for d " 0. Comparing eq. (3.34) with eq. (3.13), we observe that the BSF2 cross-section is proportional to two extra powers of the available energy to be dissipated, ω » n ` k , with respect to BSF1, as expected from the phase-space element for the second mediator d 3 P b {p2P 0 b q 9 ω 2 . This implies an extra suppression by α 4 in the regime where BSF is important, α ą v rel . We also note the suppression of BSF2 with respect to BSF1 by the numerical factor 1{p2 4 3π 2 q » 2ˆ10´3. Figure 4. Capture via emission of two scalar mediators: the contribution of the mediator quartic self-coupling to the radiative amplitude.

Amplitude and cross-section
In figs. 4 to 6, we show various contributions to BSF2 arising from the mediator selfcouplings and the mediator-DM quartic coupling for scalar DM. The diagram of fig. 4 yields the most important contribution. We shall first compute the cross-section arising from this diagram and then discuss why we neglect the other contributions. The amplitude for the λ ϕ contribution can be obtained from the second term in eq. (3.5) by replacing ρ ϕ with λ ϕ . Taking into account the discussion in section 3.1.2, we find it to be the same for both bosons and fermions, where R k,tn mu is defined in eq. (3.7). From eqs. (3.11) and (3.34) we obtain the (spinaveraged) cross-section for capture into the ground state, The cross-section (3.37) becomes comparable to the ρ ϕ contribution to BSF1 computed in section 3.1.2 for λ ϕ 4 ? 3π Comparing with the BSF1 contributions from the trilinear coupling alone (2.19) and (2.21), the λ ϕ contribution to BSF2 is more significant for The perturbativity of the couplings implies that the above condition may be meaningfully satisfied only for small coupling α À 10´4. Figure 5. Capture via emission of two scalar mediators: the contributions of the mediator trilinear self-coupling to the radiative amplitude. Figure 6. Capture via emission of two scalar mediators: the contribution of the quartic DMmediator coupling to the radiative amplitude. These diagrams exist for scalar DM only.

Other contributions
In figs. 5 and 6, we display various diagrams arising from ρ ϕ and λ jϕ . Their contributions to BSF2 are subdominant compared to the contributions of the same couplings to BSF1 for reasonable choices of the parameters, as we discuss below. Since for phenomenological purposes, we are interested mostly in the total capture rate, we do not compute these contributions in detail. Because part of the suppression arises from the phase space density of the second mediator (cf. section 3.2.1), we will consider directly the contributions of these diagrams to the cross-section rather than the amplitude, thus neglecting any cross-terms between different diagrams (except in the cases where there is a cancellation between them). It is straightforward to verify that the cross-terms do not alter our conclusions. We focus on particle-antiparticle or identical-particle pairs. As in the discussion of section 3.1.3, our comparisons refer to the regime α ą v rel where BSF is important, setting aside the ζ-dependent factors that yield the same v rel dependence between BSF1 and BSF2 in this regime.
• For particle-antiparticle or identical-particle pairs, the sum of the two diagrams on the left in fig. 5 suffers from the same cancellation as the diagrams of fig. 2  (cf. section 2.3). This introduces a α 4 suppression in the corresponding cross-section, while the phase-space density of the second emitted mediator implies an additional suppression by pµα 2 {2q 2 as discussed in section 3.2.1. These suppressions are balanced by the propagator of the off-shell mediator, which introduces a factor " pP 0 a`P 0 b q´4 9 pµα 2 {2q´4. Thus, the contribution of these diagrams to the crosssection scales as σ BSF2 v rel 9 pρ ϕ {µq 2 . This is suppressed by α with respect to the ρ ϕ contribution to BSF1, σ BSF1 v rel 9 p1{αqpρ ϕ {µq 2 [cf. eq. (3.22)], as well as by the numerical factor " 2ˆ10´3 due to the three-body phase space.
• The third and forth diagrams in fig. 5 can be estimated starting from the ρ ϕ contribution to BSF1 [cf. third diagram in fig. 3 and eq. (3.22)].
For the third diagram, the off-shell mediator yields a factor 9 pµα 2 {2q´4. Taking into account the phase-space suppression 9 pµα 2 {2q 2 , we find that the contribution to the BSF2 cross-section scales as σ BSF2 v rel 9 p1{α 5 qpρ ϕ {µq 4 . Assuming that ρ ϕ " m ϕ À µα 2 {4, effectively this yields at best σ BSF2 v rel 9 α 3 , similarly to the BSF1 cross-section (3.22) as discussed in section 3.1.3. Still, the numerical suppression of BSF2 due to the three-body phase space ensures that the BSF2 contribution is subdominant to BSF1. However, if there is a large hierarchy between m ϕ and ρ ϕ with ρ ϕ " m ϕ , then the diagram under consideration may be significant. 9 For the forth diagram, the off-shell mediator between the two emitted scalars yields a factor 9 pµαq´4. It follows that it is subdominant to the third diagram in fig. 5 and to the ρ ϕ contribution to BSF1.
• For the two diagrams on the left in fig. 6, we start from the first term of eq. (3.18), replace one factor of α with λ 2 Xϕ and introduce the phase-space suppression 9 α 4 . We obtain σ BSF2 v rel 9 α 5 λ 4 Xϕ .
• The two diagrams on the right in fig. 6 suffer from the cancellations of the BSF1 diagrams of fig. 2. Replacing one factor of α with λ 2 Xϕ results in the suppression λ 2 Xϕ α 3 for particle-antiparticle or identical-particle pairs. Introducing also the phasespace suppression yields σ BSF2 v rel 9 λ 2 Xϕ α 7 .

Conclusion
Dark matter coupled to a light scalar mediator is motivated in the context of Higgs portal models, self-interacting DM, and more generally models beyond the SM with extended scalar sectors. Light force mediators generically imply the existence of DM bound states, whose formation has important implications for the DM phenomenology. While the cosmological and astrophysical formation of bound states is significant in models with vector mediators, it is less efficient in models with scalars mediators. Indeed, considering only a trilinear DM-DM-mediator coupling, the radiative capture of particle-antiparticle or identical-particle pairs into bound states with emission of a scalar mediator is subject to cancellations that suppress the cross-section by higher powers in this coupling. While the trilinear coupling determines the long-range interaction between the DM particles, it is not necessarily the only contributor to the radiative part of the transition. Quite generically, in theories with scalar mediators the couplings in the scalar potential also contribute. In this work, we investigated the contribution of the mediator self-couplings, as well as the DM-mediator quartic coupling in the case of scalar DM, to the radiative capture into bound states. We considered capture both via one and two scalar emission. We have found that these couplings can enhance or dominate the capture rate in sizeable parts of the parameter space, particularly for small values of the trilinear DM-DM-mediator coupling.
Importantly, the new contributions are s-wave and remain significant even at very low velocities, thereby enhancing the radiative signals that can be probed in the indirect searches. This can potentially strengthen the resulting constraints. While models with light mediators and s-wave annihilation are severely constrained by the CMB and other indirect probes [10,62,63], models with p-wave annihilation -such as fermionic DM coupled to a scalar mediator -remain largely unconstrained. In this case, the formation and subsequent decay of unstable bound states offers a source of detectable indirect signals [15]. (Bremsstrahlung of dark mediators has also been invoked to lift the p-wave suppression [64].) Moreover, s-wave annihilation processes may lead to a period of reannihilation in the early universe, after DM kinetic decoupling [65].
Besides the signals produced from the decay of unstable bound states, the radiation emitted during the capture process is another source of indirect signals even in the case of stable bound states [19]. This is particularly important for asymmetric DM [66], whose annihilation signals are suppressed due to the asymmetry (albeit can still be significant due to the Sommerfeld effect [12]). Asymmetric DM, in turn, offers an excellent host of self-interacting DM, both because it allows for large couplings to light mediators, and it evades the indirect detection constraints, provided that the asymmetry is sufficiently large [13]. The enhanced radiative signals expected due to the contributions computed here can improve the prospects of probing asymmetric and self-interacting DM.
Moreover, for asymmetric DM coupled to a light scalar, the enhanced rate of formation of two-particle bound states can result in more efficient cosmological formation of multiparticle bound states [24], which in the case of scalar DM may lead to solitosynthesis [67][68][69]. We leave detailed phenomenological studies for future work.

Appendices A Non-relativistic potential for scalars and fermions
The non-relativistic static potential describing the X 1´X2 interaction is (see e.g. [31,70]) where M 2PI pkq encompasses all 2-particle-irreducible diagrams that contribute to the X 1X 2 elastic scattering. To leading order, this is the one boson exchange shown in fig. 7. k " p´p 1 Figure 7. One boson exchange diagrams that yield the leading order contribution to the nonrelativistic potential between two different paticles X 1 and X 2 (left), or a particle-antiparticle pair (right). The parentheses denote the spin of the incoming and outgoing particles, in the case of fermions.
From the Lagrangians of eqs. (2.1) to (2.3), we find for the interaction of a scalar pair X 1 X 2 , a fermionic pair X 1 X 2 and a fermion-antifermion pair XX, where we made the approximation k 0 ! |k|. with α " α sc or α " α f depending on whether the interacting particles are scalars or fermions respectively, where α sc " g 1 g 2 16π and α f " The interaction is attractive if g 1 g 2 ą 0.

B Scattering state and bound state wavefunctions
The bound state and scattering state wavefunctions obey the Schrödinger equations Coulomb limit, m ϕ Ñ 0. The ground-state and scattering state wavefunctions are where 1 F 1 is the confluent hypergeometric function of the first kind, and (B.5) φ C k prq can also be decomposed in angular momentum modes as follows Outside the Coulomb regime: Hulthen approximation. For the Yukawa potential (2.6), the solutions of the Schrödinger equation are not known in a closed form. However, it is possible to obtain analytical solutions for the " 0 modes of the wavefunctions, if we replace the Yukawa for the Hulthen potential [71,72], For m˚" m ϕ , the Hulthen potential reproduces the behavior of the Yukawa potential at short and large distances. Both potentials admit bound state solutions provided that the screening length scale is sufficiently large. The thresholds for the n-th bound level in the Hulthen potential are m˚ď 2µα{n 2 , while for the Yukawa potential (2.6), bound states exist if m ϕ ď µα{0.84 [9]. We shall pick m˚" 1.68m ϕ , (B.8) such that the conditions for the existence of the lowest bound state coincide. The solutions to the Schrödinger equation for the Hulthen potential can be expressed in terms of the two dimensionless parameters ζ " α{v rel and ξ " 2µα{m˚. (B.9) Note that ζ and ξ are always defined using the appropriate fermionic or scalar coupling α. ζ compares the average momentum transfer between two unbound particles (" µv rel ) with the relative average momentum of the particles inside the bound state (" µα), while ξ compares the Bohr momentum (κ " µα) that determines the size of the bound state, with the screening scale (m˚) that determines the range of the interaction. The interaction manifests as long-range roughly if ξ Á 1; this is the regime where non-perturbative phenomena, the Sommerfeld effect and bound states, emerge. The Coulomb limit is recovered for ξ Ñ 8.
For the ground state, the binding energy is For the scattering state, we make the partial-wave decomposition φ k prq " a S 0 pζ, ξq where we define S 0 pζ, ξq below. It is possible to find an analytical solution for S " 0 only, which suffices for our purposes, which reduces to eq. (B.5) in the limit ξ Ñ 8.

C Overlap integrals
For the capture into bound states via emission of two scalars, we need to compute the overlap integrals k pqqψn m ppq rpq´pq 2`m2 ϕ s 2 .

(C.2)
We Fourier transform the wavefunctions as follows ψ n m prq " ż d 3 p p2πq 3ψ tn mu ppq e ipr ,ψ n m ppq " In the second step in eq. (C.6), we expanded the decaying exponential inside the integral. The bound-state wavefunction implies that the integrand is significant for r À n{pµαq while the kinematic threshold (3.9) for BSF imposes m ϕ ă µα 2 {p2n 2 q; therefore m ϕ r À α ! 1. The zeroth order term in the expansion vanishes due to the orthogonality of the wavefunctions, leaving the first order term to be the dominant contribution. In the following, we will focus on the capture into the ground state, tn mu " t100u. While in the Coulomb limit (m ϕ Ñ 0) it is possible to obtain analytical expressions for the integrals eqs. (C.5) and (C.6), outside the Coulomb regime they should be evaluated numerically. However, it is possible to obtain an analytical approximation, using the Hulthen potential (cf. appendix B). Since ξ " 2µα{m˚" 1 in all the parameter space where BSF is kinematically allowed to occur [cf. eq. (3.9)], the bound state wavefunction ψ 100 prq can be well approximated by its Coulomb value. On the other hand, the scattering-state wavefunction φ k prq is close to its Coulomb limit for r À 1{m˚, up to an overall normalisation which is determined by the factor S 0 pζ, ξq of eq. (B.17). Since ξ " 1, this encompasses all the range in which ψ 100 prq and therefore the integrands in eqs. (C.5) and (C.6) are important, r À 1{pµαq. It follows that V k,t100u and R k,t100u can be well approximated by their Coulomb values even outside the Coulomb regime provided that we replace S C 0 pζq with S 0 pζ, ξq.