Bound states of WIMP dark matter in Higgs-portal models. Part II. Thermal decoupling

The Higgs doublet can mediate a long-range interaction between multi-TeV particles coupled to the Weak interactions of the Standard Model, while its emission can lead to very rapid bound-state formation processes and bound-to-bound transitions. Using the rates calculated in a companion paper, here we compute the thermal decoupling of multi-TeV WIMP dark matter coupled to the Higgs, and show that the formation of metastable dark matter bound states via Higgs-doublet emission and their decay decrease the relic density very significantly. This in turn implies that WIMP dark matter may be much heavier than previously anticipated, or conversely that for a given mass, the dark matter couplings to the Higgs may be much lower than previously predicted, thereby altering the dark matter phenomenology. While we focus on a minimal singlet-doublet model in the coannihilation regime, our calculations can be extended to larger multiplets where the effects under consideration are expected to be even more significant.


Introduction
Model-independent unitarity arguments [1], as well as model-specific calculations show that viable thermal-relic dark matter (DM) scenarios in the multi-TeV regime feature interactions that generate long-range potentials and give rise to bound states. The formation and decay of metastable bound states alters the DM decoupling in the early universe [2] and contributes to its indirect signals [1,[3][4][5][6][7][8]. In fact, in a variety of models, bound-state formation (BSF) can be faster than annihilation [9][10][11][12][13][14][15]. This is particularly striking in models where DM couples to a light scalar boson charged under a symmetry. It was recently shown that the emission of such a scalar boson by a pair of interacting particles alters their effective Hamiltonian and can result in extremely rapid monopole capture processes [13]. This may potentially have severe implications for multi-TeV DM coupled to the Higgs doublet. Moreover, it has been shown that the 125 GeV Higgs boson can mediate a significant long-range interaction between TeV-scale particles, despite being heavier than all other Standard Model (SM) force mediators [12,16].
Here we will focus on a minimal singlet-doublet realisation of this class of scenarios. In a companion paper, we have computed the non-relativistic potentials, the BSF cross-sections, the bound-state decay rates and bound-to-bound transitions [35]. In this paper, we employ the results of [35] to compute the DM thermal decoupling from the primordial plasma in the early universe. Several aspects of this calculation are examined in detail, including the effect of both radiative BSF and BSF via scattering on the relativistic thermal bath, the importance of bound-to-bound transitions and the late recoupling of the DM destruction processes due to BSF via Higgs emission. As in [35], we carry out all computations assuming electroweak symmetry; indeed, for DM heavier than 5 TeV, freeze-out begins before the EWPT, even though the complete thermal decoupling may occur much later. The validity of the approximation is discussed in detail.
We note that our calculations are important broadly for scenarios that introduce new stable species in the fashion considered here, even if these species do not account for DM. Besides specifying the parameters for which the relic density of the lightest new mass eigenstate matches the observed DM abundance, we show how the relic density is affected by the new effects within a broader parameter space. As is standard, the relic density of stable particles sets a cosmological constraint on new physics scenarios.
This work is organized as follows. In section 2, we summarise the model following [35], and briefly review its basic properties in the broken electroweak phase. Moreover, we discuss the temperature dependence of the Higgs-doublet mass in relation to the BSF processes of interest. In section 3, we lay out the formalism for computing the DM thermal decoupling, and discuss the interplay among bound-state formation, ionisation, decay and bound-to-bound transition processes. The results of this computation are presented and discussed in section 4. We conclude in section 5.

Lagrangian and mass eigenstates
We begin by specifying the model following ref. [35], before summarising the mass eigenstates and interactions in the broken electroweak phase.
We consider a gauge-singlet Majorana fermion S = (ψ α , ψ †α ) T of mass m S , and a Dirac fermion D = (ξ α , χ †α ) T of mass m D with SM gauge charges SU L (2) × U Y (1) = (2, 1/2). S and D are assumed to be 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. We take m S and m D to be real. This can 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 will focus on the regime where S and D can co-annihilate efficiently before the EWPT of the universe. This occurs if their masses are similar, within about 10%. A larger mass difference would imply that the number density of the heavier species in the non-relativistic regime becomes negligible, since it would be suppressed with respect to the lighter species by exp[−(δm/m)x] 0.1, with x ≡ m/T ∼ 25 during DM freeze-out. Since the masses must be similar for the processes considered here to be relevant, we take them to be equal for simplicity, (2. 2) It will be useful also to introduce the reduced mass of a pair of DM particles, In addition, in order to reduce the number of free parameters, we set which we take to be real. (The CP violation is not important for our purposes.) Our computations can of course be extended to more general Yukawa couplings. As is standard, we define the couplings For later convenience, following [35], we also define the couplings The model of eq. (2.1) and various aspects of its phenomenology have been considered for general parameters extensively in the past [17][18][19][20][21][22][23][24][25][26]. Here we will only 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).

