Bound states of WIMP dark matter in Higgs-portal models. Part I. Cross-sections and transition rates

We investigate the role of the Higgs \emph{doublet} in the thermal decoupling of multi-TeV dark matter coupled to the Weak interactions of the Standard Model and the Higgs. The Higgs doublet can mediate a long-range force that affects the annihilation processes and binds dark matter into bound states. More importantly, the emission of a Higgs doublet by a pair of dark matter particles can give rise to extremely rapid monopole bound-state formation processes and bound-to-bound transitions. We compute these effects in the unbroken electroweak phase. To this end, we consider the simplest renormalisable fermionic model, consisting of a singlet and a doublet under $SU_{L}(2)$ that are stabilised by a $\mathbb{Z}_2$ symmetry, in the regime where the two multiplets coannihilate. In a companion paper, we use the results to show that the formation of metastable bound states via Higgs-doublet emission and their decay decrease the relic density very significantly.


Introduction
Particles coupled to the Weak interactions of the Standard Model (SM), known as WIMPs, have been arguably the most widely considered candidates for dark matter (DM) in the past decades. Among the archetypical WIMP models are scenarios where DM is a linear combination of the neutral components of electroweak multiplets that couple to the Higgs doublet. The discovery of the Higgs boson and the measurement of its properties impel the investigation of its implications for such scenarios. While most related research in the past focused on electroweak-scale WIMP masses, m ∼ 100 GeV, the current experimental constraints motivate considering the multi-TeV mass regime more carefully. This is particularly important in view of the numerous existing and upcoming observatories probing high-energy cosmic rays, such as H.E.S.S., IceCube, CTA and KM3Net. The experimental exploration of the multi-TeV scale urges the comprehensive theoretical understanding of the dynamics and possibilities in this regime.
The hierarchy between the multi-TeV and electroweak scales implies the emergence of new effects. In particular, the Weak interactions between particles with multi-TeV mass manifest as long-range, since the interaction range l ∼ (100 GeV) −1 may be comparable or exceed the de Broglie wavelength (µv rel ) −1 and/or the Bohr radius (µα) −1 of the interacting particles, where µ = m/2 TeV, v rel and α are the reduced mass, relative velocity and coupling to the force mediators. The long-range nature of the interactions gives rise to non-perturbative phenomena, the Sommerfeld effect and the existence of bound states.
Bound states form invariably with dissipation of energy. It has been recently shown that the emission of a scalar boson charged under a symmetry alters the effective Hamiltonian between the interacting particles and gives rise to monopole transitions; this renders bound-state formation (BSF) extremely rapid even for small couplings [1]. For WIMPs, this implies that BSF via emission of a Higgs doublet may be a very significant inelastic process. Moreover, it has been shown in simplified models, that the 125 GeV Higgs boson can mediate a sizeable long-range force between TeV-scale particles, despite being heavier than all SM gauge bosons [2,3].
The phenomenological importance of the above is rather large. It has been long known that the Sommerfeld effect [4,5] -the distortion of the wavefunction of interacting particles due to a long-range force -affects the DM annihilation cross-sections at low relative velocities [6]. This, in turn, alters the DM chemical decoupling in the early universe, and consequently the predicted DM mass and couplings to other particles [7]. It also affects the radiative signals expected from the DM annihilations during CMB and inside galaxies today [8]. More recently, it has been realised that the formation and decay of metastable (e.g. particle-antiparticle) bound states in the early universe can decrease the DM density [9], and contribute to the DM indirect detection signals [10][11][12][13][14][15][16][17]. Importantly, BSF can be faster than annihilation in a variety of models [1,3,[18][19][20], and it can also produce novel indirect signals even in models where annihilation is absent or suppressed [21][22][23][24][25][26][27].
Even neglecting bound state effects during the DM freeze-out, many of these models are predicted to reside in the TeV scale (see [37] for a summary). BSF in the early universe increases the DM destruction rate and, for a given set of couplings, pushes the predicted DM mass to even higher values [9]. If DM is heavier than 5 TeV, freeze-out begins before the electroweak phase transition (EWPT), even though it may be completed much later. Here we will assume electroweak symmetry throughout; this considerably simplifies the anyhow lengthy calculations, and is a practice that has been followed in past literature [37]. We discuss in detail the extent of validity of this approximation for the DM thermal decoupling in ref. [28]. Here we simply note that, by virtue of the Goldstone boson equivalence theorem, BSF via Higgs-doublet emission remains important in the broken electroweak phase, where it occurs via emission of the Higgs boson and the longitudinal components of the Weak gauge bosons [47].
In the epoch prior to the EWPT, if the two DM multiplets are nearly mass-degenerate (perhaps as a result of a larger symmetry group), the DM freeze-out is determined by all their self-annihilation and co-annihilation processes [48]. Our computations focus on this regime. Then, in the class of models considered, the entire Higgs doublet contributes to the (long-range) potential between the DM multiplets. Moreover, bound states can form via emission of a Higgs doublet. As we shall see, the cross-sections for BSF via Higgsdoublet emission can exceed those for annihilation or BSF via vector emission by orders of magnitude, and result in a late period of DM chemical recoupling [1]. This, in turn, opens the possibility that thermal-relic WIMP DM may be much heavier than anticipated, potentially approaching the unitarity limit [28], which in the non-relativistic indeed implies (or presupposes) the presence of long-range interactions [15].
This work is organized as follows. In section 2, we introduce the DM model and derive the long-range potentials between the DM multiplets in the unbroken electroweak phase. We identify the scattering and bound eigenstates of these potentials by appropriate spin and gauge projections, before computing the DM annihilation processes and bound-state decay rates. In section 3, we calculate all the radiative BSF cross-sections, while in section 4 we consider BSF via scattering on a relativistic thermal bath. These results are employed in the companion paper [28] to compute the DM thermal decoupling in the early universe. We conclude in section 5. Several technical aspects are addressed in the appendices, including a proof of the (anti)symmetrisation of the two-particle-irreducible kernels and wavefunctions in the case of identical particles (appendix A), and the analytical computation of monopole bound-to-bound transitions in the Coulomb limit (appendix C.) 2 Long-range dynamics in the unbroken electroweak phase

The model
We introduce a gauge-singlet Majorana fermion S = (ψ α , ψ †α ) T of mass m S , as well as a Dirac fermion D = (ξ α , χ †α ) T of mass m D with SM gauge charges SU L (2)×U Y (1) = (2, 1/2). We assume that S and D are odd under a Z 2 symmetry that leaves all the SM particles unaffected. Under these assignments, the new degrees of freedom (dof) allow to extend the SM Lagrangian by the following interactions where H is the SM Higgs doublet of mass m H and hypercharge Y H = 1/2, and D L ≡ P L D = (ξ α , 0) T and D R ≡ P R D = (0, χ †α ) T , with P R,L = (1 ± γ 5 )/2 being the right-handed and left-handed projection operators. In the above, D µ ≡ ∂ µ − ig 1 Y B µ − ig 2 W a µ t a is the covariant derivative, with t a = 1 2 (σ 1 , σ 2 , σ 3 ) and σ being the Pauli matrices. The particle content of eq. (2.1) is summarised in table 1. Here, we have taken the mass parameters m S and m D to be real. This can always be achieved by rephasing ψ and either ξ or χ. Rephasing the remaining spinor eliminates the phase of one of the Yukawa couplings. Thus the free parameters of the present model are 4 real couplings (two masses and two dimensionlesss Yukawa couplings), and a phase that allows for CP violation.
We are interested in the regime in which S and D can co-annihilate efficiently before the EWPT of the universe, which occurs if their masses are similar, within about 10%. This is because the number density of the heavier species in the non-relativistic regime is suppressed with respect to that of the lighter species by a factor exp[−(δm/m)x]. Since x ≡ m/T 25 during DM freeze-out, if δm/m 10% the co-annihilations and the selfannihilations of the heavier species are typically subdominant to the self-annihilations of the lightest species, although their relative importance depends also on the corresponding cross-sections. Such a small discrepancy in the masses does not significantly affect (most of) the cross-sections that we compute in the following; we shall thus take the masses to be equal, which greatly simplifies the computations and allows to obtain analytical results, m D = m S ≡ m. (2.2) We will very often use the reduced mass of a pair of DM particles, µ ≡ m/2.
Moreover, in order to reduce the number of free parameters, we set which we take to be real. (The CP violation is anyway not important for our purposes.) Our computations can of course be extended to more general Yukawa couplings. As is standard, we define the couplings Various aspects of its phenomenology and experimental constraints of the model (2.1) have been considered extensively in the past [29][30][31][32][33][34][35][36][37][38], and we shall not review them here. We only note that after the EWPT, S and D mix to produce three neutral mass eigenstates, the lightest of which is stable and can play the role of DM. In the following, we focus on computing the long-range effects in the unbroken electroweak phase. In the companion paper [28], we briefly review the mass eigenstates and their interactions after electroweak symmetry breaking for the choice of parameters denoted in eqs. (2.2) and (2.4), before computing the DM decoupling in the early universe, which alters the predicted masscoupling relation, thereby affecting all experimental constraints.

Static potentials
The D,D and S fermions interact with each other via the B, W and H boson exchanges that give rise to long-range potentials. The kernels generating these potentials are shown in fig. 1. To compute them, we decompose the incoming and outgoing momenta as where P is the total momentum and ±p (±p ) are the momenta of incoming (outgoing) particles in the center-of-momentum (CM) frame. For low momentum transfers, we find (see e.g. [1,[18][19][20]49]) In determining the sign of each contribution in the above, we have taken into account the number of fermion permutations needed to perform the Wick contractions. This is the origin of the relative minus sign between the t-and u-channels of the DD ↔ SS and DD ↔ DD interactions. The factors 1/2 appearing in eqs. (2.7a) and (2.7c) ensure that the resummation of the kernels does not double-count the loop diagrams by exchanging identical particles in the loops; this is shown in appendix A.
To obtain the non-relativistic potentials, we must diagonalise the interactions (2.7) in momentum, spin, and gauge space.
Momentum space. The kernels of both the t-and u-channel diagrams depend only on the momentum transfer, which is however different in the two cases, K t (p − p ) and K u (p + p ). This implies that in position space, the u-channel potential depends on the orbital angular momentum mode of the state under consideration [1, appendix A]. Specifically, the static potentials generated by t-and u-channel diagrams are [1,18] Inserting eqs. (2.7) into (2.8) yields the well-known Coulomb and Yukawa potentials.
Spin diagonalisation. The factors δ s 1 s 1 δ s 2 s 2 and δ s 1 s 2 δ s 2 s 1 arising from t-and u-channel exchanges respectively, can be written in matrix form in the basis {↑↑, ↑↓, ↓↑, ↓↓} as t-channel: Clearly, the t-channel interactions conserve spin along each leg of the ladder, and the corresponding spin factor is simply the unity operator. On the other hand, the u-channel interactions conserve the total spin only. Indeed, the u-channel spin eigenvalues are {−1, 1, 1, 1}, and correspond to the eigenvectors of total spin In the following, we shall therefore project the asymptotic states of pairs of S, D,D fermions onto eigenstates of total spin.
Gauge diagonalisation. Since the interactions respect the SU L (2) symmetry, this amounts to projecting on SU L (2) representations of the incoming or outgoing pairs. For two multiplets transforming under the representations R 1 and R 2 of a gauge group, the Coulomb potential generated by the gauge-boson exchange in the configuration R ⊂ R 1 ⊗ R 2 of the pair is V R (r) = −α R /r with [50] where α is the fine structure constant of the group, and C 2 (R) is the quadratic Casimir operator of the representation R. For SU (2), 2 ⊗ 2 = 1 + 3, and C 2 (2) = 3/4, C 2 (3) = 2. Thus, the DD, DD andDD pairs appear in SU L (2) singlet and triplet configurations, with α 1 2 = 3α 2 /4 and α 3 2 = −α 2 /4 respectively. For the identical-particle pairs DD andDD, the (anti)symmetry of the SU L (2) (singlet) triplet states in gauge space generates also the factor (−1) I+1 for the u-channel diagrams, where here I = 0, 1 stands for the Weak isospin. Note that the U Y (1) potentials are of course not affected by the SU L (2) diagonalisation. 1 The remaining task is the H-mediated interaction of eq. (2.7a). This occurs only in the SU L (2) singlet state. Projecting DD on the singlet, we obtain Kernel (anti)symmetrisation. We close this discussion by noting the importance of considering properly both the t-and u-channel diagrams when identical particles are present in the initial and/or final states, as e.g. in the DD ↔ SS and DD ↔ DD interactions of the present model. Using the correct kernel, as derived in appendix A, and taking into account the above yields the overall factor that enforces the proper particle statistics. This result is valid for pairs of fermions as well as pairs of scalars.
Collecting the above, in table 2 we summarise the potentials generated by the oneboson-exchange diagrams of fig. 1.
Sign of the potential  fig. 1. and s denote the orbital angular momentum mode and the total spin respectively.

Asymptotic scattering and bound states
The potentials of table 2 determine the asymptotic states. For all gauge assignments except the singlet states, SU L (2) × U Y (1) = (1, 0), finding the corresponding wavefunctions is rather straightforward; it amounts to solving a single Schrödinger equation with the corresponding potential, and antisymmetrising the wavefunction in the case of identical particles. For the gauge-singlet states, the DD ↔ SS interaction implies that a system of coupled Schrödinger equations must be solved. We work out this case in detail in section 2.3.2, after we discuss the general properties of the wavefunctions and the underlying hierarchy of scales in section 2.3.1. All the results on the wavefunctions of the scattering and bound states are summarised in tables 3 and 4.

Wavefunctions and hierarchy of scales
The non-relativistic potentials of In the approximation (2.14), we can express all wavefunctions in terms of those for a Coulomb potential V (r) = −α/r, which we shall denote as ϕ(r; α) and we now summarise. For clarity, we denote by α S and α B the couplings of the scattering and bound states. The momentum of each particle in the CM frame in the scattering states is k ≡ µv rel , with v rel being the relative velocity. The Bohr momenta in the scattering and bound states are κ S ≡ µα S and κ B ≡ µα B . For convenience, we define the parameter ζ S ≡ κ S /k = α S /v rel (in section 3 we will also use ζ B ≡ κ B /k = α B /v rel ), and the variables x S ≡ kr and x B ≡ κ B r. The energy eigenvalues of the scattering and bound states are The scattering state wavefunctions can be decomposed in partial waves, and is the well-known s-wave Sommerfeld factor. The bound-state wavefunctions are where the normalisation of the spherical harmonics is dΩ Y m (Ω)Y * m (Ω) = δ δ mm . The emergence of the Sommerfeld effect and the existence of bound states emanate from the different scales involved in the interactions of the D,D and S particles, µv 2 rel /2 µv rel µ < m and µα 2 /(2n 2 ) µ|α|/n µ < m, (2.20) where here α = α S or α B . (Note that α S may be negative.) In computing the BSF cross-sections in sections 3 and 4, we make approximations based on these hierarchies. The hierarchies (2.20) imply that the couplings (2.5) must be evaluated at the appropriate momentum transfer in every occurrence. The average momentum transfers are (cf. section 3 for the last two) annihilation vertices: m scattering-state potentials: µv rel bound-state potentials: µα B emission vertices for BSF: (µ/2)(α 2 B /n 2 + v 2 rel ) emission vertices for bound-to-bound transitions: (µ/2)|α 2 B /n 2 − α B 2 /n 2 | Here for simplicity we will neglect the running of the couplings, although it is easy to restore the scale dependence in all of our analytical results.
In computing the DM freeze-out in [28], we adopt the values of the gauge couplings α1 and α2 at the Z pole, α1(mZ) 0.00973 and α2(mZ) 0.0339. Since αH increases with Q, we consider the quoted values of αH to correspond to the highest relevant scale, Q = m, such that the theory remains well-defined at all Q m. For large values of αH, a Landau pole appears at fairly low energies (but larger than m); however this may be cured by other new physics around those scales. For the renormalisation group equations in the present model, we refer to [38]. where the component wavefunctions are Here, S j (1,0) denotes the gauge-singlet states. Along with their wavefunctions Φ j , they carry quantum numbers that define their energy, angular momentum and spin, and that we have here left implicit. The superscript j aims to differentiate between states with the same quantum numbers but different boundary conditions. Ω here stands for the vacuum of the interacting theory, T is the time ordering operator, S, D andD are the field operators.
The resummation sketched in fig. 2 implies that Φ obey the Schrödinger equations where V (1,0) is the 2 × 2 potential matrix of the gauge-singlet states In the limit m H → 0, the system (2.23) easily decouples. We first define the couplings as well as the unitary matrix whose columns are the normalised V (1,0) eigenvectors with eigenvalues −α R and α A . We note that α R , α A 0. Then, the rotated wavefunctionsΦ ≡ P † Φ obey the Schrödinger equations with the potentialV We now seek scattering and bound state solutions to eq. (2.27).
Scattering states. The scattering state solutions of eq. (2.27) are given by eq. (2.16), albeit the normalisation of each component is allowed to vary and will be determined by the boundary conditions on Φ j . Analysing eq. (2.27) in partial waves, the solutions arê , (2.29) where at x S → ∞, the wavefunctions ϕ k, (x S ; α) obey We seek scattering-state solutions to eq. (2.23) that asymptote at r → ∞ to a pure SS or DD state. In terms of partial waves, this implies that at r → ∞, where i = SS, DD denotes the component. In eq. (2.31a), we included the antisymmetrisation factor √ 2δ +s,even due to the identical particles in the SS-like state, with s = 0 or 1 being the total spin (cf. appendix A.2).
SinceΦ |k|, (x S ) = P † Φ |k|, (x S ), the asymptotic behaviours (2.30) and (2.31) imply We recall that α R and α A depend on +s, and note that while the DD-like state should not be antisymmetrised, the SS components of the symmetric +s = odd modes do vanish due to the antisymmetrisation of the kernel (α R = 0 for + s = odd). Combining eqs. (2.29) and (2.32), we obtain the wavefunctions Φ j |k|, (r) = PΦ j |k|, (r). The results are summarised in table 3.
Bound states. It is easy to see that eq. (2.27) has only one set of bound-state solutions, Φ n m (r) = (0, ϕ n m (r; α A )) T , with E n = −µα 2 A /(2n 2 ). The unrotated wavefunctions are Φ n m (r) = PΦ n m (r). The result is summarised in table 4.

Validity of the Coulomb approximation m H → 0
Scattering states. The scattering states are listed in table 3. The DS states with + s = even are subject to an attractive Higgs-mediated potential, and the Coulomb approximation is good as long as [19] µv rel m H . (2.33) The condition becomes somewhat stronger for the DS scattering states with + s = odd, where the Higgs-mediated potential is repulsive. Moreover, it is relaxed or strengthened in the presence of an additional attractive or repulsive Coulomb potential due to B or W exchange, as is the case with the SS-like and DD-like scattering states for + s = even [ fig. 13].) The mixed SS/DD states are bound by the combined attraction of the B, W and H bosons, which removes the condition of existence and relaxes the condition for the Coulomb approximation [3, fig. 6]. Thus, the strongest condition for the Coulomb approximation to be satisfactory is We note that this condition is essentially satisfied everywhere BSF via Higgs emission is phenomenologically significant (cf. section 3.) Indeed, the energy available to be dissipated must exceed the Higgs mass, (µ/2)[(α H /n) 2 + v 2 rel ] > m H , while BSF is most significant when v rel α H /n. Since α H < 1, this yields a stronger condition.
Bound states can also form via B or W emission, in which case there is no phase-space suppression that ensures the validity of the Coulomb approximation. However, in [28] we shall see that these processes are less significant for the DM thermal decoupling in the present model.
We further discuss the Coulomb approximation for the DM freeze-out in ref. [28].

