Closing the window on WIMP Dark Matter

We study scenarios where Dark Matter is a weakly interacting particle (WIMP) embedded in an ElectroWeak multiplet. In particular, we consider real SU(2) representations with zero hypercharge, that automatically avoid direct detection constraints from tree-level Z-exchange. We compute for the first time all the calculable thermal masses for scalar and fermionic WIMPs, including Sommerfeld enhancement and bound states formation at leading order in gauge boson exchange and emission. WIMP masses of few hundred TeV are shown to be compatible both with s-wave unitarity of the annihilation cross-section, and perturbativity. We also provide theory uncertainties on the masses for all multiplets, which are shown to be significant for large SU(2) multiplets. We then outline a strategy to probe these scenarios at future experiments. Electroweak 3-plets and 5-plets have masses up to about 16 TeV and can efficiently be probed at a high energy muon collider. We study various experimental signatures, such as single and double gauge boson emission with missing energy, and disappearing tracks, and determine the collider energy and luminosity required to probe the thermal Dark Matter masses. Larger multiplets are out of reach of any realistic future collider, but can be tested in future γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}-ray telescopes and possibly in large-exposure liquid Xenon experiments.


Introduction
The possibility that Dark Matter (DM) is a new weakly interacting massive particle (WIMP), thermally produced in the early Universe and freezing out through 2 → 2 annihilations into Standard Model (SM) states, remains one of the a e-mail: salvatore.bottaro@sns.it (corresponding author) main motivations for new physics in the 10 GeV-100 TeV range. Under these simple assumptions, the lower bound on the WIMP mass comes from astrophysical constraints on DM annihilations into SM products [1], while the upper bound is a consequence of s-wave unitarity of the DM annihilation cross-section [2].
A particularly interesting possibility within this framework, because of its minimality and predictive power, is that the DM is the lightest neutral component of one electroweak (EW) multiplet. In particular, fermionic and scalar n-plets of SU(2) with odd n and zero hypercharge automatically avoid strong constraints from direct detection searches, and will be taken here as a minimal realization of the EW WIMP scenario. The lightest particle in any such representation can be made stable by enforcing a symmetry acting on the DM only (for multiplets with n ≥ 5 such a symmetry arises accidentally in the renormalizable Lagrangian). However, we shall see that in general this can require additional assumptions about the completion of the theory at some high UV scale.
The main purpose of this paper is to precisely determine the WIMP freeze-out predictions in a systematic way. For any given n-plet, computing the EW annihilation cross-section in the early Universe allows to infer the WIMP cosmological abundance. By requiring it to match the measured value of the DM abundance today, DM h 2 = 0.11933 ± 0.00091 [3], the mass of the n-plet can be univocally determined. These mass predictions are an essential input to assess if and how the future experimental program will be able to fully test the EW WIMP scenario. In contrast to previous papers on the subject [4][5][6][7], our approach here is to minimize the theory assumptions and fully classify the calculable freeze-out predictions. Because of its infrared-dominated nature, the calculability of freeze-out depends purely on the partial wave unitarity of the total annihilation cross-section  [2], which we re-analyze here for EW n-plets. All in all, demanding perturbative unitarity requires n ≤ 13 for both bosonic and fermionic DM. Approaching this boundary the theory uncertainty on the mass prediction grows as shown in Fig. 1. Stronger constraints on n can be imposed by demanding the EW interactions to remain perturbative up to scales well above the thermal DM mass.
The effects of Sommerfeld enhancement (SE) and of bound state formation (BSF) are known to significantly affect the freeze-out predictions and need to be included. The first effect has long been recognized to lead to an enhancement of the annihilation cross-section at small relative velocities [8][9][10][11]. The effects of BSF for WIMP freeze-out have been first computed in Ref. [12] for the n = 5 fermionic multiplet (see Refs. [13,14] for earlier computations in other contexts).
Here we extend their treatment to fermionic and scalar representations of arbitrary high n, up to the break-down of perturbative unitarity. At growing n, we find that bound states (BS) are more tightly bound, with their ionization rate being exponentially suppressed. At the same time, the multiplicity of accessible BS channels grows significantly. These two effects result in an increase of the annihilation cross-section compared to the estimates of Ref. [15].
The freeze-out mass predictions are summarized in Table 1 and Fig. 1 for the real n-plets considered here. With masses ranging from several TeV to tens or hundreds of TeV, most of the EW WIMP candidates are still out of reach of present experiments, but could be tested in the future, thanks to the forthcoming progress in collider physics and DM detection experiments. With the mass predictions at hand, we thus commence a systematic survey of the WIMP phenomenology: (i) at very high energy lepton colliders with 10-30 TeV center of mass energy [16,17]; (ii) at direct detection experiments with 100 tons/year of exposure like DARWIN [18,19]; (iii) at high-energy γ -ray telescopes like CTA [20][21][22][23]. We first examine the reach of a hypothetical future muon collider, studying in detail for which values of center-of-mass energy and integrated luminosity the EW 3-plets and 5-plets can be fully probed through direct production. We instead find direct production of the EW multiplets with n > 5 to be beyond the reach of any realistic future machine (this is in contrast with the results of the recent study [24] due to the increase of the thermal mass of the 7-plet with the inclusion of BSF effects). These larger n-plets are possibly within the reach of large exposure direct detection experiments, and will probably be tested more easily with future high energy γ -ray telescopes. A careful study of the expected signals in indirect detection is left for a future work [25].
This paper is organized as follows. In Sect. 2 we summarize the EW WIMP paradigm, in Sect. 3 we illustrate the main features of our freeze-out computation, and in Sect. 4 we discuss the unitarity bound assessing the theory uncertainties. These three sections provide a full explanation on the results of Table 1 and Fig. 1. In Sect. 5 we discuss the implications of our study for a future muon collider, while in Sect. 6 we briefly re-examine the reach of direct and indirect detection experiments in light of our findings. In Appendix A we give further details on the nature of next-to-leading order corrections and we detail the BS dynamics for the 7-plet. Appendix B contains further information on the collider studies.

