TeV-Scale Thermal WIMPs: Unitarity and its Consequences

TL;DR

This paper revisits unitarity bounds on thermal WIMP dark matter, showing that bound state effects and Sommerfeld enhancement significantly lower the maximum allowed mass, impacting detection prospects.

Contribution

It introduces revised unitarity bounds considering bound states and Sommerfeld effects, extending the understanding of thermal WIMP mass limits.

Findings

Unitarity bound lowered from 139 TeV to below 100 TeV for non-self-conjugate dark matter.

Largest allowed mass for composite dark matter is 1 PeV.

Bound-state formation affects late-universe annihilation signals.

Abstract

We re-examine unitarity bounds on the annihilation cross section of thermal-WIMP dark matter. For high-mass pointlike dark matter, it is generic to form WIMP bound states, which, together with Sommerfeld enhancement, affects the relic abundance. We show that these effects lower the unitarity bound from 139 TeV to below 100 TeV for non-self-conjugate dark matter and from 195 TeV (the oft-quoted value of 340 TeV assumes $\Omega_{DM} h^2 = 1$) to 140 TeV for the self-conjugate case. For composite dark matter, for which the unitarity limit on the radius was thought to be mass-independent, we show that the largest allowed mass is 1 PeV. In addition, we find important new effects for annihilation in the late universe. For example, while the production of high-energy light fermions in WIMP annihilation is suppressed by helicity, we show that bound-state formation changes this. Coupled with rapidly improving experimental sensitivity to TeV-range gamma rays, cosmic rays, and neutrinos, our results give new hope to attack the thermal-WIMP mass range from the high-mass end.

Figures

Tables

Table 1. Table 1: Relative weight factors of the Sommerfeld effect for the leading channels.

Equations

$$\frac{dY}{dz}=-\frac{\lambda ⟨{\sigma }_{\mathrm{eff}}{v}_{\mathrm{rel}}⟩(z)}{z^2}(Y^2-Y_{\mathrm{eq}}^2),$$

$$⟨{\sigma }_{\mathrm{eff}}v⟩=⟨{\sigma }_{\mathrm{ann}}{v}_{\mathrm{rel}}⟩+\sum _{I=1}^{\infty }⟨{\sigma }_{\mathrm{I}}{v}_{\mathrm{rel}}⟩BR({\mathrm{B}}_{\mathrm{I}}\to \mathrm{SM}),$$

$$\text{with }BR({\mathrm{B}}_{\mathrm{i}}\to \mathrm{SM})=\frac{{\mathrm{\Gamma }}_{\mathrm{ann}}}{{\mathrm{\Gamma }}_{\mathrm{ann}}+{\mathrm{\Gamma }}_{\mathrm{break}}}$$

$$={\left[1+\frac{⟨{\sigma }_{\mathrm{I}}{v}_{\mathrm{rel}}⟩g_{\chi }^2{M}_{\text{DM}}^3{e}^{-z{E}_{B_I}/M_{\text{DM}}}}{2g_I(4\pi z{)}^{3/2}{\Gamma }_{\mathrm{ann}}}\right]}^{-1},$$

$$Y_{\mathrm{DM}}(\infty )=\frac{1}{\lambda }{\left({\int }_{z_f}^{\infty }\frac{⟨{\sigma }_{\mathrm{eff}}{v}_{\mathrm{rel}}⟩(z)}{z^2}dz+\frac{⟨{\sigma }_{\mathrm{eff}}{v}_{\mathrm{rel}}⟩(z_f)}{z_f^2}\right)}^{-1},$$

$$z_f$$

$$(\sigma v_{\mathrm{rel}}{)}_{\mathrm{total}}^J<(\sigma v{)}_{\max }^J=\frac{4\pi (2J+1)}{M_{\text{DM}}^2{v}_{\mathrm{rel}}}.$$

$$(\sigma {)}_{\mathrm{total}}<\sum _{\ell =0}^{{\ell }_{\text{class}}}\frac{4\pi (2\ell +1)}{p_{\mathrm{DM}}^2}\approx 4\pi R_{\mathrm{DM}}^2,$$

$$\frac{\sigma v_{\mathrm{ann}}}{{\sigma }_0}=\frac{2R^4+17R^2-19}{32R},$$

$$\frac{{\sigma }_{\mathrm{S},\mathrm{ann}}}{{\sigma }_{\mathrm{ann}}}=\sum _I{f}_I{S}_{\mathrm{ann}}({\lambda }_I)\approx \sum _{I=1,3}{f}_I{S}_{\mathrm{ann}}({\lambda }_I),$$

$$(\sigma v_{\mathrm{rel}}{)}_{\ell }^n=(\sigma v_{\mathrm{rel}}{)}_{\ell }^n({\lambda }_i,{\lambda }_f)\times {\left|C_J+{\gamma }_{\ell }^n{C}_{\tau }\right|}^2.$$

$$\frac{(\sigma v_{\mathrm{rel}}{)}_{\mathrm{BSF}}^{1,0}}{\left(\frac{{\sigma }_0(2S+1)2^{11}\pi }{3g_{\chi }^2}\right)}\approx {\left|C_J+\frac{1}{{\lambda }_f}{C}_{\tau }\right|}^2\times \frac{{\lambda }_i^3{\alpha }_2}{{\lambda }_f{v}_{\mathrm{rel}}}{e}^{-4{\lambda }_i/{\lambda }_f},$$

$$C_J(1\to 3)=C_J(3\to 1)=\sqrt{\frac{R^2-1}{4}},$$

$$C_{\tau }(1\to 3)=-C_{\tau }(3\to 1)=-\sqrt{\frac{R^2-1}{4}},$$

$$C_J(3\to 5)=C_J(5\to 3)=\sqrt{\frac{R^2-4}{2}},$$

$$C_{\tau }(5\to 3)=-C_{\tau }(3\to 5)=\sqrt{2\left(R^2-1\right)}.$$

$$\frac{{\Gamma }_{1S=0}^{n\ell =0}}{M_{\text{DM}}}=\frac{{\alpha }_2^5{C}_2(R{)}^3{T}_R^2{d}_{\text{adj}}}{2\eta n^3{d}_R}=\frac{{\alpha }_2^{5}R{\left(R^2-1\right)}^5}{3\eta 2^{11}{n}^3},$$

$$\frac{{\Gamma }_{3S=1}^{n\ell =0}}{M_{\text{DM}}}=\frac{{\alpha }_2^5{\lambda }_3^3{n}_f{T}_R{T}_{\text{SM}}}{3\eta 2^2{n}^3}=\frac{n_f{\alpha }_2^{5}R\left(R^2-1\right){\left(R^2-5\right)}^3}{9\eta 2^{11}{n}^3}.$$

$$(\sigma v{)}_{\mathrm{total}}^0=(f_{\mathrm{ann}}^0{)}^2({\sigma }_{\mathrm{ann}}v)+\sum _I(f_I^0{)}^2({\sigma }_{\mathrm{BSF}}^{I}v).$$

$$({\sigma }_{n=1,\ell =0}^S{v}_{\mathrm{rel}})=\frac{32^4\left(2S+1\right){\alpha }_D^2}{\pi g_{\chi }^2}{\int }_0^{\Delta E}dk{k}^3{\left(\Delta E-k\right)}^3$$

$$\times \sum _{n=2}^{\infty }{\left|⟨{\varphi }_{\ell =0}|r|{\psi }_{n,1}⟩⟨{\psi }_{n,1}|r|{\psi }_{1,0}⟩\right|}^2$$

$${\left\{\frac{1}{-E_{1,0}+E_{n,1}-k}+\frac{1}{-E_{1,0}+E_{n,1}-(\Delta E-k)}\right\}}^2$$