Mass eigenstates in the broken electroweak phase
At electroweak symmetry breaking, the neutral component of the Higgs doublet acquires a vacuum expectation value. In terms of their SU L (2) components, the H, ξ and χ fields are with v H 246 GeV being the Higgs vacuum expectation value. We define the left-handed multiplet of the neutral statesN α ≡ (ψ α , ξ α , χ α ) T . Then, by inserting eq. (2.7) into the Lagrangian (2.1), we find the corresponding mass terms, We diagonalise eq. (2.8) by setting where U is the unitary matrix The corresponding mass eigenvalues are In addition to the neutral states, there is a charged Dirac fermion of mass m.

Interactions in the broken electroweak phase and constraints
The interactions among the neutral states, N ≡ (N 1 , N 2 , N 3 ) T , are described by the Lagrangian [19] where c W = g 2 / g 2 1 + g 2 2 . Since the coupling to Z µ is non-diagonal, with the mass splitting being always much larger than ∼ O(100 keV) for the y values we will consider here (cf. section 3), the constraints from direct detection experiments due to this interaction are evaded. On the other hand, the coupling to the Higgs boson is expected to yield sizable DM-nucleus scattering and potentially strong constraints. Existing analyses of the direct detection data for this model do not extend to the multi-TeV regime that is of interest here. Moreover, the direct detection constraints on the model (2.1) are generally significantly relaxed around the socalled blind spots where the coupling to the Higgs vanishes, roughly when y L = −y R (see e.g. [19,22].) Constraints may also arise from electroweak precision observables, and in particular from the contribution to the T parameter. While the latter scales as (y 2 L − y 2 R ) 2 and vanishes in the limit considered here, it may become important for large values of the Yukawa coupling(s) if |y L | |y R | or |y L | |y R |. We refer to [18][19][20] for related studies. While a detailed phenomenological analysis is beyond the scope of the present work, our results are important for interpreting the experimental constraints, since they imply a different relation between the DM mass and couplings in order for the observed DM density to be attained via thermal freeze-out.

Bound states
The interactions of eq. (2.1) -in particular the exchange of B, W and H bosons -generate long-range potentials among the S, D andD species that distort the wavefunctions of scattering states and give rise to bound states. The long-range dynamics in the symmetric electroweak phase has been discussed in detail in the companion paper [35], and shall not be repeated here. However, since our focus is the effect of bound states on the DM thermal decoupling, in table 2 we summarise for convenience the bound levels we consider.

Higgs doublet mass and EWPT
The cross-sections for BSF via H emission depend on the Higgs doublet mass [35]. Taking into account the finite temperature 1-loop corrections to the effective potential (see e.g. [36,37]), we estimate that before the EWPT of the universe, the Higgs doublet mass was where m h 125 GeV is the Higgs boson mass at zero temperature, λ = m 2 h /(2v 2 H ) 0.13 is the Higgs quartic coupling, V SM ⊃ −λ|H| 4 , and y t 0.994 is the top quark Yukawa coupling. The EWPT occurs as m 2 H (T ) → 0, i.e. at estimated temperature In computing the DM decoupling, we use eq. (2.14) at T T EWPT , and set m H → m h at T < T EWPT while still using the annihilation and BSF rates computed under the assumption of electroweak symmetry. We discuss this approximation in section 4.2.
We may now estimate whether or when m H (T ) implies that BSF via Higgs emission is kinematically suppressed. In a thermal distribution, the energy dissipated during BSF averages to ω = 3T /2 + |E B | (cf. eq. (3.23) and [35].) The first term suffices to provide for m H (T ) for all T > T EWPT since m H (T > T EWPT ) 0.63T , as well as after the EWPT, down to temperatures T ∼ 2m h /3 83 GeV. However, since the BSF cross-sections weigh preferentially low values of v rel , the kinematic suppression may become important at somewhat larger T than this estimate implies, unless |E B | is sufficient to provide for m H .