Which WIMP?
We summarize here the logic of our WIMP classification very much inspired by previous papers on the subject [4][5][6][7]27]. Requiring the neutral DM component to be embedded in a representation of the EW group imposes that Q = T 3 + Y , where T 3 = diag n+1 2 − i with i = 1, . . . , n, and Y is the hypercharge. At this level, we can distinguish two classes of WIMPs: (i) real EW representations with Y = 0 and odd n; (ii) complex EW representations with arbitrary n and Y = ± n+1 2 − i for i = 1, . . . , n. Here we focus on the first class of WIMPs, which is particularly interesting because the DM does not couple to the Z -boson at tree level, avoiding strong constraints from direct detection experiments. Other possibilities will be discussed elsewhere.
At the renormalizable level, the extensions of the SM that we consider are Freeze-out mass predictions for WIMP DM in real EW multiplets with Y = 0. The annihilation cross-section includes both the contribution of SE and BSF. We provide a measure of how close the DM annihilation cross-section is to the unitarity bound for s-wave annihilation (σ v) J =0 max = 4π/M 2 DM v. Approaching the unitarity bound, the error on the WIMP mass grows proportionally to the enhancement of the next-to-leading order (NLO) contributions estimated in Eq. (23). We derive the scale where EW gauge coupling will develop a Landau pole by integrating-in the WIMP multiplet at its freeze-out mass. The stability of both scalar and fermionic DM can always be enforced by requiring a Z 2 symmetry in the DM sector to forbid DM decays. This symmetry forbids the scalar and fermionic 3-plets decay at renormalizable level as indicated by the *. The value of the UV cut-off UV gives an idea of the required quality for this symmetry to make DM stable and avoid stringent bounds on decaying DM (τ DM > 10 28 s) [26]: a new physics scale lower than UV would require a Z 2 to explain DM stability, while a cut-off higher than UV would make DM stability purely accidental DM spin EW n-plet for scalars and fermions, respectively, where D μ = ∂ μ − ig 2 W a μ T a χ is the covariant derivative, and T a χ are generators in the n-th representation of SU (2). The Lagrangian for the real scalar in Eq. (1) also admits quartic self-coupling and Higgs-portal interactions at the renormalizable level. The latter is bounded from above by direct detection constraints (see Fig. 8 right) and gives a negligible contribution to the annihilation cross-section. 1 The neutral component and the component with charge Q of the EW multiplet are splitted by radiative contributions from gauge boson loops. In the limit m W M DM these contributions are non-zero and independent on M χ . This fact can be understood by computing the Coulomb energy of a charged state at distance r 1/m W or the IR mismatch (regulated by m W ) between the self-energies of the charged and neutral states. The latter can be easily computed at 1-loop [28][29][30], with the uncertainty dominated by 2-loop contributions proportional to α 2 2 m t /16π . These have been explicitly computed in Refs. [31,32] giving a precise prediction for the lifetime of the singly-charged component, which decays to the neutral 1 No other quartic coupling is allowed since χ T a χ χ identically vanishes. Indeed, (T a χ ) i j is antisymmetric in i, j, being the adjoint combination of two real representations, while χ i χ j is symmetric. one mainly by emitting a charged pion with where 2T + 1 = n. The suppression of the lifetime with the size of the EW multiplet can be understood in the M χ m W limit where the mass splitting between the charged and neutral components is independent of n while the coupling to W is controlled by √ T (T + 1)/2. As we will discuss in Sect. 5.2, the production of a singly charged DM component at colliders gives the unique opportunity of probing EW multiplets with n = 3 and n = 5 through disappearing tracks [4,24,[33][34][35].
Interestingly, the IR generated splitting from gauge boson loops is not modified substantially by UV contributions. The latter are generated only by dimension 7 (dimension 6) operators if the DM is a Majorana fermion (real scalar) and can be written as with n I = 3, 2 for I = f, s. This corresponds to a splitting M I c I v 4 / n I UV M 3−n I χ which is always negligible with respect to the residual error on the 2-loop splitting for UV 100 TeV and c I ∼ O (1).
We now move to discuss DM stability. In the case of the EW 3-plet, the renormalizable operators χ H † H and χ H L, for scalars and fermions, respectively, can induce fast DM decay. We assume these operators to be forbidden by a sym-metry (e.g. a discrete Z 2 -symmetry) acting only on the DM sector. For all the other n-plets with n ≥ 5, instead, Z 2 -odd operators are accidentally absent at renormalizable level.
Higher dimensional operators that break the Z 2 -symmetry are in general expected to be generated at the ultraviolet cutoff scale UV . We sketch here the operators of lowest dimension that can induce the decay of scalar and fermionic WIMPs for generic n: where SU(2) contractions are implicit, and the dots indicate operators of the same dimension with different combinations of W and H fields. 2 Higher-dimension operators with additional SM fields or derivatives are of course also possible. The first operators in the two equations above are just the renormalizable operators of the 3-plet case "dressed" with extra Higgs insertions. The dominant contribution to the decay width at tree-level always comes from the operator with the highest number of W insertions (namely (n − 3)/2 for fermions and 2 (n − 1)/4 for scalars). Notice that for fermionic DM, dipole-like operators with an odd number of W fields can always be constructed. In the last operator in both Eqs. (6) and (7), χ 3 is the unique isospin triplet constructed out of three SU(2) irreducible representations of odd isospin [27,36]. These operators contribute to the WIMP decay at one-loop as discrete symmetries in quantum gravity [37]. For fermionic representations, DM decay is instead induced by dimension 6 operators for n ≤ 5, and dimension 7 operators for n > 5, and the DM stability can be determined within quantum field theory.
A lower bound on UV is obtained by requiring the DM lifetime to be long enough to circumvent cosmological bounds [38,39] (τ DM 10 19 s) or astrophysical bounds on the decay products of decaying DM [26,40,41] (τ DM 10 28 s). We can then quantitatively measure the required quality of the Z 2 -symmetry by considering the ratio between the minimal UV allowed by the constraints and the WIMP freeze-out mass. A naive dimensional analysis (NDA) estimate of UV , assuming all the Wilson coefficients to be O(1), is given in Table 1 for all the relevant n-plets.
Requiring perturbativity of the EW gauge coupling above the WIMP thermal mass can provide an upper bound on the dimension of the SU(2) representation. Indeed, large SU(2) n-plets will make the EW gauge coupling run faster in the UV, eventually leading to a Landau pole. In Table 1 we provide the value of the scale Landau such that g 2 ( Landau ) = 4π . We integrate the RGE equations for the SM gauge couplings at 2loops and integrate-in the n-plet at the WIMP thermal mass. 3 Comparing Landau and UV , we see that the stability of the fermionic n-plets with n ≤ 5 only depends on physics in a regime where the EW coupling is still perturbative. Instead, the stability of n-plets with n > 5 requires specifying a UV completion for the EW gauge group that does not give rise to the dangerous operators of Eqs. (6) and (7). In this sense, the Majorana 5-plet studied in Ref. [4] is special, because it can be made accidentally stable by raising the scale UV , without any further assumption on the nature of the UV completion at Landau . Requiring UV /M χ 10 to ensure perturbativity of the theory up to well above the WIMP mass would select n ≤ 9 for fermions, and n ≤ 11 for scalars. However, requiring a large hierarchy between Landau and M χ is not necessary to ensure the calculability of thermal freeze-out, which depends only on EW processes at energies much below the DM mass. A more robust upper bound on the dimension of the SU(2) nplets will be derived in Sect. 4, analyzing the s-wave unitarity of the annihilation cross-section. This bound will require n ≤ 13 for both fermionic and scalar WIMPs.
Finally, let us comment on the EW WIMPs in complex representation of SU (2). For odd n, complex multiplets with Y = 0 are allowed. Their freeze-out dynamics shares many similarities to the one of the real multiplets discussed here and has been partially discussed in Ref. [27]. For complex representations with Y = 0, direct detection constraints can be circumvented only by introducing a splitting between the two Weyl spinors forming the Dirac pair. The required splitting can be generated via dimension 5 operators above the WIMP thermal mass leaving the freeze-out predictions unaffected. A full classification of these WIMP scenarios and their phenomenological probes is left for a future work [42].

WIMP cosmology
The determination of the DM thermal mass hinges on a careful computation of the DM annihilation cross-section in the non-relativistic regime. In particular, the potential generated by EW gauge boson exchange between DM pairs is attractive for isospins I √ 2n resulting into Bound State Formation (BSF) through the emission of an EW gauge boson in the final state. The energy of the emitted gauge boson is of the order of the Bound State (BS) binding energy α eff is the effective weak coupling defined in Eq. (17), and we neglected corrections of order m 2 W /M 2 χ . In the nonrelativistic limit, and at leading order in gauge boson emission, the BSF process is encoded in the effective dipole Hamiltonian described in Refs. [12,43] which dictates the BS dynamics and it is written for completeness in Appendix A. The BS dynamics relevant for DM freeze-out is well described by the unbroken phase of SU(2) so that the configuration of the DM pair can be decomposed into eigenstates of the isospin I of the pair where C(I I z |i j) are the Clebsch-Gordan coefficients and I is the dimension of the isospin representation. Denoting with L and S the total angular momentum and the spin, the isospin-Lorentz structure of the dipole Hamiltonian enforces the following selection rules: (i) S = 0 because the dipole Hamiltonian is spin-independent; (ii) | L| = 1 because the dipole operator transform as a vector under rotations; (iii) | I | = 2 because a single, G-parity odd weak boson is emitted.
Since we are dealing with real representations, spinstatistics imposes further restrictions on the allowed quantum numbers, depending on the fermionic or scalar nature of the wave function. In particular we have which implies that for scalars n B s (n B p) bound states, i.e. with L = 0 (L = 1), can exist only with even (odd) I −1 2 , while for fermions odd (even) I −1 2 states with L = 0 are forced to have S = 1 (S = 0).
We are now ready to describe the system of coupled Boltzmann equations for the evolution of the number densities of DM and BS. Following [12], we will discuss how this coupled system can be reduced to a single equation for the DM number density with an effective annihilation cross-section. The Boltzmann equations for DM and BS read where B I,J,... labels the different bound states, z = M χ T , s is the entropy density and Y = n s is the number density per co-moving volume.
The dynamics of a given BS B I in the plasma is described by Eq. (12b) and depends on: (i) its ionization rate B I ,break ; (ii) its annihilation rate into SM states B I ,ann ; (iii) its decay width into other bound states B I →B J . The ionization rate B I ,break ≡ n γ σ I ,break v rel encodes the probability of a photons from the plasma to break the BS B I . Assuming thermal equilibrium, detailed balance relates the cross-section for the BS breaking σ I ,break v rel to the BSF cross-section σ B I v rel where g B I and g χ count the number of degrees of freedom of the bound state B I and of the DM multiplet, respectively. If either the BS decay or the annihilation rate satisfies H , we can neglect the LHS in Eq. (12b), obtaining algebraic relations between the DM and the BS yields.
Plugging these relations into Eq. (12a), we arrive at the final form of the DM Boltzmann equation where and we defined the effective cross-section as the sum of the direct annihilation processes, S ann , and the ones which go through BSF, S B J . In particular, S ann can be written as where σ I ann is the hard cross-section for a given isospin channel I , S I E is the Sommerfeld enhancement (SE) of the Born cross-section, and v rel is the relative velocity of the two DM particles. In the limit of small relative velocity between the DM particles (but larger than m W /M χ ), the SE factor can be approximated as The finite mass effects modify the behavior of the SE at v rel m W /M χ and are included in our full computation (see Ref. [11] for explicit formulas). However, Eq. (17) will be enough to estimate the behavior of the SE at the temperatures most relevant for freeze-out. Analogously we can factorize the BSF processes as where S I,l B J is the "hard" BSF cross-section of the state B J starting from a free state with angular momentum l and isospin I multiplied by the SE factor of that particular isospin channel as defined in Eq. (17). Explicit expressions for this can be found in Refs. [12,43]. R B J gives instead the effective annihilation branching ratio into SM states which depends on the detailed BS dynamics (i.e. annihilation, ionization and decay). In particular, R B J approaches 1 once the temperature of the plasma drops below the binding energies of the bound states involved in the decay chains. In the case of a single BS, R B J takes a rather intuitive form which applies to 1s I and 2s I BS with I ≤ 5. The latter, once formed, annihilate directly into pairs of SM vectors and fermions, with rates ann α 5 eff /n 2 B M χ . These BS together make up for more of the 50% of the BSF cross-section. More complicated examples of BS dynamics will be illustrated in Appendix A.2 where we detail the case of the EW 7-plet.
While the effect of BSF has already been computed for the fermionic 5-plet in Ref. [12], here we include it for the first time for all the real WIMP candidates with n ≥ 7. For larger EW multiplets, we find the relative effect of BS dynamics on the total cross-section increases, as can be seen from Fig. 2. This is the consequence of two effects: (i) the binding energy grows at large n, suppressing the ionization rate with respect to the annihilation one; (ii) at larger n the number of attractive channels increases and thus the BS multiplicity per energy level grows linearly with n. For example, for n = 5 the attractive channels have I = 1, 3, 5, for n = 7 BS with Fig. 2 Effective cross-section for BSF normalized over the total annihilation cross-section as a function of z = M χ /T assuming vanishing ionization rates, i.e. R BS = 1 (see Eq. (18) and below). The dashed lines for the fermionic 5-plet (dark blue) and 7-plet (cyan) show the deviation of the real bound state dynamics from the approximation of vanishing ionization rates. For n > 5 the error due to the R BS = 1 is subdominant compared to the virtual and real effects at NLO in gauge boson emission I = 7, 9 can also form. The relevance of these higher isospin channels was not recognized in [15], where only the I = 1, 3 channels were included, significantly underestimating the thermal mass already for n = 7. In Appendix A.2 we show explicitly the relative contributions coming from the different isospin channels for the 7-plet. The 7-plet thermal mass was computed including all the BS up to 3s and 2 p but we checked that the contribution from 4s and 3 p BS is negligible.
As we increase the dimension of the multiplet, the bound states become more tightly bounded and the effect of the ionization rate becomes smaller. This can be explicitly seen from Eq. (13) where the binding energy controls the Boltzmann suppression of the ionization rate. For this reason, we only account for the detailed BS dynamics for n ≤ 7 while for n > 7 we set the annihilation branching ratios to 1. We assume, as explicitly checked for the 7-plet, that the formation cross sections for 4s and 3 p BS are negligible. In fact, the cross sections of BS differring only for their principal quantum number have the same parametric dependence on n, so that the hierarchy between different energy levels is independent on n. Close to the unitarity bound limit, excited states with larger angular momentum can become important. However, their long lifetimes and small binding energies limit their contributions to the thermal mass. Moreover, since the typical velocity inside the bound state is α eff /n B , relativistic corrections can also be important. We leave the discussion of these contributions to a future work.
In Appendix A.2 we estimate the error on the WIMP mass due to this approximation by comparing its effect on the thermal masses of 5-plet and the 7-plet against the full computa-tion. We find a shift in mass M DM 5 TeV for both n = 5 and n = 7 resulting in a smaller relative error for n = 7, as expected. We keep 5 TeV as an estimate of the error induced by this approximation for the larger multiplets. As we will discuss in the next section, the uncertainty for n ≥ 7 will be anyhow dominated by the next-to-leading order (NLO) contributions to the SE which are not included here.
Finally, we comment on the theory uncertainty on the mass prediction for the 5-plet. This is dominated by the approximate treatment of EW symmetry breaking effects in computation of the BSF cross-sections. The SU(2)-symmetric approximation fails once the DM de Broglie wavelength becomes of the order of m W (i.e. for z 10 4 for n ≥ 5). After the EW phase transition, Coulomb and Yukawa potentials appear at the same time so that employing either the Coulomb or the Yukawa centrifugal correction to the SE (see Ref. [11]) overestimate and underestimate, respectively, the freeze out cross-section. This gives us a rough way of determining the theory uncertainty: (i) to set the lower bound on the freeze-out mass we include BSF in σ eff until z = 10 4 with the centrifugal correction coming from the Yukawa; (ii) to set the upper bound we push the effect of BSF, neglecting the vector masses in the centrifugal correction, to arbitrary large values of z. We observe that the abundance saturates already for z ≈ 10 5 . This procedure gives the uncertainty for the 5-plet in Table 1 which is different than the one quoted in Ref. [12], where the BS contribution was switched off at z = 10 4 , underestimating the effect of BSF.

The WIMP unitarity bound
We now analyze the constraint of perturbative unitarity on the annihilation cross-section, including bound state formation. The perturbative unitarity of the S-matrix sets an upper bound on the size of each partial wave contribution to the total annihilation cross-section 4 where J = L + S is the total angular momentum. The stronger inequality comes from the s-wave channel (i.e. J = 0) which can be written as where f 0 B i selects the BS contributions that can be formed by J = 0 initial wave, which are limited by the selection rules discussed in the previous section.
For a scalar WIMP selecting the s-wave implies L = 0, and only BS in p-orbitals can contribute to the s-wave crosssection with f 0 BS = 1. The spin statistics of the wave function in Eq. (11) forces these BS to have odd (I −1)/2. In practice, the s-wave unitarity bound for scalars is determined solely by the SE. For fermionic WIMP selecting the s-wave implies the same selection rules of the scalar when S = 0. Additional contributions arise from S = 1 s-orbital states, whose isospin must be odd due to Fermi statistics. In this case, the projection onto the J = 0 wave gives f 0 BS = 1 9 . Solving the constraint in Eq. (21) we find that s-wave unitarity is violated for n ≥ 15 for both fermion and scalar WIMPs. In both cases the s-wave cross-section is largely dominate by the SE. We checked that a similar constraint can be obtained by looking at the p−wave unitarity, where the cross-section is instead dominated by the formation of 1s BS.
The selection rules that regulates the BS dynamics derive from the dipole Hamiltonian which is written for completeness in Eq. (A.1). These selection rules are only broken by NLO contributions in gauge boson emission which can be estimated as where the extra α 2 eff correctly accounts for the phase space suppression in the limit of small velocities as detailed in Appendix A.1. As a result, the LO selection rules apply all the way till the breaking of perturbative unitarity.
Interestingly, the upper bound on n from perturbative unitarity derived from Eq. (21) is significantly stronger than the one derived from the perturbative unitarity of the Born cross-section which is violated for n ≥ 38 (i.e. α eff ≥ 4π ). This suggests that because of SE, the ratio between the NLO and the LO cross-section should appreciably deviate from the NDA scaling of the Born cross-section: σ NLO Born /σ LO Born ∼ α eff /4π . Estimating the NLO correction to the potentials controlling the SE we indeed get where the NLO potential is resumming ladder diagrams like the ones in Fig. 3, and where we substituted the de Broglie length 1/Mv rel ≈ √ z/M χ as the typical length scale for the annihilation process. Our estimate above matches the explicit NLO computation of the SE for the 3-plet in Ref. [47]. Requiring this correction to be 1 across the freezeout temperatures leads to a similar upper bound on n than the one inferred from perturbative unitarity. We use the estimate above to assess the theory uncertainty on the WIMP thermal masses in Table 1. Indeed, Eq. (23) results in a correction to the Sommerfeld factor S E , which affects both S ann and S B J as introduced in Eq. (14). We find that neglecting the NLO contribution dominates the DM mass theory uncertainty for n ≥ 7. The uncertainty grows as we increase the dimensionality of the multiplet becoming as large as O(30%) for n = 13.
Finally, we compare our results to the ones obtained in Ref. [15]. Numerically, the upper bound on the WIMP mass corresponding to the saturation of the unitarity bound is roughly 500 ± 200 TeV, which is the expected thermal mass for n = 15 as can be seen from Fig. 1. The unitarity boundary was set instead to 150 TeV for n = 13 in Ref. [15] without a quoted theory uncertainty. Beside the numerical differences, our computation differ from the one in Ref. [15] in two crucial instances: (i) at large n we find that large isospin channels enhance significantly the BSF cross-section making the WIMP DM mass heavier than in Ref. [15] at fixed n; (ii) we find that including BSF does not accelerate by much the saturation of the unitarity bound because of the selection rules of the dipole Hamiltonian at LO. As we discussed above, the LO selection rules are not lifted by NLO corrections until the boundary of perturbative unitarity is reached. These two effects together push the heaviest calculable WIMP mass very close to the PeV scale appreciably enlarging the EW WIMP scenarios beyond the reach of any realistic future collider.

WIMP at high energy lepton colliders
We now look at the possible detection strategies for direct production of WIMPs at collider experiments. From the results in Table 1 one can immediately see that DM masses 50 TeV are required to achieve thermal freeze-out for EW multiplets with n > 5. Pair-production of these states would require center-of-mass energies exceeding 100 TeV, which are unlikely to be attained at any realistic future facility. On the other hand, multiplets with n ≤ 5 have thermal masses in the few TeV range, potentially within the reach of present and future colliders.
Direct reach on these dark matter candidates at hadron colliders is limited by the absence of QCD interactions for the DM candidates, which can be produced only via electro-weak interactions. As such the limits at the LHC (see e.g. [48]) are rather far from the interesting thermal mass targets and only a future pp collider may have the reach for some low-n candidates if collisions around 100 TeV can be attained [34,49,50]. Lepton colliders tend to have reach mainly through indirect effects, e.g. the modification of the angular distributions in simple ff production at center of mass energies below the threshold to produce the DM pair. The reach in this case is up to masses a factor a few above the center of mass energy [51,52].
A very-high-energy lepton collider, such as a muon collider, would be the perfect machine to hunt for these WIMPs, due to its large center-of-mass energy, relatively clean collision environment, and the capability of pair-producing weakly interacting particles up to kinematical threshold. Here we consider in particular a future muon collider with center-of-mass energy of 10 TeV or more and the baseline integrated luminosity of [16] L 10 ab −1 · √ s 10 TeV While such a machine is currently not feasible, various efforts to overcome the technological challenges are ongoing. Early developments on machine performances [53,54] found the luminosity Eq. (24) to be achievable for √ s 6 TeV, and further development to push it to larger energies is currently in progress [55].
We consider various search channels for EW 3-plets and 5plets, and determine the minimal center-of-mass energy and luminosity required to directly probe the freeze-out predictions. First, we detail in Sect. 5.1 the prospects for the observation of DM as undetected carrier of momentum recoiling against one or more SM objects. We systematically study all the "mono-V" channels, where DM is recoiling against a SM gauge boson V = γ, Z , W . We also investigate double vector boson production, that we dub "di-V" channels, where requiring a second SM gauge boson in the final state could help ameliorating the sensitivity. Second, in Sect. 5.2 we study the reach of disappearing track searches -which are robust predictions of WIMPs in real EW representations as discussed in Sect. 2 -recasting the results of [35]. Notice that our study is in principle applicable both to high-energy μ + μ − and e + e − colliders, even though soft QED radiation, beamstrahlung, and the presence of beam-induced backgrounds could affect the results in different ways.
The projections for direct production derived here have to be contrasted with similar studies in the context of future high energy proton machines [33,34] (which are limited by the partial reconstruction of the collision kinematics) or electronpositron machines [56,57] (which are limited by the moderate center-of-mass energy and hence more effective to hunt for lighter DM candidates) .
Complementary studies have also considered indirect probes of WIMPs at future high energy lepton colliders, focusing on the modifications of Drell-Yan processes [52]. Given the freeze-out masses of Table 1, EW n-plets with n > 5 are beyond the reach of any realistic future collider both directly and indirectly, even though a definitive statement about indirect observables would require further studies.

WIMPs as missing momentum
We perform a full study of the different channels to observe DM as undetected carrier of momentum. The generic strategy is to measure a hard SM particle or a set of particles X recoiling against a pair of invisible objects, Notice that we treat all the components χ i of the EW multiplet as invisible, assuming the soft decay products of the charged states to be undetected. Additional soft SM radiation is also implicit in Eq. (25). The prospects for the "monophoton" topology at a future muon collider have been already studied in [24]. Here, we want to extend this analysis by enlarging the set of SM objects recoiling against the invisible DM multiplets. Mono-V. We start by considering "mono-V" scattering processes where V = γ, Z , W is a generic EW gauge boson that accompanies the production of χ states from the n-plet, mono-Z : The main contribution to all these processes comes from initial-and final-state radiation of a vector boson, which have sizeable rates because of the large weak charge of the DM multiplet and the weak charge of the beams. 5 We sum over all components of the multiplet χ i , but the dominant signal corresponds to the production of the state with largest electric charge (i = ±n), subsequently decaying into DM plus soft SM particles.
For each of these signals, the corresponding SM background is dominated by a single process, mono-γ bkg: mono-Z bkg: mono-W bkg: where the missing transverse momentum is carried by neutrinos; the mono-W background also requires a lost charge along the beam. We simulate signal and background events with Mad-Graph5_aMC@NLO [59,60], for different DM mass hypotheses and different collider energies. The W and Z bosons are assumed to be reconstructed from all their visible decay products and are treated as single objects. We impose basic acceptance cuts on the rapidity and transverse momentum of the vectors, requiring |η V | < 2.5 and p T,V > 10 GeV. Other detector effects are neglected.
We then perform a cut-and-count analysis, estimating the significance of the signal as where S, B are the numbers of physical signal and background events, and sys parametrizes the systematic uncertainties. The signal is isolated from the background employing the kinematics of the visible object, parametrized in terms of its transverse momentum p T,V , its pseudo-rapidity η V , and the missing invariant mass (MIM) which is a function of the energy of the visible particle itself We where the p T and η selection cuts are chosen to maximize the significance for each value of M χ . The precise values of the selection cuts, together with the expected number of events and the reach of the various search channels, are given in Table 2 in the Appendix.
The background rates for mono-γ and mono-Z are very similar, with fiducial cross-sections of around 3 pb that depend weakly on the collider energy. As already pointed out in [24] for the mono-γ case, the optimal reach on M χ is obtained for low signal-to-noise ratios -in other words, systematic uncertainties could be important. For this reason, we present results for different values of sys = 0, 1 0 / 00 , 1%. We point out that in presence of larger systematic uncertain- ties, the optimal selection cuts are stronger (as can be seen in Table 2) and lead to higher values of S/B.
The mono-W differs from the other two channels. The SM background is dominated by vector boson fusion (VBF) processes, that lead to forward leptons (lost along the beam pipe) and W bosons. The signal is instead made of events where the W is radiated from the initial or final states, leading to a more central distribution. The cut on p T,W can efficiently suppress the VBF background, with a lesser impact on the signal compared to the mono-γ or mono-Z cases. As a consequence, we find that the mono-W search has the best sensitivity among the various mono-X channels. The 95% C.L. exclusion reach on M χ for a Majorana 3-plet and 5-plet is shown in Fig. 4 as a function of collider center-of-mass energy √ s and luminosity L. We also show the expected values of S/B for the excluded signal in absence of systematic errors, which are rather low also for the mono-W search.
Due to the presence of initial-state radiation, the W boson of the signal has a preference for being emitted in the forward (backward) direction, measured with respect to the flight direction of the − beam, if its charge is negative (positive). Since the charge of the W boson is potentially observable for leptonic decays, we can envisage a strategy to isolate the signal from the background using the full distribution in η W (instead of its absolute value). We thus also perform an analysis of leptonic mono-W events, where we impose the additional cut η W ± ≶ 0. We find the reach of this search to be weaker than the one of the inclusive mono-W because of the small leptonic branching ratio. However, the leptonic mono-W search possesses signal-free regions of the η W distribution which would allow for an in situ calibration of the background from the data itself, leading to possible reduction of the systematic uncertainties.
Di-V. We now consider scattering processes with multiple emission of vector bosons. While generally being suppressed by higher powers of the gauge coupling constant, these processes can be enhanced for large center-of-mass energies, and for multiplets with large weak charge. They can therefore provide very useful handles to probe WIMPs in the regimes where the mono-V searches have very low signal-to-noise ratios. Of course, a too large rate for multiple boson radiation would indicate the breakdown of the perturbative expansion, requiring the resummation of large logarithms. We have checked that for the EW 3-plet and 5-plet, and for the energies under consideration here, the fixed-order computations are still accurate.
First, we consider the di-photon process We apply the same acceptance cuts of the mono-γ analysis, and in addition we require a separation R γ γ > 0.4 between the two photons. We employ the same event selection strategy of the mono-γ case, using as variables η X , p T,X , where X is the compound γ γ system. Moreover, we require each photon to be as central as the γ γ system itself. For the 5-plet, we find that the di-γ search can be stronger than the mono-γ in presence of large systematic uncertainties, where suppressing the SM background is more important. For the 3-plet, which has a smaller EW charge, the signal yield is too much affected by the requirement of a second emission to be competitive with the mono-V. In both cases, the values of S/B for the excluded di-γ signal are much larger than for the monoγ signal, and systematic errors thus have a smaller impact. Details of the results are reported in Table 2 in the Appendix. Second, we consider the double W emission which holds a potentially very clean signature due to the two same-sign W bosons. We focus on leptonically decaying W bosons to ensure that their charge can be accurately tracked.
A potential SM background consists in events with two lost charged particles, with the leading contribution being where two W bosons of same sign are lost. This background is however negligible, as pairs of W bosons with opposite charge tend to be radiated from the same external leg and to be collinear: requiring only one of two collinear W bosons to be within detector acceptance reduces the rate to negligible levels. The other possible background is given by events with a misidentified charge, where in the second case the charged final-state leptons are lost along the beam line. Requiring p T,WW √ s/10 makes the process in Eq. (37b) subdominant with respect to the νν background Eq. (37a). On top of this p T cut, we do not apply further selection cuts, and simply require the two W bosons to be within the geometrical acceptance of the detector, |η W | < 2.5. As an estimate for the charge misidentification probability we take misid = 10 −3 .
Due to the negligible background contamination, the same-sign di-W signal has a much higher signal-to-noise ratio than the mono-V channels and even than the di-photon signal, reaching up to S/B ∼ O(1). This makes this channel very robust against systematic uncertainties, and particularly effective for large n-plets n ≥ 5 at higher energies due to their large EW charge. This signature may be one of the most robust and convincing signal of n = 5 multiplets at colliders. Further sources of background and a proper characterization of the missing (transverse) momentum in this reaction depend on detector performances, as well as on the knowledge of the initial state of the collision to be used in the computation of kinematic variables. We leave a careful evaluation of these aspects to future work.
We summarize the results of all the mono-V and di-V signatures discussed above in Fig. 5, where we show the 95% C.L. exclusion on M χ for real fermion 3-plets and 5-plets, together with the 5σ discovery potential, at two benchmark muon colliders. We also show the combined reach from all these missing mass channels. The bands with different shadings correspond to different systematic uncertainties. One can see that the inclusive mono-W yields the strongest exclusion for both the 3-plet and the 5-plet. The main effect of di-V searches is to reduce the impact of systematic uncertainties. A 14 TeV muon collider with the benchmark luminosity of Eq. (24) would be able to probe a thermally-produced Majorana 3-plet WIMP, while a center-of-mass energy of slightly above 30 TeV is needed to probe the thermal freeze-out mass with missing energy searches in the case of the 5-plet. Scalar WIMPs have lower production cross-sections. Missing mass searches do not allow to put stringent constraints on their mass, nor to probe the masses required for thermal freeze-out. We provide more details on the collider signatures, and results for real scalars in Appendix B.1.

Disappearing tracks
A second handle to tag the production of EW WIMPs at colliders is the detection of tracks from the charged states in the n-plet. As discussed in Sect. 2, the decay of χ ± → χ 0 π ± has a lifetime of roughly cτ χ + 48 cm/(n 2 − 1), which is sufficiently long-lived to give rise to reconstructed tracks of length O(cm) for n = 3, 5 that can be observable at colliders. The resulting tracks from these processes are somewhat too short for regular track reconstruction to work efficiently and they will show up as disappearing tracks (DTs), with missing hits in the outermost layers of the tracker and with little or no activity in the calorimeter and the muon chamber. States with higher electric charge in larger multiplets decay promptly to χ ± , and eventually contribute to the number of disappearing tracks.
A full-detector level study has shown that a high energy lepton collider like CLIC at √ s = 3 TeV can reconstruct them sufficiently well to separate them from other sources of look-alike short tracks [61,62]. A recent study [63] has attempted a first evaluation of the performance of this type of search at a multi-TeV muon collider. A main source of worry and a main difference with respect to e + e − machines is the abundant number of tracker hits from underlying event activity due to the muon beam decay and to the resulting secondary particles from the interactions with the machine and detector materials. These hits can accidentally become a potentially severe source of background for searches aimed at highlighting the presence of short tracks of BSM origin. We do not enter in the details of these issues here, and simply follow the analysis of [63], which is based on a simulation of beam-induced background at 1.5 TeV, and recast their results for the EW 3-plet and the 5-plet. We remind that the background from decaying muons is expected to decrease We consider mono-photon events with disappearing tracks, and search for events compatible with a WIMP signal. Following [63], we distinguish two event-selection strategies to hunt for disappearing tracks: (i) events with at least a disappearing track with p T > 300 GeV and a hard photon with E γ > 25 GeV; (ii) events with a hard photon, and two disappearing tracks originating from the same point along the beam axis. To estimate the reach we work in the cut-andcount scheme as in Eq. (32), and ignore systematic uncertainties. Further details are summarized in Appendix B.3 for completeness.
The result of our recast is shown in the last two columns of Fig. 5 for Majorana 3-plets and 5-plets at two benchmark colliders, and in Fig. 6 as a function of collider energy and luminosity. One can see that DTs are especially powerful in the case of the 3-plet, where the reach goes almost up to the kinematical threshold. In particular, an EW 3-plet WIMP of mass as predicted by thermal freeze-out can be discovered already at a 6 TeV muon collider as suggested in [24,63]. For higher n-plets DT substantially loose exclusion power because the lifetimes of the χ ± → χ 0 π ± decay become shorter. For the 5-plet the DT reach is comparable to the combined reach of the MIM searches.
As discussed in more detail in the Appendix, DT searches are particularly important to probe scalar WIMPs, since the lower production cross-sections have no significant impact on these almost background-free searches. Disappearing tracks might be the only direct signature of scalar WIMPs at collider experiments.

WIMP direct and indirect detection
In this Section we briefly summarize the opportunities of the future experimental program in direct and indirect detection in light of the mass predictions derived in Table 1.

Indirect detection
The current and upcoming ground-based Cherenkov telescopes are in a very good position to probe heavy WIMP n-plets with n > 5, which would be inaccessible otherwise. Indeed, these telescopes are designed to detect very high energy gamma-rays (i.e. E γ 100 GeV) coming from different astrophysical objects and they are therefore sensitive to the gamma-ray signal from the annihilations of EW n-plets. The typical spectrum is characterized at very high energy by gamma-ray lines, peaking at the DM mass E γ M χ , from the loop-induced annihilations into γ γ and γ Z . The cross-section in this channel is largely boosted by the SE (see e.g. [5,64,65]) and can raise above the gamma-ray continuum from the showering, hadronization and decays of the electroweak gauge bosons [66].
From the astrophysical point of view, the reach of high energy gamma lines searches depends very much on which portion of the sky the telescopes will be pointed at. In finding the optimal choice, a balance has to be found between the maximization of photon flux at Earth and the control over the systematical uncertainties. Two very well studied astrophysical targets are the Galactic Center (GC) [20,67] and the Milky Way's dwarf Spheroidal galaxies (dSphs) [20]. In the GC, the uncertainties are dominated by the importance of the baryonic physics in the inner most region of the Milky Way which comes together with the poor knowledge of the DM distribution at the center of the Milky Way [68][69][70][71]. On the contrary, dSphs stands out as very clean environments to search for high energy γ -lines only residually affected by systematics related to the determination of their astrophysical parameters in the presence of limited stellar tracers [72,73].
Motivated by the above considerations, we show a very preliminary analysis of ID signals coming from annihilations of the WIMP 7-plet. We focus on the CTA prospects by considering 50 h of observations time towards two dSph targets in the northern hemisphere: the classic dSph Draco and the ultra-faint one Triangulum II. Notice that the DM properties of Draco come from hundreds of stellar tracers, while those from Triangulum II are based on just 13 tracers, making the latter more speculative and subject to large systematics in the determination of the geometrical J -factor [74]. Hence, the reach of Draco should be taken as the baseline reach for CTA.
Our analysis is simplified because the signal shape we consider is essentially a single line at E γ M χ . Consistently we take the CTA prospects derived in Ref. [20] for a pure line. We ignore the contributions of the continuum spectrum, the extra features of the spectral shape induced by the resummation of EW radiation and the contribution of the BSF to the photon flux. While neglecting BSF is justified if we focus on very high energy photons, a careful computation of the γ + X cross-section, where X is any other final state would be needed to precisely assess the experimental sensitivity [75]. In the last decade, many different groups have investigated the impact of large Sudakov logarithms and large collinear logarithms on the indirect detection reach, focusing mainly on the case of the fermionic 3-plet [76][77][78][79][80]80]. The inclusion of these effects has been shown to increase the reach of ∼ 20 ÷ 30% for the 3-plet [20,67,81] and it is expected to be even more important for higher DM masses.
In Fig. 7 we overlay the SE annihilation cross-section for the 7-plets at v = 10 km/s against the CTA experimental reaches. In order to compute the SE in this velocity regime, we took advantage of the parametrization introduced in [65] and used the full expressions for the SE at leading order, including EW breaking effects. The SE saturate already at v 10 −3 ÷ 10 −2 far away from the resonances. As we can see, both a 50 hour observation of Triangulum II and of Draco have good chances to detect the high energy γ line in the 7-plet annihilation spectrum.
As we see from Fig. 7, given the strong mass-dependence of the features of the SE cross-section, a major source of theoretical uncertainty on the reach of indirect detection is still the determination of the 7-plet thermal mass. Therefore, a full computation of the thermal relic mass including NLO effects is required together with a careful computation of the γ + X cross-section along the lines of Ref.s [76][77][78][79][80]80] to careful assess the indirect detection reach for the 7-plet.
Independently on our current inability of making a conclusive statement because of the large theory uncertainties, it is clear that large n-plets are a perfect target for future Cherenkov telescopes which deserves further theoretical study. A complementary open phenomenological question is if the low energies gamma lines at E γ E B associated to BSF can be actually disentangled from the continuum (see [12,82] for preliminary work in this direction). An analogous question can be asked for monocromatic neutrinos from BS annihilations. Fig. 7 Expected CTA sensitivities (dashed black lines) with 68% and 95% CL intervals derived as in Ref. [20] assuming 50 h observation time towards Draco (green) and Triangulum II (magenta). We show the SE annihilation cross-section into the channels that contribute to the monocromatic gamma line signal (i.e. γ γ an γ Z ) for a scalar 7plet (blue) and a fermionic 7-plet (red). The vertical bands show the predicted thermal masses for the scalar 7-plet (blue) and the fermionic 7-plet (red), where the theory uncertainty is dominated by the neglected NLO contributions (see Table 1)