Annihilation
The S, D andD fermions annihilate into SM particles via various processes that we list in table 5 together with their tree-level cross-sections and Sommerfeld factors. We consider s-wave contributions only. Because the non-relativistic potentials between the annihilating particles depend in many cases on their spin and gauge representations, we project the initial states on eigenstates of total spin and Weak isospin. With the help of table 3, it is straightforward to obtain the Sommerfeld factors for all states except the spin-0 (1, 0) ones, which involve mixing of the SS and DD channels and we discuss in detail below. The Sommerfeld factors are expressed in term of the S 0 (ζ) function defined in eq. (2.18), and the variables where the couplings α 1 , α 2 , α H , α A and α R have been defined in eqs. (2.5) and (2.25). In fig. 3, we present the total s-wave 2-to-2 annihilation cross-section, averaged over the dof of the incoming particles.

Mixed (1, 0) spin-0 states
The annihilation amplitudes of the SS-like and DD-like states (denoted by M) are related to the perturbative amplitudes (denoted by A) as follows where f stands for the final state, and [φ k (k )] j i are the momentum-space wavefuntions, with the j and i indices denoting the state and the component respectively, as in section 2.3. Since the perturbative s-wave SS annihilation vanishes, in our approximation there is no interference between SS and DD channels. Taking into account that the perturbative s-wave annihilation amplitudes are independent of the momentum to lowest order in v rel , we obtain the standard result, where [φ k (r)] j i are the position-space wavefunctions computed in section 2.3. Using the results quoted in table 3, we find All channels with Higgs potential without Higgs potential Figure 3. The s-wave annihilation cross-sections, by initial state (left) and total (right), averaged over the dof of the incoming particles, with and without the Higgs-mediated potential (cf . table 5). Because the processes affected by the latter have either low multiplicity or small perturbative crosssections, the Sommerfeld effect at low velocities arises mostly due to the B µ and W µ gauge bosons.
For α H 0.01, the perturbative annihilation cross-section is also not substantially affected by the coupling to the Higgs (right panel, blue line). Note that we have weighted the contribution of each annihilation channel with the number of DM particles eliminated in each process as estimated upon thermal averaging (cf. ref. [28].) Note that eq. (2.38a) includes the symmetry factor of the spin-0 SS-like state, and that [φ k (0)] SS DD 2 spin-1 = 0. We are ultimately interested in the reduction of the DM density via the various annihilation processes. For the spin-0 (1, 0) states, the rate is where n and · denote densities and thermal averages, and the factor 2 in the second term appears because D andD are not identical (cf. ref. [28].) In the limit m S = m D , the densities are (n S n S ) spin-0 (1,0) = (n D nD) spin-0 (1,0) , thus the DM density destruction rate can be computed by regarding that the spin-0 (1, 0) DD perturbative annihilation cross-sections are enhanced by the factor We quote this result in table 5, but emphasise that (n S n S ) spin-0 (1,0) and (n D nD) spin-0 (1,0) depend exponentially on the corresponding masses, thus eq. (2.40) ceases to be a good approximation already for fairly small mass differences |m D − m S |. Table 5. Annihilation processes, their tree-level s-wave velocity-weighted cross-sections σ 0 v rel and Sommerfeld factors. All σ 0 v rel are averaged over the degrees of freedom of the corresponding projected scattering state (5th column). For DD, σ 0 v rel includes the symmetry factor due to the identical initial-state particles. For the gauge-singlet spin-0 DD channels, see text for discussion. S 0 (ζ) and the various ζ parameters are defined in eqs. (2.18) and (2.35).

