{'fraction_non_alphanumeric': 0.04483526560555447, 'fraction_numerical': 0.015278050697583022, 'mean_word_length': 4.267747778832054, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Higgs-mediated bound states in dark-matter models

It has been recently demonstrated that the 125 GeV Higgs boson can mediate a long-range force between TeV-scale particles, that can impact considerably their annihilation due to the Sommerfeld effect, and hence the density of thermal relic dark matter. In the presence of long-range interactions, the formation and decay of particle-antiparticle bound states can also deplete dark matter significantly. We consider the Higgs boson as mediator in the formation of bound states, and compute the effect on the dark matter abundance. To this end, we consider a simplified model in which dark matter co-annihilates with coloured particles that have a sizeable coupling to the Higgs. The Higgs-mediated force affects the dark matter depletion via bound state formation in several ways. It enhances the capture cross-sections due to the attraction it mediates between the incoming particles, it increases the binding energy of the bound states, hence rendering their ionisation inefficient sooner in the early universe, and for large enough couplings, it can overcome the gluon repulsion of certain colour representations and give rise to additional bound states. Because it alters the momentum exchange in the bound states, the Higgs-mediated force also affects the gluon-mediated potential via the running of the strong coupling. We comment on the experimental implications and conclude that the Higgs-mediated potential must be taken into account when circumscribing the viable parameter space of related models.


Introduction
The nature of dark matter (DM) is one of the biggest open questions in modern particle physics. An intriguing possibility is that the Higgs boson participates in the physics of the dark sector. The discovery of the Higgs and the measurement of its mass impel the thorough investigation of its implications for various DM scenarios. One of the most generic scenarios for the cosmological production of DM asserts that DM arose from the thermal freeze-out of non-relativistic particles that were previously in equilibrium with the Standard Model (SM) plasma in the early universe. This possibility necessitates a sizeable coupling between the SM and DM, and is thus being probed by a variety of direct, indirect and collider experiments. The increasingly more stringent constraints have pushed the viable DM mass range in this class of models toward and beyond the TeV regime, and have given impetus to co-annihilation scenarios [1][2][3][4][5], which provide more flexibility in reproducing the observed DM density via freeze-out.
If the 125 GeV Higgs couples to a (multi-)TeV dark sector, then the difference between the two mass scales may give rise to long-range effects. Indeed, in models where DM or its co-annihilating partners couple to a much lighter force mediator, the resulting longrange force distorts the wavepackets of the interacting pairs, giving rise to the well-known Sommerfeld effect [6,7] that alters their annihilation rate [8,9]. Long-range interactions imply also the existence of bound states; the formation and decay of particle-antiparticle bound states opens an additional two-step annihilation channel that can deplete the DM abundance [10,11] and contribute to the DM indirect detection signals [12][13][14][15][16]. While such effects emanating from hidden-sector mediators, the electroweak gauge bosons, or the gluons in co-annihilation scenarios, have been considered , only recently the longrange effect of the Higgs was pointed out [40].
The reasoning for disregarding the latter appears to be twofold: the Higgs boson was considered too heavy to yield a long-range force, and the coupling of the Higgs to DM was assumed to be smaller than (or at best comparable to) the SM gauge couplings. However, in co-annihilation scenarios, the mass of DM or its co-annihilating partners can be in the multi-TeV regime and their coupling to the Higgs can be sizeable, such that the range of the Higgs-mediated interaction becomes comparable to or exceeds the Bohr radius of the interacting pair of particles. Then, the interaction manifests as long-range. Reference [40] demonstrated the significant impact of the Higgs-generated potential on the annihilation cross-sections due to the Sommerfeld effect, and consequently on the DM density. In the present paper, we investigate the impact of the Higgs-mediated force on the existence and formation of bound states, and the associated effect on the DM abundance.
For definiteness, we shall consider the simplified model of ref. [40], where DM is assumed to co-annihilate with a coloured scalar that transforms under the fundamental of SU (3) c and possesses a sizeable coupling to the Higgs. This scenario is encountered in the Minimal Supersymmetric SM (MSSM) -with the coloured particle being typically the lightest stop eigenstate -and constitutes one of the last refuges of neutralino DM; as such, it has received a lot of attention recently [2,41,42], with the effect of radiative bound-state formation (BSF) due to gluon exchange computed in [37]. We shall focus on the DM mass range (0. 5 -5) TeV. The lower bound of this interval largely ensures viability against current experimental constraints, while the upper bound implies that DM freezes-out after the electroweak symmetry breaking (EWSB), when the neutral component of the Higgs doublet becomes the mediator that couples to the DM co-annihilating partners. Importantly, this mass interval encompasses the range probed at the LHC and the range of highest interest for the resolution of the hierarchy problem. In fact, retaining the supersymmetry-breaking scale relatively low such that the hierarchy problem is better addressed, and reproducing the measured Higgs mass at the same time, is possible if the coupling of the stops to the Higgs is sizeable. Altogether, this scenario contains the necessary ingredients for the Higgs-mediated force to have an important effect. Similar features are present also in various Higgs-portal models (see e.g. [3,43]), for which we expect our findings to have important implications.
In particular, the impact of the Higgs-mediated potential on the DM abundance affects the prediction of the DM couplings and mass. This in turn changes the interpretation of the experimental results and affects the viability of specific scenarios. The viability of models with long-range interactions with respect to indirect detection constraints is rather sensitive to the predicted DM mass, since the expected γ-ray spectrum from annihilations inside galaxies exhibits sharp parametric resonances with varying DM mass. Moreover, thermal-relic DM models are probed at precision collider experiments. The measurement of the DM abundance to an unprecedented precision, Ωh 2 = 0.120 ± 0.001, by the Planck satellite [44] yields a powerful constraint that can meaningfully complement those from colliders and indirect detection. This however requires that the theoretical computation of the relic density is also sufficiently accurate. During the last couple of years, important progress has been made along this direction. While state-of-the-art numerical tools [45][46][47] have employed so far mainly leading order calculations, it has been demonstrated that higher order corrections to the annihilation processes can impact the DM density up to 20% [48][49][50][51]. The Sommerfeld effect and BSF due to the SM gauge bosons can have an even greater impact [37]. As already shown in [40] and we shall see again in the following, the Higgs-mediated force can have a comparable or larger effect that far exceeds the experimental precision in the measurement of the DM density.
The paper is organized as follows. In section 2, we describe the simplified model, first introduced in [40], and review the calculation of the DM relic density, including the formation and decay of bound states. We then discuss the scattering and bound states in a mixed Coulomb and Yukawa potential in section 3. In section 4, we compute the BSF cross-sections. Since several different momentum scales enter the calculation, in section 4.4 we discuss the running of the strong coupling and explore its effect on the cross-sections. In section 5, we present our results. We discuss the impact of the Higgs boson on the existence and formation of bound states, and on the effective annihilation cross-section, before demonstrating its effect on the predicted DM density. We conclude in section 6. For easy reference, we summarize the notation used throughout the paper in table 1. s, [8] in annihilation and BSF vertices (singlet, octet) α ann s and α BSF s, [1] , α BSF s, [8] Gluon-generated Coulomb potential coupling in the scattering states (singlet, octet)

Description
Dimensionless radial space coordinates for bound and scattering states x B ≡ κr and x S ≡ kr Dimensionless time parameters for freeze-out x ≡ m χ /T,x ≡ m X /T 2 Simplified model and relic density

Simplified model
We assume that DM is a Majorana fermion χ with mass m χ , that co-annihilates with a complex scalar X with mass m X . X is a triplet under SU (3) c and couples to a real scalar h of mass m h = 125 GeV, that aims to resemble the SM Higgs boson. Both χ and X are odd under a Z 2 symmetry, and assumed to be the lightest (LP) and next-to-lightest particles (NLP). The Lagrangian is [40] is the covariant derivative with G a µ being the gluon fields and T a the corresponding generators. We denote the LP-NLP relative mass splitting as and define the total and the reduced mass of an XX † pair, We also define the couplings that will appear in the non-relativistic potential in section 2.3 Let us briefly comment on the origin of the trilinear term h|X| 2 . In the corresponding dimensionful coupling, we have factored out m X for convenience. This factorization does not imply that the coupling originates solely from the quartic term |H| 2 |X| 2 ⊃ |H 0 | 2 |X| 2 , with H and H 0 being the Higgs doublet and its neutral component, whose radial excitation after the EWSB is h. The h|X| 2 term may be also sourced by a gauge-invariant trilinear coupling HX † 1 X 2 , where X 1 and X 2 have different SU (2) L × U (1) Y charges and mix after the EWSB to generate the mass eigenstate X. This dynamics is in fact encountered in the MSSM, where X 1 and X 2 would be the two stop fields that couple to the Higgs via a SUSY-breaking A term. Moreover, the X, X 1 , X 2 fields -being scalars -may also have SU (2) L × U (1) Y invariant mass terms. In order to remain agnostic about the underlying UV complete model, we shall treat g h and m X as independent parameters that are not delimited by one another, or by the Higgs vacuum expectation value, v H 246 GeV.
Here, we are interested in the long-range effect of the Higgs, which emanates from the h|X| 2 coupling. Thus for simplicity, we do not include quartic terms in the scalar potential, or the couplings of the Higgs to the SM particles.

Boltzmann equation for the relic density
If ∆ 1, the DM abundance is determined not only by the χ − χ annihilations, but also by the χ − X, χ − X † and X − X † (co)annihilation processes. Following [52], we defineỸ to be the sum over the yields of all (co)annihilating particles, Figure 1. The gluon and Higgs exchange are the leading order contributions to the 2-particleirreducible (2PI) kernel that gives rise to the long-range interaction between X and X † .
where Y i is the ratio of the number density n i over the entropy density of the universe s ≡ (2π 2 /45) g * S T 3 . We denote by g * S and g * the entropy and energy degrees of freedom respectively, and define The evolution of the DM density is described by the Boltzmann equation As usual, we have defined the dimensionless time parameter x ≡ m χ /T . The yields in equilibrium of χ, X and X † are (g χ = 2, g X = 3) While σ eff v rel in general denotes the effective, thermally averaged cross-section that includes all annihilation and co-annihilation processes, each weighted by the number densities of the interacting species, we assume for the purpose of our study that the NLP self-annihilation is the dominant contribution. Then, the effective cross-section becomes includes all processes that deplete the XX † pairs. It comprises not only of the direct annihilation channels, but also of the BSF processes, whose cross-sections must be weighted by the fraction of bound states that decay rather than being ionised, as we discuss below.

Non-relativistic potential
In the presence of a long-range interaction between two particles, the asymptotic states are not well approximated by plane waves. In order to determine their properties, we must resum the two-particle-irreducible (2PI) diagrams (see e.g. [11]). In our model, the X, X † pairs interact via gluon and Higgs exchange, both of which may result in a sizeable longrange effect [40]. In the non-relativistic and weak coupling regime, the resummation of the gluon and Higgs exchange diagrams shown in fig. 1 amounts to solving the Schrödinger equation with the mixed Coulomb and Yukawa potential where α h has been defined in eq. (2.5) and derived in ref. [15]. The coupling α g is related to α s defined in eq. (2.4) but depends also on the colour representation of the XX † state. The product of a colour triplet and an anti-triplet decomposes into a singlet and an octet state, 3 ⊗3 = 1 ⊕ 8, such that the coupling α g amounts to [53] α g, [1] ≡ α s, [1] × C F = (4/3)α s, [1] , for the singlet and the octet configurations, where C F = 4/3 and C A = 3 are the quadratic Casimir invariants of the fundamental and adjoint representations of SU (3). In eqs. (2.13), we have differentiated between the strong coupling in the singlet and octet states, α s, [1] and α s, [8] respectively, to accommodate the possibly different momentum transfer along the 2PI diagrams, which affects the value of α s due its running. For the massive mediator h, the interaction manifests as long-range if the range of the potential m −1 h is comparable to or larger than the corresponding Bohr radius (µα h ) −1 . To parametrise the range of the interaction we use the dimensionless ratio (2.14) The Sommerfeld effect and BSF become significant when the average momentum transfer between the scattering particles, µv rel , is comparable to or lower than the Bohr momentum, µα. To compare the relevant scales in our set-up, we introduce the parameters where in eq. (2.15a), R = 1, 8 denotes the colour representation. Because the average momentum transfer along the 2PI diagrams is different in the scattering and the bound states, we shall also discern between ζ S g, [R] and ζ B g, [R] , for which we evaluate the strong coupling at the respective scale, as we shall see in the following (cf. section 4.4 for a summary). However, since the running of α h depends on the underlying UV model, we neglect it in the following. For convenience, we define Finally, for the bound states, it will be useful to define the parameter The dimensionless variables defined in eqs. (2.14) to (2.17) suffice to parametrise the pieces of the annihilation and BSF cross-sections that need to be evaluated numerically.
Generically, the Sommerfeld effect and the capture into bound states are expected to be significant for |ζ g | O(1) and ζ h , d h O(1), for the Coulomb and Yukawa potentials, respectively. However, in ref. [40] it was shown that the Coulomb potential affects how long-range the Yukawa potential manifests. In section 3, we revisit this point and determine the conditions for the existence of bound states and their properties.
Thermal effects. The coupling of the mediators -the gluons and the Higgs -to the relativistic SM plasma modifies the zero-temperature potential (2.12) in several ways. On one hand, it generates thermal contributions to the mediator masses that screen the longrange effect. On the other hand, it gives rise to frequent scatterings between X, X † and the relativistic plasma that allow for non-radiative dissipation of energy from an XX † pair, and thus ensure equilibrium between bound and scattering states. In fact, the XX † spectral function morphs into a continuum that encompasses energies both above and below 2m X [27,28]. These thermal effects have been considered in DM freeze-out computations using linear response theory, first in refs. [29][30][31], where radiative capture processes were neglected, and subsequently in ref. [32] which incorporated BSF via gluon emission. Recently, ref. [54] derived ab initio a formalism for the description of DM long-range interactions in a plasma background, starting from non-equilibrium quantum field theory.
In our analysis, we neglect thermal effects. In the scattering states, the screening of the gluon-mediated force due to the thermal gluon masses is negligible, because of the fairly large average velocity of the X, X † particles during DM freeze-out, which ensures that the average momentum transfer in an XX † pair is larger than the thermal gluon mass. The radiative BSF processes we consider are sufficiently rapid to bring bound states in equilibrium at early times. In equilibrium, and at temperatures much higher than the binding energy, the DM depletion via bound-state decay becomes independent of the BSF rate, whether this is due to radiative and/or scattering processes; we return to this point in section 2.5. More importantly, we find that the impact of the bound states on the DM relic density occurs mostly at later times, when the plasma temperature approaches or falls below the binding energy. In this regime, the thermal corrections become less important, and the radiative capture is typically the dominant BSF process.

Annihilation
The dominant contributions to the direct annihilation of XX † pairs arise from the XX † → gg and XX † → hh channels, shown in fig. 2. We neglect the p-wave suppressed contributions XX † → qq, gh. (Similarly we neglect the p-wave suppressed contribution of the s-channel diagram of XX † → gg). To leading order, the colour-averaged perturbative cross-sections are [40] (σv rel ) perturb where α ann s ≡ α s (Q = m X ) arises from the gluon emission vertices in the hard annihilation process and must be evaluated at momentum transfer Q = m X . Both the colour-singlet and Figure 2. Tree-level diagrams contributing to the X − X † annihilation into gluons and into Higgs bosons. The u-channel diagram that match the corresponding t-channel diagram (first diagram) is not shown. The s-channel annihilation into gluons (third ) yields a p-wave contribution such that we neglect this process in our computations. octet scattering states contribute in eq. (2.18a), while only the singlet can annihilate into hh. Due to their different gluon-mediated potential the two states are affected differently by the Sommerfeld effect. The full cross-sections are [40] where S 0, [1] and S 0, [8] indicate the s-wave Sommerfeld factors for the colour-singlet and octet states respectively. For the mixed Coulomb and Yukawa potential (2.12), the s-wave Sommerfeld factor depends on the parameters defined in eqs. (2.14) and (2.15), i.e. S 0 = S 0 (ζ S g , ζ h , d h ), as we discuss in section 3.1. Thus, taking eqs. (2.13) into account, [8] denotes the strong coupling evaluated at the average momentum transfer in the scattering states, Q = µv rel , which is independent of the colour representation. The thermally-averaged total annihilation cross-section is

Bound-state formation, ionisation and decay
Formation. Attractive long-range interactions imply the existence of bound states that may form from the scattering states via dissipation of energy. In the non-relativistic regime, Figure 4.
The amplitude for the radiative capture into bound states consists of the (nonperturbative) initial and final state wavefunctions, and the perturbative 5-point function that includes the radiative vertices. The wavefunctions (cf. section 3) are determined by the resummation of the 2-particle-irreducible (2PI) diagrams. The 2PI kernel is shown in fig. 1. the energy available to be dissipated is the difference between the initial and final states, i.e. the kinetic energy of the scattering state in the centre-of-momentum (CM) frame is the generalised Bohr momentum of the bound state and γ n parametrises the departure of E n from its Coulomb limit d h → ∞ (cf. section 3.2). For capture into the ground state, the dissipated energy is The capture into bound states via emission of an ultrasoft gluon is depicted in figs. 4 and 5. The decompositions 3 ⊗3 = 1 ⊕ 8 and 8 ⊗ 8 = 1 S + 8 A + 8 S + 10 A + 10 A + 27 S imply that the allowed transitions are (X + X † ) [8] → B(XX † ) [1] + g [8] , In (2.23b) and (2.23c), the indices outside the curly brackets denote the colour representation of the entire final state (bound state plus gluon). When considering the gluon potential solely, only the colour-singlet XX † combination supports a bound state, thus (2.23a) is the sole capture process via gluon emission. However, if the coupling to the Higgs is sufficiently strong, the Higgs exchange can overcome the gluon-mediated repulsion in the octet state and give rise to octet bound states. Then, the processes (2.23b) and (2.23c) may be possible. In section 3.2, we determine the condition for the existence of bound states in a mixed repulsive Coulomb plus attractive Yukawa potential, and in section 4.4 we adapt this condition to the colour-octet states in our model, taking into account the α s running.
If the energy (2.22) available to be dissipated exceeds the Higgs mass, then bound states may form also via Higgs emission, according to However, the capture of particle-antiparticle pairs via scalar emission is subject to cancellations that suppress the processes (2.24) by higher orders in α s and α h with respect to (2.19) and (2.23) [15, eqs. (4.10)]. Although other couplings in the scalar potential, such as the quartic h 2 |X| 2 term, can enhance and even dominate the cross-sections for the processes (2.24) [55], these cross-sections remain subdominant to the direct annihilation and/or BSF via gluon emission. In the following, we thus neglect BSF via Higgs emission. The thermally-averaged BSF cross-section is where f g (ω) = 1/(e ω/T − 1) is the gluon occupation number, with ω being the energy of the emitted gluon, given by eq. (2.22). The factor 1+f g (ω) accounts for the Bose enhancement due to the final-state gluon, and is necessary to ensure the detailed balance between the bound-state formation and ionisation processes at T ω, which encompasses a significant temperature range that is relevant to the DM freeze-out [10]. 1 Once bound states form, they may either be ionised back into their constituents by the ambient radiation, or decay into radiation, as we now discuss.
Ionisation. The ionisation cross-section is related to the BSF cross-section via the Milne relation (cf. e.g. [37, appendix D]), where R S , R denote the colour representations of the scattering and bound XX † states, with g [1] = 1 and g [8] = 8 being the colour-singlet and octet degrees of freedom respectively, and g g = 8 being the gluon degrees of freedom. In eq. (2.26), we have made explicit the dependence of the gluon energy (2.22) on the colour representation of the bound state. Note that for the ionisation cross-section of the colour-octet bound states -if they existwe must include the contributions from both processes (2.23b) and (2.23c), i.e. sum over R S = 1, 8 as indicated in eq. (2.26). The ionisation rate of a bound state is then As is standard, we have omitted the Bose enhancement factor in the thermal averaging of the annihilation cross-sections (2.21), since the energy of the annihilation products, E = mX, far exceeds the temperature of the relativistic bath at the temperature range relevant for freeze-out, T mX/20.
where ω min [R] is the minimum energy required for ionisation, i.e. the binding energy, recovered from eq. (2.22) for v rel = 0. Using eqs. (2.22) and (2.26), we arrive at Note that in eq. (2.27), γ 1,0 also depends on R even though it is not explicitly indicated. From eqs. (2.25) and (2.27), we can explicitly verify the principle of detailed balance, where n eq X and n eq B, [R] are the equilibrium densities of X and of XX † bound states of colour representation [R], respectively. Note that eq. (2.28) is more general than our derivation here may suggest; it holds true independently of what processes -radiative or scattering -may contribute to the capture into and the ionisation of bound states. It is clear from (2.27) and (2.28), that Γ ion decreases exponentially when T drops below the binding energy and most particles in the plasma do not have enough energy to dissociate the bound states.
Decay. The decay rate of the = 0 bound states is related to the perturbative s-wave annihilation cross-section times relative velocity as (see e.g. [11]) where ψ [R] n00 (0) indicates the n-th level s-wave bound-state wavefunction of the colour representation R, evaluated at the origin, and will be discussed in section 3.2. The cross-section σ s−wave ann, [R] v rel corresponds to the colour configuration of the bound state and should be averaged over the bound-state colour degrees of freedom only, rather than those of an unbound XX † pair. Using the perturbative annihilation cross-sections (2.18) and taking into account the colour decomposition in eq. (2.19a), we obtain for the ground states (2.30b) Effective bound-state formation cross-section. The effect of unstable bound states on the DM relic abundance is described by a system of coupled Boltzmann equations for the bound and unbound particles that describe the interplay between bound-state formation, ionisation and decay [10]. Only the bound states that decay into radiation before being ionised contribute to the depletion of the DM abundance. It is thus possible to describe the impact of bound states using a single Boltzmann equation for the unbound particles and an effective BSF cross-section that incorporates the branching ratio of the bound states that decay rather than being ionised. For our model, this is [8] . (2.31) The decay becomes faster than ionisation at temperatures comparable to or smaller than the binding energy, when most particles in the thermal bath do not possess enough energy to disassociate the bound states and Γ ion decreases exponentially with decreasing temperature. The effective BSF cross-section (2.31) together with the thermally averaged annihilation cross-sections (2.19) yield the total effective annihilation cross-section (2.11) that determines the DM relic density according to the Boltzmann eq. (2.8).
Ionisation equilibrium. If the BSF and ionisation cross-sections are sufficiently large to set Γ ion Γ dec at high temperatures, then the effective BSF cross-section (2.31) becomes independent of the actual BSF cross-section during that time. Instead, it is proportional to the annihilation cross-section, which determines the bound-state decay rate. Indeed, using the detailed balance eq. (2.28) and the decay rate (2.29), we find To assess the importance of eq. (2.33), it is informative to consider a simple attractive Coulomb potential of strength α, for which |ψ 1,0,0 (0)| 2 = (µα) 3 /π. Then, eq. (2.33) yields π −5/2 [m X α 2 /(4T )] 3/2 1, where the square bracket contains the binding-energy-totemperature ratio, |E 10 |/T . This contribution can be compared with the thermally averaged Sommerfeld factor from the scattering states, which in this regime is Clearly, the enhancement due to the bound states is subdominant to that due to the scattering states. Including excited states in eq. (2.32) would somewhat increase the contribution of the bound levels, but not change this conclusion.
The depletion of DM via BSF becomes more significant as the temperature of the universe lowers and approaches the binding energy [10]; this is in part manifested by the exponential factor in eq. (2.32). In this regime, the bound-state decay becomes comparable to or faster than ionisation, which is now exponentially suppressed, [cf. eq. (2.27)], and σ BSF v rel eff rapidly saturates to σ BSF v rel [cf. eq. (2.31)]. It is important to note that if σ BSF is comparable to exceeds σ ann , then it is possible that the DM depletion via BSF contributes significantly to the effective annihilation rate even starting from temperatures that are larger than the binding energy by a factor of a few (cf. refs. [10,37] and section 5).

Scattering states
The non-relativistic potential due to gluon and Higgs exchange for the scattering states is where we will consider both positive and negative α S g , but only α h 0.
Wavefunctions. The scattering states are described by a wavefunction φ k (r), parametrised by the continuous quantum number where k and v rel are the expectation values of the momenta of the interacting particles in the CM frame, and of their relative velocity. The strong coupling α S s that determines α S g in eq. (3.1) is evaluated at Q = |k|. The scattering state wavefunction obeys the Schrödinger equation They are normalised according to To solve the Schrödinger equation numerically, we perform the separation of variables and introduce the dimensionless radial space coordinate such that the radial Schrödinger equation reads where we have used the parameters ζ S g ≡ α S g /v rel , ζ h ≡ α h /v rel and d h ≡ µα h /m h defined in eqs. (2.14) and (2.15). At x S → 0, and for > 0, the second term of eq. (3.8) is dominated by the centrifugal contribution. In this region, the two independent solutions of eq. (3.8) scale as x +1 S (regular) and x − S (irregular). We are interested in the regular solutions, which imply the boundary condition lim x→0 χ |k|, (x S ) = ( + 1) lim (3.9) The condition (3.9) will be valid also for = 0. To fully specify the wavefunction χ |k|, , we shall also use the asymptotic behaviour at x S → ∞. At large x S , the wavefunction χ |k|, behaves as (see e.g. ref. [56, chapter 7]) where the phase shifts δ depend on ζ S g , ζ h and d h . This implies that (3.10) Coulomb limit. In the limit d h → ∞, the Schrödinger eq. (3.8) depends only on ζ S = ζ S g + ζ h , and can be solved analytically, (3.11a) where 1 F 1 is the confluent hypergeometric function of the first kind, and (3.11b) The sum over the modes in eq. (3.6) can be expressed in closed form such that the total scattering-state wave function in the Coulomb limit reads (3.11c) Note that the limit d h → 0 (pure gluon-mediated potential) can be also obtained from the above with the substitution ζ S → ζ S g .
Sommerfeld factor S 0 for s-wave annihilation. This is (see e.g. [57]) (3.12) In the Coulomb limits d h → 0 and d h → ∞, the Sommerfeld factor reduces to with S C 0 (ζ) given in eq. (3.11b). The features of S 0 and the impact of Higgs enhancement on the relic density have been discussed in detail in ref. [40] (see in particular figs. 2 and 3 therein).

Bound states
The non-relativistic potential for the bound states, due to gluon and Higgs exchange, is Note the difference in α g with respect to eq. (3.1). While α S g in the scattering potential is to be evaluated at the scale Q = µv rel , α B g in the bound state potential is to be taken at the Bohr momentum scale, as will be discussed below. 16) where {n m} are the standard principal and angular momentum quantum numbers, and the wavefunctions are normalised as