Direct detection
For Y = 0 the elastic scattering of DM with the nuclei is induced by EW loop diagrams first computed in [87,88]. After EW gauge bosons are integrated out, the structure of the UV effective Lagrangian describing the DM interactions reads where we focus on the DM spin independent (SI) interactions with quarks and gluons [89]. The quark twist-2 operator is defined as O q μν ≡ i 2q D μ γ ν + D ν γ μ − g μν / D/2 q. The Wilson coefficients of the operators for general EW nplets with Y = 0 have been computed in Ref. [90] and at the leading order in M χ /m W,h 1 read where m h = 125 GeV is the SM Higgs mass, q ∈ (c, b, t) and κ c = 1.32, κ b = 1.19, κ t = 1. Following Ref. [89], starting from the UV DM interactions we derive the IR interaction of DM with the nucleons. All in all, the SI elastic cross-section per nucleon in the limit M χ m N reads where m N is the nucleon mass and k EW N is defined as with the dimensionless nucleon form factors defined as (2)), where q(2) andq(2) are the second moments of the parton distribution functions for a quark or antiquark in the nucleon taken from [90]. Notice that we choose a different set of values for the nucleon form factors with respect to previous studies [91] which explain the difference in our results. In particular, we take the FLAG average of the lattice computations in the case of N f = 2 + 1 + 1 dynamical quarks [92][93][94].
By propagating LQCD uncertainties on the elastic crosssection Eq. (41), we obtain the vertical uncertainties on the SI cross-section predictions in Fig. 8. We find the partial accidental cancellation between the one loop and the two loop contribution to reduce the elastic cross-section up to 30%. The horizontal bars represent the uncertainties coming from the computation of the thermal masses through the relic abundance. As shown in the plot, while all the WIMP crosssections lie above the Xenon neutrino floor as computed in [86] but only a very large exposure experiment like DARWIN [19] would be able to probe the heavy thermal WIMPs.
Spin dependent (SD) interactions of DM with the nuclei are also induced by EW loops where the Wilson coefficient was computed in Ref. [90] and we expanded it at zeroth order in M χ /m h 1. The corresponding SD cross-section is too small to be probed even at a very large exposure experiment like DARWIN.
Finally, we comment on the new opportunities for direct detection that arise for scalar DM. Here, a non-zero Higgs portal quartic in Eq. (2) leads to a new contribution to the SI DM scattering cross-section with the nuclei, which again in the M χ m N limit reads where Fig. 8 In dark green we show the present constraints from XENON-1T [83] and PandaX-4T [84], the green dashed line shows the reach of LZ [85] and the brown green dot-dashed line the ultimate reach of DARWIN [19]. The light gray region show the neutrino floor for 200 ton/year exposure derived in Ref. [86]. Left: expected spin independent (SI) direct detection cross-section for Majorana n-plets (red) and for real scalar n-plets (blue) (assuming the Higgs portal coupling λ H = 0).
The vertical error bands correspond to LQCD uncertainties on the elastic cross-section in Eq. (41) while the horizontal error band comes from the theory determination of the WIMP freeze out mass. Right: current and future reach on the Higgs portal quartic λ H defined in Eq. (1) for scalar DM. In the shaded dark red region the quartic modifies the freeze-out cross-section by O(1) or more. The dashed red contours indicate smaller ratios of the Higgs-portal and the EW annihilation cross-sections with f N 0.31 obtained from lattice QCD results (see [95] for a more detailed discussion on the scalar triplet). In the right panel of Fig. 8 we show the regions of parameterspace where the Higgs-portal interaction can be tested in direct detection. The requirement of not significantly affecting the freeze-out dynamics bounds the annihilation crosssection induced by the Higgs portal to be smaller than the EW cross-section, σ H ann /σ EW ann 1, which results in an upper bound on the quartic coupling λ H shown by the red shading in Fig. 8. An estimate for this bound can be obtained by comparing the hard annihilation cross-sections, and reads λ 2 H (n 2 − 3)(n 2 − 1)g 4 2 /8. Interestingly, XENON1T and PANDAX-4T already exclude a large part of the region where the Higgs portal induces O(1) modifications of the freeze-out predictions, while LZ will completely exclude this possibility.