$$\approx \frac{\tilde{\kappa }\left(2S+1\right){\bar{\sigma }}_0{\alpha }_D^3}{g_{\chi }^2}\frac{{\alpha }_D}{v_{\mathrm{rel}}},$$

$$\frac{E_{1s_3}^{\gamma }}{M_{\text{DM}}}=\frac{{\lambda }_3^2{\alpha }_2^2}{4}=\frac{{\alpha }_2^2(R^2-5{)}^2}{64},$$

$$\frac{E_{2p_3}^{\gamma }}{M_{\text{DM}}}=\frac{{\lambda }_3^2{\alpha }_2^2}{16}=\frac{{\alpha }_2^2(R^2-5{)}^2}{256},$$

$$\frac{E_{2p_3-1s_1}^{\gamma }}{M_{\text{DM}}}=\frac{{\alpha }_2^2}{256}\left(3R^4+2R^2-21\right).$$

$$n_{\max }={\ell }_{\max }+1={\left(\frac{{\alpha }_{\mathrm{eff}}^3}{12{\lambda }_R}\right)}^{1/4}\sqrt{M_{\mathcal{Q}}/\Lambda }$$

$$=\frac{{\alpha }^{3/4}{\lambda }_R^{1/2}}{12^{1/4}}\sqrt{M_{\mathcal{Q}}/\Lambda }.$$

$${\ell }_{\mathrm{class}}=p_Q\times R=v_{\mathrm{rel}}{M}_{\mathcal{Q}}/\Lambda \approx 4\sqrt{M_{\mathcal{Q}}T}/\Lambda ,$$

$${\sigma }_R^{\mathrm{ann}}=p_{R}4\pi \sum _{\ell =0}^{L_{\max }(T)}\frac{(2\ell +1)\sin ({\delta }_{\ell }{)}^2}{M_{\mathcal{Q}}^2{v}_{\mathrm{rel}}^2}{\mathbb{P}}_{\ell }^{\mathrm{ann}},$$

