Higgs-mediated bound states in dark-matter models

TL;DR

This paper investigates how the Higgs boson can mediate long-range forces and bound state formation between dark matter particles, significantly affecting their relic abundance and constraining model parameters.

Contribution

It introduces a simplified model analyzing Higgs-mediated bound states and their impact on dark matter depletion, highlighting effects overlooked in previous studies.

Findings

Higgs-mediated forces enhance bound state formation and dark matter depletion.

Higgs coupling influences the binding energy and ionization efficiency of bound states.

The Higgs potential must be considered when constraining dark matter models.

Abstract

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.

Figures

Tables

Table 1. Table 1: Summary of notation. For the running of the strong coupling, see table   2 .

Description	Symbol
DM particle	$\chi$
DM mass	$m_{\chi }$
Coloured coannihilating partners of DM	$X$ and $X^{†}$ [scalars, $\mathrm{𝟑}$ and $\overline{\mathrm{𝟑}}$ under $SU{\left(3\right)}_c$]
Mass of coannihilating partners	$m_X$
Total mass of pair of coannihilating partners	$M\equiv 2m_X$
Reduced mass of pair of coannihilating partners	$\mu \equiv m_X/2$
Relative mass difference between DM and coannihilating partners	$\Delta \equiv \left(m_X-m_{\chi }\right)/m_{\chi }$
Mass of Higgs-like boson $h$	$m_h=125\mathrm{GeV}$
$X,X^{†}$ coupling to the Higgs	${\alpha }_h=g_h^2/\left(16\pi \right)$
Strong coupling  [$Q$: momentum transfer]	${\alpha }_s\left(Q\right)\equiv g_s^2\left(Q\right)/\left(4\pi \right)$
in scattering and bound states (singlet, octet)	${\alpha }_s^S$   and   ${\alpha }_{s,\left[\mathrm{𝟏}\right]}^B,{\alpha }_{s,\left[\mathrm{𝟖}\right]}^B$
in annihilation and BSF vertices (singlet, octet)	${\alpha }_s^{\mathrm{ann}}$   and   ${\alpha }_{s,\left[\mathrm{𝟏}\right]}^{\mathrm{BSF}},{\alpha }_{s,\left[\mathrm{𝟖}\right]}^{\mathrm{BSF}}$
Gluon-generated Coulomb potential coupling
in the scattering states (singlet, octet)	${\alpha }_{g,\left[\mathrm{𝟏}\right]}^S\equiv \left(4/3\right){\alpha }_s^S,{\alpha }_{g,\left[\mathrm{𝟖}\right]}^S\equiv -{\alpha }_s^S/6$
in the bound states (singlet, octet)	${\alpha }_{g,\left[\mathrm{𝟏}\right]}^B\equiv \left(4/3\right){\alpha }_{s,\left[\mathrm{𝟏}\right]}^B,{\alpha }_{g,\left[\mathrm{𝟖}\right]}^B\equiv -{\alpha }_{s,\left[\mathrm{𝟖}\right]}^B/6$
Bohr momentum of the bound states (R =1,8)	${\kappa }_{\left[𝐑\right]}\equiv \mu \left({\alpha }_{g,\left[𝐑\right]}^B+{\alpha }_h\right)$
Expectation value of relative velocity of interacting particles in the scattering state	${\text{v}}_{\mathrm{rel}},v_{\mathrm{rel}}=\left|{\text{v}}_{\mathrm{rel}}\right|$
Expectation value of momentum in the scattering state in the CM frame	$\text{k}\equiv \mu {\text{v}}_{\mathrm{rel}},k\equiv \left|\text{k}\right|$
Dimensionless variables that parametrise the long-range effects and the parts of the annihilation and BSF cross-sections that are evaluated numerically (colour representation R =1,8)	$d_h\equiv \mu {\alpha }_h/m_h$
	${\lambda }_{\left[𝐑\right]}\equiv {\alpha }_{g,\left[𝐑\right]}^B/{\alpha }_h$
	${\zeta }_h\equiv {\alpha }_h/v_{\mathrm{rel}}$
	${\zeta }_{g,\left[𝐑\right]}^S\equiv {\alpha }_{g,\left[𝐑\right]}^S/v_{\mathrm{rel}}$
	${\zeta }_{g,\left[𝐑\right]}^B\equiv {\alpha }_{g,\left[𝐑\right]}^B/v_{\mathrm{rel}}$
	${\zeta }_{S,\left[𝐑\right]}\equiv \left({\alpha }_h+{\alpha }_{g,\left[𝐑\right]}^S\right)/v_{\mathrm{rel}}$
	${\zeta }_{B,\left[𝐑\right]}\equiv \left({\alpha }_h+{\alpha }_{g,\left[𝐑\right]}^B\right)/v_{\mathrm{rel}}$
Binding energy of $n\ell m$ bound state	${ℰ}_{n\ell }=-{\gamma }_{n\ell }^2(\lambda ,d_h)\times {\kappa }^2/\left(2\mu \right)$
Kinetic energy of scattering state in CM frame	${ℰ}_{\text{k}}={\text{k}}^2/\left(2\mu \right)$
Wavefunction of $n\ell m$ bound state	${\psi }_{n\ell m}\left(\text{r}\right)={\kappa }^{3/2}\left[\frac{{\chi }_{n\ell }\left(\kappa r\right)}{\kappa r}\right]Y_{\ell m}\left({\Omega }_{\text{r}}\right)$
Wavefunction of scattering state	${\varphi }_{\text{k}}\left(\text{r}\right)=\sum _{\ell =0}^{\infty }(2\ell +1)[\frac{{\chi }_{\left|\text{k}\right|,\ell }\left(kr\right)}{kr}]P_{\ell }\left(\hat{\text{k}}\cdot \hat{\text{r}}\right)$
Dimensionless radial space coordinates for bound and scattering states	$x_B\equiv \kappa r$  and  $x_S\equiv kr$
Dimensionless time parameters for freeze-out	$x\equiv m_{\chi }/T,\tilde{x}\equiv m_X/T$

Table 2. Table 2: The momentum transfer Q 𝑄 Q at which the strong coupling α s ​ ( Q ) subscript 𝛼 𝑠 𝑄 \alpha_{s}(Q) is evaluated. For the bound states, the functions γ n ​ ℓ ​ ( λ , d h ) subscript 𝛾 𝑛 ℓ 𝜆 subscript 𝑑 ℎ \gamma_{n\ell}(\lambda,d_{h}) are computed numerically (cf.  sections   3.2 and  6 ). We recall that λ [ 𝐑 ] ≡ α g , [ 𝐑 ] / α h subscript 𝜆 delimited-[] 𝐑 subscript 𝛼 𝑔 delimited-[] 𝐑 subscript 𝛼 ℎ \lambda_{{\bf[R]}}\equiv\alpha_{g,{\bf[R]}}/\alpha_{h} with 𝐑 = 𝟏 , 𝟖 𝐑 1 8 {\bf R=1,8} , and d h = μ ​ α h / m h subscript 𝑑 ℎ 𝜇 subscript 𝛼 ℎ subscript 𝑚 ℎ d_{h}=\mu\alpha_{h}/m_{h} .

Vertices	${𝜶}_{𝒔}$	${𝜶}_{𝒈}$	Average momentum transfer $𝑸$
Annihilation: gluon emission	${\alpha }_s^{\mathrm{ann}}$		$m_X$
Scattering state wavefunction (ladder)	${\alpha }_s^S$	Colour-singlet ${\alpha }_{g,\left[\mathrm{𝟏}\right]}^S=\frac{4{\alpha }_s^S}{3}$	$\mu v_{\mathrm{rel}}$
		Colour-octet ${\alpha }_{g,\left[\mathrm{𝟖}\right]}^S=-\frac{{\alpha }_s^S}{6}$