Conclusions
After many years of hard experimental and theoretical work, the possibility that Dark Matter is part of an EW multiplet is still open and deserves theoretical attention in view of the future plans for experimental searches. In this paper we made a first step in sharpening the theoretical predictions computing all the calculable thermal WIMP masses for real EW representations with vanishing hypercharge. We included both Sommerfeld enhancement and bound-state-formation effects at LO in gauge boson exchange and emission. Our results are summarized in Table 1.
We find that the largest calculable SU(2) n-plet at LO is the 13-plet, which is as heavy as 350 TeV. Stronger requirements about the perturbativity of the EW sector up at high scales can further lower the number of viable candidates. We consistently assign a theory error to our predictions by estimating the NLO corrections to the SE. The latter dominate the theory uncertainty for n ≥ 7, while for n = 5 the error is dominated by the approximate treatment of EW symmetry-breaking effects in the computation of the BSF cross-sections.
Given the updated mass predictions from thermal freezeout, we re-examined various phenomenological probes of WIMP DM.
High energy lepton colliders in the 10-30 TeV range, such as a future muon collider, can directly produce EW multiplets with n ≤ 5. In order to probe a Majorana fermion with n = 3 (n = 5) with missing-mass searches, a collider with at least √ s ∼ 12 TeV ( √ s ∼ 35 TeV) and the baseline integrated luminosity of Eq. (24) would be required. The highest mass reach is obtained by means of an inclusive mono-W search.
Interestingly, disappearing tracks originating from the decay of the singly-charged state into the neutral one are robust predictions of real EW multiplets with Y = 0, and ameliorate the sensitivity for the 3-plet compared to missingmass searches. For the 5-plet we find the expected sensitivity of disappearing tracks to be very similar to the one of missing-mass searches due to the shorter average lifetime of the tracks. Scalar WIMPs can not be probed through missing-mass searches, due to their smaller production cross-section. However, disappearing tracks searches are very powerful tests even for scalar multiplets, thanks to their very low background contamination. This signature is therefore a crucial ingredient to fully explore the parameter space of thermally produced WIMP Dark Matter at future colliders.
Heavy EW WIMPs with n > 5 are too heavy to be produced at colliders. However, they are perfect targets for indirect detection at upcoming ground-based Cherenkov telescopes like CTA. More theoretical work is necessary to make a robust forecast both on the determination of the photon spectrum for large n-plets and on improved precision predictions for the freeze-out masses.
Finally, large-exposure liquid Xenon experiments like DARWIN can in principle probe all the relevant EW WIMPs through their weak interaction with nuclei. Scalar WIMPs can further be tested through their Higgs-portal quartic interaction. Interestingly, O(1) modification of the thermal freeze-out masses due to the Higgs portal are already partially excluded by the XENON1T and PANDAX-4T results, and will be completely excluded by LZ.
A natural continuation of the work done here would be to consider complex EW multiplets. For vanishing hypercharge both the cosmology and the phenomenology will be very similar to the ones discussed here. The suppression of the annihilation cross-section, resulting in lower thermal masses, together with the enhancement of the production cross-section at colliders will favour the direct exploration of complex multiplets at colliders. More interestingly, EW multiplets with nonzero hypercharge, such as the Higgsino in supersymmetric models, are also phenomenologically viable if the DM elastic scattering with nucleons is suppressed or kinematically forbidden. Classifying the predictions of this class of models would give a complete picture on EW DM multiplets. We hope to come back to this open issues in the near future.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study and no experimental data has been listed.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendic A: Bound states dynamics
In this Appendix we discuss in detail the dynamics of Bound State Formation (BSF). First we discuss in Sect. 1 the general features of BSF at leading order (LO) in gauge boson emission and at next-to-leading order (NLO) in gauge boson emission. Then we detail in Sect. 1 the 7-plet BS dynamics, focusing on the differences with the 5-plet case.
A.1 Bound state formation at LO and at NLO At leading order, bound states form through the emission of a single vector boson V a : χ i + χ j → BS i j + V a . The non-relativistic limit of the amplitude can be recast in the form of an effective interaction Hamiltonian, such that the full amplitude can be obtained from its matrix element with the wave function of the initial and final two-particle states (reconstructed from the resummation of the ladder diagrams). The leading order contribution to this effective hamiltonian comes in the form of electric dipole interaction terms [12,43]: where the first to terms are a simple generalization of the standard QED dipole interaction while the last one is a purely non-abelian term which arises from vector boson emission from a vector line. The computation of the transition amplitudes from Eq. (A.1) simplifies if we assume the SU(2) L -invariant limit. This approximation applies when the DM (BS) de Broglie wavelength is much smaller than the range of the Yukawa interaction 1/m W and therefore for z ≤ (M χ /m W ) 2 . In this regime the Yukawa potential is well approximated by the Coulomb one which turns out to be a good approximation to describe WIMP freeze-out. The BS dynamics can then be understood by using isospin selection rules while the main consequence of having finite vector masses is to provide an energy threshold to the emission of a single massive boson in the formation or the decay of a BS.
Since α eff ∼ n 2 , increasing the dimensionality of the DM multiplet enhances next to leading order (NLO) processes in gauge boson emission such as χ i + χ j → BS i j + V a + V b . These could be in principle relevant for both the computation of the thermal mass and the saturation of the perturbative unitarity bound.
The main NLO contributions to BSF come from diagrams like the ones in Fig. 9 and are essentially of two types: (i) the first diagram is essentially the second order Born approximation of the LO Hamiltonian, with the intermediate state being a free or a BS; (ii) the second diagram, where the two Fig. 9 Examples of diagrams controlling the BS effective Hamiltonian at next-to-leading order in gauge boson emission The first diagram corresponds to the second order Born approximation for the dipole operators in Eq. (A.1), where the resummation of the vector boson insertions between the two emission reconstructs the wave function of an intermediate BS or scattering state. The second diagram, instead, is obtained from the O( A 2 ) terms in the interaction Hamiltonian, at leading order in the Born approximation emitted vectors come from the same vertex, is generated by the effective Hamiltonian at order O(A 2 ). The latter contains terms of the form where we focus here on the abelian part of the hamiltonian, postponing a full study for a future work. Given the above Hamiltonian and the LO one in Eq. (A.1) we can estimate the corresponding contribution to the double emission BSF cross-section as: where g χ = 1 for Majorana fermions (g χ = 2) for real scalars. In the LO estimate, a factor 2 α eff M χ 2πα eff v rel comes from the overlap integral while a factor E 8π from the two-body phase space. Similarly, in the NLO estimate a factor 1 2 E 3 256π 3 comes from the 3-body phase space, taking into account the two identical final vectors, and 2 α eff M χ 3 2πα eff v rel from the overlap integrals between the wave functions. From the above formula we derive the scaling of the NLO corrections in Eq. (22).
We now discuss the contributions from second order Born expansion whose general expression is given by where we defined with where E B is a typical binding energy. This quick estimate, supported by the full numerical computation, suggests that C BS contribution is fully captured in the Narrow Width Approximation (NWA) for the intermediate BS. Therefore, neglecting the interference terms, one gets which is exactly the single emission result. To estimate the contribution from C free we need to estimate I q k which encodes the contribution from intermediate continuum states. For simplicity, we stick to the abelian contribution which reads The integral above can be split into small and large r regions, roughly separated by the Bohr radius a 0 = 1 which plugged into Eq. (A.6) gives an estimate to C free . All in all, plugging these estimates in Eq. (A.5) and replacing q = M χ v rel we get that the contribution from NLO exchange In conclusion, NLO corrections to BSF are suppressed by ∼ α 3 eff /64π with respect to the LO ones. As a consequence, the leading NLO contributions to the total annihilation crosssection are the ones correcting the LO SE. The latter are logenhanced as detailed in Eq. (23) and first computed in [47] for the fermionic 3-plet. It would be interesting to extend these computation to higher EW n-plets.