Coupled Boltzmann equations
Let Y j ≡ n j /s and Y B ≡ n B /s be the number-density-to-entropy-density ratios of the free species j and the bound state B respectively. In our model, j = S, D,D and B = SS/DD, DD, DD, DS (cf. table 2.) We are ultimately interested in the total DM yield (3.1) Note that the bound states are metastable and their abundance becomes eventually negligible, so we do not include them in eq. (3.1). As is standard, we will use the time parameter The entropy density of the universe is s = (2π 2 /45)g * S T 3 = (2π 2 /45)g * S m 3 /x 3 . We denote by g * S and g * the entropy and energy dof respectively, and define The evolution of Y j and Y B is governed by the coupled Boltzmann equations , (3.5) and the equilibrium densities in the non-relativistic regime are where g i are the spin and SU L (2) dof of the free species, with g S = 2, g D = gD = 4. The dof g B and the binding energies E B of the bound states we consider are listed in table 2.
For later convenience, we also define the total DM dof g DM ≡ g S + g D + gD = 10, and the equilibrium density of (3.1) (3.7) In the above, Γ dec B , Γ ion B→ij and Γ trans B→B are respectively the rates of B decay into radiation, ionisation (a.k.a. dissociation) to ij, and transition into the bound level B . The rates Γ j→i describe the transitions between free particles, due to decays, inverse decays and/or scatterings on the thermal bath; overall, these processes do not change the DM number density, but retain equilibrium among the dark species. Note that in eqs. (3.4) we must use the thermally averaged rates, Γ . The thermal average introduces Lorentz dilation factors for decay processes -which however are insignificant in the non-relativistic regime -as well as Bose-enhancement factors in the case of transitions and capture or ionisation processes. We discuss this in more detail in section 3.3. The thermally-averaged rates and cross-sections of inverse processes are related via detailed balance that we have already employed in writing eqs. (3.4), The fractional relic DM density is where Y ∞ is the final yield, and we have included the mass shift of the lightest state that arises after the electroweak symmetry breaking (cf. eq. (2.12)); this is significant only for the lower end of the mass range we consider and for large couplings α H . In eq. (3.9), s 0 2839.5 cm −3 and ρ c 4.78 · 10 −6 GeV cm −3 are the entropy and critical energy densities of the universe today [38].

Effective Boltzmann equation
The system of coupled Boltzmann eqs. (3.4) is numerically difficult to solve. We shall thus adopt an effective method that reduces eqs. (3.4) to one equation for the DM yield (3.1).
For convenience, we first define the total formation cross-section, ionisation rate and transition rate of every bound state B, We begin by assuming that the i ↔ j interactions are sufficiently rapid to ensure chemical equilibrium among the free species, such that Y i /Y eq i = w, where w is the same for all i = S, D,D. Due to their rapid decays, inverse decays and transitions to other bound levels, the bound states are typically close to equilibrium, thus dY B /dx 0. Under this assumption, eqs. (3.4b) yield a system of linear equations for Y B that can be solved and re-employed in eq. (3.4a) [39]. For bound states that do not participate in any boundto-bound transitions, this simplifies to (3.11) In the model under consideration and within our approximations [35], the spin-1 DD and DD bound states do not participate in any bound-to-bound transitions, while the spin-0 SS/DD and DS bound states can rapidly transit into each other via H absorption/emission (cf. [35, table 6].) For the latter, eqs. (3.4b) read where we set Γ dec DS 0 [35], and we recall that Γ trans SS/DD→DS = Γ trans DS→SS/DD (Y eq DS /Y eq SS/DD ), due to detailed balance eq. (3.8a). The factors 2 in the first row account for transitions to and from the two conjugate bound states DS andDS.
Next, we use eq. (3.11) for the DD and DD yields and the solution of eq. (3.12) for the SS/DD and DS yields, in the Boltzmann eq. (3.4a). Summing over all free particle species, we find that the evolution of Y is governed by the Boltzmann equation where the equilibrium density Y eq is given in eq. (3.7). The DM destruction cross-section σv rel eff receives contributions from direct annihilation and BSF processes, and where the factors 2 in the DD and DS terms account also for the formation of the conjugate bound states. The individual contributions are found as follows. For the bound-states that do not participate in any bound-to-bound transitions, while for the coupled bound states In eqs.