Wavefunctions. The Schrödinger equation for bound states is
We define the generalised Bohr momentum such that the discrete binding energy levels of the system are where the factor γ n is determined numerically (see below) and parametrises the departure from the Coulomb limit d h → ∞ where γ n = 1/n. We perform the separation of variables Defining the dimensionless space coordinate we find that χ n is normalised as and obeys the radial Schrödinger equation We recall that the dimensionless parameters λ ≡ α B g /α h and d h ≡ µα h /m h have been defined in eqs. (2.14) and (2.17). It follows that the bound-state wavefunctions χ n and the (normalised) energy eigenvalues γ n depend on λ and d h only. Assuming α h > 0, the necessary condition κ > 0 for bound states to exist implies λ > −1. This condition is necessary but not sufficient, as we shall see below.
As in the case of the scattering states, we are interested in the regular solutions of eq. (3.23), which scale as (3.24) The discrete spectrum of eigenvalues γ n = γ n (λ, d h ) is then determined by requiring that χ n vanish at infinity, lim We present some numerical results in fig. 6 and discuss them next. Coulomb limits. In the left panel of fig. 6, the two plateaus indicate the two Coulomb limits, d h → 0 (pure gluon-mediated potential) and d h → ∞ (Coulomb limit of the Higgsmediated potential), which allow for analytic solutions of eq. (3.23).
• In the limit d h → ∞, the spectrum of eigenvalues and the radial wavefunction are where L a n are the generalised Laguerre polynomials of degree n, and we assume the normalisation condition • In the limit d h → 0, bound states exist only for λ > 0 (α B g > 0). In this case, the binding energy takes the Coulomb value for the gluon-generated potential only, where κ g ≡ µα B g . This implies . (3.29) Existence condition. If the Coulomb component of the potential is not attractive, λ 0 (equivalently, α B g 0), the Yukawa potential has to be sufficiently long-range for bound states to exist, i.e. d h has to be sufficiently large. In the right panel of fig. 6, we show the critical value d h,crit (λ) such that for d h d h,crit (λ) at least one bound level exists. We find that in the range −0.96 λ 0, d h,crit can be fitted by the analytical expression