Colour-singlet bound-state wavefunction (ladder)	${\alpha }_{s,\left[\mathrm{𝟏}\right]}^B$	${\alpha }_{g,\left[\mathrm{𝟏}\right]}^B=\frac{4{\alpha }_{s,\left[\mathrm{𝟏}\right]}^B}{3}$	$\begin{array}{c}{\kappa }_{\left[\mathrm{𝟏}\right]}{\gamma }_{n\ell }({\lambda }_{\left[\mathrm{𝟏}\right]},d_h)=\\ \mu \left({\alpha }_h+\frac{4{\alpha }_{s,\left[\mathrm{𝟏}\right]}^B}{3}\right)\times {\gamma }_{n\ell }(\frac{4{\alpha }_{s,\left[\mathrm{𝟏}\right]}^B}{3{\alpha }_h},d_h)\end{array}$ (4.10)
Colour-octet bound state wavefunction (ladder)	${\alpha }_{s,\left[\mathrm{𝟖}\right]}^B$	${\alpha }_{g,\left[\mathrm{𝟖}\right]}^B=-\frac{{\alpha }_{s,\left[\mathrm{𝟖}\right]}^B}{6}$	$\begin{array}{c}{\kappa }_{\left[\mathrm{𝟖}\right]}{\gamma }_{n\ell }({\lambda }_{\left[\mathrm{𝟖}\right]},d_h)=\\ \mu \left({\alpha }_h-\frac{{\alpha }_{s,\left[\mathrm{𝟖}\right]}^B}{6}\right)\times {\gamma }_{n\ell }(-\frac{{\alpha }_{s,\left[\mathrm{𝟖}\right]}^B}{6{\alpha }_h},d_h)\end{array}$ (4.11)
Formation of colour-singlet bound states: gluon emission	${\alpha }_{s,\left[\mathrm{𝟏}\right]}^{\mathrm{BSF}}$		$\frac{\mu }{2}\left[v_{\mathrm{rel}}^2+{\left({\alpha }_h+\frac{4{\alpha }_{s,\left[\mathrm{𝟏}\right]}^B}{3}\right)}^2{\gamma }_{n\ell }^2(\frac{4{\alpha }_{s,\left[\mathrm{𝟏}\right]}^B}{3{\alpha }_h},d_h)\right]$
Formation of colour-octet bound states: gluon emission	${\alpha }_{s,\left[\mathrm{𝟖}\right]}^{\mathrm{BSF}}$		$\frac{\mu }{2}\left[v_{\mathrm{rel}}^2+{\left({\alpha }_h-\frac{{\alpha }_{s,\left[\mathrm{𝟖}\right]}^B}{6}\right)}^2{\gamma }_{n\ell }^2(-\frac{{\alpha }_{s,\left[\mathrm{𝟖}\right]}^B}{6{\alpha }_h},d_h)\right]$

Equations

$$\Delta \equiv (m_X-m_{\chi })/m_{\chi },$$

$$M=2m_X\text{and}\mu =m_X/2.$$

$$\Gamma (X+\mathrm{SM}\leftrightarrow \chi +\mathrm{SM})\gg H,$$

$$\sigma (X+\mathrm{SM}\leftrightarrow \chi +\mathrm{SM})\gg \frac{17\text{leavevmode}{x}_{\mathrm{dec}}}{m_{\chi }{M}_{\mathrm{Pl}}}\sim 6\times 10^{-11}\text{leavevmode}\mathrm{pb}\left(\frac{\mathrm{TeV}}{m_{\chi }}\right),$$

$$\delta \mathcal{L}=(D_{\mu ,ij}{X}_j{)}^{†}(D_{i{j}^{\prime }}^{\mu }{X}_{j^{\prime }})-m_X^2{X}_j^{†}{X}_j+\frac{1}{2}({\partial }_{\mu }h)({\partial }^{\mu }h)-\frac{1}{2}{m}_h^2{h}^2-g_h{m}_{X}h{X}_j^{†}{X}_j,$$

$${\alpha }_s$$

$${\alpha }_h$$

$$\tilde{Y}\equiv Y_{\chi }+Y_X+Y_{X^{†}}=Y_{\chi }+2Y_X,$$

$$g_{\ast ,\mathrm{eff}}^{1/2}$$

$$\frac{d\tilde{Y}}{dx}=-\sqrt{\frac{\pi }{45}}\frac{m_{\mathrm{Pl}}{m}_{\chi }{g}_{\ast ,\mathrm{eff}}^{1/2}}{x^2}⟨{\sigma }_{\mathrm{eff}}{v}_{\mathrm{rel}}⟩({\tilde{Y}}^2-{\tilde{Y}}_{\mathrm{eq}}^2).$$

$$Y_{\chi }^{\mathrm{eq}}$$

$$Y_X^{\mathrm{eq}}=Y_{X^{†}}^{\mathrm{eq}}$$

$$⟨{\sigma }_{\mathrm{eff}}{v}_{\mathrm{rel}}⟩=\frac{2Y_X^{\mathrm{eq}}{Y}_{X^{†}}^{\mathrm{eq}}⟨{\sigma }_{X{X}^{†}}{v}_{\mathrm{rel}}⟩}{{\tilde{Y}}_{\mathrm{eq}}^2}=⟨{\sigma }_{X{X}^{†}}{v}_{\mathrm{rel}}⟩\left(\frac{2g_X^2(1+\Delta {)}^3{e}^{-2x\Delta }}{{\left[g_{\chi }+2g_X(1+\Delta {)}^{3/2}{e}^{-x\Delta }\right]}^2}\right),$$

$$⟨{\sigma }_{X{X}^{†}}{v}_{\mathrm{rel}}⟩\equiv ⟨{\sigma }_{\mathrm{ann}}{v}_{\mathrm{rel}}⟩+⟨{\sigma }_{\mathrm{BSF}}{v}_{\mathrm{rel}}{⟩}_{\mathrm{eff}}$$

$$V(r)=-\frac{{\alpha }_g}{r}-\frac{{\alpha }_h}{r}{e}^{-m_{h}r},$$

$${\alpha }_{g,[1]}$$

$${\alpha }_{g,[8]}$$

$$d_h\equiv \mu {\alpha }_h/m_h.$$

$${\zeta }_{g,[\mathbf{R}]}$$

$${\zeta }_h$$

$${\zeta }_{S,[\mathbf{R}]}$$

$${\zeta }_{B,[\mathbf{R}]}$$

$${\lambda }_{[\mathbf{R}]}$$

$$(\sigma v_{\mathrm{rel}}{)}_{X{X}^{†}\to gg}^{\mathrm{perturb}}$$

$$(\sigma v_{\mathrm{rel}}{)}_{X{X}^{†}\to hh}^{\mathrm{perturb}}$$

$$(\sigma v_{\mathrm{rel}}{)}_{X{X}^{†}\to gg}$$

$$(\sigma v_{\mathrm{rel}}{)}_{X{X}^{†}\to hh}$$

$$S_{0,[1]}$$

$$S_{0,[8]}$$

$$⟨{\sigma }_{\mathrm{ann}}{v}_{\mathrm{rel}}⟩={\left(\frac{\mu }{2\pi T}\right)}^{3/2}\int d^3{v}_{\mathrm{rel}}\exp \left(-\frac{\mu v_{\mathrm{rel}}^2}{2T}\right)\left[(\sigma v_{\mathrm{rel}}{)}_{X{X}^{†}\to gg}+(\sigma v_{\mathrm{rel}}{)}_{X{X}^{†}\to hh}\right].$$