Effective cross-section
We now consider in more detail the contributions to the effective DM destruction crosssection in our model, based on the computations of ref. [35]. We begin with direct annihilation in section 3.3.1, and then discuss BSF in section 3.3.2. In figs. 1 and 2 we illustrate the contributions to BSF, while in fig. 3 we compare all contributions to the DM destruction cross-section for a chosen set of parameters, showcasing the effect of the Higgs potential and of BSF via Higgs emission.

Annihilation
In our model, the total annihilation cross-section is where the indices denote the two-particle scattering states, with dof g SS = 4, g DD = 16, g DD = 16, g DS = 8. The DD and DS contributions carry factors of 2 to account also for the annihilation of the conjugate states, and DD and DS carry factors of 2 to account for the two distinguishable particles annihilated in each process.
In the above, the s-wave Sommerfeld factor is where ζ S ≡ α S /v rel , with α S being the strength of the Coulomb potential of the scattering state. The various ζ S appearing in eqs. (3.18) are Each of the terms in eqs. (3.18) is the product of the dof, the perturbative annihilation cross-section and the Sommerfeld factor of a spin-and gauge-projected state. The DD cross-section includes contributions from both the SU L (2) singlet and triplet projections, which are characterised by different non-relativistic potentials and thus have different Sommerfeld factors. For the singlet states, the potential depends also on the spin. Indeed, the spin-0 SU L (2) singlet SS and DD states mix due to the Higgs mediated potential. Since the perturbative s-wave annihilation of the SS component vanishes, the contribution from the annihilation of the SS-like state has been included in the DD-like state, for simplicity. We refer to [35] for more details.

Bound state formation
The radiative BSF cross-sections have been summarised in [35, tables 7-10], and we shall denote them here as σ rBSF B v rel [xx] with xx and B being the scattering and bound states. 1 As is well known, for pairs of identical particles (here SS, DD,DD), this factor is canceled upon thermal averaging by the factor 1/2 needed to avoid double-counting of the initial particle states [40].
Bound states can also form through scattering on the thermal bath, via exchange of an off-shell mediator; the corresponding cross-section factorise in their radiative counterparts and temperature-dependent functions [35,41,42]. Collecting these results, the velocityweighted cross-sections σ BSF B v rel for the formation of the various bound-state species B receive the following contributions from the individual channels (cf. [35, tables 7-10]) In each term above, the first factor accounts for the number of DM particles destroyed (upon thermal averaging), as well as the capture of the conjugate scattering state if applicable, in analogy to eq.
is the energy dissipated in the capture process [9,35], with α B being the strength of the potential of the corresponding bound state (cf . table 2 and [35, table 6].) R H depends also on m H /ω, and essentially replaces the phase-space suppression due to on-shell H ( †) emission that is included in the radiative cross-sections. The R H , R B , R W factors can be found in [35]. Next, we must thermally average eqs. (3.22). In BSF, the emitted boson carries away a small amount of energy that can be comparable to the temperature of the primordial plasma during the DM decoupling. The Bose enhancement due to the final state boson can thus be significant, and must be included in thermal averaging the BSF cross-sections to ensure that detailed balance holds [2], 2 As seen in eq. (3.16), the contributions of each bound level to the effective DM destruction cross-section (3.14) have to be waited by the appropriate branching fractions that account for the portion of bound states that decay into radiation thereby reducing the DM density. The bound-state decay and transition rates needed to compute these branching fractions can be found in [35, table 6]. In thermally averaging these rates, we may neglect the Lorentz dilation factor that is 1 in the non-relativistic regime. However, the low-energy boson emitted in bound-to-bound transitions implies a Bose enhancement that must be included to ensure detailed balance at temperatures higher than the dissipated energy. So, where the dissipated energy here is ω = m(α 2 A − α 2 H )/4, and h H (ω) is the phase-space suppression defined in eq. (3.24). Finally, the bound-state ionisation rates are computed using the detailed balance eq. (3.8b), and summing over all ionised states as in eqs. (3.10a) and (3.10b); this yields