Overlap integrals
In order to evaluate the amplitudes for the BSF processes of interest, we need the following overlap integrals of the scattering and bound state wavefunctions [37] where we used the following Fourier transforms For the numerical evaluation of the integrals (3.31), we express them in terms of the radial wavefunctions χ k, and χ n introduced in sections 3.1 and 3.2. The integrals J k,{n m} have been previously analysed in ref. [15, appendix B]. The leading order contribution to J k,{n m} (b) with b ∝ P g is independent of b [11,15] such that in the following we denote J k,{n m} = J k,{n m} (b = 0). Adapting the result to our notation for capture into the ground state, {n m} = {100}, we arrive at 2 where we recall from eq. (2.16b) that ζ B = κ/k. Similarly, the integral (3.31b) becomes The angular integral above is (see e.g. ref. [15, eq. (B.7b) with b → 0]) dΩr P (k ·r) = δ ,1 (4π/3)k.
Thus, we arrive at the final expression We shall use eqs. (3.33) to evaluate the BSF cross-sections in section 4.
Coulomb limit. In the limit d h → ∞, the scattering and bound state wavefunctions, (3.11c) and (3.27) are 11b)], and 1 F 1 is the confluent hypergeometric function of the first kind. Following refs. [11,15,37,58], we derive the Coulomb limit of the overlap integrals using the identity [59] and we used again κ/k = ζ B . As noted before, it suffices to evaluate J k,{100} (b) at b = 0 [11,15] such that Y k,{100} and J k,{100} are related as with the squared amplitude of J C k,{100} being In section 4.3, we shall use eqs. (3.37) and (3.38) to obtain analytical expressions for the BSF cross-sections in the Coulomb limit.