$$\omega \simeq {\mathcal{E}}_{\mathbf{k}}-{\mathcal{E}}_{10}=\frac{\mu }{2}[v_{\mathrm{rel}}^2+({\alpha }_g^B+{\alpha }_h{)}^2{\gamma }_{10}^2].$$

$$\displaystyle(X+X^{\dagger})_{\bf[8]}\$$

$$\displaystyle(X+X^{\dagger})_{\bf[1]}\$$

$$\displaystyle(X+X^{\dagger})_{\bf[8]}\$$

$$\displaystyle(X+X^{\dagger})_{\bf[1]}\$$

Full text

††preprint:

TUM-HEP-1186-19

Nikhef-2019-004

††arxiv: 1901.10030

1 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 Griest:1990kh ; Ibarra:2015nca ; Freitas:2015hsa ; Baker:2015qna ; ElHedri:2018atj , 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 long-range force distorts the wavepackets of the interacting pairs, giving rise to the well-known Sommerfeld effect Sommerfeld:1931 ; Sakharov:1948yq that alters their annihilation rate Hisano:2002fk ; Hisano:2003ec . 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 vonHarling:2014kha ; Petraki:2015hla and contribute to the DM indirect detection signals Pospelov:2008jd ; An:2016gad ; An:2016kie ; Petraki:2016cnz ; Cirelli:2016rnw . While such effects emanating from hidden-sector mediators, the electroweak gauge bosons, or the gluons in co-annihilation scenarios, have been considered Hisano:2002fk ; Hisano:2003ec ; Pospelov:2008jd ; Beneke:2014hja ; Beneke:2014gja ; Hryczuk:2010zi ; Hryczuk:2011tq ; Harz:2014gaa ; vonHarling:2014kha ; Petraki:2015hla ; Cirelli:2015bda ; An:2016kie ; Petraki:2016cnz ; ElHedri:2016onc ; Liew:2016hqo ; Asadi:2016ybp ; An:2016gad ; Kouvaris:2016ltf ; Cirelli:2016rnw ; Kim:2016kxt ; Kim:2016zyy ; Biondini:2017ufr ; Biondini:2018pwp ; Biondini:2018xor ; Biondini:2018ovz ; Baldes:2017gzw ; Baldes:2017gzu ; Pearce:2015zca ; Ellis:2018jyl ; Geller:2018biy ; Harz:2018csl ; Cirelli:2018iax ; Bhattacharya:2018ooj ; Schmiemann:2019czm , only recently the long-range effect of the Higgs was pointed out Harz:2017dlj .

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 Harz:2017dlj 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. Harz:2017dlj , 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 Keung:2017kot ; Ibarra:2015nca ; Pierce:2017suq , with the effect of radiative bound-state formation (BSF) due to gluon exchange computed in Harz:2018csl . 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. Freitas:2015hsa ; Lopez-Honorez:2017ora ), 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 $\gamma$-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, $\Omega h^{2}=0.120\pm 0.001$, by the Planck satellite Aghanim:2018eyx 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 Belanger:2013oya ; Bringmann:2018lay ; Ambrogi:2018jqj 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% Harz:2012fz ; Harz:2014tma ; Harz:2016dql ; Baro:2007em ; Schmiemann:2019czm . The Sommerfeld effect and BSF due to the SM gauge bosons can have an even greater impact Harz:2018csl . As already shown in Harz:2017dlj 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 Harz:2017dlj , 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.

2 Simplified model and relic density

2.1 Simplified model

We assume that DM is a Majorana fermion $\chi$ with mass $m_{\chi}$, 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\leavevmode\nobreak\ {\rm GeV}$, that aims to resemble the SM Higgs boson. Both $\chi$ and $X$ are odd under a $\mathbb{Z}_{2}$ symmetry, and assumed to be the lightest (LP) and next-to-lightest particles (NLP). We denote the LP–NLP relative mass splitting as

[TABLE]

and define the total and the reduced mass of an $XX^{\dagger}$ pair,

[TABLE]

If $\Delta\ll 1$, the DM abundance is determined not only by the $\chi-\chi$ annihilations, but also by the $\chi-X$, $\chi-X^{\dagger}$ and $X-X^{\dagger}$ (co)annihilation processes. We shall assume that the $X-X^{\dagger}$ annihilation processes dominate the total depletion rate of the $\mathbb{Z}_{2}$-odd particles, and neglect all other contributions. This is a reasonable approximation, since the former involve the strong coupling and the coupling to the Higgs, which can be quite significant, while the latter involve the Weak coupling and are typically less efficient. Moreover, our goal here is to investigate, within a self-consistent framework, the long-range effect of the Higgs on the formation of bound states by particles charged under $SU(3)_{c}$. Including the $\chi-\chi$, $\chi-X$ and $\chi-X^{\dagger}$ processes in the computation of the DM density, would necessitate further model specifications and result in loss of generality. We note though, that a scenario where the depletion of DM is dominated by the annihilation of its coloured co-annihilating partners is realistic. For instance, within the MSSM it is possible to find viable stop-neutralino co-annihilation scenarios where, for small relative mass differences, the $\tilde{t}-\tilde{t}^{\dagger}$ annihilation contributes about $75\%-80\%$ to the neutralino depletion rate Harz:2014gaa ; Schmiemann:2019czm .

Implicit in our computations is the hypothesis that the LP and the NLP stay in chemical equilibrium, i.e. their interconversion rate exceeds the Hubble expansion rate,

[TABLE]

at least until the decoupling of the $X-X^{\dagger}$ annihilation processes. This implies a lower limit on the interconversion cross-section, $\sigma(X+\mathrm{SM}\leftrightarrow\chi+{\rm SM})\sim\Gamma(X+\mathrm{SM}\leftrightarrow\chi+{\rm SM})/n_{{\rm SM}}^{\mathrm{eq}}$, where $n_{{\rm SM}}^{\rm eq}$ stands for the number density of the SM particles participating in these processes. With the latter being relativistic during the DM freeze-out, $n_{{\rm SM}}^{\mathrm{eq}}\sim T^{3}$, and $H\sim 17\leavevmode\nobreak\ T^{2}/M_{\mathrm{Pl}}$, we find the condition

[TABLE]

where we set $x_{\rm dec}\equiv{m_{\chi}}/T_{\rm dec}\sim 100$ to encompass even $X-X^{\dagger}$ annihilations that occur somewhat after the DM freeze-out and can be significant in scenarios with light force mediators. Evidently, the required cross-section is many orders of magnitude below the Weak scale. How constraining the assumption of LP-NLP chemical equilibrium in reality is, depends on the model and the specific processes under consideration. For instance, in the MSSM, the interconversion between stops and neutralinos is mediated by SM particles, and the rate comfortably satisfies the condition (2.4). On the other hand, the gluino-neutralino interconversion is usually mediated by a heavy squark, and the condition of chemical equilibrium sets an upper limit on its mass, of the order of 100 TeV Nagata:2015hha . We discuss this and other caveats of our calculations further in section 5.3.

Given the above considerations, the Lagrangian terms relevant for our purposes are Harz:2017dlj

[TABLE]

where $D_{\mu,ij}=\delta_{ij}\partial_{\mu}+ig_{s}\,G_{\mu}^{a}T^{a}_{ij}$ is the covariant derivative with $G_{\mu}^{a}$ being the gluon fields and $T^{a}$ the corresponding generators. We define the couplings that will appear in the non-relativistic potential in section 2.3, as