A.2 The 7-plet bound states in detail
The 7-plet has a richer bound states dynamics with respect to the 5-plet, essentially because of the additional layers of isospin and energy levels. As we will discuss here, keeping track of this dynamics is crucial to compute correctly the relic abundance.
The left panel of Fig. 10 shows the relative importance of the different BS and the SE to the effective cross-section at fixed T = 10 −3 M χ . As we can see, BSF accounts for most of the total cross-section. Compared to the 5-plet case, the new attractive isospin channels with I = 7 give a sizeable contribution to the 7-plet cross-section as well as the 2 p states which were instead irrelevant for the 5-plet. In the right panel of Fig. 10 we show how the details of the bound state dynamics are especially important at temperatures around the freeze-out (i.e. z = 10 2 ) where the effects of BS breaking due to interactions with the plasma are non negligible. This can be seen by comparing the behavior of the full computation of the effective cross-section (solid lines) against the BSF cross-section with zero ionization rate (i.e. R BS = 1 in the notation of Eq. (18)). In particular taking R BS = 1 yields an overestimate of the final thermal mass of about 6 TeV. Interestingly, we see that for z = 10 3 all the BSF rates approach the R BS = 1 limit, signalling that the ionization rate is already heavily Boltzmann suppressed.
We now illustrate the details of the BS dynamics for the 7-plet. The general computation outlined around Eq. (12) is in general cumbersome, but it simplifies singling out the specific features of each BS. These are summarized in Fig. 11. We now discuss them in turn, going from the largest to the smallest binding energy. M χ . Since their decay rate can be neglected, the effective cross-section can easily be obtained from Eq. (19). (ii) The 1s 7 BS cannot decay directly into SM pairs because of its large isospin so that its annihilation rate arises at NLO in gauge boson emission. Similarly, the decay to lower 1s states can only go through NLO processes or velocity-suppressed magnetic transitions. As a consequence, this BS can only be excited to 2 p 5 at LO, and its effective cross-section can be written in terms of the one of the 2 p 5 : where the excitation rate can be written in terms of the decay rate I →J g I /g J J →I e − E T times the prob- (iv) The annihilation rates into SM state of the 2 p I BS are suppressed by α 2 eff compared to the ones of the 2s I BS. Their dynamics is then dominated by the decay (excitation) rates into lower (higher) s−orbital BS which scale as dec ∼ α 5 eff M χ . A simple example of this dynamics is provided by the two-state system 2 p 1 − 1s 3 where 2 p 1 dominantly decays to 1s 3 , which promptly annihilates to SM. The effective cross-section of 2 p 1 reads as we would intuitively expected. The other 2 p states have more intricated chains, which involve also excitations 3s states.
(v) We also include 3s I BS which annihilate directly to SM for I ≤ 5 and decay into p−orbitals states.
Finally, we checked that p states with n > 2, s states with n > 3, and BS with I = 9 have a negligible impact on the cosmological evolution.