Bound-state formation cross-sections
Having outlined how to evaluate the scattering and bound state wavefunctions, we have the tools at hand to calculate the cross-sections for the BSF processes discussed in section 2.5.

Amplitudes
The radiative capture via gluon emission has been computed in ref. [37] for general groups and representations, in terms of the overlap integrals (3.31). Adapting the results to our model, we find the following colour-averaged squared amplitudes: where α BSF s ≡ α s (Q = |P g |) is the strong coupling evaluated at momentum transfer equal to the momentum of the emitted gluon; this factor arises from the gluon emission vertices in the diagrams of fig. 5. |P g | is found from energy-momentum conservation to be |P g | = E k − E n = (µ/2)[(α h +α B g ) 2 γ 2 n +v 2 rel ]. We recall that in eqs. (4.1), the superscripts on the overlap integrals denote, in order, the scattering-state and bound-state colour representations.

Coulomb limit
In the limit d h → ∞, we may use the analytical expressions (3.37) and (3.38) for the overlap integrals. We find the colour-averaged BSF cross-sections to be where f c is a numerical factor that depends on the transition, s, [8] α h − α B s, [8]  and (4.5c) We emphasise that, even when not explicitly denoted above, the couplings α BSF s , α B s , α B g , α S g , and therefore also ζ S ≡ (α S g + α h )/v rel and ζ B ≡ (α B g + α h )/v rel , depend on the colour representations of the initial and final states, and are thus different for every transition.
The function S C BSF (ζ S , ζ B ) encapsulates all the velocity dependence of σ C BSF v rel . The first two factors in eq. (4.5c) arise solely from the scattering state wavefunction and correspond to the Sommerfeld effect on p-wave processes in the Coulomb limit,  [8] , ζ B, [1]  (4.7)