Ionisation equilibrium
Equation ( v rel is sufficiently large. If so, the system reaches a state of ionisation equilibrium, where the effective BSF cross-sections (3.16) become independent of the actual ones [43], where for the DS bound state whose direct decay into radiation is suppressed, we must use the effective decay rate (cf. eq. (3.16d)) (3.28) Since the bound-state decay rates are proportional to the annihilation cross-sections of the corresponding scattering states (cf. e.g. ref. [35]), eq. (3.27) implies that at high temperatures and while ionisation equilibrium holds, the BSF contribution to the DM destruction rate is negligible in comparison to that of direct annihilation (cf. e.g. [13, eq. (3.20)].) However, as T approaches or drops below |E B |, the ionisation rates become exponentially suppressed and are overcome by the bound-state decay and/or bound-to-bound transition rates. For the uncoupled bound states DD and DD, this implies that the effective BSF cross-sections increase exponentially until they saturate to their maximum values, the actual BSF cross-sections. For the SS/DD and DS coupled system, Γ trans DS→SS/DD > Γ ion DS , occurs before the decay rates surpass the ionisation rates; in this interval, the effective BSF cross-sections (3.16) together with the detailed balance eq. (3.26c), imply that ionisation equilibrium holds for the sum of the SS/DD and DS contributions, where again we neglected the SS/DD decay against ionisation rate. At even lower temperatures, when ionisation becomes slower than decay, the effective BSF cross-sections reach their actual values. We illustrate the above in figs. 1 and 2, where we also compare radiative BSF and BSF via scattering. Two observations are useful more generally for calculations of freeze-out with bound states: • In some (but not all) cases, BSF via scattering dominates at early times; BSF via Higgs exchange may also dominate at late-times over on-shell emission due to the phase-space suppression of the latter. Nevertheless BSF via scattering does not change significantly the effective BSF cross-section with respect to considering radiative BSF only, because overall it becomes subdominant while the system is still in ionisation equilibrium, or around the time it exits it.
• For the DD, DD bound states, ionisation equilibrium ceases at T > |E B | (cf. fig. 1.) In contrast, the bound-to-bound transitions prevent the SS/DD and DS coupled system to reach ionisation equilibrium. However, it closely tracks it until much lower temperatures, T |E B |, due to the largeness of the BSF cross-sections (cf. fig. 2.) We also note here that the computation of the DM thermal decoupling (cf. section 4) shows that much of the BSF effect on the relic density arises after the system exits ionisation equilibrium. (This was also found in ref. [2].) The above imply that it is not safe to estimate the BSF effect by assuming ionisation equilibrium until T ∼ |E B | and neglecting any effect thereafter, an approach previously adopted in refs. [44][45][46][47][48]. Considering instead the BSF cross-sections is necessary for an accurate computation.

Timeline and relic density
Collecting all the above, we are now ready to compute the DM decoupling and relic density. We consider and compare the cases described in table 3, and recall that our calculations always assume electroweak symmetry. We discuss this approximation in section 4.2.   fig. 3), we thus do not present these cases separately.