B.1 The scalar WIMPs
Probing scalar WIMPs with typical missing mass searches is quite hard. This is due to multiple reasons: (i) the scalar production cross-sections are roughly one order of magnitude smaller than for fermions with same n, as shown on the left of Fig. 13. A factor of 4 suppression comes from the lower number of degrees of freedom for scalar final states, while the remaining suppression comes from a velocity suppressed production cross-section compared to the fermionic case. Since the reach is a very slow function of the mass of the WIMP M χ , as shown in the right panel of Fig. 13, a reduction of the signal cross-section implies a drastic change in the reach. (ii) The scalar WIMPs have typically larger freezeout masses compared to fermionic WIMPs with same EW charge n. All in all, scalar WIMPs give dimmer signals at colliders and are generically heavier than fermionic WIMP. It is thus not surprising that the results expected from collider searches of scalar WIMPs, shown in Fig. 12, are far less exciting than those for fermions in Fig. 5. The overall picture in the landscape of possible beam energy and luminosity options for a future very high energy lepton collider is displayed in Fig. 14. At variance with the fermionic case presented in Fig. 4, the potential to probe scalar WIMPs with mono-X signals is very limited. Details on the optimized analyses we carried out are given in Table 3.
We stress that our results are based purely on Drell-Yan production of χ , which accounts perfectly for the total production rate of WIMPs of mass comparable with √ s. For significantly lighter WIMPs it is possible to add further production modes and discovery channels, such as production by vector boson fusion and mono-muon channels studied for lighter fermionic WIMPs [24], which may result in a bound for light enough scalar WIMPs. In Fig. 15 we plotted the cross-sections for scalar χχ production in W -fusion (as a representative for VBF modes) and Drell-Yan as a function of M χ . It can be seen that the VBF cross-section decreases quickly, while DY remains almost constant except near the kinematic threshold. In particular, for the real scalar 5-plet at √ s = 30 TeV our DY 2σ reaches can be trusted, as the VBF contribution is smaller than 10% of the DY one. For the scalar triplet at √ s = 14 TeV, the inclusion of VBF modes is not expected to improve the reach for masses 1 TeV.
It is remarkable that for real scalars the mass splitting between charged and neutral states in the n-plet is dominated by EW interactions. Indeed, no splitting term with the Higgs can be written at the quartic level, due to the antisymmetry of the SU(2) contraction. By hypercharge conservation, and assuming the scalar does not get any extra VEV, the leading terms contributing to the mass splitting are dimension 6 in the SM. Therefore the stub-track prediction is robust and does not depend on peculiar UV completions of the model. Results for searches of scalar WIMPs from stub-track analyses are reported in Fig. 16. Details on the recasting of results contained in Ref. [63] to obtain our results are given in the following Appendix B.3.