$$(\sigma^{\rm ann})\leq\frac{4\pi\,\ell_{\rm max}(T)^{2}}{M_{\mathcal{Q}}^{2}v_{\rm rel}^{2}}=\sigma_{\rm geom.}\left\{\begin{array}[]{ll}4T_{c}/M_{\mathcal{Q}}v_{\rm rel}^{2}&\mathrm{;}\,T>T_{c},\\ 1&\mathrm{;}\,T<T_{c}\,,\end{array}\right.$$

$$\langle\sigma^{\rm ann}v_{\rm rel}\rangle\leq\sigma_{\rm geom.}\,\sqrt{\frac{16T}{\pi\,M_{\mathcal{Q}}}}\left\{\begin{array}[]{ll}T_{c}/T&\mathrm{for}\quad T>T_{c},\\ 1&\mathrm{for}\quad T<T_{c}\,.\end{array}\right.$$

$$Y_B(\infty )=\sqrt{\frac{45}{g_{\mathrm{SM}}\pi }}\frac{1}{M_{\mathcal{Q}}{M}_{\mathrm{Pl}}}$$

$${\left(\frac{1}{Y_B(z_{\Lambda })}+{\int }_{z=M_{\mathcal{Q}}/\Lambda }^{\infty }\frac{⟨{\sigma }^{\mathrm{ann}}{v}_{\mathrm{rel}}⟩}{z^2}dz\right)}^{-1}.$$

$$Y_B(\infty )=\frac{9\sqrt{5M_{\mathcal{Q}}/g_{\mathrm{SM}}}}{8M_{\mathrm{Pl}}{\sigma }_{geom.}\left(3\sqrt{\Lambda /T_c}-2\right)T_c^{3/2}}$$

$$\approx \frac{3\sqrt{5/g_{\mathrm{SM}}}}{8M_{\mathrm{Pl}}{\sigma }_{geom.}{T}_c}\sqrt{\frac{M_{\mathcal{Q}}}{\Lambda }}\approx \frac{7.610^{-21}}{\mathrm{GeV}{\lambda }_{\mathrm{R}}}\sqrt{\frac{M_{\mathcal{Q}}\Lambda }{{\alpha }^3}}.$$

Full text

I Introduction

The unknown particle nature of dark matter has inspired a plethora of imaginative models Ge:2019voa ; Berges:2019dgr ; Blennow:2019fhy ; Arcadi:2019lka ; Cirelli:2018iax ; 1705.05567 . Among them, one well-motivated model is unique in its simplicity and specificity, and that is a thermal-relic WIMP that annihilates to Standard-Model particles 0404175 ; Steigman:1984ac ; Arcadi:2017kky ; Roszkowski:2017nbc . While this may not be the correct description of nature, it is essential that this hypothesis be fully tested. This is challenging but possible.

In this model, the early universe annihilation rate factor is determined from the dark matter relic abundance as $\langle\sigma v\rangle=(2.2\times 10^{-26}\,\text{cm}^{3}\,\text{s}^{-1})(0.12/\Omega_{DM}h^{2})$, where this is the total cross section to all final states SteigmanDasguptaBeacom . (We quote the value at large dark-matter masses; at smaller masses, it is larger.) If annihilation proceeds through $s$-wave scattering, as is well motivated, then the late-universe annihilation rate factor is the same. Given the density distribution of dark matter, determined by gravitational probes Aghanim:2018eyx , upper limits on the fluxes of energetic particles (gamma rays, cosmic rays, and neutrinos) then determine upper limits on $\langle\sigma v\rangle_{i}$, where $i$ denotes particular final states. Importantly, in a generic thermal-WIMP model, only the total annihilation rate factor is model-independent; those to particular final states — as well as production rates at colliders DM@LHC ; Baltz:2006fm ; Han:2018wus and elastic-scattering rates in underground detectors Goodman:1984dc ; Cushman:2013zza ; Schumann:2019eaa — are model-dependent.

Fully testing the thermal-WIMP hypothesis requires reaching a cross section sensitivity below this prediction, and doing so for all WIMP masses. Moreover, though searches for annihilation products test the partial cross sections to particular final states, one must combine those to determine the upper limit on the cross section to all Standard-Model final states UpperBoundCross ; LeaneSlatyerBeacom . If one excludes neutrinos from being dominant final states, then the lower limit is $\simeq 20$ GeV LeaneSlatyerBeacom , to be contrasted with the limits claimed assuming favorable final states, which approach $\simeq 100$ GeV FermiDwarf . (WIMPs that annihilate only to neutrinos are excluded below $\sim 10$ MeV by big bang nucleosynthesis NollettSteigman1 ; NollettSteigman2 ; if there is any annihilation to other channels, the mass limit from the cosmic microwave background is much larger LeaneSlatyerBeacom ; Galli:2013dna ; Padmanabhan:2005es .) It is well known that improvements in the sensitivity of existing searches are possible, and that, if no signals are found, the thermal-relic mass range will be progressively attacked from the low-mass end 1408.4131 .

But what about the high-mass end? The largest allowed mass of a thermal WIMP is then determined only by the theoretical bound from $s$-wave unitarity, and is $\simeq 340$ TeV (for $\Omega_{DM}h^{2}=1$; see below for $\Omega_{DM}h^{2}<1$) for pointlike dark matter GriestKamionkowski . For composite dark matter, the lower limit on the radius is $7.5\times 10^{-7}\,\rm fm$, independent of mass GriestKamionkowski . In 1990, when these limits were set, the experimental sensitivity to high-mass dark matter annihilation was vastly inadequate. Today, it is much better but still inadequate, though that will change due to new generations of experiments and better understandings of astrophysical backgrounds. This prompts a re-examination of thermal WIMPs at the largest masses. New theoretical developments Nima ; Cassel ; Petraki ; Wise ; CosmoBS make this especially relevant.

Saturating unitarity requires a large cross section. In the minimal scenario, we consider, where the only new particle is the WIMP, the options to achieve this are limited. (For example, it is not possible to invoke resonances Pavel ; DiracDM , as these require additional particles to be the mediators.) A key insight is that when the WIMP mass is very large compared to all Standard-Model particles, there are only light mediators, which induce long-range potentials and enhance cross sections through the Sommerfeld effect Nima ; Cassel ; Hisano:2003ec ; Hisano:2004ds ; Hisano:2006nn ; 0812.0559 . Including near threshold effects the Sommerfeld effect was further investigated in 1609.00474 ; 1706.01894 . (The dark matter itself can not serve as a mediator, as this topology would unavoidably open a decay channel to Standard Model particles.) Further, based on more recent work, there are also generically bound-state effects Petraki ; CosmoBS . Only by taking all of these effects into account, as we do in this paper, can accurate results be obtained. Further, these effects change the prospects for the detection of late-universe annihilation products.

This paper is structured as follows. In Sec. II, we review bound-state effects on WIMP freeze-out and how these affect the unitarity bound. In Sec. III, we calculate the bound-state effects on the unitarity bound for pointlike dark matter and calculate the largest allowed WIMP masses. In Sec. IV, we extend our calculations to composite dark-matter states of finite size. In Sec. V, we summarize our results, emphasizing the path towards experimental sensitivity at high masses overtaking the unitarity bounds for the first time. The ultimate goal is to test thermal WIMPs over the full mass range by attacking from both the low-mass and high-mass ends.

II Overview of Bound-State Effects for Heavy Dark Matter

In this section, we briefly review previous studies of bound-state effects on annihilation, define in what way the present framework is different, and note the implications for freeze-out.

II.1 Bound-State Effects on Freezeout

Bound-state effects on dark matter (DM) annihilation have been discussed only recently, and only in certain cases. In the context of a $U(1)$ model, they were found in Ref. Petraki to be important, though Ref. Wise argued the opposite. A follow-up study PetrakiNew discusses bound-state formation as a possibility to reach the unitarity bound in a perturbative abelian model, with a new hypothetical force carrier and a coupling strength close to its non-perturbative value. Bound-state effects in a weakly interacting DM system were studied in Refs. Slatyer ; 0901.2125 ; 0812.0559 . However, only late time annihilation was considered and not the freeze-out process. A consistent framework for an effective theory approach to the weak triplet (Wino) scenario was developed in Refs. Braaten1 ; Braaten2 ; Braaten3 . While models like these, or those with a new force carrier, such as a new neutral gauge boson suggested in Ref. PranNath , are certainly appealing, there are ways to introduce a more minimal dark sector with only one additional particle.

In this work, we consider the simplest set-up for a dark sector, with one new particle charged under the weak force ($SU(2)_{L}$). Despite the existence of multiple components of different charge in the multiplet, there are only two free parameters, the particle mass and its representation under $SU(2)_{L}$. As will become evident, the representation of the particle is its effective charge and determines its coupling strength to the weak gauge bosons. We show below that larger representations are extremely relevant for DM systems close to the unitarity bound. An important aspect arises from the fact that the weak force is non-abelian and thus, given two particles in a representation $R$ and $\bar{R}$, there is always an attractive (singlet) channel that makes Sommerfeld enhancement and bound-state formation absolutely generic phenomena.

The claim of Ref. Wise (see their Sec. 4), that bound state dissociation makes the net effect on the freeze-out irrelevant, was refuted in Refs. CosmoBS ; PetrakiNew . This was done through developing an analytic asymptotic solution to the system of Boltzmann equations arising from all allowed bound states Ellis:2015vaa ; CosmoBS . We use that mathematical framework to show the unavoidable effects of bound-state formation on the unitarity bound.

Freeze-out in the presence of bound states is described by a set of coupled Boltzmann equations, one for the dark-matter number density and one for each allowed bound state $I$ with the formation cross section $\langle\sigma_{I}v_{\rm rel}\rangle$ and annihilation rate $\Gamma_{I}$. The analysis relies on the fact that the DM bound state is not stable, as for example positronium, and can annihilate to SM particles or dissociate into WIMP constituents. The crucial simplification follows from the fact that for the WIMP bound states we study, their lifetimes are much less than the Hubble time. Then all the Boltzmann equations for bound states can be treated algebraically, neglecting the time derivatives, and inserted in the Boltzmann equation for the DM abundance CosmoBS . This leads to an effective cross section for DM, where the bound states lead to new channels with temperature-dependent branching ratios into SM particles. Intuitively, this effective cross section takes into account the fact that a bound state can be broken by the plasma before it annihilates, and thus not affect the total DM number density. The breaking rate is related to the formation rate in equilibrium by the Milne relation $\Gamma_{\rm break}\,n_{I}=\langle\sigma_{I}v_{\rm rel}\rangle n_{\rm DM}^{2}$. Once the temperature drops far below the bound-state binding energy, the breaking is strongly suppressed.

Quantitatively, this treatment leads to a single Boltzmann equation for $Y=Y_{\rm DM}=n_{\rm DM}/s$:

[TABLE]

where $\lambda=\sqrt{\frac{g_{\rm SM}\pi}{45}}M_{\rm Pl}M_{\text{DM}}$, $z=M_{\rm DM}/T$ and

[TABLE]

[TABLE]

where the rate for breaking of bound states follows from the Milne relation.

This equation can be easily integrated numerically, but has an analytic asymptotic solution, which agrees very well with the numerical treatment,

[TABLE]

with the inverse temperature at freeze-out $M_{\text{DM}}/T_{f}=z_{f}$ given by the transcendental equation

[TABLE]

For multi-TeV DM, $z_{f}\approx 25$ is typical.

In Fig. 1, we demonstrate the effects of the Sommerfeld enhancement and bound-state formation on freeze-out. Including the Sommerfeld effect leads to additional attraction among WIMPs and enhances the annihilation rate, which in turn reduces the relic abundance by $\mathcal{O}(10)$. The consideration of bound states is an additional effective annihilation channel and leads to a further $\mathcal{O}(10)$ reduction. This is not surprising since it is known that in the SM non-relativistic $e^{+}$ $e^{-}$ annihilation is dominated by positronium formation and its successive annihilation. Additionally, the importance of the decay width of the considered bound state is highlighted. The typical annihilation width scales as $\alpha^{5}M_{\rm DM}$, where $\alpha$ is the coupling strength of the interaction considered, and thus a typical width in a perturbative model would be of the order $10^{-5}M_{\rm DM}$ or smaller. The observation we want to stress is, that while a bound state can be a reaction product of dark-matter interactions, its effect on the relic density strongly depends on its binding energy and decay width to SM particles.

II.2 Effects on the Unitarity Bound

As discussed in the classic paper of Griest and Kamionkowski GriestKamionkowski , conservation of probability limits the reaction cross section of DM annihilating to any final state for each partial wave by

[TABLE]

Note the scaling of the bound with $v_{\rm rel}^{-1}$, which is not expected from contact type interactions, but is generic in the presence of long range forces. To understand the physical implications of the above inequality, we first discuss the cross sections that are relevant for the physical system. In the following, $(\sigma v_{\rm rel})$ denotes non-averaged cross sections and $\langle\sigma v\rangle$ denotes thermally averaged cross sections. The total (inelastic) reaction cross section is $(\sigma v_{\rm rel})_{\rm total}=\sum_{J}(\sigma v_{\rm rel})_{\rm total}^{J}$. The total reaction cross section is composed of an annihilation part and the bound-state formation cross section $(\sigma v_{\rm rel})_{\rm total}=(\sigma v_{\rm rel})_{\rm ann}+\sum_{I}(\sigma_{I}v_{\rm rel})_{\rm BSF}$. The relevant quantity for the freeze-out, as we have shown, is $(\sigma v_{\rm rel})_{\rm eff}=(\sigma v_{\rm rel})_{\rm ann}+\sum_{I}(\sigma_{I}v_{\rm rel})_{\rm BSF}\hbox{BR}(\text{B}_{I}\rightarrow\text{SM})\leq(\sigma v_{\rm rel})_{\rm total}$. The equality saturates only at zero temperature, otherwise the inequality holds, due to the fraction of bound states broken by ambient plasma quanta.

In Ref. GriestKamionkowski , the total reaction cross section is approximated as $(\sigma v_{\rm rel})_{\rm total}\approx(\sigma v_{\rm rel})_{\rm ann}$ and taken for the freeze-out computation, not considering the bound-state effects. The scaling, with the inverse velocity of this cross section, is however only possible in the presence of light mediators, which unavoidably lead to bound state formation 1703.00478 . Thus, in any perturbative physical system, saturating the unitarity bound on $(\sigma v_{\rm rel})_{\rm total}^{J}$, the inequality $(\sigma v_{\rm rel})_{\rm eff}\leq(\sigma v_{\rm rel})_{\rm total}$ leads to a lower maximally attainable DM mass than expected from considering only annihilation. This is one of the main findings of our paper and will be made quantitative in the coming sections.

The second case considered in Ref. GriestKamionkowski is the annihilation of extended DM objects. Then the total cross section becomes geometrical as a result of the sum over all partial waves allowed in the annihilation process,

[TABLE]

since by momentum conservation $\ell_{\text{class}}\approx p_{\rm DM}\,R$. Note that in the case of the geometrical interaction the cross section $\sigma$ is constant and the quantity relevant for the freeze-out scales as $\sigma v_{\rm rel}\propto v_{\rm rel}$. In principle, the required freeze-out cross section now sets a lower bound on the size of the dark object. This is, however, a zero-temperature computation and does not take into account that some angular momentum eigenstates will lose energy less efficiently and be broken by ambient plasma particles prior to annihilation. The situation is then similar to the point-like particle case since higher partial waves contribute to the total cross section at zero temperature, they will not efficiently contribute to annihilation in the hot plasma and thus change the annihilation efficiency. In Sec. IV, we demonstrate how this finite-temperature effect severely alters the obtained bounds.

III Unitarity for Pointlike WIMPs

In this section, we describe the non-perturbative effects in the simplest model of weakly interacting DM. The basic assumption is that DM interacts with SM particles and that, in the case of heavy DM, those particles induce a long-range force if, at the freeze-out temperature, we have $M_{\text{DM}}>M_{\rm SM}/v_{\rm rel}$ and $M_{\text{DM}}>M_{\rm SM}/\alpha$, where $M_{\rm SM}$ is the mass of a weak scale particle and $\alpha$ the strength of the interaction. Explicit examples have been studied in Refs. HiggsMediatedBS ; CosmoBS .

To parametrize the explicit model, we use $SU(2)$ symmetry as the guiding principle. We consider a particle in the $SU(2)_{\rm L}$ representation $R$. To have a neutral component, the hypercharge has to be [math] for odd $R$ and $1/2$ for even $R$. For the heavy freezeout computation, the hypercharge gives just a subleading effect and $R\geq 2$ will be treated as a continuous parameter. For concreteness, we consider real representations of $SU(2)$ and interpolate in between. Furthermore, DM is assumed to be a Majorana fermion, such that the effective number of degrees of freedom $g_{\chi}=2R$. This model, even though specific in terms of $SU(2)_{L}$ quantum numbers, serves as an excellent template to parametrize a generic model which has a coupling to weak gauge bosons and we expect the general statements to hold universally. The tree level annihilation cross section is

[TABLE]

with $\sigma_{0}=\pi\alpha_{2}^{2}/M_{\text{DM}}^{2}$ MinimalDM .

III.1 Sommerfeld Enhancement

The annihilation rate is enhanced by the Sommerfeld effect, which decomposes into isospin channels as

[TABLE]

where the weight factor $f_{I}$ can be found by explicit computation and is listed in Table 1, the $S_{\rm ann}(\lambda)$ are the Sommerfeld factors; see Ref. Feng:2010zp for their explicit form. The $\lambda_{I}$ is the effective coupling of each channel and controls the interaction strength, this coupling can be expressed in terms of quadratic Casimir operators, which leads to $\lambda_{I}=1/2[2C_{2}(R)-C_{2}(I)]$. For any $R$ the combination contains $R\otimes\bar{R}\supset 1\oplus 3\oplus...2R-1$. Since among these representations the singlet and adjoint channels have the smallest quadratic Casimirs the corresponding channels are the most attractive ones, explicitly $\lambda_{1}=(R^{2}-1)/4$ and $\lambda_{3}=(R^{2}-5)/4$ and, sub-dominantly (for $R\geq 3$), $\lambda_{5}=(R^{2}-12)/4$. The larger representations are less attractive, or even repulsive with an effective Sommerfeld suppression factor, similar to Ref. 1611.04606 . We point out that since the computation is performed in the $SU(2)$ multiplet formalism, we stop our integration of the Boltzmann equations at $z\approx 10^{4}$, because at lower temperatures mass-splitting effects would become relevant MinimalDM ; 1411.0752 . However, for large DM masses this is a very good approximation and errors remain below the percent level.

III.2 Bound-State Effects

For the estimate of the bound-state formation effect, it is sufficient to take into account the bound states in the 1 (singlet) and the 3 (adjoint) channel, as they correspond to the most attractive channels. As discussed in Ref. CosmoBS , the wave function symmetry leads to selection rules, such that the singlet bound state has spin-[math] for $\ell$ even and spin-$1$ for $\ell$ odd and the adjoint bound state has spin-[math] for $\ell$ odd and spin-$1$ for $\ell$ even. In the $SU(2)$-symmetric limit, the singlet bound state can only be formed from an adjoint initial configuration, but the adjoint state can be formed from the singlet initial configuration and, if $R\geq 3$, from the $5$ configuration.

III.2.1 Formation Rates

The total formation cross section can be written in a similar way to the Coulomb cross section, but taking into account different potentials for the initial and final states and a factored-out non-abelian structure,

[TABLE]

For example, in the non-relativistic limit, the capture cross section to the ground state can be written as

[TABLE]

where $S$ is the spin of the bound state and $C_{J}$ and $C_{\tau}$ are group theory factors for the abelian and non-abelian emission, respectively. The full expressions, including the excited states, can be found in Ref. CosmoBS .

The group-theory factors depend on the transition and have in general a complex structure, see Refs. CosmoBS ; 1805.01200 for the full expressions. In the cases of the transitions we are interested in, however, they can be rewritten in a general simple way:

[TABLE]

With these factors and the $\lambda_{I}$ for $I=\{1,2,3\}$, the cross sections can be readily computed. Note that capture to the adjoint state consists of two transitions from the singlet and (if $R\geq 3$) the $5$ state, which makes it relevant for larger representations. As discussed in Ref. CosmoBS , the effects of massive mediators can be taken into account by a kinematic factor, which will suppress bound-state formation, once the binding energy is not sufficient to emit the gauge boson.

A further aspect to take into account is the formation of excited bound states, which can be taken into account by means of Kramer’s formula; see Ref. Kramer . The resummation leads to a logarithmic enhancement of the capture rate with $\approx\left(1+\kappa\ln{\left(\alpha_{\rm eff}/2v_{\rm rel}\right)}\right)$.

III.2.2 Annihilation and Decay Rates

From Fig. 1, it becomes clear that the annihilation rates are crucial for judging the impact of a bound state on the final DM relic abundance. As discussed in Ref. CosmoBS , the excited $\ell>1$ states have strongly suppressed annihilation rates and decay to the ground states in order to annihilate. We work with selection rules for real $SU(2)$ representations, which have $Y=0$, as pseudo-real representations are only viable as admixtures, as discussed in Ref. Tytgat . Given the selection rules for identical fermions, the singlet ground-state bound states will have spin-[math] and annihilate to two gauge bosons with a rate of

[TABLE]

where $T_{R}$ is the representation index. On the other hand, the selection rules dictate that spin-$1$ states are triplets of $SU(2)$ and their annihilation to gauge bosons is suppressed. They will annihilate democratically into SM fermions with the rate

[TABLE]

The total rate contains the SM fermion multiplicity $n_{f}=3(3+1)$, $T_{\rm SM}=1/2$ and $\eta=2$ for self-conjugrate DM and $\eta=1$ for non-self-conjugrate DM.

In Fig. 2, we show the DM mass that leads to the correct DM relic abundance given a representation $R$. We see that our estimates agree well with the full computation for the $R=\{2,3,5\}$ representation. In addition, we show that for $R\geq 4$, the bound-state formation effects lead to substantial corrections on top of the Sommerfeld corrections.

III.2.3 Partial-Wave Contributions

As already discussed in Ref. GriestKamionkowski , the annihilation cross section $(\sigma_{\rm ann}v)$ is dominated by the s-wave annihilation process. We perform analogously the partial wave decomposition of the amplitudes for the bound-state formation cross sections. This is done by obtaining the different $J$ contributions by projecting the transition amplitude on the Legendre polynomials. For example the $(\sigma_{\rm BSF}v)^{n=2,\ell=1}$ contains a significant $s$-wave fraction, in agreement with 1611.01394 . By explicit computation, it can be verified that the amplitudes for higher excited states, $(\sigma_{\rm BSF}v)^{n\ell}$, when summed over $\ell$, contain the same s-wave contributions. Those states are the ones relevant for the derivation of Kramer’s resummation formula. We find that the s-wave processes contribution to the total cross section is thus given by

[TABLE]

This cross section will be used to find the value of the coupling strength, controlled by $R$, at which the $s$-wave contribution saturates the unitarity bound.

III.3 Unitarity Bounds for Pointlike WIMPs

Each value of the $R$ parameter corresponds to a DM model where the lightest new particle can be a superposition of fermions but has a dominant contribution, which transforms in the $R$ representation of $SU(2)_{L}$. Therefore, the only free parameter for such a model is the mass of the DM candidate, and we determine it to satisfy the observed relic-density value.

In Fig.  3, we show the relic abundance as a function of the DM mass given that its cross section saturates the unitarity bound. The computation is performed in the following way: we demand that $(\sigma v)_{\rm total}^{0}<(\sigma v)_{\rm max}^{0}$, which is set by partial wave unitarity, Eq. (6), and extract for each given DM mass a value for $R_{\rm max}$. With this maximal coupling and the DM mass, we compute the relic abundance using now the full bound-state formation cross sections, not limiting them to the $s$-wave processes. The allowed DM mass range is found to be $M_{\text{DM}}\leq 144\,\rm TeV$ ($M_{\text{DM}}\leq 96\,\rm TeV$ for a non-self-conjugate particle). We have to stress that the obtained mass bound is subject to significant uncertainty since the effective coupling values are large. Thus higher-order processes, such as multiple gauge boson emission, are likely to become relevant. The DM mass bound is an important result of this paper and applies to asymmetric DM as well. Note that without taking into account the bound-state formation, which helps to saturate the unitarity bound, but does not fully contribute to the annihilation cross section, the allowed mass range would be significantly higher, $M_{\text{DM}}\leq 200\,\rm TeV$.

This allows an interesting argument for new physics at the 100-TeV scale. As a historical example, the breakdown of perturbative unitarity in the gauge boson scattering cross section near 1 TeV indicated the existence of the light SM Higgs boson HiggsUnitarity . Here, the upper bound on the DM annihilation cross section is a consequence of probability conservation only. In the case of a pointlike WIMP, this bound induces an upper bound on the WIMP mass and thus renders the WIMP parameter space finite. Furthermore, as we have demonstrated, the scale of new physics, containing a least the WIMP particle is below 100 TeV. This mass scale is low enough to be explored by future collider experiments. This conclusion is reached without appealing to naturalness or other aesthetic considerations.

III.4 The Abelian WIMP

Previous papers investigated the freezeout of fermionic DM interacting via a dark Abelian force, also in the case of asymmetric self-interacting DM 1712.07489 . While Ref. Petraki concluded that taking into account bound-state formation reduced the unitarity limit on the DM mass to $139$ TeV, a later study suggested that this was incorrect 1611.01394 .

Here, we consider the impact on freezeout of bound states with spin-1 or spin-0 with $n=1$. It is argued that at leading order the selection rule $\Delta\ell=1$ leads to $p$-wave dominance of the bound-state formation in the ground state. Therefore, the bound-state formation does not contribute to the saturation of the $s$-wave unitarity bound and the maximal dark gauge coupling allowed by unitarity remains $\alpha_{D}<0.86$, which is determined only from the Sommerfeld enhanced annihilation cross section. This is the reason, why in 1611.01394 no DM mass reduction was found.

The argument based on the leading order selection rule is, however, problematic given the size of the gauge coupling involved. We therefore investigate the effect of two-photon emission adopting the formalism of Ref. Bhanot:1979vb . The $\Delta\ell=0$ contribution to the bound-state formation in the ground state is $s$-wave dominated and given by

[TABLE]

where $\phi$ is the radial part of the free two-particle state wave-function, $\psi$ are the radial parts of bound-state wave-functions, the numerical factor $\tilde{\kappa}\approx 20$ and $\bar{\sigma}_{0}=\pi\alpha^{2}_{D}/M_{\text{DM}}^{2}$. Taking this contribution to the ground-state formation into account reduces the $s$-wave unitarity bound on the dark gauge coupling to $\alpha_{D}<0.6$ and leads furthermore to a reduced maximally attainable DM mass of the order of $M_{\text{DM}}<150\,\rm TeV$. (Compared to the bound of $M_{\text{DM}}<200\,\rm TeV$ if only the Sommerfeld effect was considered). This finding is in agreement with the expectations that we elaborated on earlier. Because bound-state formation contributes to the total inelastic cross section, it reduces the maximum viable coupling strength, but because it is not always efficient in reducing the DM abundance, the bound on the maximal DM mass is decreased. This example shows at the same time that close to the unitarity bound, the large coupling strength makes higher-order processes relevant and therefore the computed cross sections are expected to have a substantial theoretical uncertainty.

III.5 Observable Signatures

The effective WIMP framework makes concrete predictions for various channels and connects different observables. We show that several unexpected new signatures arise. This predictive power is the result of $SU(2)_{L}$ gauge symmetry, which we used to construct the effective WIMP model.

A new exciting possibility for indirect detection is the observation of gamma rays emitted in the formation of the dark-matter bound states. Those capture gamma rays are nearly monochromatic and carry away the binding energy.

A few generic statements can be made about this new process. It is generally true that the neutral component of the weak multiplet, which is the lightest particle in the spectrum and the natural DM candidate, has no triplet component in the wave-function. Therefore, no bound state can be formed in the singlet configuration at zero temperature.

The largest formation cross sections will thus be attributed to the adjoint, i.e., the triplet configuration. The most probable main quantum numbers will be the ground state, with $n=1$, i.e., the $1s_{3}$ state and the lowest excited state, with $n=2$ and $p=1$, i.e., the $2p_{3}$ state. The production cross section for the latter is sizable, since it is produced from the s-wave configuration and is not suppressed at low velocities. Furthermore, each excited state will decay to the ground state before annihilation, in particular there will be the $2p_{3}$ to $1s_{1}$ transition. We can predict, for each representation $R$, the relative position of the dominant gamma lines in the Coulomb limit

[TABLE]

In Fig. 4 (a), the annihilation cross section at freeze out in the early universe is shown. The cross section, in the presence of long-range forces scales as $1/v_{\rm rel.}$. Since the average velocities at freeze-out are of order $v_{\rm rel.}\sim\mathcal{O}(0.3)$, the cross section is close to the well known value of $2.2\times 10^{-26}\,\text{cm}^{3}\,\text{s}^{-1}$. In contrast, the late time annihilation is significantly enhanced, since the average DM velocity in galaxies today is $v_{\rm rel.}\sim\mathcal{O}(10^{-3})$. This is indicated in Fig. 4 (b) for all allowed final states.

In Fig. 4 (b) ($\gamma-$line), the capture gamma line cross sections are shown as a function of DM mass. Measuring the ratio of those three line energies provides enough spectroscopic data to determine the $R$ parameter and thus, the DM gauge charge. This seems to be a unique possibility, in which indirect detection experiments can probe not only the magnitude of the annihilation cross section to test the freeze-out hypothesis but also directly measure the DM connection to the SM. Since given the large DM mass the capture photons emitted will be also hard gamma rays, resummation techniques can become important for the precise spectrum prediction 1712.07656 ; 1808.08956 . Because at late times DM entirely consists of the neutral $\chi_{0}$ component, we have to project out its $SU(2)$ multiplet components to compute the late time annihilation signals. In our analysis this is done by use of the Clebsch-Gordan coefficients. This procedure was already used in Ref. CosmoBS and a direct comparison of 5-plet spectra to previous exact results 0706.4071 ; 1507.05519 ; 1608.00786 shows that it is a good estimator for the order of magnitude signal strength.

In Fig. 5 panel (a), we show the capture photon spectra for a $14$ TeV DM candidate with $R=5$ on top of the gamma-ray continuum. The figure demonstrates that given a $10\%$ energy resolution, the lines will be visible. Panel (b) shows the high energy line signal from direct annihilation, which is a particularly good target for the H.E.S.S observatory 1810.00995 . We point out that because the signals are strongly separated in their respective energy range, the experimental searches for them will be complementary.

In Fig. 6, the spin-0 bound state annihilation in shown, which contributes to the same final states as the dominant direct WIMP annihilation channel, and thus enhances the expected signals.

In Fig. 7, we show another new avenue, which the bound states open up: reactions with light final states. In many models, annihilations to neutrinos or electrons are suppressed due to the smallness of their mass (chirality suppression, e.g., Ref. LeptophillicDM ) or due to the smallness of their Yukawa couplings to the Higgs, e.g., Ref. DirectVSIndirect . On the contrary, weak WIMP bound states with spin-1 cannot decay to two gauge bosons, and have suppressed annihilation channels to three gauge bosons. The dominant annihilation channel for those states will thus have two SM fermions in the final state. The $SU(2)_{L}$ invariance dictates that $25\%$ of those will be leptons and that half of those leptons will be neutrinos. In contrast to spin-1 bound-state annihilation, the direct WIMP annihilation to leptonic final states is subdominant for large representations MinimalDM . Given the current low energy resolution of neutrino observatories, the fact that the produced neutrinos are monochromatic is less relevant than the overall enhancement of the high-energy neutrino flux.

In Fig. 4 (b), we show the predicted cross sections for those nearly monochromatic, high-energy neutrinos. These deserve particular attention since they will be easier to distinguish from power-law astrophysical backgrounds. An upper bound on the DM annihilation cross section using neutrino signals was already derived in Refs. UpperBoundCross ; NuBoundsYuksel . It is intriguing that at DM masses approaching the unitarity bound, non-perturbative effects unavoidably lead to an enhancement, raising the value of the annihilation cross section above the naively expected one.

Finally, it is possible to predict the spin-independent scattering cross section with nucleons, needed for the direct detection experiments, as discussed in Ref. RunningJap . The basic assumption is that the DM candidate does not couple directly to the Z boson. This is always true for representations with odd $R$, while representations with even $R$ can be viewed as variants of two or more representations, which mutually mix due to Higgs Yukawa interactions. An example is a Higgsino-like neutralino in supersymmetric models. The general framework for such models has been discussed in Ref. Tytgat .

In Fig. 8, we show the predictions for the spin-independent cross section as function of the representation parameter $R$. It is encouraging that almost the entire parameter space of the weakly coupled DM models predicts direct detection cross sections above the neutrino floor. This is an interesting feature of the suggested framework since the available parameter space is finite due to theoretical considerations. It is even more intriguing, since upcoming direct detection experiments, as the DARWIN experiment DARWIN , will cover it entirely.

III.6 Theory for Large Representations

For DM candidates close to the unitarity limit, there is a question of theoretical consistency because, at large coupling constants of large representations, low-lying Landau poles could render the models nonviable; see Ref. MinimalDM . Since in that case, a DM theory needs a low-lying cut off scale, DM stability might well be compromised. However, recent developments in quantum field theory show that the paradigm of so-called triviality could be challenged.

As discussed in Ref. BubbleSum , new resummation techniques indicate that it is possible that a theory that exhibits a strong renormalization flow running enters an asymptotically safe fixed point. In particular, it has been demonstrated that in an explicit construction with vector-like fermions, the weak coupling constant has a safe fixed-point. The computation relies on an expansion in the small parameter $1/N_{f}$ , where $N_{f}$ is the number of fermion flavors; see Ref. Holdom for the resummation technique.

A different argument, presented in 1512.05353 , shows that introducing a millicharge $\epsilon$ instead of the charge neutrality requirement for DM leads to a stable DM candidate without the accidental symmetry assumption. Thus, even in an EFT with sub-planckian Landau pole electroweak multiplets in larger representations are viable DM candidates.

Those considerations open up an entirely new realm of models that have been previously ignored. Concretely, it implies that fermions in $SU(2)_{L}$ representations larger than the quintuplet are good, minimal DM candidates. Conversely, the discovery of a DM particle in a large representation would be a smoking-gun signal of an asymptotically safe model.

IV Unitarity Bounds for Extended Objects

In this section, we show that bound states are also crucial for the unitarity bound for composite DM systems. As discussed in Ref. GriestKamionkowski , the unitarity argument can be also applied to the annihilation of extended composite objects, such as DM atoms, see for example Refs. Khlopov:2011tn ; Khlopov:2014bia . These particles can be seen as WIMPs in the wider sense since they are massive, have initial thermal equilibrium with the SM plasma and a cross section, weak enough such that a significant relic abundance survives asymptotically. The intuitive reason why their cross section is big is that all partial waves up to $l_{\rm class}\approx p\times R$ contribute, where $R$ is the size of the object and $p$ its momentum. For example, the atomic cross section for hydrogen is large, since the Bohr radius $R\propto 1/\alpha_{\rm em}m_{e}$ is controlled by the mass of the electron, which is much smaller than the mass of the entire atom.

IV.1 The Simplest Dark Atom

In the case of heavy, multi-TeV DM, the low number density makes hydrogen-like recombination very inefficient Tesi . Thus if heavy DM is an atom, it is formed by a non-abelian confining gauge force. The simplest model for this is an $SU(N)$ gauge theory with a fermion in the adjoint representation. Such a set-up is present in any non-abelian supersymmetric Yang-Mills theory. The minimal example is an $SU(3)$ theory with a fermion ${\mathcal{Q}}$ in the $8$ representation, as introduced in Refs. gluino ; gluequark .

The fermion can have an explicit vector-like mass. Below this mass scale, the running of the gauge coupling $\alpha$ will lead to a strong coupling regime at a scale $\Lambda$ given by $\Lambda\approx M_{\mathcal{Q}}\exp{\left(-6\pi/11N\alpha(M_{\mathcal{Q}})\right)}$. That means that at energies below $\Lambda$ the fermion will be forced to form bound states.

In the relevant mass regime for DM, the fermion number density is low. Therefore, at the phase transition string breaking will lead to the formation of bound states of adjoint fermions and dark gluons $B={\mathcal{Q}}g$. The stability of this state is guaranteed by an accidental global symmetry. The size of the bound state is set by the confinement scale $R\approx 1/\Lambda$, but its mass is controlled by $M_{\mathcal{Q}}\gg\Lambda$ and its de-Broglie wavelength is much shorter than its size, as in the case of hydrogen.

IV.2 The Atomic Relic Abundance

Annihilation happens via ${\mathcal{Q}}g+{\mathcal{Q}}g\rightarrow{\mathcal{Q}}{\mathcal{Q}}+gg$, where the $G=gg$ is the glueball, the lightest degree of freedom in this theory, and the ${\mathcal{Q}}{\mathcal{Q}}$ state can self-annihilate efficiently once it has fallen to the ground state. Given that $M_{\mathcal{Q}}\gg\Lambda$, the deeply bound states of ${\mathcal{Q}}{\mathcal{Q}}$ can be treated perturbatively, as heavy quarkonium.

The fermions in the adjoint representation can arrange into a singlet or an adjoint bound state. Higher bound-state representations will be less bound or even repulsive, as is the case for $SU(3)$. Here the two octets can arrange in the di-quark ${\mathcal{Q}}{\mathcal{Q}}$ state in the $1\otimes 8_{S}\otimes 8_{A}$ representations, arrangements in higher representations are repulsive. Pure combinatorics suggests that the adjoint configurations are formed with a probability $p_{\rm adj.}=1/(N^{2}-1)$ and the singlet with $p_{1}=1/(N^{2}-1)^{2}$. Furthermore, the adjoint bound state can lose energy efficiently by radiating a gluon in the adjoint representation, while the singlet has to emit at least two gauge bosons, which is a higher order process. Therefore, it is sufficient to focus on the reactions of the adjoint configuration.

Now we use the principles of energy and momentum conservation to derive a reaction cross section bound for the extended object. The unitarity limit on the dark-matter mass will be a direct consequence.

The bound state is formed in all partial waves allowed by the energy conservation and momentum conservation. The energy balance of the reaction is $2E_{B}^{{\mathcal{Q}}g}+E_{\rm kin}^{\rm CMS}\geq E_{B}^{{\mathcal{Q}}{\mathcal{Q}}}+M_{G}$. Assuming that for the binding energies of the non-perturbative states, the energy of the gluon field in the octet configuration of ${\mathcal{Q}}{\mathcal{Q}}$ and the glueball mass $M_{G}$ are of the order of $\Lambda$, we have the condition for the binding energy between the heavy quarks $E_{B}^{{\mathcal{Q}}{\mathcal{Q}}}<E_{\rm kin}^{\rm CMS}\approx T<\Lambda$, as we consider the post-confinement regime.

Heavy quarkonium can be well described by a Coulomb potential with a linear potential as a perturbation, relevant for states close to the continuum. Their energy is approximately given by $E_{n,\ell}=E_{\text{Coulomb}}+E_{\text{Linear}}\approx\alpha_{\rm eff}^{2}M_{\mathcal{Q}}/4(-1/n^{2}+12\epsilon n^{2})$, with $\epsilon=\sigma/(\alpha_{\rm eff}^{3}M_{\mathcal{Q}}^{2})\approx\lambda_{R}\Lambda^{2}/(\alpha_{\rm eff}^{3}M_{\mathcal{Q}}^{2})$ Potential . Here the potential part, dominated by the string tension $\sigma$, is approximated by $\lambda_{R}\Lambda^{2}$. The combination of the representations’ quadratic Casimir operators $\lambda_{R}=1/2\left(2C_{\rm adj.}-C_{R}\right)$ defines the strength of the channel.

The above energy-conservation condition ($E_{n,\ell}<\Lambda$) demands that the correction from the linear potential is not dominant and it thus limits the maximal partial wave $\ell$ by

[TABLE]

After deriving the general bound $\ell_{\text{max}}$ from energy conservation, we come now to the momentum conservation condition. The maximal $\ell$ allowed by classical momentum conservation is

[TABLE]

where $1/\Lambda$ is the size of the ${\mathcal{Q}}g$ bound state. Since $\ell_{\rm class}$ scales with the temperature, there is a $T_{c}$ below which $\ell_{\rm class}<\ell_{\rm max}$, implying that all classically allowed partial waves contribute to the cross section, leading to a geometric interaction. Solving for this temperature, we get $T_{c}=\Lambda\,\lambda_{R}\,\alpha^{3/2}/8\sqrt{3}$. This leads to a compact expression for the maximal partial wave allowed by energy conservation, $\ell_{\rm max}^{2}=4\,T_{c}M_{\mathcal{Q}}/\Lambda^{2}$.

The cross section for the reaction for the two quark state in representation R decomposed in partial waves reads

[TABLE]

where we define $L_{\rm max}(T)=\text{min}[\ell_{\rm max},\ell_{\rm class}]$ which is set by either energy or momentum conservation respectively. $\mathbb{P}_{\ell}^{\rm ann}$ is the probability that a bound state with given $\ell$ loses energy fast enough, falls to the ground state and self-annihilates.

For this cross section, unitarity is saturated if all phases have $\sin(\delta_{\ell})=1$. To derive the unitarity limit, it is sufficient to assume that all states that are allowed by energy and momentum conservation can self-annihilate fast enough and thus $\mathbb{P}_{\ell}^{\rm ann}=1$ for all of them. In fact, in Ref. gluino , it has been explicitly shown that this is close to the numerical result for the gluino. Furthermore, to obtain an upper bound on the cross section, the probability $p_{R}$ stemming from combinatorics can be set to unity as well.

Summing the resulting cross sections up to the maximally allowed partial wave at a given temperature leads to the following maximal cross section compatible with unitarity for $T<\Lambda$:

[TABLE]

with $\sigma_{\rm geom.}=4\pi/\Lambda^{2}$. This implies for the thermally averaged cross section,

[TABLE]

It is obvious that the result of Ref. GriestKamionkowski is recovered at zero temperature, however, as in the perturbative calculation in the first part of this paper, we find that the effect of the thermal plasma reduces the cross section above the critical temperature $T_{c}$.

The relic abundance is given by

[TABLE]

Integrating over both scaling regimes we get the asymptotic relic abundance

[TABLE]

In fact, the relic density is dominantly determined at temperatures close to the phase transition, such that the geometrical cross section is less relevant for its asymptotic value than the cross section, which scales with temperature. This is a new result, obtained by taking into account the finite temperature effect and the basic principle of energy and momentum conservation.

IV.3 Implications of Geometrical Unitarity

In Fig. 9, the results are summarized, with our bound coinciding with the zero temperature bound of Ref. GriestKamionkowski only if $M_{\mathcal{Q}}\approx\Lambda$. In this regime however, the dark atom mass is set by the scale $M_{\rm DM}\approx\Lambda$ and thus $R_{\rm DM}\approx 1/\Lambda\approx 1/M_{\rm DM}$, which implies that $\ell_{\rm max}=0$. The reaction is therefore just an s-wave process and not geometrical because of the contributions of many partial waves.

On the contrary, if many partial waves contribute and $\ell_{\rm max}\gg 1$, we are in the regime with $M_{\mathcal{Q}}\gg\Lambda$ and $R_{\rm DM}\approx 1/\Lambda\gg 1/M_{\rm DM}$. We find that in this case the unitarity bound depends on the mass $M_{\mathcal{Q}}$ and results in a larger DM radius corresponding to a monotonically decreasing $\Lambda$ with growing $M_{\mathcal{Q}}$.

We demand that the scale $\Lambda>\Lambda_{\rm QCD}$, as otherwise big bang nucleosynthesis will be affected. The observation is that the existence of dark atoms would require a new confinement scale below approximately $60\,\rm TeV$. This contrasts with the finding of Ref. GriestKamionkowski that the limit on the DM size is independent of its mass. So, in the case of the lowest possible confinement scale $\Lambda\approx\Lambda_{\rm QCD}$, the DM mass can be as high as a PeV, but not much above that. Note that this upper bound on the mass also applies to a scenario where DM is co-annihilating with a partner particle, which has a non-abelian charge. For a detailed discussion, see Ref. gluino .

In Fig. 10, the late time and low-temperature geometrical annihilation cross section of the dark atoms is shown. It exceeds by far the cross section expected from $s$-wave unitarity in the regime where $\Lambda\gg M_{\mathcal{Q}}$ and many partial waves contribute. This can be easily understood, as due to the effects of the plasma, the annihilation in the early universe is less efficient and the effective cross section setting the relic abundance is smaller given the same model parameters.

When $\Lambda\approx M_{\mathcal{Q}}$ the cross section is $s$-wave dominated, but since the first step before the annihilation is a rearrangement process which is not strongly exothermic, we have $\sigma\approx\rm constant$ instead of $\sigma v\approx\rm constant$, and the annihilation is suppressed at late times. We note at this point that in this regime the quantum numbers of the DM particles with respect to the SM are important. When the geometrical annihilation cross section drops, they will become dominant and match the composite DM regime to the regime with elementary particles with dominant $s$-wave annihilation, as discussed in the first part of the paper.

In Fig. 10, we compare the total annihilation cross section to the HAWC experiment bounds from Ref. HAWC assuming conservatively a dominant annihilation into neutrino final states, as suggested in Refs. UpperBoundCross ; NeutrinoConstraints ; VHDM , and a Burkert DM profile, which features a central core. The HAWC experiment can set limits on the annihilation cross section into neutrino final states, as at the large energies involved electro-weak corrections lead to unavoidable photon emission. Those corrections are taken into account by the PYTHIA software 1410.3012 , which was used in Ref. HAWC . Note that the predicted large annihilation cross sections at large DM masses are compatible with the excess of high energy neutrino events found in IceCube IceCubeexcess . We point out that composite dark matter, which saturates the unitarity can have intriguing direct detection signatures, as discussed in 1311.6386 ; 1505.02772 .

V Conclusions

In searching to discover the particle nature of dark matter, the minimal thermal-relic WIMP model must be fully tested. Because this model is defined by its annihilation cross section at freezeout in the early universe, the decisive test is to use searches for late-universe annihilation products to probe down to below the corresponding cross section scale. Once this sensitivity is reached for a given mass, which requires allowing any combination of Standard-Model final states, then the minimal model is ruled out at that mass. Importantly, the mass range is finite. At the low-mass end, the thermal-WIMP hypothesis is being strongly tested (excluding neutrinos as final states, the current lower limit is about 20 GeV LeaneSlatyerBeacom ), and that great progress is expected from new experiments.

But, to fully test thermal WIMPs, this hypothesis must also be strongly tested at high-mass end. Here there are so far no experimental limits that approach the required sensitivity, so we must fall back on the unitarity bound, which intersects the thermal-relic cross section line in the range of hundreds of TeV.

Here we re-examine the physics of heavy WIMPs, focusing on a minimal model where the only new particle is the heavy WIMP itself so that its interactions must be mediated by Standard-Model particles. Because those are relatively lighter than the heavy WIMP, they act as light mediators, which brings in Sommerfeld-enhancement effects, which are well studied, and bound-state effects, which are much less so. We take advantage of the new formalism and results of Refs. Petraki ; Wise ; PetrakiNew ; CosmoBS ; Tytgat , which allow us to set up a formalism to study the freezeout of a particle in a generic $SU(2)_{L}$ representation, thus a true WIMP. As the unitarity bound is approached, the growing interaction strength renders a tree-level approach inadequate. Crucially, we use resummation techniques to take into account the non-perturbative effects of Sommerfeld enhancement and bound states.

For pointlike WIMPs, we find the following:

Taking into account bound-state formation effects in the early universe, the upper mass bound is reduced. In a UV-complete model of elementary dark-matter particles linked to the Standard Model by weak gauge bosons, the expected maximal dark-matter mass is lowered by about $30\%$. We find an upper bound of 144 TeV for self-conjugate dark matter and about 96 TeV for non-self-conjugate dark matter. Those mass scales, while high, are conceivably within reach of future collider experiments, more so then the oft-quoted value of $340$ [TeV]. 2. 2.

The late-time annihilation of WIMPs is significantly enhanced due to Sommerfeld and bound-state formation effects. Furthermore, bound-state formation leads to new signals, such as capture photon emission and enhanced annihilation to monochromatic neutrinos.

We find that within the finite parameter space of the WIMP particle carrying $SU(2)_{L}$ charge, the direct detection signals are above the neutrino floor. Thus the entire parameter space will be tested by the next generation of underground direct-detection experiments.

For composite WIMPs, we find the following:

Taking into account the energy conservation of different partial waves contributing to geometrical annihilation, the resulting cross section cannot be constant at high temperatures. In particular, this finding leads to an upper bound on the dark matter radius, which depends on the dark matter mass, in sharp contrast to previous results GriestKamionkowski . Further, demanding that the physics of the big bang nucleosynthesis remains unchanged, leads to an upper bound on the composite dark matter mass of about 1 PeV. 2. 2.

The newly found temperature dependence of the annihilation cross section is such that late-time annihilation is greatly enhanced. Therefore, despite the large mass, thermally produced extended dark matter particles appear to be testable by indirect detection experiments. Additionally, we find that in order to have a dark matter candidate near the heavy end of the allowed spectrum, the size of the composite dark matter grows, indicating a confining dynamics below the 100 TeV scale.

The phenomenological consequences of our framework are quite intriguing. We found that the bound-state formation processes at late times open up two new observational opportunities. On the one hand capture photons, which have the binding energy of the formed bound state and are nearly monochromatic are a new exciting signature for gamma-ray searches. On the other hand spin$-1$ bound states in our set-up have by $SU(2)_{L}$ symmetry significant annihilation rates to neutrinos and thus are good candidate processes to search for with experiments like IceCube, HAWC and Antares; see Refs. IceCubeexcess ; HAWC ; Antares . In particular, the monochromatic nature of the neutrino signals allows a powerful background rejection of astrophysical sources and can become a unique tool for heavy dark matter searches. Concluding, we emphasise that the gamma-ray (HAWC 1804.00628 , LHAASO Ma:2018mzp , CTA 1709.07997 ), cosmic-ray (AMS Zimmermann:2017abs , CALET Motz:2015qna , DAMPE TheDAMPE:2017dtc ), and high-energy neutrino experiments (IceCube Flis:2017uay , Antares Antares ) provide a very promising route to test annihilations of TeV scale elementary and composite dark matter.

acknowledgments

We thank Eric Braaten, Marc Kamionkowski, Mark Wise, Yue Zhang, and especially Iason Baldes, Simone Biondini, Marco Cirelli, Basudeb Dasgupta, Nicola Dondi, Ranjan Laha, Rebecca Leane, Kallia Petraki, Michele Redi, and Bei Zhou for helpful comments on our manuscript.

J.S. is thankful for the hospitality of CCAPP, where this project was initiated. The work of J.S. was supported by the INFN division in Florence and by the CP3-Origins centre, which is partially funded by the Danish National Research Foundation, grant number DNRF90. The work of J.F.B. was supported by NSF grant PHY-1714479.