The running of α s
As it has been pointed out throughout this paper, the various vertices where the strong coupling appears in the annihilation and BSF diagrams are characterised by different momentum transfers Q. 3 Since the strong coupling is sensitive to the momentum transfer, α s = α s (Q), it is important to carefully account for its running. An overview of the various scales is given in table 2. 4 In the occasions that the momentum transfer itself depends on the strong coupling, Q = Q(α s ), we solve numerically the following equation forα, The effect of the running of α s is depicted in figs. 7 to 10, with related discussion in the captions. A few remarks are in order here.
Scattering states. The average momentum transfer Q = µv rel , implies that α S s becomes suppressed for large m X . Consequently, for large m X , the attraction due to the Higgs can significantly ameliorate or even overcome the repulsion due to the gluons in the colour-octet state. Due to the large multiplicity of the latter, (cf. eq. (2.19a)), this enhances significantly the total annihilation rate, even for moderate values of α h . 5 This is depicted in fig. 7.
Bound states. Because the momentum transfer depends on the strong coupling, we solve numerically eq. (4.8). Generally, α h increases the average momentum transfer, thus suppressing α B s , as shown in fig. 8. Overall though, the Bohr momentum κ ≡ µ(α B g + α h ) increases with increasing α h . If α h is sufficiently strong, it overcomes the gluon-mediated repulsion in the octet states, and gives rise to bound levels. In the octet ground state, where we recall that d h = µα h /m h , and the function γ 1,0 (λ, d h ) has been determined numerically and is depicted in fig. 6. For this momentum transfer, we see graphically in the left panel of fig. 9 that eq. (4.8) has a solution when the coloured lines of fixed α h intersect the horizontal black line. Let us examine the scaling of α s (Q(α))/α withα.
• For sufficiently smallα (α 6α h ), the momentum transfer (4.9) is independent of α, i.e. Q µα h × γ 1,0 (0, d h ). Since α s ∼ 0.1 for the m X values of interest, the ratio α s (Q(α))/α is simply 0.1/α, i.e. starts from values 1 and decreases withα. 3 In fact, the smallness of the momentum transfer along the mediators in the ladder diagrams is responsible for the emergence of the non-perturbative phenomena, the Sommerfeld effect and bound states. 4 An improved treatment would incorporate the running of the coupling in the Schrödinger equations for the scattering and bound states, by setting Q = 1/r. The approximations of table 2 are sufficient for our purposes. We refer to [60][61][62] for the computation of the αs running in NRQCD and pNRQCD. 5 This holds provided that v rel is not too low. See related discussion in section 5.1.
The minimum value of the Higgs coupling, α h, crit , required for the octet bound states to exist is shown in the right panel of fig. 9.
Gluon-emission vertices in BSF. The momentum transfer -which depends here as well on the strong coupling -is the softer scale that enters our calculations, and thus yields the largest values of α s (for fixed values of the other parameters). Since α h increases the binding energy, it suppresses α BSF s . This is seen in fig. 10.