B.2 Details of the missing momentum analyses
In Tables 2 and 3 we provide the results of the optimized cuts for all the considered mono-V and di-V channels, for the case of a Majorana n-plet, or a real scalar, respectively. The optimization was carried out in an equally spaced 25×12 Table 2 95% C.L. reach on the mass of a Majorana 3-plet and 5-plet from the various mono-X channels. The excluded number of signal events S 95% and the relative precision S 95% /B are also given, together with the values of the optimal event selection cuts on η X and p T,X , where X is either the single vector boson or the compound diboson system for Di-W and Di-γ .   We also report the expected number of signal events, the signal-to-noise ratio, and the value of the mass that can be excluded at 95% C.L. We provide results for muon colliders with √ s = 3, 14, 30 TeV with integrated luminosity as in Eq. (24), and for systematic uncertainties sys = 0, 1 0 / 00 , 1%.
Among all the channels considered, the only background that needs some careful treatment is the mono-W one. We split this background in two contributions. For pseudorapidities of the final state lost muon η μ > η match (computed with respect to the direction of the initial state muon with the same charge), we compute the cross-section of the process γ μ ∓ → W ∓ ν, using the improved Weizsäcker-Williams approximation [96]. For 2.5 < η μ < η match , we compute the full hard process μ − μ + → W ∓ ν ± . The values used for η match are 5.4, 7.0, 7.5 for √ s = 3, 14, 30 TeV, respectively. These values are such that the two back-ground contributions are the same in the pseudorapidity region (η match , η match + 0.2) for the lost muon.