Ground-level bound states and their decay rates
Besides annihilating directly into radiation, the S, D andD fermions can form unstable bound states that decay rapidly into radiation, thereby enhancing the DM destruction rate. However, the DM annihilation via BSF is impeded by the inverse (ionisation) processes. The efficiency of BSF in reducing the DM density depends on the interplay of bound-state ionisation and decay processes and bound-to-bound transitions [9]. The ionisation processes become inefficient as the temperature drops around or below the binding energy. This occurs earlier for the most deeply bound states, which also have the largest decay rates. Thus, the tightest bound states have the greatest effect on the DM density, and for this reasons we shall consider the ground level of each bound state species only, {n m} = {100}, which has the largest binding energy. 3 We list the ground-level bound states in table 6.
The BSF cross-sections and bound-to-bound transition rate are computed in sections 3 and 4. Here we discuss the bound-state decay into radiation and the relevance of boundto-bound transitions.

Decay into radiation
The decay rate of bound states with zero angular momentum into radiation is where X 1 X 2 represent the constituent fields of the bound state and f stands for the final state particles. (σ 0 v rel ) X 1 X 2 →f is the s-wave annihilation cross-section of an X 1 X 2 scattering state with the same quantum numbers (spin, gauge and global) as the bound state, averaged over the dof that correspond to those quantum numbers only (rather than all the dof of an X 1 X 2 scattering state.) For an attractive Coulomb potential of strength α, the squared ground-state wavefunction evaluated at the origin is |ψ n00 (0)| 2 = κ 3 B /π = µ 3 α 3 /π, where κ B = µα is the Bohr momentum.
Taking into account the s-wave annihilation processes of table 5, we compute the total decay rates of the ground states and list them in table 6. The decay of the mixed SS/DD bound state occurs via its DD component, and the rate is computed analogously to the annihilation of the mixed scattering states described in section 2.4. For the DD bound state, a factor 2 due to the antisymmetrisation of the wavefunction has already been included in the corresponding σ 0 v rel in table 5 and should not be included again when computing the decay rate of the bound state.
As seen in table 5, the DS bound state cannot decay into two particles. It may decay instead into three bosons, however the corresponding rates are suppressed by higher powers of the couplings, O(α 2 , as well as the three-body final-state phase Negligible. Transition to B(SS/DD): eq. (3.26), eqs. (4.12) and (4.19) to (4.21) Table 6. The ground-level bound states, {n m} = {100}, and their rates of decay into radiation, in the limit m H → 0. The decay rate of the DS bound state is suppressed, and we reference instead the formulae for its transition rate to the SS/DD bound state. In the first row, α A and α R are found from eqs. (2.25) for = s = 0. The binding energy of each bound state is space. The DS bound state instead decays much faster into the tighter SS/DD bound state, as we shall now see.

Transitions into deeper bound levels
Besides decaying directly into radiation, bound states may transition into lower-lying bound levels via dissipation of energy. The bound-to-bound transition rates are computed in a similar fashion to BSF processes, as we shall see in sections 3 and 4. Here we note that in radiative transitions (either scattering-to-bound or bound-to-bound), spin is conserved at leading order in the non-relativistic regime. Spin-flipping transitions can occur via emission of a vector boson, but rely on spin-orbit interaction and are suppressed by higher powers of the couplings. Consequently, they are subdominant to the direct bound-state decay rates into two relativistic species.
Considering the bound states of table 6, there are only two spin-conserving transitions, of which only one may occur with emission of a single boson contained in the theory. Noting that α A α H (cf. eq. (2.25b)), this transition is In fact, much like BSF, bound-to-bound transitions may occur either radiatively or via scattering on the relativistic thermal bath. In sections 3 and 4, we compute the corresponding rates for the transition (2.42) and reference the results in table 6. Figure 4. Radiative bound-state formation and bound-to-bound transitions.

Radiative bound-state formation and bound-to-bound transitions
We now compute the cross-sections for the radiative formation of the ground-level bound states of table 6. We first outline the elements of the computation, explain the approximations involved, and define some useful quantities. Then we proceed with the computation of the amplitudes and cross-sections. The final results are listed in tables 7 to 10 and illustrated in fig. 9.

Preliminaries
Full amplitudes M k→n m and M n m →n m . As depicted schematically in fig. 4, the full amplitudes consist of the radiative transition part A T computed perturbatively and convoluted with the initial and final-state wavefunctions [18], Equations (3.1) can accommodate the possibility that the incoming and/or outgoing states are superpositions of different Fock states, as is the case with the mixed SS/DD states discussed in section 2.3. Then, the wavefunctions are vectors and A T becomes a matrix.
Transition amplitudes A T . The radiative parts of the BSF and transition diagrams are shown in figs. 5 to 8. Following refs. [1,18,20], we compute the amplitudes A T , applying the standard approximations due to the hierarchy of scales and retaining only the leading order terms. Among else, we shall use the following approximate spinor identities valid for low momentum changes, p p , Moreover, we emphasise that the signs arising from the fermion permutations needed to perform the Wick contractions must be carefully taken into account, as they often differ among various diagrams contributing to the same amplitude. An example calculation is presented in detail in appendix B.
Overlap integrals. The scattering and bound state wavefunctions are listed in tables 3 and 4. To express M k→n m and M n m →n m compactly, we define the Coulombic overlap integrals [1,[18][19][20] and entail the coupling strengths α S and α B of the potential in the scattering and bound states respectively, and we have defined the functions [1,9,[18][19][20] Note that S scl and S vec include the s-and p-wave Sommerfeld factors, S 0 (ζ S ) ≡ 2πζ S /(1 + e −2πζ S ) and S 1 (ζ S ) = S 0 (ζ S )(1 + ζ 2 S ), respectively. Indeed, eq. (3.4a) arises from the S = 0 and eqs. (3.4b) to (3.4d) arise from the S = 1 modes of the scattering state wavefunctions. In eq. (3.4d) we took the limit m H → 0 to be consistent with our approximations, although it is easy to obtain an analytical result for m H = 0 (but using the Coulomb wavefunctions.) BSF cross-sections. The cross-sections for BSF via emission of a massless vector (B or W ) or scalar (H or H † ) boson are, respectively [18,19] where P V and P H are the momenta of the emitted bosons, which dissipate the kinetic energy of the relative motion and the binding energy, In eq. (3.7a), we have summed over polarisations of the emitted vector. As we shall see in the following, to working order, the amplitudes for BSF via vector emission are M k→100 ∝ k, while the amplitudes for BSF via scalar emission are independent of the solid angle Ω. Then, eqs. (3.7) simplify to [18,19] where h H ≡ |P H |/E H is the phase-space suppression due to the Higgs mass, Bound-to-bound transition rates. Similarly to the above, the rates for the radiative de-excitation of bound states are dΓ V n m →n m where |P V | and |P H | are determined again by the amount of dissipated energy, which is now the difference between the binding energies of the two states, For monopole transitions via H or H † emission, the amplitude is independent of the P H direction, and eq. (3.11b) simplifies to where the phase-space suppression h H is defined in eq. (3.10). The BSF processes are listed in table 7, and the radiative part of the diagrams contributing to these processes are shown in fig. 5. We project the bound-state fields on the spin-0 state via [U −1 spin-0 ] r 2 r 1 = r 1 r 2 / √ 2, and the DD component on the SU L (2) singlet via δ i j / √ 2.

SS/DD
The perturbative parts of the amplitude are In eq. (3.14a), the fermion permutations introduced factors (−1) and (+1) for the t-and u-channel diagrams. The projection on the antisymmetric spin-0 eigenstate alloted another factor (−1) to the u-channel. The full amplitude (3.1a) is where only the S = 1 component of the scattering state wavefunction is meant to be kept, and here α A and α R should be evaluated from eqs. (2.25) for = s = 0. This becomes Note that we have neglected the ±P B /2 terms inside the δ-functions of eq. (3.14b) that give rise to higher order corrections [18,19]. Squaring and summing over the initial-state gauge indices and spins selects the (1, 0) spin-0 DD scattering state, which has one dof. Using the overlap integrals (3.4), we find The cross-section is obtained from eq. (3.9a) setting α B → α A , and is shown in table 7.

DD → B(SS/DD) + W
The perturbative parts of the amplitude are where the signs of the t-and u-channel diagrams in eq. (3.18a) are as in DD → SS + B above, and the first factor 2 in the first term of eq. (3.18b) is the quadratic Casimir of SU L (2) (see ref. [20] for details of this computation.) The full amplitude (3.1a) is where again α A and α R should be evaluated from eqs. (2.25) for = s = 0, and we have neglected the ±P W /2 terms inside the δ-functions in eq. (3.18b). Squaring, summing over the initial and final state gauge indices and spins selects the (3, 0) spin-0 DD state, which has three dof. Using eqs. (3.4) and (3.6), we find The cross-section is obtained from eq. (3.9a) setting α B → α A , and is shown in table 7.

DS → B(SS/DD) + H
The perturbative parts of the amplitude are where now the fermion permutations introduced factors (+1) and (−1) for the t-and uchannel DS → SS + H diagrams respectively, and (−1) for the DS → DD + H diagram. We note that there are two diagrams where an off-shell vector boson (B or W ) emitted from one leg and an off-shell Higgs emitted from the other leg fuse to produce the final-state Higgs. In appendix D, we show that these diagrams are of higher order, thus we do not consider them here. The full amplitude (3.1a) is where we have neglected the ±P H /2 terms inside the δ-functions in eqs. (3.21). Squaring and summing over the initial and final state gauge indices and spins selects the (2, 1/2) spin-0 DS state, which has two dof. Using the overlap integrals (3.4), we find (3.23) The cross-section is obtained from eq. (3.9b) setting α B → α A , and is shown in table 7.

Bound-to-bound transition B(DS) → B(SS/DD) + H
Using the perturbative amplitudes (3.21), we may now compute the rate of bound-to-bound transition (2.42). Projecting on the spin-0 DS state, and taking into account the boundstate wavefuntions of table 4, we find that the full amplitude is, analogously to eq. (3.22), given by Squaring and averaging over the initial and final state gauge indices, and using the overlap integral (3.4e), we find (3.25) The transition rate is found from eq. (3.13) to be    The BSF processes are listed in table 8, and the radiative parts of the diagrams contributing to these processes are shown in fig. 6. For all processes below, we project the bound-state fields on the SU L (2) singlet via δ i j / √ 2. Moreover, since spin is conserved at working order, as already seen in section 3.2, we project both the scattering and the bound states on the spin-1 configuration by contracting the spin indices with [U σ spin-1 ] s 1 s 2 [U ρ † spin-1 ] r 2 r 1 , where the indices σ, ρ = −1, 0, 1 run through the three states of the spin-1 multiplets. While the Clebsh-Gordan coefficients contained in U spin-1 are well-known, we shall not need them explicitly. We will instead only invoke that the operators [U σ † spin-1 ] s 1 s 2 are symmetric in s 1 , Besides the projections mentioned above, here we also project the DD component of the scattering state on the SU L (2) singlet configuration via δ ij / √ 2. The perturbative parts of the amplitude are In eq. (3.27a), the fermion permutations introduced factors (−1) and (+1) for the t-and u-channel diagrams. Upon projection on the symmetric spin-1 eigenstate, their relative sign does not change. The factor √ 2 appearing in the beginning of eq. (3.27a) arises from the projection onto the SU L (2) singlet final DD state. Using the wavefunctions listed in tables 3 and 4, we find the full amplitude (3.1a) to be where here α A and α R should be evaluated from eqs. (2.25) for S = s = 1 (scattering state). As before, we have neglected the ±P B /2 terms inside the δ-functions in eq. (3.27b). Next, we square, sum over the spins, and average over the three dof of the incoming spin-1 state. Using the overlap integrals eq. (3.4), we find where here α B = (α 1 + 3α 2 )/4 and correspondingly ζ B = (ζ 1 + 3ζ 2 )/4. The cross-section is obtained from eq. (3.9a), and is shown in table 8.

DD → B(DD) + W
The perturbative part of the amplitude is (cf. eq. (3.18b)) Using the wavefunctions listed in tables 3 and 4, and anticipating that the scattering state will be an SU L (2) triplet, we find the full amplitude (3.1a), Squaring and summing over gauge and spin indices, projects the scattering state on the spin-1 SU L (2) triplet that has 9 dof. Using the overlap integrals (3.4), we find The cross-section is obtained from eq. (3.9a), and is shown in table 8. It agrees with the results of refs. [3,20] appropriately adjusted.

DS → B(DD) + H
The perturbative part of the amplitude is Using the wavefunctions listed in tables 3 and 4, we find the full amplitude (3.1a), and taking into account the overlap integral (3.4a),

DD bound states: (1, 1), spin 1, n m = {100}
The BSF processes are listed in table 9, and the radiative part of the diagrams contributing to these processes are shown in fig. 7. In all processes below, we project the bound-state fields on the SU L (2) singlet via i j / √ 2. Moreover, as in section 3.3, we project both the scattering and the bound states on the spin-1 configuration via [U σ spin-1 ] s 1 s 2 [U ρ † spin-1 ] r 2 r 1 , and invoke that [U σ † spin-1 ] s 1 s 2 are symmetric in s 1 , s 2 , and [U σ spin-1 ] s 1 s 2 [U ρ † spin-1 ] s 1 s 2 = δ σρ .

DD → B(DD) + W
The perturbative part of the amplitude is where the last line accounts for the u-channel diagrams. The different number of fermion permutations in the t-and u-channel diagrams introduces a relative (−1) factor between the two, while the projection on the symmetric spin-1 state does not. Here we have neglected the ±P W /2 terms inside the δ-functions already at the level of the perturbative amplitude. It is easy to check that the gauge factors in eq. (3.38) are symmetric in i ↔ j, as expected, since the scattering state must be an SU L (2) triplet. Convoluting with the scattering state wavefunction and setting k → −k for the u channel renders the latter identical to the t-channel up to the extra factor −(−1) S = +1 since S = 1. Thus, the t and u channels add up, and we find the full amplitude (3.1a) to be where we also took into account the symmetry factors of the scattering and bound state wavefunctions, as stated in tables 3 and 4. Considering the relation (3.4c) between the overlap integrals, the above simplifies to where here α B = (−α 1 + 3α 2 )/4, and G a ij is the gauge factor (3.42) Squaring eq. (3.40) and summing over gauge indices and spins projects the scattering state on the spin-1 SU L (2) triplet configuration that has nine dof. Considering the overlap integral (3.4b), we find The cross-section is obtained from eq. (3.9a) and is shown in table 9.

DS → B(DD) + H †
The perturbative part of the amplitude is where the fermion permutations alloted factors (+1) and (−1) to the t and u channels respectively. As seen from eq. Here α B = (−α 1 + 3α 2 )/4, and we have included the symmetry factors of the scattering and bound state wavefunctions. Squaring and summing over the gauge indices and spins, and using the overlap integral eq. (3.4a), we find where we averaged over the six dof of the spin-1 SU L (2) doublet scattering state. The cross-section is obtained from eq. (3.9b) and is shown in table 9.

DS → B(DS) + B
The perturbative part of the amplitude is where the fermion permutations alloted factors (+1) and (−1) to the t and u channels. Upon projection on the spin-0 states, the u channel acquired another factor (−1). Using the wavefunctions of tables 3 and 4, we find the full amplitude (3.1a) to be  Next we square, sum over the gauge indices and average over the two dof of the spin-0 scattering state. Using the overlap integrals (3.4), we obtain The cross-section is found from eq. (3.9a) for α B = α H and is shown in table 10.

DS → B(DS) + W
The perturbative part of the amplitude is where we used Tr(t a t a ) = 3/2. The cross-section is found from eq. (3.9a) for α B = α H and is shown in table 10.

SS-like → B(DS) + H †
The perturbative parts of the amplitude are In eq. (3.53a), the fermion permutations alloted signs (+1) and (−1) factors to the t and u channels. Upon projection on the spin-0 state, the u channel acquired another factor (−1). We now project the DD component of the scattering state in eq. (3.53b), on the SU L (2) singlet configuration via δ ij / √ 2, Considering the wavefunctions of tables 3 and 4, the full amplitude (3.1a) is where α A and α R should be evaluated from eq. (2.25) for S = s = 0. Using the overlap integrals (3.4), The cross-section is found from eq. (3.9b) for α B = α H and is shown in table 10.

DD-like → B(DS) + H †
Starting from the perturbative parts (3.53a) and (3.53c), and considering the wavefunctions of tables 3 and 4, the full amplitude (3.1a) is where again α A and α R should be evaluated from eq. (2.25) for S = s = 0. Using the overlap integrals (3.4), The cross-section is found from eq. (3.9b) for α B = α H and is shown in table 10.

DD → B(DS) + H †
The perturbative part of the amplitude is given in eq. (3.53b). We project it the DD scattering state on the SU L (2) triplet configuration via t a ji / C(2) = √ 2t a ji , where C(2) = 1/2 is the Casimir of the SU L (2) doublet representation, and obtain Considering the wavefunctions of tables 3 and 4, the full amplitude (3.1a) is Squaring, summing over the final state gauge indices, and averaging over the three dof of the scattering state, we obtain where we used Tr( √ 2t a √ 2t a ) = 3. The cross-section is found from eq. (3.9b) for α B = α H and is shown in table 10.

DD → B(DS) + H
The perturbative part of the amplitude is where the fermion permutations alloted signs (+1) and (−1) factors to the t and u channels. Upon projection on the spin-0 state, the u channel acquired another factor (−1). We now project the DD scattering state on the SU L (2) triplet configuration via the symbolic

Unitarity and BSF via Higgs emission
The unitarity of the S matrix implies an upper limit on the partial-wave inelastic crosssections, σ inel σ uni = (2 + 1)π/k 2 , where is the partial wave and k is the momentum of either of the interacting particles in the CM frame. In the non-relativistic regime, k = µv rel with µ being the reduced mass, thus As already discussed in ref. [1], the high efficiency of BSF via charged scalar emission (here, the Higgs doublet) implies that the unitarity limit (3.65) may be saturated already for rather small values of α H . If the incoming particles interact via an attractive long-range force, then this occurs for the continuum of velocities v rel α B , otherwise only for a finite range or discrete values of v rel (cf. fig. 11.) The apparent violation of unitarity at larger α H by the computations of sections 3.2 to 3.5, whether it occurs for an infinite or finite range of v rel , indicates that these computations must be amended. At small values of α H , higher order corrections to the perturbative transition amplitudes A T are expected to be insignificant. Restoring unitarity necessitates instead that the two-particle interactions at infinity are resummed [15].
The results of sections 3.2 to 3.5 already include the resummation of the (leading-order) long-range interaction between the incoming particles, computed in section 2.2. However, according to the optical theorem, all elastic and inelastic processes to which the incoming state may participate contribute to its self-energy. Typically, contact-type interactions can be neglected, as they do not distort significantly the wavefuctions of the interacting particles. Nevertheless, if a contact-type interaction is very strong -as is the case when the corresponding cross-section approaches (or even appears to exceed) the unitarity limit (3.65) -then its contribution to the two-particle self-energy may be significant.
The contributions to the 2PI kernels arising from inelastic processes that involve Higgs emission are shown in fig. 10. These diagrams include both scattering and bound intermediate states of the S, D andD particles, and therefore include bremsstrahlung, BSF and bound-to-bound transitions. Note that in fig. 10 we have not included the resummation of the long-range kernels of section 2.2 in the incoming and outgoing pairs, as this would result in double-counting; only 2PI diagrams must be included in the kernels that determine the potential.
The proper resummation of the diagrams of fig. 10 requires developing suitable formalism, and is beyond the scope of the present work. To ensure that our cross-sections are consistent with partial-wave unitarity, we shall instead adapt the result of ref. [51] that resummed the box diagrams arising from the perturbative part of s-wave annihilation into radiation (hard scattering), to compute the effect on the scattering-state wavefunctions and ultimately on the full cross-sections for s-wave annihilation into radiation. This procedure iK ⊃ iA T G (4) H iA T Figure 10. The contributions to the 2PI kernels arising from inelastic processes that involve Higgs emission. The solid lines stand for any of the S, D orD particles, and G (4) includes their scattering and bound states. The dashed line represents the Higgs doublet. A T are the perturbative transition amplitudes with H ( †) emission, computed in sections 3.2 to 3.5 regulates the annihilation cross-sections as follows [51, eq. (40)], where we have generalised the result of ref. [51] to multiple annihilation channels. Equation (3.66) ensures that the unitarity limit is respected by the total s-wave inelastic crosssection, since it implies r reg = r/(1 + r/4) 2 We emphasise that the assumptions made in deriving eq. (3.66) are not strictly satisfied in our case, for at least two reasons: (i) Reference [51] assumed that for the resummed inelastic processes (hard scattering), σv rel is independent of v rel . For BSF via Higgs emission, the corresponding cross-sections can be found from tables 7 to 10 by setting ζ S → 0; eq. (3.6a) then shows that they depend on v rel for v rel α B . (ii) In the present case, the new contributions to the kernel may affect both the initial (scattering) and final (bound) state wavefunctions, while only the former is relevant for annihilation into radiation in the analysis of ref. [51]. Nevertheless, we shall adopt eq. (3.66) as a perscription that regulates the inelastic cross-sections in the velocity range where the base calculation violates unitarity, while leaving them essentially unaffected outside that range. We leave a more precise treatment for future work.
In fig. 11, we show how the prescription (3.66) adjusts the inelastic cross-sections for the scattering states that may participate in BSF via Higgs emission. We include the total inelastic cross-section (BSF plus annihilation) in the resummation, although only BSF is significant. However, both BSF and annihilation are regulated by the same factor, which implies that the annihilation cross-sections are affected significantly even while themselves being well below the unitarity limit.

Bound-state formation and bound-to-bound transitions via scattering
The dissipation of energy necessary for the capture into bound states or transitions between bound levels, may occur via exchange of an off-shell mediator with particles of the thermal bath [52][53][54][55][56][57][58]. 5 References [57,58] showed, in the context of a U (1) model, that the crosssections for BSF via scattering factorise into the radiative ones (with any phase-space suppression due to the mass of the emitted vector removed), and a factor that depends on the thermal bath and the interaction that mediates the scattering. In the following, we consider BSF and bound-to-bound transitions via off-shell Higgs exchange. We derive a similar factorisation and then compute the BSF cross-sections and transition rates. We also adapt the results of refs. [57,58] to our model, for BSF and transitions via off-shell B and W exchange.

Factorisation of the effective BSF cross-sections and transition rates
For simplicity, we lay out the discussion in terms of the BSF cross-sections only. The derivation for bound-to-bound transition rates is analogous.
The thermally averaged rate per unit volume for BSF via off-shell Higgs exchange is 5 Other rearragement processes have been considered in ref. [59].
where we defined Here, the indices 1,2 correspond to the two incoming DM fields, while q i and q f denote the initial and final bath particle momenta. We consider scattering only on relativistic species, whose density in a thermal environment is large, and set q 0 i = |q i | and q 0 f = |q f |. f ± (E) = (e E/T ±1) −1 are the phase-space occupation numbers for fermions (+) and bosons (−), and in eq. (4.2) the upper and lower signs correspond to scattering on fermionic or bosonic dof respectively. ω k→n m is the energy dissipated by the DM fields in the capture process, given by eq. (3.8). When the DM particles are non-relativistic, it depends only on the relative momentum of the incoming DM particles and the binding energy of the bound state, but is independent of the momenta of the bath particles. (We note that this does not hold for the 3-momentum exchange |q f − q i | along the off-shell mediator.) The factor [1 + f − (ω k→n m )] in eq. (4.1) is compensated by the inverse factor in eq. (4.2); this definition ensures that the thermal averaging of the cross-section (4.2) is the same as that of its radiative counterpart, which includes a Bose-enhancement factor for the radiated boson [9] (cf. ref. [28].) Next, we conjecture that the amplitude for off-shell H exchange, M H * -BSF k→n m , can be factorised into the corresponding amplitude with on-shell H emission, M H-BSF k→n m , and a function of the momenta of the bath particles, as follows 6 where the sum on the left side runs over the dof (spin and gauge) of the initial and final bath particles. M H-BSF k→n m may depend only on the momentum exchange q ≡ q f − q i rather than on q i and q f separately. R 0 must be Lorentz invariant since the amplitudes are. It may thus depend only on the 4-vector products q 2 i = q 2 f = 0 and q i · q f ; in eq. (4.3) we have denoted its dependence on the latter. The R 0 factors will be specified in section 4.1.2 for the processes shown in fig. 12. Switching the integration from q f to q, eq. (4.2) gives where The second line of eq. (4.4) would form the BSF cross-section via on-shell emission, except for two complications: (i) The dispersion relation of the radiated momentum is given by eq. (4.5) rather than the on-shell condition of the radiated boson, and depends on the variable q i . (ii) The last line of eq. (4.4) depends on q. These complications are resolved within the non-relativistic approximation, where the cross-section of BSF via scattering can be shown to be proportional to that of radiative BSF, as we shall now see.

Non-relativistic approximation
Similarly to radiative BSF, if the incoming particles are non-relativistic and the final state particles are weakly bound, we may neglect the recoil of the bound state. Then, the energymomentum conservation implies q 0 ω with ω given by eqs.
with τ ≡q i ·q. It is also easy to show that q i · q f = (q 2 − ω 2 )/2. Using the δ 4 -function to perform the integration over d 3 P d|q| (as is standard in the computation of 2-to-2 crosssections), eq. (4.4) gives with P 0 2m and |q| given by eq. (4.6). Note that the second line of eq. (4.7) differs from (σ H-BSF k→n m v rel ) by the factor ω + |q i | |q| + |q i |τ ω |q| due to the different dispersion relation of the radiated momentum q, here given by eq. (4.5).
The radiative BSF amplitudes are typically computed by expanding in powers of the radiated momentum q [19]. 7 As seen in section 3, the dominant contribution to the various 7 The expansion is in effect in the dimensionless combination |q|/ κ 2 /n 2 + k 2 = |q|/ √ 2µω. For BSF via on-shell emission, the radiated momentum is limited by the available energy, |q| ω, thus the expansion parameter is always ω/(2µ) 1. However, for BSF via scattering, |q| can be comparable to or larger than √ 2µω, particularly at T κ/n, which puts in question the validity of the expansion. Nevertheless, for the purposes of DM freeze-out, BSF typically reaches ionisation equilibrium at high T , where the DM destruction rate via BSF is independent of the BSF cross-sections [60] (cf. ref. [28].) At lower T , where the magnitude of the BSF cross-sections matters, the |q| expansion is a valid approximation.

M H-BSF
k→n m amplitudes arises from the zeroth order term [1], i.e. M H-BSF k→n m are independent of q (and therefore τ ) at leading order. Reshuffling the various factors, eq. (4.7) becomes The integration over d 3 q i in the second line of eq. (4.8) eliminates any dependence of the integrand on the orientation of the vector q, allowing us to identify the first line as the cross-section for on-shell Higgs emission with the phase-space suppression removed (cf. eqs. (3.7b) and (3.9b).) Thus, the effective cross-section for off-shell Higgs exchange can be factorised at leading order as follows 8 where h H (ω) is the phase-space suppression (3.10) of the on-shell emission due to the Higgs mass, with ω k→n m being the dissipated energy (3.8), and we restored the indices for concreteness. The dimensionless factor R H (ω) is (4.10) We recall that |q| is given by eq. (4.6). Since the entire integrand is rotationally invariant (recall that R 0 is Lorentz invariant), we perform the d 3 q i integration by setting q on the z axis. Then the azimuthal angle is parametrised by τ defined above. Changing integration variables from |q i | and τ to u ≡ |q i |/ω and z ≡ q 2 /ω 2 − 1, eq. (4.10) simplifies to We compute R H next. The final result can be found in eqs. (4.19) to (4.21).
Following the same steps, we find that the bound-to-bound transition rate via off-shell Higgs exchange is related to the radiative one via where ω n m →n m is the dissipated energy (3.12). 8 Note that eq. (4.9) holds for the effective cross-section of BSF via off-shell Higgs exchange as defined in eq. (4.2). Since M H-BSF k→n m is presumed to be independent of q and therefore q · k, it has not been necessary to integrate over the angular variables of k in order to obtain a factorised expression for the cross-section. This is in contrast to the case of vector exchange, where the amplitude does not factorise as in eq. (4.3) and in fact depends on qi · k and q f · k. A factorisation similar to eq. (4.9) is obtained only for the cross-section averaged over the k solid angle [57,58] (cf. section 4.2.)

Amplitudes
Similarly to their radiative analogues, the amplitudes for BSF and bound-to-bound transitions via scattering consist of the perturbative transition amplitudes that encode the scattering on the bath particles, convoluted with the initial and final state wavefunctions. Focusing again on BSF (cf. eq. (3.1a)), We now compute A H * -BSF T (k , p) for the scattering processes shown in fig. 12, and deduce from eq. (4.11) the corresponding R factors.

Scattering on fermions
The Higgs couples to the SM fermions via the operators where the a, b superscripts indicate the SU L (2) contractions, while the family indices are suppressed. These couplings give rise to the scattering processes shown in fig. 12 (top). The corresponding BSF perturbative transition amplitudes (projected on the desired spin and gauge, scattering and bound states) are 9 where h, h are the SU L (2) indices of the left-handed SM fermion field and the exchanged Higgs, respectively. The scattering and bound states gauge indices, if any, are left implicit. We use y F to denote collectively the SM Yukawa couplings of eq. (4.14). Inserting eq. (4.15) into (4.13), squaring and summing over the bath particle spin and gauge dof, we arrive at the R 0 factors (cf. eq. (4.3)), where we introduced a factor 2 to account for the partner process controlled by the same coupling, where the initial (final) fermion becomes the final (initial) antifermion.

Scattering on bosons
The perturbative transition amplitude for scattering on gauge bosons ( fig. 12, bottom left), projected on the desired spin and gauge, scattering and bound states, is In eq. (4.15), the sign of the γ5 term and whether the Yukawa coupling should be yF or y * F depend on the exact process we are considering. For scattering on antifermions, the spinors ui, u f become vi, v f . In addition, for a scattering involving an up-type right-handed (anti)quark, δ hh should be replaced by hh . However, all these differences do not affect the R0 factors.
where h, h and a are the SU L (2) indices of the outgoing and exchanged Higgs bosons and the incoming gauge boson respectively. T a and g stand for the generators and the gauge coupling of the gauge group under consideration. Inserting eq. (4.17) into (4.13), squaring and summing over the bath particle polarisations and gauge dof, we find the R 0 factors where, as before, we introduced a factor 2 to account for the partner process where H † is the incoming bath particle ( fig. 12, bottom right). C 2 (R H ) is the quadratic Casimir of the Higgs representation under the gauge group considered; here, C 2 (R H ) = Y 2 H = 1/4 for Hypercharge and C 2 (R H ) = 3/4 for SU L (2).

BSF cross-sections and transition rates
Both eqs. (4.16) and (4.18) depend only on q i · q f , as presumed in eq. (4.3), and in fact in the same fashion. Inserting them into eq. (4.11), and carrying out the integration over z, we find the contributions of scattering on fermions and bosons to BSF, where R ± are dimensionless functions of two parameters, ω/T and m H /ω, (4.20) The factor 1/2 in eq. (4.19a) is due to the SM fermions being chiral. The R H factor that determines the BSF via scattering cross-section (4.9) is Among the SM fermions, the top quark yields the largest contribution as long as it remains relativistic. The R ± factors (4.20) diverge at m H → 0, which during the DM thermal decoupling around the EWPT (cf. ref. [28].) This divergence can be removed by a full next-to-leadingorder calculation, as done in ref. [58] in the context of a U (1) gauge theory. Performing such a computation for the model considered here is beyond the scope of this work. However, comparing the results of refs. [57] and [58] for a massive and massless vector mediator respectively, we find that, upon thermal averaging, the former approximates well the latter at temperatures higher than the binding energy if the screening scale (i.e. the mediator mass) is set to 0.74 ω. 10 Considering this, in eq. (4.20) we shall do the replacement m H → max(m H , ω). (4.22) We present R ± and R H in fig. 13. It is clear that they are more significant for ω/T 1. This implies that for bound-to-bound transitions, they enhance the rates at T ω n →n = |E n − E n |. For BSF via scattering, the R ± factors weigh preferentially the contribution of DM pairs with low relative velocity. We note that even though in a thermal bath ω k→n m = |E n | + (3/2)T > T (cf. eq. (3.8)), lower values of ω k→n m may still incur in a sizeable portion of the DM collisions while T |E n |.
Even when the R H factor (4.21) is less than 1, BSF via scattering may potentially be (i) faster than radiative BSF, which is suppressed by the h H (ω k→n m ) phase-space factor (3.10), becoming entirely inaccessible for m H /ω k→n m > 1, and (ii) significant with respect to direct annihilation and BSF via vector emission, since these processes are suppressed by one (two) extra power(s) of couplings compared to BSF via H ( †) off-shell exchange (on-shell emission.) Analogously, bound-to-bound transitions via scattering may dominate over their radiative counterparts and/or the direct bound-state decay into radiation.
To assess realistically the impact of BSF and bound-to-bound transitions via scattering we must thermally average the cross-sections and rates of eqs. (4.9) and (4.12), and consider the interplay of bound-state formation, decay, ionisation and transition processes in the thermal bath. This is done in ref. [28]. Here we only note that considering BSF via scattering does not increase the DM destruction rate proportionally, since at early times a state of ionisation equilibrium is typically reached where the DM destruction due to BSF is independent of the actual BSF rate provided that the latter is sufficiently large [60].

B and W exchange
References [57,58] showed in the context of a U (1) gauge theory that the effective crosssection for BSF via off-shell vector exchange, defined via the thermally averaged rate per volume (cf. footnote 8) is, at leading order in the non-relativistic regime, proportional to the cross-section for onshell emission, 1) . As in section 4.1, the factor 2 accounts for the partner processes related via exchanging the initial (final) bath particle with the final (initial) bath antiparticle, and α is the fine structure constant of the group. The factor R U (1) depends only on ω/T provided that V is massless, and has been derived in [58] via a next-to-leading order calculation where the colinear and infrared the divergences are cancelled. For a massive V , a simple analytical formula that depends on 10 Reference [58] found that using the binding energy as the minimum screening scale provides a good approximation. While this is indeed so, the above prescription ensures that the R factor depends only on ω/T as predicted by the full computation, besides being a somewhat better approximation. ω/T and m V /ω has been computed in [57]. We define R U (1) to correspond to scattering on one species of relativistic Dirac fermions with charge unity, and use the results of [58]. 11 Adapting the result to the present model, the cross-sections for BSF via off-shell B and W exchange are related to those of on-shell emission as follows The factors c B and c W account for scattering on the relativistic SM fermions. The contribution of a chiral fermion F transforming under the representation R F of a gauge group is c F = C(R F )/2, where C is the Casimir operator. When all the SM fermions are relativistic, c B = (1/2) F Y 2 F = 5 and c W = (1/2)C(2) × 12 = 3, where C(2) = 1/2 for SU L (2). Note that eq. (4.25b) includes only scattering on SM fermions via off-shell W exchange. However, non-Abelian gauge bosons may also scatter on themselves due to the trilinear gauge coupling. Estimating this effect necessitates a dedicated next-to-leading order computation that is beyond the scope of this work. We shall thus neglect this contribution.
Formulae analogous to eqs. (4.25) hold for bound-to-bound transitions via off-shell B and W exchange, however no such transition is of interest here.
We present the R B and R W factors in fig. 13.

Conclusion
Renormalisable scenarios in which DM is the lightest mass eigenstate arising from the mixing of two electroweak multiplets that couple to the Higgs, are among the archetypical WIMP DM models. Here, we have considered the role of the Higgs doublet in the nonperturbative phenomena -the Sommerfeld effect and the formation of bound states -that take place during the thermal decoupling of multi-TeV DM from the primordial plasma. We have shown that the effect of the Higgs doublet is two-fold: (i) it mediates a long-range interaction that affects the wavefunctions of both scattering and bound states, and (ii) its emission precipitates extremely fast monopole transitions, including capture into bound states and transitions between bound levels. In a companion paper [28], we show that the above effects can reduce the relic density of stable species very significantly, thereby altering experimental constraints.
These results build on the work of ref. [1] that showed the importance of bound-state formation via emission of a scalar charged under a symmetry (see also [47]), as well as the work of refs. [2,3] that demonstrated the long-range effect of the Higgs boson between TeV-scale particles. In the present first computation of such effects involving the Higgs doublet, we have focused on the simplest model, comprised by two (nearly) mass degenerate fermionic SU L (2) multiplets, a singlet and a doublet. Our calculations can of course be extended to other models of WIMP DM coupled to the Higgs doublet.
We have considered transitions -both BSF and bound-to-bound transitions -via radiative emission of an on-shell Higgs doublet, as well as via scattering on the thermal bath through off-shell Higgs-doublet exchange. We showed that the rates for the latter factorise into the former times a temperature-dependent function. This parallels the results of refs. [57,58] that considered capture via exchange of an off-shell gauge boson and found a similar factorisation. However, the temperature-dependent function depends on the multiple mode of the transition, and is thus different for the monopole transitions occurring via Higgs-doublet (or generally charged-scalar) exchange and the dipole transitions occuring via gauge-boson exchange.
We finish by summarising the main assumptions and approximations involved in our computations. (a) We have assumed that the coannihilating multiplets have equal masses, in order to obtain analytical results even for the processes that involve mixing of different states (here SS/DD; for processes that do not involve mixing, this is not necessary, as the wavefunctions can be expressed in terms of the reduced mass.) (b) We have neglected the Higgs mass in the Higgs-mediated potential in order again to obtain analytical results; this approximation is discussed in section 2.3.3. (c) The calculations are done in the symmetric electroweak phase. (d) While not considered here, the capture into excited states can be significant (even if not dominant) due to the monopole transitions occurring in this class of models [1]. (e) Monopole transitions can lead to the apparent violation of unitarity even for small couplings; here we have adopted an prescription to treat this issue, however a dedicated study is required, as discussed in section 3.6. Among the above, we believe that (d) and (e), in particular, merit further work in the future. The suitability of these approximations for the computation of the DM relic density is discussed in detail in [28].

A.1 The 2PI kernel
For identical particle (IP) pairs, t-channel diagrams have u-channel counterparts. However, adding them up and resumming them double-counts the loop diagrams because it corresponds to exchanging identical particles in the loops. The proper resummation necessitates using an (anti-)symmetrised kernel, as we will now show. Note that this holds not only for tree-level 2PI diagrams, but more generally for 2PI diagrams involving loops in t-and u-type configurations. We consider the 4-point function of a pair of identical particles, G IP p, s A , s B p , s A , s B , where the momentum and spin assignments are shown in fig. 14. To ease the notation, the dependence of G IP on the total momentum P is left implicit. Let A be a function of the same variables that stands for the sum of either the t-or u-type 2PI diagrams. Clearly, Then, the sum of complementary 2PI diagrams (u-or t-type respectively) is where f = 0 or 1 if the interacting particles are bosons or fermions. This factor arises from the different number of fermion permutations needed in the t-and u-type cases, in order to perform the Wick contractions. The 4-point function G IP includes the two ladders shown in the two columns of fig. 14. These ladders are related by exchanging the momenta and spins of the initial (or final) state particles, thus we may write In eq. (A.4), we can change the integration variable q → −q and switch r A ↔ r B . Adding up the resulting equation with eq. (A.4), we obtain Evidently, eq. (A.7) is the average of the t-and u-type 2PI diagrams, The factor 1/2 ensures that the loop diagrams are not double-counted. From eq. (A.7), we can also deduce the following relation which we use below in the discussion on the (anti-)symmetrisation of the wavefunctions.
Finally, we note that if the interacting particles carry additional conserved numbers, e.g. non-Abelian (gauge) charges, then appropriate factors ensuring their conservation may appear in the 0th order terms of eq. (A.6), as well as inside A and consequently K. However, eq. (A.8) remains generally valid as is.

A.2 Wavefunctions
The 0th order terms of the Dyson-Schwinger equations determine the normalisation of the wavefunctions (see e.g. [18,62].) The two contributions appearing in the second line of eq. (A.6) ensure that the wavefunctions of identical particles are properly (anti)symmetrised, as we will now show. Instead of deriving the normalisation conditions from eq. (A.6), we shall deduce them by comparing to the case of distinguishable particles (DP), whose wavefunctions are normalised as standard [18,62].
For DP with equal masses, and incoming and outgoing momenta and spins as in fig. 14, the Dyson-Schwinger eq. for the four-point function G DP is (compare with eq. (A.6)) G DP p, s A , s B p , s A , s B = S(P/2 + p ) S(P/2 − p )× whereΨ DP n,s (p) are the solutions to the DP Dyson-Schwinger eq. (A.11), assuming the kernel is the same as that of eq. (A.13). Plugging eq. (A.17) into (A. 16), and considering (A.12), we re-express G IP n,s as G IP n,s (p, p ) = 1 2 G DP n,s (p, p ) + G DP n,s (−p, −p ) + (−1) s 2 G DP n,s (p, −p ) + G DP n,s (−p, p ) .

B Perturbative transition amplitudes: an example
We demonstrate the calculation of diagrams contributing to the perturbative part of the transition amplitudes of section 3. We will work out in detail the diagram shown in fig. 15. Figure 15.
Example of diagram contributing to the perturbative parts of the amplitudes of various transition processes considered in section 3.

C Overlap integral for monopole bound-to-bound transitions
We want to compute the overlap integrals defined in eq.

D Scalar emission via vector-scalar fusion
In many of the BSF processes considered in this work, the radiative parts of the amplitudes receive contributions from diagrams where off-shell vector and Higgs bosons fuse to produce the on-shell radiated Higgs boson; one such diagram is pictured in fig. 16. These diagrams resemble the ones where an on-shell vector is emitted from an off-shell vector or scalar mediator exchanged between the interacting particles (cf. figs. 5 to 8.) This suggests that they may be significant. Here we show that the diagrams of the type of fig. 16 are of higher order than those featuring emission of a vector from an off-shell mediator. Moreover, BSF via vector emission is of higher order than BSF via emission of a charged scalar [1]. Thus, the Higgs emission diagrams of the type of fig. 16 are very subdominant. The contribution from the diagram of fig. 16 is iA =ū(P/2 + p, r 1 ) ig 2 γ µ t a i i u(K/2 + k , s 1 ) Applying the standard approximations due to the scale hierarchies, the above becomes . (D. 2) The 4-vector product in the numerator is of order ∼ m 2 (α 2 B + v 2 rel ), which renders eq. (D.2) of higher order than eqs. (3.14) and (3.18), and even more so that eqs. (3.21).