Vertices
s, [8] 6α h , d h Table 2. The momentum transfer Q at which the strong coupling α s (Q) is evaluated. For the bound states, the functions γ n (λ, d h ) are computed numerically (cf. section 3.2 and fig. 6). We recall that λ [R] ≡ α g,[R] /α h with R = 1, 8, and αh=0.05  Upper left: α S s evaluated at the average momentum exchange in the scattering state, Q = µv rel ; for the colour-singlet state α S g, [1] = (4/3)α S s , while for the colour-octet state α S g, [8] = −α S s /6 within the indicative velocity range 0.2 < v rel < 0.4 that is typical during the DM freeze-out. For comparison, we show α ann s occurring in the perturbative annihilation processes taken at a momentum transfer Q = m X (red diamonds).
Upper right: The parameter ζ g = α S g /v rel that determines the Sommerfeld effect due to gluon exchange only, for the singlet and octet configurations.
Lower left: The coloured bands show the thermally-averaged Sommerfeld factor for s-wave annihilation processes,S 0 , evaluated within the indicative range 20 <x ≡ m X /T < 50, during which the DM relic density is mostly determined. The Higgs exchange enhances the attraction / reduces the repulsion between the scattering-state particles (blue bands: α h = 0.05) with respect to gluon-only exchange (gray bands: α h = 0). The thin solid and dotted lines mark the Coulomb limit, m h → 0, which becomes a good approximation for large enough m X .
Lower right: The thermally-averaged total Sommerfeld factor for annihilation into gluons. The Higgs exchange leads to a significant enhancement (blue band: α h = 0.05, purple band: α h = 0.2) with respect to gluon-only exchange (gray band: α h = 0).  Figure 8.
The strong couplings α B s, [1] and α B s, [8] and the corresponding α B g, [1] = 4α B s, [1] /3 and α B g, [8] = −α B s, [8] /6 that determine the colour-singlet (left) and the colour-octet (right) bound-state wavefunctions. The coupling to the Higgs, α h , increases the average momentum transfer within the bound states, thereby suppressing α s . While the effect is modest in case of the colour singlet, a large enough α h is required for the existence of the octet bound state.   Figure 10. The strong coupling α BSF s, [1] and α BSF s, [8] at the gluon emission vertex during the radiative capture of a colour-singlet (left panel ) or colour-octet (right panel ) bound state. The emitted gluon carries away the binding energy of the bound state plus the kinetic energy of the scattering state (see table 2) that depends on the relative velocity v rel . As before, we show α BSF s, [1] and α BSF s, [8] in the typical velocity range 0.2 v rel 0.4 during freeze-out. As v rel decreases, the momentum transfer at the vertex drops and α BSF s increases. In the right panel, the cut-off of the bands for low m X reflects the fact that colour-octet bound states exist only for large enough m X and α h (cf. fig. 9).

Binding energy [GeV]
Colour-octet bound states Figure 11. The binding energy of the singlet and octet bound states, B [1] and B [8] , for various values of the coupling to the Higgs, α h (solid lines). The dotted red lines show the approximate temperature of freeze-out, T FO ≈ m X /30.