B.2 Recasting the disappearing tracks
We recast the two search strategies discussed in Ref. [63] that exploit the presence of a single short reconstructed disappearing track or a two-track analysis that require at least one of them to be a short disappearing track, in addition to a trigger photon. The requirements are summarized in Table 4 from Ref. [63].
Single-track search. For the single-track analysis we take the background cross-section quoted in [63]. This rate is mainly determined by the combinatorial of track reconstruction induced by beam-induced backgrounds. 6 To determine the rate of the single-track events, we compute the monophoton cross-section doubly differential in the polar angles of the charged particles χ 1 , χ 2 . This dσ/dθ 1 dθ 2 is obtained at LO in perturbation theory with MadGraph5_aMC@NLO and is further reweighted to take into account angular and distance sensitivity to stub-tracks reported in Ref. [63]. Let P(θ 1 ) be the probability that the particle χ 1 is reconstructed as a track: P(θ, r min , r max ) = r max r min dr rec (r, θ) cτβγ sin θ e −r/(cτβγ sin θ) , (A.14) where r is the transverse radius and rec (r, θ) is the probability to reconstruct as a track a particle travelling at an angle θ that decayed at a transverse radius r given in Fig. 11 of Ref. [63]. For single tracks rec (r, θ) is 0 outside the interval r ∈ [50 mm, 127 mm], and outside π/6 < θ < 5π/6. The radial condition reflects the fact that tracks can only be reconstructed if the particles make at least 4 hits in the vertex detector, which for the considered geometry means that the particle must travel at least a minimum distance of 50 mm in the detector, while the upper limit stems from the disappearing condition of the track. The latter condition will be relaxed in the 2-tracks search. With the knowledge of rec the integral in Eq. (A.14) can be performed numerically. As per Table 4, the hard cross-section σ S,γ is subject to trigger requirements: the leading observed track is required to have p T > 300 GeV (B.1) to help discriminate it against fake tracks, and it must lie within the cone 2π 9 < θ < 7π 9 .
(B.2) 6 As acknowledged in [63], this estimate of the background is quite conservative because it is based on detailed beam dynamics simulation for √ s = 1.5 TeV. Due to the relativistic dilution of muon decays, we expect smaller background cross-section at higher √ s.
In our recast, due to lack of a detailed tracking and detector simulation, these cuts are implemented at parton level on the DM particles momenta, which leads us to overestimates the number of events that pass the selection. To account for this effect we assume that only a fraction tran of the events with parton p T > 300 GeV gives a track whose p T fulfils the same conditions. The transfer factor tran ≈ 0.5 is estimated from the p T distribution of χ obtained at generator level, and track p T distribution given in Ref. [63]. We assume that tracks with p T > 300 GeV can only come from χ with p T > 300 GeV. To properly avoid over-counting events with two reconstructed tracks, we divide the final state phase space into two non-overlapping regions that require different reconstruction constraints: (i) Both χ fulfil the conditions to be considered as leading track (Eqs. (B.1) and (B.2)). In this case both tracks are subject to the detection and reconstruction efficiencies tran and rec (θ, r ). These events may give rise to zero, one, or two reconstructed stub-tracks. We count events with at least one stub-track.
(ii) Exactly one χ fulfils the conditions to be considered as leading track. Only events in which this track is reconstructed according to detection and reconstruction efficiencies tran and rec (θ, r ) are counted. The fate of the sub-leading χ (if any) is irrelevant.
The largest contribution to the single-track cross-section comes from events in region (i), where both DM particles satisfy the p T and θ requirements to be considered as a leading track. The preference for this configuration reflects the approximate 2-body kinematics of the mono-γ events with small p T . In order to understand the nature of signal we can split it into two further sub-categories with: (a) exactly one reconstructed track which fulfils the conditions Eq. (B.1) and Eq. (B.2); (b) exactly 2 reconstructed stub-tracks, of which at least one fulfils the same conditions. The respective rates are given by: S,γ d cos θ 1 d cos θ 2 · tran 2P(θ 1 )(1 − P(θ 2 )) 1 track, 1 − (1 − tran ) 2 P(θ 1 )P(θ 2 ) 2 tracks, where the hard cross-section σ 1T S,γ is restricted to the phasespace region where both χ particles fulfil the requirements of Eqs. (B.1) and (B.2). The boost factor βγ and the angular distribution are both taken from a MC sample with cuts only on the photon at generator level.
The resulting number of events is used to compute the reach on the DM mass reported in Fig. 5, according to Eq. (32) with sys = 0.
Interestingly, the results obtained from the MC sample can also be understood semi-analytically thanks to the simple kinematics of the mono-photon process. Given that the photon tends to be soft, the kinematics of the three body process is not too different from direct production of a pair of oppositely charged DM particles without the photon. Therefore a very good analytic approximation of the above results can be obtained, with the χ boost factor and flight directions approximated by the ones for pair-produced DM particles with energy √ s/2, βγ ≈ s 4M 2 χ − 1, θ 1 = π + θ 2 . (B. 3) The angular distribution can also be computed analytically in the 2-body limit,  We additionally require the two tracks to originate from points that are close to each other along the direction of the beam axis, z < 0.1 mm (see Table 4). This effectively reduces the background to negligible levels. In this limit, we use 4 signal events as a conservative estimate of the 95% C.L. exclusion for a Poissonian counting. The angular cuts on the tracks are the same as in the single track case, while the p T cuts are much milder: p T > 10, 20 GeV for the sub-leading and leading tracks, respectively. In this case the mismatch between the p T of the reconstructed track and the p T of the charged χ obtained at generator level is negligible. The additional cuts do not affect significantly the signal events. Note that, following Ref. [63], the disappearing condition is required on at least one track, i.e. this analysis includes in the signal all events in which the second track extends up to a transverse radius of r = 1153 mm. Following Ref. [63], we assumed for such long tracks a reconstruction efficiency equal to the tracks decaying between 101 mm < r < 127 mm. Also for double tracks, the result obtained using the MC sample βγ and θ distributions are in agreement with the ones computed analytically in the 2-body limit.
We remark that for SU(2) triplets the double track analysis has a higher exclusion power than the single track analysis, whereas for n ≥ 5 it has a lower reach. This is due to the shorter life-time τ χ ∝ 1/n 2 of larger multiplets, that suppresses the exponential decay factor of Eq. (A.14) twice in the double-track rate.