[TABLE]

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}\supset|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}^{\dagger}X_{2}$, where $X_{1}$ and $X_{2}$ have different $SU(2)_{L}\times 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}\times 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}}\simeq 246\leavevmode\nobreak\ {\rm 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.

2.2 Boltzmann equation for the relic density

Following Edsjo:1997bg , we define $\tilde{Y}$ to be the sum over the yields of all (co)annihilating particles,

[TABLE]

where $Y_{i}$ is the ratio of the number density $n_{i}$ over the entropy density of the universe $s\equiv(2\pi^{2}/45)\,g_{*{S}}\,T^{3}$. We denote by $g_{*{S}}$ and $g_{*}$ the entropy and energy degrees of freedom respectively, and define

[TABLE]

Under the assumption of chemical equilibrium between the LP and NLP, the evolution of the DM density can be described by a single Boltzmann equation,

[TABLE]

As usual, we have defined the dimensionless time parameter $x\equiv m_{\chi}/T$. The yields in equilibrium of $\chi$, $X$ and $X^{\dagger}$ are ($g_{\chi}=2$, $g_{{X}}=3$)

[TABLE]

While $\langle\sigma_{\rm eff}\,v_{\rm rel}\rangle$ 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

[TABLE]

where

[TABLE]

includes all processes that deplete the $XX^{\dagger}$ 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.

2.3 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. Petraki:2015hla ). In our model, the $X,X^{\dagger}$ pairs interact via gluon and Higgs exchange, both of which may result in a sizeable long-range effect Harz:2017dlj . 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

[TABLE]

where $\alpha_{h}$ has been defined in eq. 2.7 and derived in ref. Petraki:2016cnz . The coupling $\alpha_{g}$ is related to $\alpha_{s}$ defined in eq. 2.6 but depends also on the colour representation of the $XX^{\dagger}$ state. The product of a colour triplet and an anti-triplet decomposes into a singlet and an octet state, ${\bf 3\otimes\bar{3}}={\bf 1\oplus 8}$, such that the coupling $\alpha_{g}$ amounts to Kats:2009bv

[TABLE]

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 eq. 2.15, we have differentiated between the strong coupling in the singlet and octet states, $\alpha_{s,{\bf[1]}}$ and $\alpha_{s,{\bf[8]}}$ respectively, to accommodate the possibly different momentum transfer along the 2PI diagrams, which affects the value of $\alpha_{s}$ due its running.

For the massive mediator $h$, the interaction manifests as long-range if the range of the potential $m_{h}^{-1}$ is comparable to or larger than the corresponding Bohr radius $(\mu\alpha_{h})^{-1}$. To parametrise the range of the interaction we use the dimensionless ratio

[TABLE]

The Sommerfeld effect and BSF become significant when the average momentum transfer between the scattering particles, $\mu v_{\rm rel}$, is comparable to or lower than the Bohr momentum, $\mu\alpha$. To compare the relevant scales in our set-up, we introduce the parameters

[TABLE]

where in eq. 2.17a, ${\bf 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 $\zeta_{g,{\bf[R]}}^{{S}}$ and $\zeta_{g,{\bf[R]}}^{{B}}$, 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 $\alpha_{h}$ depends on the underlying UV model, we neglect it in the following. For convenience, we define

[TABLE]

Finally, for the bound states, it will be useful to define the parameter

[TABLE]

The dimensionless variables defined in eqs. 2.16, 2.17, 2.18 and 2.19 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 $|\zeta_{g}|\gtrsim{\cal O}(1)$ and $\zeta_{h},d_{h}\gtrsim{\cal O}(1)$, for the Coulomb and Yukawa potentials, respectively. However, in ref. Harz:2017dlj 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.14) in several ways. On one hand, it generates thermal contributions to the mediator masses that screen the long-range effect. On the other hand, it gives rise to frequent scatterings between $X,X^{\dagger}$ and the relativistic plasma that allow for non-radiative dissipation of energy from an $XX^{\dagger}$ pair, and thus ensure equilibrium between bound and scattering states. In fact, the $XX^{\dagger}$ spectral function morphs into a continuum that encompasses energies both above and below $2m_{{X}}$ Kim:2016zyy ; Kim:2016kxt . These thermal effects have been considered in DM freeze-out computations using linear response theory, first in refs. Biondini:2017ufr ; Biondini:2018pwp ; Biondini:2018xor , where radiative capture processes were neglected, and subsequently in ref. Biondini:2018ovz which incorporated BSF via gluon emission. Recently, ref. Binder:2018znk 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^{\dagger}$ particles during DM freeze-out, which ensures that the average momentum transfer in an $XX^{\dagger}$ 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.

2.4 Annihilation

The dominant contributions to the direct annihilation of $XX^{\dagger}$ pairs arise from the $XX^{\dagger}\rightarrow gg$ and $XX^{\dagger}\rightarrow hh$ channels, shown in fig. 3. We neglect the $p$-wave suppressed contributions $XX^{\dagger}\rightarrow q\bar{q},gh$. (Similarly we neglect the $p$-wave suppressed contribution of the $s$-channel diagram of $XX^{\dagger}\rightarrow gg$). To leading order, the colour-averaged perturbative cross-sections are Harz:2017dlj

[TABLE]

where $\alpha_{s}^{\rm ann}\equiv\alpha_{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 octet scattering states contribute in eq. 2.20a, 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 Harz:2017dlj

[TABLE]

where $S_{0,{\bf[1]}}$ and $S_{0,{\bf[8]}}$ indicate the $s$-wave Sommerfeld factors for the colour-singlet and octet states respectively. For the mixed Coulomb and Yukawa potential (2.14), the $s$-wave Sommerfeld factor depends on the parameters defined in eqs. 2.17 and 2.16, i.e. $S_{0}=S_{0}(\zeta_{g}^{{S}},\zeta_{h},d_{h})$, as we discuss in section 3.1. Thus, taking eq. 2.15 into account,

[TABLE]

Here, $\alpha_{s}^{{S}}\equiv\alpha_{s}(Q=\mu v_{\rm rel})=\alpha_{s,{\bf[1]}}^{{S}}=\alpha_{s,{\bf[8]}}^{{S}}$ denotes the strong coupling evaluated at the average momentum transfer in the scattering states, $Q=\mu v_{\rm rel}$, which is independent of the colour representation. The thermally-averaged total annihilation cross-section is

[TABLE]

2.5 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, 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 ${\cal E}_{\bf k}={\bf k}^{2}/(2\mu)=\mu v_{\rm rel}^{2}/2$ minus the binding energy of the bound state ${\cal E}_{n\ell}=-\gamma^{2}_{n\ell}\times\kappa^{2}/(2\mu)$, where $\kappa\equiv\mu(\alpha_{g}^{{B}}+\alpha_{h})$ is the generalised Bohr momentum of the bound state and $\gamma_{n\ell}$ parametrises the departure of ${\cal E}_{n\ell}$ from its Coulomb limit $d_{h}\to\infty$ (cf. section 3.2). For capture into the ground state, the dissipated energy is

[TABLE]

The capture into bound states via emission of an ultrasoft gluon is depicted in figs. 5 and 5. The decompositions ${\bf 3\otimes\bar{3}}={\bf 1\oplus 8}$ and ${\bf 8}\otimes{\bf 8}={\bf 1}_{S}+{\bf 8}_{A}+{\bf 8}_{S}+{\bf 10}_{A}+\overline{\bf 10}_{A}+{\bf 27}_{S}$ imply that the allowed transitions are

[TABLE]

In (2.25b) and (2.25c), 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^{\dagger}$ combination supports a bound state, thus (2.25a) 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.25b) and (2.25c) 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 $\alpha_{s}$ running.

If the energy (2.24) available to be dissipated exceeds the Higgs mass, then bound states may form also via Higgs emission, according to