Scattering states, bound states, annihilation and BSF rates
Based on the computations of previous sections, we now present and summarise the impact of the Higgs-mediated force on the properties of the scattering and bound states, and on the various cross-sections of interest.
Scattering states and Sommerfeld effect. The Higgs-mediated potential enhances the attraction in the colour-singlet state, and suppresses or overcomes the repulsion in the octet state. Notably, the effect of the Higgs potential is influenced by the presence of the gluon-mediated potential, as was already shown in ref. [40]. In particular: • In the singlet state, the Higgs enhancement becomes sizeable -and potentially reaches its Coulomb limit for a given coupling strength -for lower masses (more generally, for lower d h ≡ µα h /m h ) than in the absence of the attractive gluon-mediated force [40, fig. 2].
• Conversely, in the octet state, a larger d h is required for the Higgs to affect the longrange interaction. However, the gluon-mediated repulsion is suppressed by a colour factor with respect to the full strength of the strong coupling, α g, [8] = −α s, [8] /6, suggesting that even a modest α h can potentially counteract it. This holds so long as v rel is not too low; at low v rel , the gluon-induced repulsion becomes exponential [cf. the Sommerfeld factor in the Coulomb limit (3.11b) with ζ = ζ S g, [8] −1], and cannot be surmounted by a finite-range attractive force (cf. ref. [40, fig. 3, left panel]). However, the DM freeze-out occurs while v rel is fairly large and the gluon-induced repulsion is milder; this implies that the Higgs-mediated force has a significant effect.
Tighter bound states. The Higgs attraction increases the absolute value of the binding energies, as shown in fig. 11. This renders the bound-state dissociation processes inefficient earlier, when the DM density is greater, thus enhancing the efficacy of BSF in depleting DM. In fig. 11, we compare the binding energy with the typical temperature around DM freeze-out, T FO ≈ m X /30. If the binding energy equals or exceeds T FO , then the dissociation of bound states is inhibited already at freeze-out, and the efficiency of BSF in depleting DM is maximal. While this occurs only for very large values of α h , we emphasise that the effective BSF rate can become comparable to or exceed the annihilation rate even at temperatures that are a factor of a few larger than the binding energy, even if it is not maximal (cf. fig. 13).
Additional bound states. As already discussed in section 4.4 and shown in fig. 9, for sufficiently large masses and couplings, the Higgs-mediated attractive force implies the existence of colour-octet bound states. However, these bound states are considerably looser than the singlet states, as seen in fig. 11, and have a limited effect on the DM density (except perhaps for α h 0.3).
Impact on the strong coupling. The Higgs-mediated force increases the momentum transfer in the bound states and in the gluon-emission vertices of the radiative capture processes. Because of the running of the strong coupling, this suppresses the corresponding values of α s , as discussed in section 4.4 and depicted in figs. 7, 8 and 10. However, the Bohr momentum and the binding energy increase with α h , as seen in fig. 11.
Cross-sections. In fig. 12, we illustrate the velocity dependence of the BSF crosssections, and compare them with the direct annihilation processes.
• For large v rel , the BSF cross-sections are rather suppressed and subdominant to the annihilation. This is due to the small overlap of the wavefunctions: the average momentum in the scattering states, k = µv rel , is much greater than that in the bound states, which in the Coulomb approximation is the Bohr momentum, κ = µ(α B g +α h ). In this regime, the BSF cross-sections scale as g + α h , the overlap of the wavefunctions is nearly maximal (although the precise value of v rel at which this occurs depends also on the Bohr momentum of the scattering state, hence on α S g +α h ). In this regime, the cross-section for capture into the tightest bound state (colour singlet) dominates over annihilation. Note that since α s ∼ 0.1, this velocity range is relevant for freeze-out.
• At v rel α S g + α h , the Sommerfeld effect becomes important. Because of the gluon-mediated repulsion in the octet scattering states, the [8] → [1] and [8] → [8] capture processes become suppressed with decreasing v rel . The coupling to the Higgs ameliorates this suppression, and large α h can even make these BSF cross-sections increase temporarily as the velocity drops; this can be seen in the right panel of fig. 12. However, at sufficiently low v rel , the exponential Coulomb repulsion cannot be overcome by the finite-range Higgs-mediated attraction.
On the other hand, the [1] → [8] cross-section increases steadily with decreasing v rel , due to the gluon and Higgs-mediated attractive forces in the colour-singlet scattering state. Similarly, the annihilation cross-section becomes dominated by the coloursinglet contribution and increases with decreasing v rel .
Note that in fig. 12, we show the various cross-sections normalized to the perturbative annihilation cross-section, such that the direct dependence on m X cancels out. Nevertheless, m X affects the scale of α s and thus has indirect impact on the various lines.
Let us now focus on the effect of the Higgs. As expected, the annihilation cross-section increases with larger α h . This is both due to the XX † → hh channel, and because the Higgs-mediated attractive potential enhances the annihilation cross-sections for all final states. In part because of the increase in the (perturbative) annihilation cross-section, the relative strength of BSF appears to diminish with increasing α h . The suppression of α BSF s at larger α h due to the running of the strong coupling (cf. fig. 10) contributes to this trend. However, the larger Bohr momentum of the bound states implies that BSF peaks at a larger v rel . This is very important for the DM depletion via BSF, as we discuss next.
Thermally averaged cross-sections. In fig. 13, we depict the impact of the Higgs on the thermally averaged annihilation and BSF rates. The features of σ ann v rel follow from the discussion above. For BSF, we present both the actual and the effective thermally averaged cross-sections, σ BSF v rel and σ BSF v rel eff . Several comments are in order.
• Despite the suppression of σ BSF v rel at large v rel discussed above, σ BSF v rel appears to be large even at high T ; this is due to the Bose-enhancement factor [cf. eq. (2.25)].
• Because of the rapid ionisation processes, σ BSF v rel eff is suppressed at high T , and in fact becomes independent of σ BSF v rel , as discussed in section 2.5. However, σ BSF v rel eff increases as T drops.
For capture from an octet state, σ BSF v rel eff peaks at T ∼ binding energy. The reason is three-fold: (i) This condition is the thermally-averaged equivalent of k ∼ κ, which maximises the overlap of the bound and scattering-state wavefunctions, as discussed above. (ii) At and below this temperature, the ionisation rate becomes exponentially suppressed, therefore σ BSF v rel eff saturates to σ BSF v rel . (iii) At even lower temperatures, the gluon repulsion in the octet state suppresses σ [8]→ [1] BSF v rel and σ [8]→ [8] BSF v rel .
• At least for 0.02 BSF v rel eff overcomes σ ann v rel , even before it saturates to σ [8]→ [1] BSF v rel , i.e. at temperatures higher than the binding energy of the singlet bound state. Overall, BSF yields a sizeable contribution to the total effective annihilation cross-section σ XX † v rel eff for the range of couplings considered.

Dark matter relic density
We now turn to the impact of the Higgs-mediated force on the relic density. We perform the calculation as described in section 2, using the BSF cross-sections of section 4. In the left panels of figs. 14 and 15, we show, for different values of α h , the NLP-LP mass difference ∆m versus the DM mass m χ that give rise to the observed DM abundance. The width of the bands reflects the uncertainty in the measurement of the DM abundance by Planck. We compare different computations that incorporate combinations of the various nonperturbative effects discussed in this work: the Sommerfeld effect on the direct annihilation processes and BSF, due to the gluon and Higgs mediated potentials. We present the impact of these effect on the relic density in the right panels of figs. 14 and 15.
For α h = 0.02, the Higgs enhancement of the direct annihilation gives only a moderate correction on the predicted relic density. However, the full calculation that includes annihilation and BSF with gluon and Higgs exchange, changes the predicted mass gap by up to 10 GeV, with respect to considering gluon exchange only. Given that the usual tools for computing the relic density include a perturbative calculation only, their estimation can be off by a factor of 2.5 − 8.5 depending on the DM mass, even for a small value of α h .
With increasing α h , the effect of BSF with gluon exchange only appears to become less significant. This is because it has remained unaffected by α h , while the direct annihilation into two Higgs bosons has become very rapid. However, including the Higgs-mediated potential gives rise to a very sizeable effect. This is true when considering direct annihilation only, and even more so when including BSF. We point out indicatively that for α h = 0.1 and α h = 0.2, ∆m is predicted to be larger by up to 40 GeV and 70 GeV respectively.
The above clearly demonstrate that in coloured co-annihilation scenarios, both the Sommerfeld effect and BSF, as well as both the gluon-exchange and Higgs-exchange potentials must be considered in order to obtain a reliable estimation of the relic density.