AnnS
In fig. 4, we present an example of the time evolution of the effective cross-section and the DM density. For the parameters chosen, the exponential increase of σv rel eff due to BSF when the ionisation processes cease, gives rise to a second period of DM destruction that decreases the DM density by two orders of magnitude! In fig. 5, we show the timeline of the DM thermal decoupling. We define the recoupling period of DM destruction due to BSF as the interval between the two occurrences when d 2 (ln Y )/d(ln x) 2 = 0, and the chemical decoupling as the latest time when d(ln Y )/d(ln x) = 10%. In the same plot, we also mark the EWPT, as well as the time beyond which the finite Higgs mass affects its long-range effect. Since in part of the parameter space, the recoupling occurs after the EWPT and the full chemical decoupling occurs even later, the effect of BSF via Higgs emission is most important for the range of DM masses where the binding energies exceed the Higgs boson mass, m h 125 GeV. These ranges are also marked in fig. 5. (We discuss the validity of various approximations, including that of electroweak symmetry, in section 4.2.) In fig. 6, we show the values of α H vs. m that reproduce the observed DM density, as well as the impact of the various processes on the relic density. As already seen in fig. 4, at m few TeV, BSF via emission of a Higgs doublet is estimated to decrease the relic density by up to two orders of magnitude. The implications are twofold. For a fixed mass m, the coupling α H is predicted to be almost up to an order of magnitude smaller than when neglecting BSF via Higgs emission. This should be expected to change (relax) experimental constraints very significantly. Conversely, for a given coupling, a much larger m is anticipated. In fact, DM masses almost up to the unitarity limit can be attained for α H < 1. (We discuss the unitarity limit in more detail in section 4.3.) This motivates experimental searches at very high masses.
To understand better the effect of the various bound states, in fig. 7 we show the α H −m relation determined by considering direct annihilation plus each of the four bound states separately. The spin-1 DD and DD bound states have only a small effect because their binding energy is independent of α H and somewhat small. This implies that ionisations inhibit the DM destruction via their formation until late, when BSF via Higgs emission is    as the grey dotted lines in fig. 7 indicate.
Even away from the correlation of parameters that reproduces the observed DM density, the BSF effect on the relic abundance of the stable species can be very large as seen in fig. 8. The parameter space where the relic density is cosmologically insignificant is greatly enlarged. This is important for scenarios that do not aspire to explain the DM density, but nevertheless predict the existence of stable particles.