[TABLE]

However, the capture of particle-antiparticle pairs via scalar emission is subject to cancellations that suppress the processes (2.26) by higher orders in $\alpha_{s}$ and $\alpha_{h}$ with respect to (2.21) and (2.25) (Petraki:2016cnz, , 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.26) Oncala:2018bvl , 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

[TABLE]

where $f_{g}(\omega)=1/(e^{\omega/T}-1)$ is the gluon occupation number, with $\omega$ being the energy of the emitted gluon, given by eq. 2.24. The factor $1+f_{g}(\omega)$ 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\gtrsim\omega$, which encompasses a significant temperature range that is relevant to the DM freeze-out vonHarling:2014kha .111 As is standard, we have omitted the Bose enhancement factor in the thermal averaging of the annihilation cross-sections (2.23), since the energy of the annihilation products, $E=m_{{X}}$, far exceeds the temperature of the relativistic bath at the temperature range relevant for freeze-out, $T\lesssim m_{{X}}/20$.

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. (Harz:2018csl, , appendix D)),

[TABLE]

where ${\bf R}_{S}$, ${\bf R}$ denote the colour representations of the scattering and bound $XX^{\dagger}$ states, with $g_{{\bf[1]}}=1$ and $g_{{\bf[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.28, we have made explicit the dependence of the gluon energy (2.24) on the colour representation of the bound state. Note that for the ionisation cross-section of the colour-octet bound states – if they exist – we must include the contributions from both processes (2.25b) and (2.25c), i.e. sum over ${\bf R}_{S}={\bf 1,8}$ as indicated in eq. 2.28. The ionisation rate of a bound state is then

[TABLE]

where $\omega_{{\bf[R]}}^{\min}$ is the minimum energy required for ionisation, i.e. the binding energy, recovered from eq. 2.24 for $v_{\rm rel}=0$. Using eqs. 2.28 and 2.24, we arrive at

[TABLE]

Note that in eq. 2.29, $\gamma_{1,0}$ also depends on ${\bf R}$ even though it is not explicitly indicated. From eqs. 2.27 and 2.29, we can explicitly verify the principle of detailed balance,

[TABLE]

where $n_{{X}}^{{\rm eq}}$ and $n_{{{\cal B},\bf[R]}}^{{\rm eq}}$ are the equilibrium densities of $X$ and of $XX^{\dagger}$ bound states of colour representation ${\bf[R]}$, respectively. Note that eq. 2.30 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.29) and (2.30), that $\Gamma_{{\rm 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 $\ell=0$ bound states is related to the perturbative $s$-wave annihilation cross-section times relative velocity as (see e.g. Petraki:2015hla )

[TABLE]

where $\psi_{n00}^{{\bf[R]}}(0)$ indicates the $n$-th level $s$-wave bound-state wavefunction of the colour representation ${\bf R}$, evaluated at the origin, and will be discussed in section 3.2. The cross-section $\sigma_{{\rm ann},{\bf[R]}}^{s-\rm wave}v_{\rm 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^{\dagger}$ pair. Using the perturbative annihilation cross-sections (2.20) and taking into account the colour decomposition in eq. 2.21a, we obtain for the ground states

[TABLE]

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 vonHarling:2014kha . 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

[TABLE]

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 $\Gamma_{{\rm ion}}$ decreases exponentially with decreasing temperature.

The effective BSF cross-section (2.33) together with the thermally averaged annihilation cross-sections (2.21) yield the total effective annihilation cross-section (2.13) that determines the DM relic density according to the Boltzmann eq. 2.10.

Ionisation equilibrium.

If the BSF and ionisation cross-sections are sufficiently large to set $\Gamma_{{\rm ion}}\gg\Gamma_{{\rm dec}}$ at high temperatures, then the effective BSF cross-section (2.33) 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.30 and the decay rate (2.31), we find

[TABLE]

where we included the ground-state contributions only, and recall that ${\cal E}_{10}^{{\bf[R]}}$ are their binding energies. Equation 2.34 implies that at $T\leavevmode\nobreak\ \gg\leavevmode\nobreak\ |{\cal E}_{10}^{{\bf[R]}}|$, the bound states contribute to the Sommerfeld factor of the colour configuration ${\bf[R]}$, by

[TABLE]

To assess the importance of eq. 2.35, it is informative to consider a simple attractive Coulomb potential of strength $\alpha$, for which $|\psi_{1,0,0}(0)|^{2}=(\mu\alpha)^{3}/\pi$. Then, eq. 2.35 yields $\pi^{-5/2}[m_{{X}}\alpha^{2}/(4T)]^{3/2}\ll 1$, where the square bracket contains the binding-energy-to-temperature ratio, $|{\cal E}_{10}|/T$. This contribution can be compared with the thermally averaged Sommerfeld factor from the scattering states, which in this regime is $1+(4\pi)^{1/2}[m_{{X}}\alpha^{2}/(4T)]^{1/2}$. Clearly, the enhancement due to the bound states is subdominant to that due to the scattering states. Including excited states in eq. 2.34 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 vonHarling:2014kha ; this is in part manifested by the exponential factor in eq. 2.34. In this regime, the bound-state decay becomes comparable to or faster than ionisation, which is now exponentially suppressed, [cf. eq. 2.29], and $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle_{{\rm eff}}$ rapidly saturates to $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle$ [cf. eq. 2.33]. It is important to note that if $\sigma_{{\rm BSF}}$ is comparable to or exceeds $\sigma_{{\rm 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. vonHarling:2014kha ; Harz:2018csl and section 5).

3 Wavefunctions and overlap integrals

3.1 Scattering states

The non-relativistic potential due to gluon and Higgs exchange for the scattering states is

[TABLE]

where we will consider both positive and negative $\alpha_{g}^{{S}}$, but only $\alpha_{h}\geqslant 0$.

Wavefunctions.

The scattering states are described by a wavefunction $\phi_{\textbf{k}}(\textbf{r})$, parametrised by the continuous quantum number

[TABLE]

where k and $\textbf{v}_{\rm rel}$ are the expectation values of the momenta of the interacting particles in the CM frame, and of their relative velocity. The strong coupling $\alpha_{s}^{{S}}$ that determines $\alpha_{g}^{{S}}$ in eq. 3.1 is evaluated at $Q=|{\bf k}|$. The scattering state wavefunction obeys the Schrödinger equation

[TABLE]

with

[TABLE]

They are normalised according to

[TABLE]

To solve the Schrödinger equation numerically, we perform the separation of variables

[TABLE]

and introduce the dimensionless radial space coordinate

[TABLE]

such that the radial Schrödinger equation reads

[TABLE]

where we have used the parameters $\zeta_{g}^{{S}}\equiv\alpha_{g}^{{S}}/v_{\rm rel}$, $\zeta_{h}\equiv\alpha_{h}/v_{\rm rel}$ and $d_{h}\equiv\mu\alpha_{h}/m_{h}$ defined in eqs. 2.17 and 2.16. At $x_{{S}}\to 0$, and for $\ell>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_{{S}}^{\ell+1}$ (regular) and $x_{{S}}^{-\ell}$ (irregular). We are interested in the regular solutions, which imply the boundary condition

[TABLE]

The condition (3.9) will be valid also for $\ell=0$. To fully specify the wavefunction $\chi_{|\textbf{k}|,\ell}$, we shall also use the asymptotic behaviour at $x_{{S}}\to\infty$. At large $x_{{S}}$, the wavefunction $\chi_{|\textbf{k}|,\ell}$ behaves as (see e.g. ref. (Sakurai_QMbook, , chapter 7))

[TABLE]

where the phase shifts $\delta_{\ell}$ depend on $\zeta_{g}^{{S}}$, $\zeta_{h}$ and $d_{h}$. This implies that

[TABLE]

Coulomb limit.

In the limit $d_{h}\to\infty$, the Schrödinger eq. 3.8 depends only on $\zeta_{{S}}=\zeta_{g}^{{S}}+\zeta_{h}$, and can be solved analytically,

[TABLE]

Note that the limit $d_{h}\to 0$ (pure gluon-mediated potential) can be also obtained from the above with the substitution $\zeta_{{S}}\to\zeta_{g}^{{S}}$.

Sommerfeld factor $\bm{S_{0}}$ for $\bm{s}$-wave annihilation.

This is (see e.g. Cassel:2009wt )

[TABLE]

In the Coulomb limits $d_{h}\to 0$ and $d_{h}\to\infty$, the Sommerfeld factor reduces to

[TABLE]

with $S_{0}^{C}(\zeta)$ 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. Harz:2017dlj (see in particular figs. 2 and 3 therein).

3.2 Bound states

The non-relativistic potential for the bound states, due to gluon and Higgs exchange, is

[TABLE]

Note the difference in $\alpha_{g}$ with respect to eq. 3.1. While $\alpha_{g}^{{S}}$ in the scattering potential is to be evaluated at the scale $Q=\mu v_{\rm rel}$, $\alpha_{g}^{{B}}$ in the bound state potential is to be taken at the Bohr momentum scale, as will be discussed below.

Wavefunctions.

The Schrödinger equation for bound states is

[TABLE]

where $\{n\ell m\}$ are the standard principal and angular momentum quantum numbers, and the wavefunctions are normalised as

[TABLE]

We define the generalised Bohr momentum

[TABLE]

such that the discrete binding energy levels of the system are

[TABLE]

where the factor $\gamma_{n\ell}$ is determined numerically (see below) and parametrises the departure from the Coulomb limit $d_{h}\to\infty$ where $\gamma_{n\ell}=1/n$. We perform the separation of variables

[TABLE]

Defining the dimensionless space coordinate

[TABLE]

we find that $\chi_{n\ell}$ is normalised as

[TABLE]

and obeys the radial Schrödinger equation

[TABLE]

We recall that the dimensionless parameters $\lambda\equiv\alpha_{g}^{{B}}/\alpha_{h}$ and $d_{h}\equiv\mu\alpha_{h}/m_{h}$ have been defined in eqs. 2.16 and 2.19. It follows that the bound-state wavefunctions $\chi_{n\ell}$ and the (normalised) energy eigenvalues $\gamma_{n\ell}$ depend on $\lambda$ and $d_{h}$ only. Assuming $\alpha_{h}>0$, the necessary condition $\kappa>0$ for bound states to exist implies $\lambda>-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

[TABLE]

The discrete spectrum of eigenvalues $\gamma_{n\ell}=\gamma_{n\ell}(\lambda,d_{h})$ is then determined by requiring that $\chi_{n\ell}$ vanish at infinity,

[TABLE]

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}\to 0$ (pure gluon-mediated potential) and $d_{h}\to\infty$ (Coulomb limit of the Higgs-mediated potential), which allow for analytic solutions of eq. 3.23.

• •

In the limit $d_{h}\to\infty$, the spectrum of eigenvalues and the radial wavefunction are

[TABLE]

[TABLE]

where $L_{n}^{a}$ are the generalised Laguerre polynomials of degree $n$, and we assume the normalisation condition $\int_{0}^{\infty}z^{a}e^{-z}L_{n}^{(a)}(z)L_{m}^{(a)}(z)dz=[\Gamma(n+a+1)/n!]\,\delta_{n,m}$.

• •

In the limit $d_{h}\to 0$, bound states exist only for $\lambda>0$ ($\alpha_{g}^{{B}}>0$). In this case, the binding energy takes the Coulomb value for the gluon-generated potential only,

[TABLE]

where $\kappa_{g}\equiv\mu\alpha_{g}^{{B}}$. This implies

[TABLE]

Existence condition.

If the Coulomb component of the potential is not attractive, $\lambda\leqslant 0$ (equivalently, $\alpha_{g}^{{B}}\leqslant 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,\rm crit}(\lambda)$ such that for $d_{h}\geqslant d_{h,\rm crit}(\lambda)$ at least one bound level exists. We find that in the range $-0.96\leqslant\lambda\leqslant 0$, $d_{h,\rm crit}$ can be fitted by the analytical expression

[TABLE]

with an accuracy of better than $3\%$. The value $d_{0}=d_{h,\rm crit}(\lambda=0)\simeq 0.84$ corresponds to a Yukawa potential only and is in agreement with the results of ref. (Petraki:2016cnz, , appendix A2).

3.3 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 Harz:2018csl

[TABLE]

where we used the following Fourier transforms

[TABLE]

Note that the integrals (3.31) depend on the representations of the scattering and bound states. For simplicity we leave this dependence implicit in this section, but will denote it explicitly in section 4, as $\bm{{\cal J}}_{\textbf{k},\{n\ell m\}}^{{\bf[R_{\textsl{S}},R_{\textsl{B}}]}}$ and $\bm{{\cal Y}}_{\textbf{k},\{n\ell m\}}^{{\bf[R_{\textsl{S}},R_{\textsl{B}}]}}$, where ${\bf R_{\textsl{S}}}$ and ${\bf R_{\textsl{B}}}$ stand for the scattering-state and bound-state colour representations respectively.

For the numerical evaluation of the integrals (3.31), we express them in terms of the radial wavefunctions $\chi_{\textbf{k},\ell}$ and $\chi_{n\ell}$ introduced in sections 3.1 and 3.2. The integrals $\bm{{\cal J}}_{\textbf{k},\{n\ell m\}}$ have been previously analysed in ref. (Petraki:2016cnz, , appendix B). The leading order contribution to $\bm{{\cal J}}_{\textbf{k},\{n\ell m\}}(\textbf{b})$ with $\textbf{b}\propto\textbf{P}_{g}$ is independent of b Petraki:2015hla ; Petraki:2016cnz such that in the following we denote $\bm{{\cal J}}_{\textbf{k},\{n\ell m\}}=\bm{{\cal J}}_{\textbf{k},\{n\ell m\}}({\bf b}=0)$. Adapting the result to our notation for capture into the ground state, $\{n\ell m\}=\{100\}$, we arrive at222 Note that in ref. Petraki:2016cnz , a factor of $(-i)$ has been missed in going from eq. (B.8b) to eq. (B.9b). While this factor is not important for the computations of ref. Petraki:2016cnz , it is important here, since it determines the relative phase of $\bm{{\cal J}}_{\textbf{k},\{100\}}$ and $\bm{{\cal Y}}_{\textbf{k},\{100\}}$.

[TABLE]

We shall use eq. 3.33 to evaluate the BSF cross-sections in section 4.

Coulomb limit.

In the limit $d_{h}\to\infty$, the scattering and bound state wavefunctions, (3.11c) and (3.27) are

[TABLE]

where $\zeta_{{S}}\equiv\zeta_{g}^{{S}}+\zeta_{h}$ [cf. eq. 2.18a], $S_{0}^{C}(\zeta_{{S}})\equiv 2\pi\zeta_{{S}}/(1-e^{-2\pi\zeta_{{S}}})$ [cf. eq. 3.11b], and ${}_{1}F_{1}$ is the confluent hypergeometric function of the first kind. Following refs. Pearce:2013ola ; Petraki:2015hla ; Petraki:2016cnz ; Harz:2018csl , we derive the Coulomb limit of the overlap integrals using the identity AkhiezerMerenkov_sigmaHydrogen

[TABLE]

Equation 3.31 become (Harz:2018csl, , appendix B)

[TABLE]

and we used again $\kappa/k=\zeta_{{B}}$. As noted before, it suffices to evaluate $\bm{{\cal J}}_{\textbf{k},\{100\}}(\textbf{b})$ at $\textbf{b}=\textbf{0}$ Petraki:2015hla ; Petraki:2016cnz such that $\bm{{\cal Y}}_{\textbf{k},\{100\}}$ and $\bm{{\cal J}}_{\textbf{k},\{100\}}$ are related as

[TABLE]

with the squared amplitude of $\bm{{\cal J}}_{\textbf{k},\{100\}}^{C}$ being

[TABLE]

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.

4 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.

4.1 Amplitudes

The radiative capture via gluon emission has been computed in ref. Harz:2018csl 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:

[TABLE]

where $\alpha_{s}^{{\rm BSF}}\equiv\alpha_{s}(Q=|\textbf{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. $|\textbf{P}_{g}|$ is found from energy-momentum conservation to be $|\textbf{P}_{g}|={\cal E}_{\textbf{k}}-{\cal E}_{n\ell}=(\mu/2)[(\alpha_{h}+\alpha_{g}^{{B}})^{2}\gamma_{n\ell}^{2}+v_{\rm rel}^{2}]$. We recall that in eq. 4.1, the superscripts on the overlap integrals denote, in order, the scattering-state and bound-state colour representations.

4.2 Cross-sections for capture into the ground state

The differential cross-section for the radiative capture into bound states is

[TABLE]

For capture into the ground state, the leading-order amplitude is $\bm{{\cal M}}_{\textbf{k}\to\{100\}}\propto\textbf{k}$, as seen from the expressions (3.33) for the overlap integrals. Thus, carrying out the angular integration, we arrive at

[TABLE]

where we recall that $\zeta_{{B}}\equiv(\alpha_{h}+\alpha_{g}^{{B}})/v_{\rm rel}$ [cf. eq. 2.18b]. Using the amplitudes (4.1) and the overlap integrals (3.33), we find from eq. 4.3 the coloured-averaged radiative capture cross-sections to be

[TABLE]

We use eq. 4.4 to evaluate numerically the BSF cross-sections.

4.3 Coulomb limit

In the limit $d_{h}\to\infty$, 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

[TABLE]

We emphasise that, even when not explicitly denoted above, the couplings $\alpha_{s}^{{\rm BSF}}$, $\alpha_{s}^{{B}}$, $\alpha_{g}^{{B}}$, $\alpha_{g}^{{S}}$, and therefore also $\zeta_{{S}}\equiv(\alpha_{g}^{{S}}+\alpha_{h})/v_{\rm rel}$ and $\zeta_{{B}}\equiv(\alpha_{g}^{{B}}+\alpha_{h})/v_{\rm rel}$, depend on the colour representations of the initial and final states, and are thus different for every transition.

The function $S_{{\rm BSF}}^{C}(\zeta_{{S}},\zeta_{{B}})$ encapsulates all the velocity dependence of $\sigma_{{\rm BSF}}^{C}v_{\rm 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, $S_{1}^{C}(\zeta_{{S}})=[2\pi\zeta_{{S}}/(1-e^{-2\pi\zeta_{{S}}})](1+\zeta_{{S}}^{2})$. The factors inside the square brackets in eq. 4.5c arise from the convolution of the scattering and bound state wavefunctions with the radiative vertices. Comparing eqs. 4.4 and 4.5, yields $S_{{\rm BSF}}^{C}(\zeta_{{S}},\zeta_{{B}})=[(1+\zeta_{{B}}^{2})/(16f_{c})]\,\lim_{d_{h}\to\infty}{\cal H}$.

Two limiting cases of eq. 4.5 are of particular interest:

• •

For $\alpha_{h}\ll\alpha_{s}^{{B}}$, only colour-singlet bound states exist, hence the transition ${\bf[8]\to[1]}$ is the only possible, with $f_{c}=289/64$. The corresponding cross-section is

[TABLE]

• •

For $\alpha_{h}\gg\alpha_{s}^{{B}}$, the colour factors (4.5b) reduce to

[TABLE]

4.4 The running of $\bm{\alpha_{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$.333 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. Since the strong coupling is sensitive to the momentum transfer, $\alpha_{s}=\alpha_{s}(Q)$, it is important to carefully account for its running. An overview of the various scales is given in table 2.444 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 Pineda:2011aw ; Pineda:2001ra ; Manohar:2000kr for the computation of the $\alpha_{s}$ running in NRQCD and pNRQCD. In the occasions that the momentum transfer itself depends on the strong coupling, $Q=Q(\alpha_{s})$, we solve numerically the following equation for $\tilde{\alpha}$,

[TABLE]

The effect of the running of $\alpha_{s}$ is depicted in figs. 7, 9, 9 and 10, with related discussion in the captions. A few remarks are in order here.

Scattering states.

The average momentum transfer $Q=\mu v_{\rm rel}$, implies that $\alpha_{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.21a), this enhances significantly the total annihilation rate, even for moderate values of $\alpha_{h}$.555 This holds provided that $v_{\rm rel}$ is not too low. See related discussion in section 5.1. This is depicted in fig. 7.

Bound states.

Because the momentum transfer depends on the strong coupling, we solve numerically eq. 4.8. Generally, $\alpha_{h}$ increases the average momentum transfer, thus suppressing $\alpha_{s}^{{B}}$, as shown in fig. 9. Overall though, the Bohr momentum $\kappa\equiv\mu(\alpha_{g}^{{B}}+\alpha_{h})$ increases with increasing $\alpha_{h}$. If $\alpha_{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,

[TABLE]

where we recall that $d_{h}=\mu\alpha_{h}/m_{h}$, and the function $\gamma_{1,0}(\lambda,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 $\alpha_{h}$ intersect the horizontal black line. Let us examine the scaling of $\alpha_{s}(Q(\tilde{\alpha}))/\tilde{\alpha}$ with $\tilde{\alpha}$.

• •

For sufficiently small $\tilde{\alpha}$ ($\tilde{\alpha}\ll 6\alpha_{h}$), the momentum transfer (4.9) is independent of $\tilde{\alpha}$, i.e. $Q\simeq\mu\alpha_{h}\times\gamma_{1,0}(0,d_{h})$. Since $\alpha_{s}\sim 0.1$ for the $m_{{X}}$ values of interest, the ratio $\alpha_{s}(Q(\tilde{\alpha}))/\tilde{\alpha}$ is simply $0.1/\tilde{\alpha}$, i.e. starts from values $\gg 1$ and decreases with $\tilde{\alpha}$.

• •

At larger $\tilde{\alpha}$ (but still $\tilde{\alpha}<6\alpha_{h}$), this scaling changes. The factor $\gamma_{1,0}$ tends to 0, i.e. $d_{h}$ approaches $d_{h,\,\rm crit}(\lambda=-\tilde{\alpha}/(6\alpha_{h}))$ [cf. sections 3.2 and 6]. This, in turn, drives $Q$ to zero and consequently $\alpha_{s}(Q)$ to infinity.

• •

Between these two limits, $\alpha_{s}(Q(\tilde{\alpha}))/\tilde{\alpha}$ reaches the value 1 if $\alpha_{h}$ is sufficiently large.

The minimum value of the Higgs coupling, $\alpha_{h,\,\rm 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 $\alpha_{s}$ (for fixed values of the other parameters). Since $\alpha_{h}$ increases the binding energy, it suppresses $\alpha_{s}^{{\rm BSF}}$. This is seen in fig. 10.

5 The long-range effect of the Higgs

5.1 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. Harz:2017dlj . 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}\equiv\mu\alpha_{h}/m_{h}$) than in the absence of the attractive gluon-mediated force (Harz:2017dlj, , fig. 2).

• •

Conversely, in the octet state, a larger $d_{h}$ is required for the Higgs to affect the long-range interaction. However, the gluon-mediated repulsion is suppressed by a colour factor with respect to the full strength of the strong coupling, $\alpha_{g,{\bf[8]}}=-\alpha_{s,{\bf[8]}}/6$, suggesting that even a modest $\alpha_{h}$ can potentially counteract it. This holds so long as $v_{\rm rel}$ is not too low; at low $v_{\rm rel}$, the gluon-induced repulsion becomes exponential [cf. the Sommerfeld factor in the Coulomb limit (3.11b) with $\zeta=\zeta_{g,{\bf[8]}}^{{S}}\lesssim-1$], and cannot be surmounted by a finite-range attractive force (cf. ref. (Harz:2017dlj, , fig. 3, left panel)). However, the DM freeze-out occurs while $v_{\rm 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_{{{\rm FO}}}\approx m_{{X}}/30$. If the binding energy equals or exceeds $T_{{{\rm 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 $\alpha_{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 $\alpha_{h}\gtrsim 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 $\alpha_{s}$, as discussed in section 4.4 and depicted in figs. 7, 9 and 10. However, the Bohr momentum and the binding energy increase with $\alpha_{h}$, as seen in fig. 11.

Cross-sections.

In fig. 13, we illustrate the velocity dependence of the BSF cross-sections, and compare them with the direct annihilation processes.

• •

For large $v_{\rm 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=\mu v_{\rm rel}$, is much greater than that in the bound states, which in the Coulomb approximation is the Bohr momentum, $\kappa=\mu(\alpha_{g}^{{B}}+\alpha_{h})$. In this regime, the BSF cross-sections scale as $\sigma_{{\rm BSF}}v_{\rm rel}\propto(\kappa/k)^{4}\approx[(\alpha_{g}^{{B}}+\alpha_{h})/v_{\rm rel}]^{4}$.

• •

At lower velocities, when $k\sim\kappa$ or $v_{\rm rel}\sim\alpha_{g}^{{B}}+\alpha_{h}$, the overlap of the wavefunctions is nearly maximal (although the precise value of $v_{\rm rel}$ at which this occurs depends also on the Bohr momentum of the scattering state, hence on $\alpha_{g}^{{S}}+\alpha_{h}$). In this regime, the cross-section for capture into the tightest bound state (colour singlet) dominates over annihilation. Note that since $\alpha_{s}\sim 0.1$, this velocity range is relevant for freeze-out.

• •

At $v_{\rm rel}\lesssim\alpha_{g}^{{S}}+\alpha_{h}$, the Sommerfeld effect becomes important.

Because of the gluon-mediated repulsion in the octet scattering states, the ${\bf[8]\to[1]}$ and ${\bf[8]\to[8]}$ capture processes become suppressed with decreasing $v_{\rm rel}$. The coupling to the Higgs ameliorates this suppression, and large $\alpha_{h}$ can even make these BSF cross-sections increase temporarily as the velocity drops; this can be seen in the right panel of fig. 13. However, at sufficiently low $v_{\rm rel}$, the exponential Coulomb repulsion cannot be overcome by the finite-range Higgs-mediated attraction.

On the other hand, the ${\bf[1]\to[8]}$ cross-section increases steadily with decreasing $v_{\rm 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 colour-singlet contribution and increases with decreasing $v_{\rm rel}$.

Note that in fig. 13, 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 $\alpha_{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 $\alpha_{h}$. This is both due to the $XX^{\dagger}\to 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 $\alpha_{h}$. The suppression of $\alpha_{s}^{{\rm BSF}}$ at larger $\alpha_{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_{\rm 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 $\langle\sigma_{{\rm ann}}v_{\rm rel}\rangle$ follow from the discussion above. For BSF, we present both the actual and the effective thermally averaged cross-sections, $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle$ and $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle_{{\rm eff}}$. Several comments are in order.

• •

Despite the suppression of $\sigma_{{\rm BSF}}v_{\rm rel}$ at large $v_{\rm rel}$ discussed above, $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle$ appears to be large even at high $T$; this is due to the Bose-enhancement factor [cf. eq. 2.27].

• •

Because of the rapid ionisation processes, $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle_{{\rm eff}}$ is suppressed at high $T$, and in fact becomes independent of $\sigma_{{\rm BSF}}v_{\rm rel}$, as discussed in section 2.5. However, $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle_{{\rm eff}}$ increases as $T$ drops.

For capture from an octet state, $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle_{{\rm eff}}$ peaks at $T\sim$ binding energy. The reason is three-fold: (i) This condition is the thermally-averaged equivalent of $k\sim\kappa$, 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 $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle_{{\rm eff}}$ saturates to $\langle\sigma_{{\rm BSF}}v_{\rm rel}\rangle$. (iii) At even lower temperatures, the gluon repulsion in the octet state suppresses $\langle\sigma_{{\rm BSF}}^{{\bf[8]\to[1]}}v_{\rm rel}\rangle$ and $\langle\sigma_{{\rm BSF}}^{{\bf[8]\to[8]}}v_{\rm rel}\rangle$.

• •

At least for $0.02\lesssim\alpha_{h}\lesssim 0.1$, $\langle\sigma_{{\rm BSF}}^{{\bf[8]\to[1]}}v_{\rm rel}\rangle_{{\rm eff}}$ overcomes $\langle\sigma_{{\rm ann}}v_{\rm rel}\rangle$, even before it saturates to $\langle\sigma_{{\rm BSF}}^{{\bf[8]\to[1]}}v_{\rm rel}\rangle$, 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 $\langle\sigma_{{XX^{\dagger}}}v_{\rm rel}\rangle_{{\rm eff}}$ for the range of couplings considered.

5.2 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 $\alpha_{h}$, the NLP-LP mass difference $\Delta m$ versus the DM mass $m_{\chi}$ 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 non-perturbative 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 $\alpha_{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\leavevmode\nobreak\ \mathrm{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 $\alpha_{h}$.

With increasing $\alpha_{h}$, the effect of BSF with gluon exchange only appears to become less significant. This is because it has remained unaffected by $\alpha_{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 $\alpha_{h}=0.1$ and $\alpha_{h}=0.2$, $\Delta 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.

5.3 Caveats

The present work along with refs. Harz:2017dlj ; Harz:2018csl 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.8, 2.9, 2.10, 2.11, 2.12 and 2.13. This in turn would shift the contours of figs. 14 and 15 to lower values of $\Delta m$ and $m_{\chi}$.

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 $\alpha_{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 $\alpha_{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. Binder:2018znk , 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 Giudice:1998dj .

6 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 Harz:2017dlj , and by the formation and decay of bound states due to gluon exchange Harz:2018csl . 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 $\alpha_{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 mass — which 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.

Acknowledgements

We thank Tobias Binder, Simone Biondini, Sacha Davidson, Adam Falkowski, Mikko Laine, Pasquale Serpico, Pietro Slavich and Michel Tytgat for useful discussions. J.H. was supported by the Labex ILP (reference ANR-10-LABX-63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche (ANR), as part of the programme Investissements d’avenir under the reference ANR-11-IDEX-0004-02. J.H. was further supported by the DFG Emmy Noether Grant No. HA 8555/1-1. K.P. was supported by the ANR ACHN 2015 grant (“TheIntricateDark” project), and by the NWO Vidi grant “Self-interacting asymmetric dark matter”.