TeV-Scale Thermal WIMPs: Unitarity and its Consequences

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.


I. INTRODUCTION
The unknown particle nature of dark matter has inspired a plethora of imaginative models [1][2][3][4][5][6]. 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 [7][8][9][10]. 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 σv = (2.2×10 −26 cm 3 s −1 )(0.12/Ω DM h 2 ), where this is the total cross section to all final states [11]. (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 lateuniverse annihilation rate factor is the same. Given the density distribution of dark matter, determined by gravitational probes [12], upper limits on the fluxes of energetic particles (gamma rays, cosmic rays, and neutrinos) then determine upper limits on σv 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 [13][14][15] and elasticscattering rates in underground detectors [16][17][18] -are model-dependent.
Fully testing the thermal-WIMP hypothesis requires reaching a cross section sensitivity below this prediction, * Email: smirnov@cp3.sdu.dk; ORCID: 0000-0002-3082-0929 † Email: beacom.7@osu.edu; ORCID: 0000-0002-0005-2631 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 [19,20]. If one excludes neutrinos from being dominant final states, then the lower limit is 20 GeV [20], to be contrasted with the limits claimed assuming favorable final states, which approach 100 GeV [21]. (WIMPs that annihilate only to neutrinos are excluded below ∼ 10 MeV by big bang nucleosynthesis [22,23]; if there is any annihilation to other channels, the mass limit from the cosmic microwave background is much larger [20,24,25].) It is well known that improvements in the sensitivity of existing searches are possible, and that, if no signals are found, the thermalrelic mass range will be progressively attacked from the low-mass end [26].
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 340 TeV (for Ω DM h 2 = 1; see below for Ω DM h 2 < 1) for pointlike dark matter [27]. For composite dark matter, the lower limit on the radius is 7.5 × 10 −7 fm, independent of mass [27]. 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 reexamination of thermal WIMPs at the largest masses. New theoretical developments [28][29][30][31][32] 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 arXiv:1904.11503v1 [hep-ph] 25 Apr 2019 limited. (For example, it is not possible to invoke resonances [33,34], 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 longrange potentials and enhance cross sections through the Sommerfeld effect [28,29]. (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 [30,32]. 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.

A. 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. [30] to be important, though Ref. [31] argued the opposite. A follow-up study [35] discusses boundstate 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 pure weak doublet (Wino) system were studied in Refs. [36][37][38]. However, only late time annihilation was considered and not the freeze-out process. A consistent framework for an effective theory approach to the Wino scenario was developed in Refs. [39][40][41]. While models like these, or those with a new force carrier, such as a new neutral gauge boson suggested in Ref. [42], 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 andR, there is always an attractive (singlet) channel that makes Sommerfeld enhancement and bound-state formation absolutely generic phenomena.
The claim of Ref. [31] (see their Sec. 4), that bound state dissociation makes the net effect on the freeze-out irrelevant, was refuted in Refs. [32,35]. This was done through developing an analytic asymptotic solution to the system of Boltzmann equations arising from all allowed bound states [32,43]. 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 σ I v rel and annihilation rate Γ 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 [32]. 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 Γ break n I = σ I v rel n 2 DM . 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 DM = n DM /s: where λ = gSMπ 45 M Pl M DM , z = M DM /T and with BR(B i → SM) = Γ ann Γ ann + Γ break 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, with the inverse temperature at freeze-out M DM /T f = z f given by the transcendental equation For multi-TeV DM, z f ≈ 25 is typical. In Fig. 1, we demonstrate the effects of the Sommerfeld enhancement and bound-state formation on freezeout. Including the Sommerfeld effect leads to additional attraction among WIMPs and enhances the annihilation rate, which in turn reduces the relic abundance by O(10). The consideration of bound states is an additional effective annihilation channel and leads to a further 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 α 5 M DM , where α 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 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.

B. Effects on the Unitarity Bound
As discussed in the classic paper of Griest and Kamionkowski [27], conservation of probability limits the reaction cross section of DM annihilating to any final state for each partial wave by Note the scaling of the bound with v −1 rel , 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 Effects on freeze-out due to the Sommerfeld effect alone and the additional effects of bound-state formation. The inset shows the qualitative behavior at the time of deviation from the thermal DM abundance. Note in particular, that the DM depletion due to bound-state formation (green lines) sets in at later times than the Sommerfeld enhanced freeze-out. In particular in the case indicated by the dot-dashed green line, where the smaller bound-state annihilation rate of Γ ann ≈ 10 −7 M DM leads to a belated annihilation. This is a direct consequence of the branching ratio introduced in Eq. (3). first discuss the cross sections that are relevant for the physical system. In the following, (σv rel ) denotes nonaveraged cross sections and σv denotes thermally averaged cross sections. The total (inelastic) reaction cross section is (σv rel ) total = J (σv rel ) J total . The total reaction cross section is composed of an annihilation part and the bound-state formation cross section (σv rel ) total = (σv rel ) ann + I (σ I v rel ) BSF . The relevant quantity for the freeze-out, as we have shown, is (σv rel ) eff = (σv rel ) ann + I (σ I v rel ) BSF BR(B I → SM) ≤ (σv rel ) 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. [27], the total reaction cross section is approximated as (σv rel ) total ≈ (σv rel ) ann and taken for the freezeout 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. Thus, in any perturbative physical system, saturating the unitarity bound on (σv rel ) J total , the inequality (σv rel ) eff ≤ (σv rel ) 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. [27] is the anni-hilation 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, since by momentum conservation class ≈ p DM R. Note that in the case of the geometrical interaction the cross section σ is constant and the quantity relevant for the freeze-out scales as σv rel ∝ v 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 zerotemperature 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 poin-tlike 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 DM > M SM /v rel and M DM > M SM /α, where M SM is the mass of a weak scale particle and α the strength of the interaction. Explicit examples have been studied in Refs. [32,44].
To parametrize the explicit model, we use SU (2) symmetry as the guiding principle. We consider a particle in the SU (2) L representation R. To have a neutral component, the hypercharge has to be 0 for odd R and 1/2 for even R. For the heavy freezeout computation, the hypercharge gives just a subleading effect and R ≥ 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 χ = 2R. This model, even though specific in terms of SU (2) L quantum numbers, serves as an excellent template to parametrise 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 with σ 0 = πα 2 2 /M 2 DM [45].

A. Sommerfeld Enhancement
The annihilation rate is enhanced by the Sommerfeld effect, which decomposes into isospin channels as σ S,ann σ ann = where the weight factor f I can be found by explicit computation and is listed in Table I, the S ann (λ) are the Sommerfeld factors; see Ref. [46] for their explicit form. The λ 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 Since among these representations the singlet and adjoint channels have the smallest quadratic Casimirs the corresponding channels are the most attractive ones, explicitly λ 1 = (R 2 − 1)/4 and λ 3 = (R 2 − 5)/4 and, sub-dominantly (for R ≥ 3), λ 5 = (R 2 − 12)/4. The larger representations are less attractive, or even repulsive with an effective Sommerfeld suppression factor, similar to Ref. [47].

B. 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. [32], the wave function symmetry leads to selection rules, such that the singlet bound state has spin-0 for even and spin-1 for odd and the adjoint bound state has spin-0 for odd and spin-1 for 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 ≥ 3, from the 5 configuration.

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, For example, in the non-relativistic limit, the capture cross section to the ground state can be written as where S is the spin of the bound state and C J and C τ are group theory factors for the abelian and non-abelian emission, respectively. The full expressions, including the excited states, can be found in Ref. [32]. The group-theory factors depend on the transition and have in general a complex structure, see Refs. [32,48] for the full expressions. In the cases of the transitions we are interested in, however, they can be rewritten in a general simple way: With these factors and the λ 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 ≥ 3) the 5 state, which makes it relevant for larger representations. As discussed in Ref. [32], 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. [49]. The resummation leads to a logarithmic enhancement of the capture rate with ≈ (1 + κ ln (α eff /2v rel )).

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. [32], the excited > 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 pseudoreal representations are only viable as admixtures, as discussed in Ref. [50]. Given the selection rules for identical fermions, the singlet ground-state bound states will have spin-0 and annihilate to two gauge bosons with a rate of 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 The total rate contains the SM fermion multiplicity n f = 3(3 + 1), T SM = 1/2 and η = 2 for self-conjugrate DM and η = 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 ≥ 4, the bound-state formation effects lead to substantial corrections on top of the Sommerfeld corrections.

Partial-Wave Contributions
As already discussed in Ref. [27], the annihilation cross section (σ 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 The green solid line shows the predicted relic density under the assumption that the s-wave cross section saturates the full unitarity bound, for a selfconjugate DM candidate. The green dashed line is the same for the non-self-conjugate case. The red line shows this prediction when bound-state formation is neglected.
contains a significant s-wave fraction, in agreement with [51]. By explicit computation, it can be verified that the amplitudes for higher excited states, (σ BSF v) n , when summed over , 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 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.

C. 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 (σv) 0 total < (σv) 0 max , which is set by partial wave unitarity, Eq. (6), and extract for each given DM mass a value for R 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 DM ≤ 144 TeV (M DM ≤ 96 TeV for a non-selfconjugate particle). We have to stress that the obtained mass bound is subject to a 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. 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 DM ≤ 200 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 [52]. 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.

D. The Abelian WIMP
Previous papers investigated the freezeout of fermionic DM interacting via a dark Abelian force. While Ref. [30] 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 [51].
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 ∆ = 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 α D < 0.86, which is determined only from the Sommerfeld enhanced annihilation cross section. This is the reason, why in [51] 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 twophoton emission adopting the formalism of Ref. [54]. The ∆ = 0 contribution to the bound-state formation in the The allowed cross section range in the early universe. Boundstate formation is suppressed below 5 TeV, since the bound-state binding energy is too low to emit W and Z bosons and the process is kinematically blocked. Thus the J = 0 unitarity bound coincides with the tree-level result. Note that the unitarity limit on the cross section is inversely proportional to the DM velocity and is thus more stringent in the early universe. The total DM annihilation cross section, denoted by the blue line, is also lower in the early universe. Since the cross section contains the 1/v rel factor from the Sommerfeld and bound-state formation effects, the unitarity bounds on the maximally allowed cross section strongly differ between the early universe in (a) and late time annihilation in (b). Today, the annihilation cross section can be significantly enhanced.
ground state is s-wave dominated and given by 1 where φ is the radial part of the free two-particle state wave-function, ψ are the radial parts of bound-state wave-functions, the numerical factorκ ≈ 20 andσ 0 = πα 2 D /M 2 DM . Taking this contribution to the groundstate formation into account reduces the s-wave unitarity bound on the dark gauge coupling to α D < 0.6 and leads furthermore to a reduced maximally attainable DM mass of the order of M DM < 150 TeV. (Compared to the bound of M DM < 200 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.

E. 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 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 rel. . Since the average velocities at freeze-out are of order v rel. ∼ O(0.3), the cross section is close to the well known value of 2.2 × 10 −26 cm 3 s −1 . In contrast, the late time annihilation is significantly enhanced, since the average DM velocity in galaxies today is v rel. ∼ O(10 −3 ). This is indicated in Fig. 4 (b) for all allowed final states.
In Fig. 4 (b) (γ−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 [55,56].
In Fig. 5, 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. 6 are suppressed due the smallness of their mass (chirality suppression, e.g., Ref. [57]) or due to the smallness of their Yukawa couplings to the Higgs, e.g., Ref. [58]. 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 [45].
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. [19,59]. 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. [60]. 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 the Higgsino-like neutralino in supersymmetric models. The general framework for such models has been discussed in Ref. [50]. In Fig. 7, we show the predictions for the spinindependent 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 [61], will cover it entirely.

F. 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. [45]. Since in that case a DM theory needs a lowlying 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. [63], 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 demon-strated that in an explicit construction with vectorlike 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. [64] for the resummation technique.
This opens 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. [27], the unitarity argument can be also applied to the annihilation of extended composite objects, such as DM atoms, see for example Refs. [65,66]. The intuitive reason why their cross section is big is that all partial waves up to l class ≈ p × 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 ∝ 1/α em m e is controlled by the mass of the electron, which is much smaller than the mass of the entire atom.

A. The Simplest Dark Atom
In the case of heavy, multi-TeV DM, the low number density makes hydrogen-like recombination very inefficient [67]. 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 supersymmertic Yang-Mills theory. The minimal example is an SU (3) theory with a fermion Q in the 8 representation, as introduced in Refs. [68,69].
The fermion can have an explicit vector-like mass. Below this mass scale, the running of the gauge coupling α will lead to a strong coupling regime at a scale Λ given by Λ ≈ M Q exp (−6π/11N α(M Q )). That means that at energies below Λ 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 formation of bound states of adjoint fermions and dark gluons B = Qg. 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 ≈ 1/Λ, but its mass is controlled by M Q Λ and its de-Broglie wavelength is much shorter than its size, as in the case of hydrogen.

B. The Atomic Relic Abundance
Annihilation happens via Qg + Qg → QQ + gg, where the G = gg is the glueball, the lightest degree of freedom in this theory, and the QQ state can self-annihilate efficiently once it has fallen to the ground state. Given that M Q Λ, the deeply bound states of QQ can be treated perturbatively, as heavy quarkonium.
The fermions in the adjoint representation can arrange into a singlet or an adjoint bound state. Higher boundstate 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 QQ state in the 1 ⊗ 8 S ⊗ 8 A representations, arrangements in higher representations are repulsive. Pure combinatorics suggests that the adjoint configurations are formed with a probability p 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 darkmatter 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 Qg Assuming that for the binding energies of the non-perturbative states, the energy of the gluon field in the octet configuration of QQ and the glueball mass M G are of the order of Λ, we have the condition for the binding energy between the heavy quarks E QQ B < E CMS kin ≈ T < Λ, as we consider the postconfinement 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, = E Coulomb + E Linear ≈ α 2 eff M Q /4(−1/n 2 + 12 n 2 ), with [70]. Here the potential part, dominated by the string tension σ, is approximated by λ R Λ 2 . The combination of the representations' quadratic Casimir operators λ R = 1/2 (2C adj. − C R ) defines the strength of the channel.
The above energy-conservation condition (E n, < Λ) demands that the correction from the linear potential is not dominant and it thus limits the maximal partial wave by After deriving the general bound max from energy conservation, we come now to the momentum conservation condition. The maximal allowed by classical momentum conservation is where 1/Λ is the size of the Qg bound state. Since class scales with the temperature, there is a T c below which class < 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 = Λ λ R α 3/2 /8 √ 3. This leads to a compact expression for the maximal partial wave allowed by energy conservation, 2 max = 4 T c M Q /Λ 2 . The cross section for the reaction for the two quark state in representation R decomposed in partial waves reads where we define L max (T ) = min[ max , class ] which is set by either energy or momentum conservation respectively. P ann is the probability that a bound state with given loses energy fast enough, falls to the ground state and self-annihilates. For this cross section, unitarity is saturated if all phases have sin(δ ) = 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 P ann = 1 for all of them. In fact, in Ref. [68], 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 < Λ: with σ geom. = 4π/Λ 2 . This implies for the thermally averaged cross section, It is obvious that the result of Ref. [27] 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 Integrating over both scaling regimes we get the asymptotic relic abundance 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.

C. Implications of Geometrical Unitarity
In Fig. 8 only if M Q ≈ Λ. In this regime however, the dark atom mass is set by the scale M DM ≈ Λ and thus R DM ≈ 1/Λ ≈ 1/M DM , which implies that 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 max 1, we are in the regime with M Q Λ and R DM ≈ 1/Λ 1/M DM . We find that in this case the unitarity bound depends on the mass M Q and results in a larger DM radius corresponding to a monotonically decreasing Λ with growing M Q .
We demand that the scale Λ > Λ 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 TeV. This contrasts with the finding of Ref. [27] that the limit on the DM size is independent of its mass. So, in the case of the lowest possible confinement scale Λ ≈ Λ 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 non-abelian charge. For a detailed discussion, see Ref. [68].
In Fig. 9, 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 Λ M 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 set-ting the relic abundance is smaller given the same model parameters.
When Λ ≈ M 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 σ ≈ constant instead of σv ≈ 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. 9, we compare the total annihilation cross section to the HAWC experiment bounds from Ref. [71] assuming conservatively a dominant annihilation into neutrino final states, as suggested in Refs. [19,72,73] 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 [74], which was used in Ref. [71]. Note that the predicted large annihilation cross sections at large DM masses are compatible with the excess of high energy neutrino events found in Ice-Cube [75]. We point out that composite dark matter, which saturates the unitarity can have intriguing direct detection signatures, as discussed in [76,77].

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 [20]), 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. Be-cause those are relatively lighter than the heavy WIMP, they act as light mediators, which brings in Sommerfeldenhancement effects, which are well studied, and boundstate effects, which are much less so. We take advantage of the new formalism and results of Refs. [30-32, 35, 50], 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: 1. 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 darkmatter 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. The late-time annihilation of WIMPs is significantly enhanced due to Sommerfeld and boundstate 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: 1. 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 [27]. 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. 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. [71,75,78]. 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 [79], LHAASO [80], CTA [81]), cosmic-ray (AMS [82], CALET [83], DAMPE [84]), and high-energy neutrino experiments (IceCube [85], Antares [78]) provide a very promising route to test annihilations of TeV scale elementary and composite dark matter.