Caveats
The present work along with refs. [37,40] aimed to demonstrate the long-range effect of the Higgs boson and the impact of BSF on the DM abundance. In order to not obfuscate these goals, various issues were not addressed, but would need to be properly considered in a comprehensive analysis. We briefly mention the main ones.
Our simplified model does not specify the LP-NLP interactions that are hypothesised to establish equilibrium between the two in the early universe. Departure from equilibrium during the NLP freeze-out would imply that the NLP annihilation (direct or via BSF) depletes DM less efficiently than implicitly assumed by the set of eqs. (2.6) to (2.11). This in turn would shift the contours of figs. 14 and 15 to lower values of ∆m and m χ .
The LP-NLP interactions, along with their mass difference, determine also the NLP decay rate. The latter must be smaller than the decay rate of the NLP-NLP bound states in order for BSF to deplete DM. This potentially implies that BSF cannot affect the DM density if the LP-NLP mass difference exceeds a threshold that depends on the strength of the LP-NLP interactions. Such a condition can be specified within complete models, e.g. MSSM scenarios, where the couplings of the UV theory determine on the one hand the NLP decay rate, and on the other hand the entirety of the annihilation processes that control the bound-state decay rate. Moreover, these couplings determine the value of α h .
Including all significant annihilation channels in the analysis of course affects also the total effective annihilation rate -both directly, and via the bound-state decay rate that determines the efficiency of BSF in depleting DM. This could have a considerable impact on the results shown in figs. 14 and 15, and therefore on the interpretation of the experimental constraints. For very large α h , the radiative capture via Higgs emission that was neglected here, may also become significant, as noted in section 2.5.
Ultimately, the precise computation of the DM density requires taking into account the thermal effects that were briefly discussed in section 2.3. This could be potentially done using the formalism of ref. [54], which can account for the transition of the BSF processes from in-to out-of equilibrium. However, this formalism must be first extended to encompass radiative processes that involve ultrasoft modes, which are currently neglected.
Finally we note that sizeable couplings of the coloured particles to the Higgs endanger the stability of the SU (3) c symmetric vacuum. The determination of vacuum stability requires a dedicated study for any specific scenario under consideration. However, it has also been argued that this kind of dynamics could imply the emergence of a new phase in the MSSM, where the standard treatment of the vacuum stability does not apply [63].      Ω / Ω DM Figure 14.
Left panels: The mass difference ∆m versus the DM mass m χ that reproduce the observed DM density when taking into account: perturbative annihilation only (grey dashed ), annihilation with Sommerfeld effect due to gluon exchange only (grey dotted ), annihilation with Sommerfeld effect due to gluon and Higgs exchange (green), annihilation with Sommerfeld effect and bound state formation with gluon exchange only (blue), annihilation with Sommerfeld effect and bound state formation with gluon and Higgs exchange (red ). The width of the bands corresponds to the 3σ experimental uncertainty on the DM abundance, using the Planck 2018 results.
Right panels: For the values of ∆m predicted by the full calculation (red bands on the left panels), we show the ratio of the relic density from each partial calculation to the observed DM abundance. The colour coding is the same as on the left panels.

Conclusions
While our searches for DM coupled to the electroweak gauge interactions have returned null results until now, it is plausible that the Higgs constitutes the portal to the dark sector. The discovery of the Higgs boson and the measurement of its properties in collider experiments urge the comprehensive assessment of its implications for DM.
If TeV-scale DM or its co-annihilating partners couple to the 125 GeV Higgs, then this interaction may have a sizeable long-range effect. In recent work, we showed that the DM relic abundance in coloured co-annihilation scenarios is affected significantly by the enhancement of the direct annihilation processes due to the Higgs-mediated force [40], and by the formation and decay of bound states due to gluon exchange [37]. Here we considered the effect of the Higgs-generated potential on the formation of bound states, and demonstrated the impact on the DM abundance. With respect to available tools that include perturbative calculations only, we have found that the predicted DM density may differ by up to almost one order of magnitude. While we focused on a particular class of models, we expect that the long-range effect of the Higgs boson is important in a variety of scenarios where the Higgs constitutes the portal to the dark sector.
Altogether, in the set-up we considered, we have shown that the Higgs-mediated force affects the DM density in a variety of ways that exhibit some salient features. In summary: • The attractive interaction in the scattering states due to the Higgs exchange enhances both the direct annihilation cross-sections, as well as the capture into bound states.
• The Higgs-mediated potential increases the binding energy of the colour-singlet bound states. This implies that the capture rate becomes maximal at larger velocities, and consequently at higher temperatures in the early universe. Moreover, the bound-state dissociation becomes insignificant starting at higher temperatures. As a result, BSF begins to deplete DM efficiently at earlier times, when the DM density is higher, and therefore has a larger impact on the relic abundance.
• The Higgs counteracts the gluon-mediated repulsion in the octet states. For large enough mass of the interacting particles and coupling to the Higgs, colour-octet bound states exist. Their formation and decay contributes to the DM depletion rate in the early universe, even if only modestly due to their small binding energy.
• The Higgs exchange increases the momentum transfer in the bound states and at the gluon emission vertex in the capture process, hence driving the strong coupling to smaller values. This somewhat quells the strength of the BSF processes at large α h .
• The interference of the (gluon-generated) Coulomb potential and the (Higgs-generated) Yukawa potential influences the long-range effect of the latter, both in the scattering and the bound states.
The enhancement of the effective DM annihilation rate due to the Higgs-mediated force, implies that the LP-NLP mass difference and/or the DM mass must be larger than previously predicted, in order for DM to attain the observed density via freeze-out. This is particularly important for collider searches, since a larger mass gap yields harder jets that are easier to probe. Similarly, the increased values of the predicted DM masswhich can lie even beyond the range considered here -strengthen the motivation for indirect searches in the multi-TeV regime. However, the accurate interpretation of the experimental constraints necessitates that the effects discussed here are considered within complete models, where various relevant technicalities can be treated properly, as discussed in section 5.3.