Major approximations and their validity
We now summarise the main approximations made in our analysis and comment on their potential effect on the estimated relic density.
(i) Considered only ground-level bound states.
BSF via vector or neutral scalar emission is dominated by dipole and quadrapole moments respectively. In these cases, the capture into the ground state is the dominant BSF process [2,[9][10][11]49], the reason being twofold: it is the most exothermic process, and the overlap of scattering and bound state wavefunctions is larger.
In contrast, BSF via emission of a charged scalar is a monopole transition, and the capture into excited states can be comparable to or faster than the capture into the ground state, despite the latter being more exothermic [13]. This suggests that in the present model, capture into excited states via Higgs emission may be important.
Nevertheless, independently of the BSF cross-sections, the relative effect of the excited states on the relic density is moderated by their smaller binding energy that renders their ionisation efficient until later. We thus anticipate that in the present model excited states may have a significant albeit not dominant effect that would further diminish the relic density and alter the coupling-mass relation along the direction found here. This is worth pursuing in more detail in the future.
(ii) Regularisation of inelastic cross-sections in parametric regimes where BSF via Higgs emission approaches or appears to exceed the unitarity limit.
The efficiency of BSF via emission of a charged scalar implies that the computed cross-sections may reach or violate the upper limit on inelastic cross-sections implied by unitarity (cf. eq. (4.4)) even for small or moderate values of the coupling of the interacting particles to the scalar [13]. The restoration of unitarity implies that resummed higher order corrections, i.e. higher-order contributions to the non-relativistic potential, must be considered [1]. In the context of the present model, the problem was discussed in ref. [35, section 3.6], where the regularisation scheme of ref. [50] was adopted as an effective method to ensure that unitarity is not violated. However, as discussed in ref. [35], this scheme is not entirely suitable for the case of interest.
In fig. 6 we compare the α H − m relations with and without regularisation. Clearly, at large m the effect is significant; the regularisation of the cross-sections ensures that m does not exceed the unitarity limit on the mass of thermal relic DM [1,51], which we discuss in section 4.3. This suggests that working out a more accurate regularisation scheme that would address the issues discussed in ref. [35] may be important in order to obtain more accurate results. We leave this for future work.
(iii) Neglected the Higgs mass in the Higgs-mediated potential.
The validity of the Coulomb approximation for the Higgs-mediated potential has been discussed in [35, section 2.3.3], where the relevant conditions have been put forward.
Here we discuss their validity during the DM thermal decoupling.
Scattering states: In a thermal bath, the condition µv rel > m H for the validity of the Coulomb approximation implies 3mT /2 m H . Considering the Higgs doublet mass (2.14) before the EWPT, this becomes T 3m, which clearly covers all of the range of interest for the DM freeze-out (T m/25.) Below the EWPT, where m H → m h 125 GeV, the condition is satisfied until after the DM chemical decoupling, as shown in fig. 5. Therefore, the Coulomb approximation does not pose any problem.
Bound states: The condition µα H > few×m H , becomes xα H > few before the EWPT, with x = m/T . This is satisfied for all relevant x and α H for which BSF has an effect (x 50 and α H 0.1, cf. figs. 5 and 6.) It is easy to check that this condition is also satisfied below the EWPT for all relevant DM masses and couplings (m > 5 TeV and α H 0.1, cf. fig. 6.) Note that this is not coincidental; BSF via Higgs emission does not have a significant effect for lower α H values because of the phase-space suppression due to the Higgs mass (cf. fig. 7.) The estimation here thus confirms the argument of ref. [35] that bound states are nearly Coulombic in the parameter space where their formation is kinematically allowed and significant.
In fig. 5, we see that the DM destruction via BSF may be efficient after the EWPT. The breaking of the electroweak symmetry has several important implications that we now discuss.
a. The Goldstone modes of the Higgs doublet are absorbed by the Z, W ± bosons.
• BSF via emission of a Higgs doublet in the unbroken electroweak phase corresponds to BSF via emission of h or the longitudinal modes of the Z, W ± bosons in the broken electroweak phase. The Goldstone boson equivalence theorem implies that the amplitudes for BSF via emission of the longitudinal Z, W ± components are the same as those for the corresponding processes in the unbroken electroweak phase, in the limit that the energy of the emitted vector boson is much larger than its mass, m 2 Z,W /ω 2 1. In our computation, the phase-space suppression sets m H → m h 125 GeV after the EWPT, ensuring that m h /ω < 1 or equivalently m 2 Z,W /ω 2 < 0.5. We thus regard the approximation as acceptable, especially in the parameter space away from the phase-space thresholds (cf. fig. 7.) The importance of monopole BSF processes in a broken gauge phase due to the Goldstone boson equivalence theorem was previously pointed out in ref. [52]. • The potential mediated by the Higgs doublet in the unbroken electroweak phase is mediated by h and the longitudinal Z, W ± components in the broken phase. To compute the non-relativistic potential generated by the latter, we need their contribution to the vector boson propagators, where q and m V = g V v H /2 denote the vector boson momentum and mass, for V = Z, W ± , with g Z = g 2 1 + g 2 2 and g W = g 2 . In general, the exchange of where q = p 1 − p 1 = p 2 − p 2 , and we used the Dirac equation / pu(p) = mu(p). Considering the mass splittings ∼ yv H , this becomes The non-zero Z, W ± masses curtail the range of the potentials generated by the exchange of both their transverse and longitudinal components, and introduce phase-space suppression to the BSF processes occurring via their emission. The validity of the Coulomb approximation for the Z, W ± bosons can be assessed as in the preceding discussion for the Higgs. However, in the present model, the B, W -generated potentials and BSF via B or W emission do not have a significant effect, due to the fact that one of the dark multiplets is a gauge singlet and the other belongs to a small representation. We thus do not consider the transverse Z, W ± components further. The effect of the longitudinal Z, W ± components was discussed above.
c. The components of the DM multiplets acquire different masses.
After acquiring a mass splitting, the various pairs of Z 2 -odd particles can oscillate into each other according to the non-relativistic potentials computed in [35, section 2] provided that the kinetic energy of their relative motion exceeds their mass difference. This necessitates mv 2 rel /4 > 2yv H , which, upon thermal averaging, becomes T > (4/3) √ 4πα H v H . This condition is not satisfied below the EWPT for the α H values of interest (α H 0.1.) We thus expect that the rates of some of the processes below the EWPT will be lower than estimated here. This is probably the most severe limitation of our computation. To assess its impact, in fig. 6 we include the coupling-mass relation obtained by integrating the Boltzmann equations only up to the EWPT. Clearly, a proper treatment would result in an α H −m relation between our this and the result obtained by integrating until late times. We see that even when the integration stops at the EWPT, the Higgs effect is still very significant, even if it appears only for larger α H values. The impact on the relic density reaches up to a factor of a few.

Unitarity limit on the dark matter mass
The unitarity of the S matrix sets an upper limit on the partial-wave inelastic cross-sections, σ inel σ uni = (2 + 1)π k 2 (2 + 1)π µ 2 v 2 rel , (4.4) where is the partial wave and k is the momentum of either of the interacting particles in the CM frame. The last approximation in eq. (4.4) concerns the non-relativistic regime, where k = µv rel with µ being the reduced mass. The upper limit (4.4) suggests that for very large masses, annihilations in the early universe may not suffice to reduce the density of thermalised particles to the observed DM value. It thus sets an upper bound on the mass of thermal relic DM annihilating predominantly via a finite number of partial waves in the early universe [51]. For selfconjugate DM in thermal equilibrium with the SM plasma, this is [1,2]   Equation (4.5) is modified by 1/ √ 2 in the case of non-self-conjugate DM. The parametric dependence of σ uni on µ and v rel implies that the limit (4.4) can be attained down to arbitrarily low velocities -thus the upper limit (4.5) on the mass of thermal-relic DM can be reached -only if there is an attractive long-range force between the interacting particles, and provided of course that the relevant couplings are sufficiently large [1]. Attractive long-range interactions imply also the existence of bound states, whose formation and subsequent decay may decrease the DM abundance more efficiently than direct annihilation [2]. This means that BSF may essentially be the dominant process that saturates the unitarity limit (4.4), and that additional partial waves to those dominating in the annihilation processes may become important, thereby increasing the upper limit on the DM mass [1].
In general, the Weak interactions of the Standard Model are not sufficiently strong to generate cross-sections that approach the unitarity limit (4.4), unless perhaps the interacting particles belong to very large SU L (2) representations. However, BSF via emission of a scalar charged under a symmetry can be very efficient even for small couplings [13]. Here, we have seen that BSF via Higgs emission can raise the predicted WIMP mass very significantly, bringing it potentially close to the unitarity limit.

Conclusion
Our DM searches are currently at the onset of the exploration of the multi-TeV regime with a variety of existing and upcoming telescopes observing high-energy cosmic rays. In this mass regime, within the thermal-relic scenario, the DM interactions are expected to manifest as long-range and give rise to non-perturbative effects, in particular the Sommerfeld effect and the formation of bound states. These effects may operate in the early universe during the DM thermal decoupling, as well as inside DM haloes today, and significantly alter the DM phenomenology.
In the present work, consisting of this and a companion paper [35], we have considered the role of the Higgs doublet as a light force mediator, in the thermal decoupling of multi-TeV WIMP DM. We have shown that the Higgs-doublet-mediated potential between DM particles and the formation of DM bound states via Higgs-doublet emission can dramatically change the predicted relic density. This, in turn, alters the coupling-mass relation that reproduces the observed DM abundance. Moreover, it greatly expands the parameter space where the stable relics do not overclose the universe, even if they are a subdominant component of DM. In the former case, the modified coupling-mass relation implies that on one hand, for a given DM mass, existing constraints may be significantly relaxed, and on the other hand, DM may be much heavier than previously anticipated, potentially approaching the unitarity limit.
While the amplitude for BSF via Higgs-doublet emission can be quite large even for small couplings of the DM multiplets to the Higgs, the Higgs-doublet mass introduces a kinematic suppression to the cross-section that renders this effect relevant for larger DM masses and/or couplings to the Higgs. In the specific singlet-doublet scenario considered here, we found that the effect is significant for m 5 TeV and α H 0.1. However, in models involving larger SU L (2) representations, the gauge interactions contribute more significantly to the binding energy of the bound states, thereby rendering the phase-space suppression less significant. We thus expect that the effect on the relic density will be important even for lower couplings.
Finally, we note that the capture into excited bound levels, which we neglected here, may also have a sizeable effect due to the monopole nature of the transitions occurring via Higgs-doublet emission. On the other hand, we have found that including BSF through scattering on the relativistic thermal bath via an off-shell Higgs doublet does not affect the relic density significantly.