The last Complex WIMPs standing

We continue the study of weakly interacting massive particles (WIMP) started in [arXiv:2107.09688], focusing on a single complex electroweak $n$-plet with non-zero hypercharge added to the Standard Model. The minimal splitting between the Dark Matter and its electroweak neutral partner required to circumvent direct detection constraints allows only multiplets with hypercharge smaller or equal to 1. We compute for the first time all the calculable WIMP masses up to the largest multiplet allowed by perturbative unitarity. For the minimal allowed splitting, most of these multiplets can be fully probed at future large-exposure direct detection experiments, with the notable exception of the doublet with hypercharge 1/2. We show how a future muon collider can fully explore the parameter space of the complex doublet combining missing mass, displaced track and long-lived track searches. In the same spirit, we study how a future muon collider can probe the parameter space of complex WIMPs in regions where the direct detection cross section drops below the neutrino floor. Finally, we comment on how precision observables can provide additional constraints on complex WIMPs.


I. INTRODUCTION
Despite the intense experimental activity in the past 50 years, the possibility that Dark Matter (DM) is a new weakly interacting massive particle (WIMP) remains open today. Within this framework, the most minimal and predictive possibility is that the DM is the lightest neutral component of a single electroweak (EW) multiplet. In this simplified setup, the experimentally viable DM candidates can be fully classified and their thermal masses predicted once the relic abundance from freeze-out is required to match the observed DM density today Ω DM h 2 = 0.11933 ± 0.00091 [2].
Experimentally, direct detection (DD) experiments robustly excluded the possibility that DM has treelevel Z-mediated interactions with the nuclei. This experimental milestone left only two possible scenarios for EW WIMPs wide open: i) WIMPs with zero hypercharge (Y = 0), interacting with nucleons only trough loop-induced interactions, ii) inelastic WIMPs with non-zero hypercharge (Y = 0), where the elastic scattering between DM and the nucleons is forbidden kinematically. These two possibilites have been historically incarnated by the Wino and the Higgsino DM scenarios in the Minimal Supersymmetric Standard Model, but a full classification of the viable EW WIMPs from an effective field theory standpoint is still lacking.
Attempting to fill up this gap in the literature, we studied EW WIMPs with Y = 0 which are naturally embedded in real EW representations in Ref. [1]. Here we focus instead on inelastic WIMPs which are embedded in complex EW representations with Y = 0. 1 1 Complex EW representations include also multiplets with Y = 0. Some of these were discussed in Ref. [3], their Similarly to our previous study, our goal here is to precisely determine and fully classify all the calculable freeze-out predictions for complex inelastic WIMPs, including Sommerfeld enhancement (SE) and Bound State Formation (BSF). The final purpose is to use these predictions to establish if and how the future experimental program will be able to test all the viable WIMP scenarios. In particular, we will show how future large exposure direct detection experiments together with a future high energy muon collider will probe complementary portions of the complex WIMP parameter space being potentially able to fully test the theoretically allowed WIMP scenarios.
Our results are summarized in Table I, whose logic can be explained as follows. Once the DM mass is fixed from the freeze-out predictions, the phenomenology of complex WIMPs depends essentially on two parameters: i) the "inelastic" splitting between thenext-to lightest neutral component and the DM; ii) the "charged" splitting between the charged components and the DM. In our setup these splittings are generated by its (non-renormalizable) interactions with the SM Higgs (generated by uspecified UV dynamics).
The inelastic splitting δm 0 is bounded from below by DD constraints [6,7] and BBN constraints on the decay of the next to lightest neutral component. Interestingly, this requrement alone selects a limited number of complex WIMPs: i) scalar and fermionic WIMPs with Y = 1/2 and even n up to the unitarity bound of the freeze-out annihilation cross section [4]; ii) scalar and fermionic WIMPs with Y = 1 and n = 3, 5. At the same time, the inelastic splitting is bounded from above by DD constraints on phenomenology is similar to the one discussed in Ref. [1]. For completeness we give the full predictions for the thermal masses up to the limit of perturbative unitarity in Appendix II. The upper bound on n for even multiplets comes from the perturbative unitarity bound, as can be seen from the (σv) J=0 tot /(σv) J=0 max , where (σv) J=0 max is the maximal allowed annihilation cross section [4]. The loss of perturbativity is also signaled by the Landau pole Λ Landau progressively approaching the DM mass. The upper bound on odd n with Y = 1 comes from the perturbativity of the higher dimensional operators generating δm0 . For multiplets with n > 5 the largest UV cutoff Λ max UV required to generate the minimal viable splitting is smaller than 10MDM. For each candidate we provide the allowed range for the mass splittings. The lower limit on δm0 comes from strongest bound between direct detection and BBN as shown in Fig. 1. The upper bound from the most stringent condition between DD constraint from PandaX-4T [5] and the perturbativity of the coupling of O0 in Eq. 1. Similarly, the lower limit on δmQ M comes from the BBN bound on the charged state decay rate, while the upper limit from the strongest limit between DD and the perturbativity of the coupling of O+ in Eq. 1.
Higgs-mediated nuclear recoils. This leaves a finite window for the inelastic splitting of every multiplet which we report in Table I. This window will be further probed by large exposure DD experiments such as LZ [8], Xenon-nT [9], and ultimately by DAR-WIN/G3 [10,11].
The natural value of the charged splitting δm Q is fixed by the radiative EW contributions [12][13][14] but (non-renormalizable) interactions with the SM Higgs can induce large deviation from this value. In particular, for all the n-plets with non-maximal hypercharge, these interactions are required to make the DM stable.
The allowed range of the two splittings above controls the hierarchy of the states within the EW multiplet. In this parameter space one can map out the expected signals in a future hypothetical muon collider [15]. Depending on the lifetime of the charged states we can have different signatures at colliders: i) long lived charged tracks; ii) disappearing tracks (DT); iii) missing energy accompanied by an EW bosons. While the first two searches rely on the macroscopic decay length of the charged states, the last is directly related to the DM pair production recoiling against one (or more) EW boson. In general, a future muon collider could complement the large exposure DD experiments in probing complex WIMPs in regions of the parameter space where the DD drops below the neutrino floor of xenon experiments [16].
The doublet with Y = 1/2 (2 1/2 ) and the triplet with Y = 1 (3 1 ) stand out as the only two complex candidates for which the lightest component of the multiplet is guaranteed to be electrically neutral, which is a necessary requirement for DM. It is well known that the 2 1/2 DD cross section lies well below the neutrino floor because of an accidental cancellation between the different contribution to the elastic cross section with the nuclei [17]. We find that a muon collider of √ s = 3 TeV would only be able to exclude the fermionic EW doublet with a large luminosity from around L = 5 ab −1 to L = 50 ab −1 , depending on the assumed systematics, by a combination of charged tracks, disappearing, tracks and missing mass searches. Somewhat less high luminosity may suffice at a √ s = 6 TeV collider, which can exclude the 2 1/2 with 4 ab −1 and make a discovery with a more demanding 20 ab −1 . 2 The DD cross section of the 3 1 can only be probed with large exposure Xenon experiment like DARWIN, while a muon collider of √ s = 10 TeV would be able to make a discovery at the thermal mass for the fermionic 3 1 using a luminosity around 10 ab −1 , that is considered as a benchmark for currently discussed project for this type of collider [20,21]. This is clearly a major theoretical motivation for a high energy muon collider which could be in the unique position of testing the few WIMP scenarios whose DD cross section lies well below (or in proximity) of the neutrino floor.
This paper is organized as follows: in Sec. II we summarize the main features of EW WIMP paradigm focusing on fermionic DM and relegating many details to Appendix A 1 and analogous formulas for scalar DM to Appendix A 2. This section provides a full explanation on the results of Table I while the more technical aspects of the freeze-out computation are given  in Appendix B. In Sec. III we show how large exposure DD experiments can probe most of the complex  WIMPs with the minimal inelastic and charged split-tings. In Sec. IV we study the full parameter space of the complex WIMPs, in the charged vs neutral splitting plane. We take as examples the doublet and 4plet with Y = 1/2 and the triplet and 5-plet with Y = 1. In Sec. V we summarize the indirect probes of complex WIMPs in electroweak observables and the electron dipole moment (EDM). In Appendix C we discuss the freeze-out prediction for complex WIMPs with Y = 0 and in Appendix D we illustrate the collider reach for all the WIMPs.

II. WHICH COMPLEX WIMP?
We recall here the logic of our WIMP classification as described in [1] (see also [3,[22][23][24][25] for previous work on the subject). Requiring the neutral DM component to be embedded in a representation of the EW group imposes that . . , 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. Within class ii), we can further distinguish two subclasses of complex WIMPs: a) complex representations Y = 0 and odd n; b) complex representation with Y = 0 with even n (odd n) for half-integer Y (integer Y ). The first subclass is a straightforward generalization of the models analyzed in [1] where the stability of the DM is guaranteed by an unbroken dark fermion number which can be gauged as was first done in Ref. [3]. For completeness, we give the freeze-out prediction of these scenarios in Appendix C.
In this paper we focus on complex WIMPs with Y = 0 whose phenomenology differs substantially from the one with Y = 0. We focus here on the fermionic case and leave the discussion about the scalar WIMPs to the Appendix A 2. The minimal Lagrangian for a fermionic complex WIMP with Y = 0 is: where T a is a SU(2) L generator in the DM representation. The main difference with respect to real WIMPs is that the renormalizable Lagrangian is no longer sufficient to make the DM model viable. In Eq. 1 we write only the minimal amount of UV operators required to make the DM model viable. As we discuss in Appendix A 1 these operators are also unique and possibly accompanied by their axial counterparts. These are obtained from the ones in Eq. 1 by adding a γ 5 inside the DM bilinear. We now illustrate the physical consequence of O 0 and O + in turn in Sec. II A and Sec II B and derive the implication for complex WIMPs in Sec. II C. We comment on DM stability in Sec. II D.

A. Inelastic splitting
The non-renormalizable operator O 0 is required to remove the sizeable coupling to the Z boson of the neutral component χ N of the EW multiplet This coupling would lead to an elastic cross section with nuclei already excluded by many orders of magnitude by present DD experiments [26]. After the EWSB, O 0 induces a mixing between χ N and χ c N . Replacing the Higgs with its VEV, (H c † ) σ a 2 H is nonzero only if we pick σ a = σ + , so that the new (pseudo Dirac) mass term in the Lagrangian reads (3) The Z-mediated scattering of DM onto nucleons is no longer elastic and the process is kinematically forbidden if the kinetic energy of the DM-nucleus system in the center-of-mass frame is smaller than the mass splitting where m N is the mass of the nucleus, µ is the reduced mass and v rel is DM-nucleus relative velocity. In particular, given the upper bound on the relative velocity v rel < v E + v esc , where v E = 240 km/sec is the Earth's velocity and v esc = 600 km/sec is the assumed escape velocity of DM in the Milky Way, the largest testable mass splitting is δm max 0 = 1/2µ(v E + v esc ) 2 which for xenon nuclei gives δm max 0 450 keV. The splitting for a given recoil energy is which explain why the maximal constrained splitting experimentally is δm max,exp 0 240 keV as shown in Fig. 1, given that Xenon1T [27] analyzed data only for E R < 40 keV. Extending the range of Xenon1T to higher recoil energies would be enough to probe splitting up to δm max 0 as already noticed in Ref. [6,7]. In principle, larger mass splittings can be reached using heavier recoil targets than xenon such as iodine in PICO-60 [28], tungsten in CRESST-II [29], CaWO 4 [30], PbWO 4 [31], 180 Ta [32], Hf [33] and Os [34]. However, these experiments currently do not have enough exposure to probe EW cross-sections.
A complementary bound on δm 0 comes from requiring that the decay χ 0 → χ DM +SM happens well before BBN. The leading decay channels are χ 0 → χ DM γ, χ 0 → χ DMν ν and χ 0 → χ DMē e with decay widths The first process is induced by a dipole operator generated at 1-loop for fermionic DM as computed in [35]. The three body decays are instead induced at tree-level by the EW interactions both for fermionic and scalar DM. For fermionic DM, the dipole induce decay dominates the width in the mass range of interest. In order for these processes not to spoil BBN, we have to impose the following condition on the decay rate of χ 0 : where τ −1 BBN = 6.58 × 10 −25 GeV. The lower bounds on the neutral mass splitting for fermions are shown in Fig. 1 together with those for scalars computed in Appendix A 2. The main difference between scalars and fermions is that the former are typically more long lived due to the suppression of χ 0 → χ DM γ. As a consequence the BBN bouns are stronger for scalar WIMPs.

B. Charged splitting
The operator O + in Eq. 1 is necessary to make the DM the lightest state in the EW multiplet for all the n-plets where the hypercharge is not maximal. Indeed, EW interactions induce at 1-loop mass splittings between the charged and the neutral components of the EW multiplet which in the limit m W M χ are [12][13][14] ∆M EW where δ g = (167 ± 4) MeV and Q = T 3 + Y . This implies that negatively charged states with Q = −Y are pushed to be lighter than the neutral ones by EW interactions. Notable exceptions are odd-n multiplets with Y = 0 and all the multiplets with maximal hypercharge |Y max | = (n − 1)/2 where negatively charged states are not present. For these multiplets, having y + = 0 would be the minimal and phenomenologically viable choice. Including the contribution of O + , the final splittings δm Q = M Q − M DM between the DM and the charged components read In the 31 case, no mixing between the charged components of the multiplet can occur and δm+ is a monotonic function of δm2+. In the 4 1/2 case, instead, because of the mixing induced by O0, the positively charged mass eigenstate χ + is always heavier than χ − . The splitting between χ + and χ − has a minimum of order δm0, which was taken 50 MeV in this plot for display purposes.
where sgn(Q) in Eq. 11 accounts for the presence of opposite charge states that are not related by charge conjugation, as implied by to the non-zero hypercharge of our WIMPs. The second term inside the square root comes from the mixing between the charged gauge eigenstates χ Q and χ c −Q induced by O 0 . This obviosuly vanishes for Q > (n − 1)/2 − Y .
The different charged-neutral mass splittings can all be written in terms of two independent splittings, which we choose to be δm 0 and δm Q M , where Q M ≡ Y + (n − 1)/2 is the largest electric charge in the multiplet. Since c nY Q M = 0, δm Q M is a monotonic function of y + and Eq. 11 can be inverted. In Fig. 2 we show as an example the mass splittings of 3 1 and 4 1/2 as a function of δm 2+ . In the former case, no mixing occur within the components of the multiplet and δm + is a monotonic function of δm 2+ . For the 4 1/2 , instead, the mixing induced by O 0 between the components with Q = ±1 makes the positively charged mass eigenstate χ + heavier than χ − .
In Sec. IV we will explore the parameter space spanned by δm 0 and δm Q M , fixing the thermal DM mass of every EW multiplet as shown in Table I. Crucially, the operators inducing the splitting in Eq. 1 also generate new higgs-exchange contributions to the spin-independent scattering cross-section of DM on nucleons. Therefore, the current best upper limit on the DM elestic cross section onto nucleons set by PandaX-4T [5] translates into upper bounds on the neutral and charged splittings as reported in Table I.
Charged-neutral splittings smaller than the EW one in Eq. 10 require a certain amount of fine-tuning between UV operators and the EW contribution. To quantify this we define the Fine Tuning (F.T.) where the index I runs over the three contributions in the definition of δm Q in Eq. 11. Large values of F.T. imply a significant amount of cancellation between two or more parameters.

C. Viable complex WIMPs
The EFT approach used to write Eq. 1 is meaningful only if the UV physics generating O + and O 0 is sufficiently decoupled from DM. When Λ UV approaches the DM mass the cosmological evolution of the DM multiplet cannot be studied in isolation, since the heavy degrees of freedom populates the thermal bath at T M DM and are likely to modify our freezeout predictions. To avoid these difficulties we restrict ourselves to Λ UV ≥ 10M DM .
This condition, together with the required inelastic splittings in Fig. 1, can be used to select the viable complex WIMPs. Starting from Eq. 3 and imposing δm 0 > δm min 0 , we derive the viable window for Λ UV : We are now interested in estimating for which multiplets the viable window shrinks to zero. Setting y 0 = (4π) 4Y in Eq. 13 that is the largest value allowed by Naive Dimensional Analysis (NDA) we derive the values of n and Y having a non zero cutoff window in Eq. 13. These are for both scalar and fermionic WIMPs n 1/2 multiplets with n ≤ 12 together with the 3 1 and the 5 1 mupltiplets. This result can be understood as follows. The upper bounds on n for Y = 1/2 multiplets come from the perturbative unitarity of the annihilation cross section as discussed in Appendix B 2. The maximal cutoff required to obtain the phenomenologically viable splittings is of order ∼ 10 7 TeV as shown in the fifth column of Table I. This is many orders of magnitude larger than the DM masses allowed by freeze-out so that Eq. 13 results in a wide range of allowed Λ UV .
For Y = 1 multiplets the maximal required cutoff is of order ∼ 10 2 TeV so that the n-dependence of the allowed window in Eq. 13 becomes relevant. Given that the DM mass grows with n as M DM ∼ n 5/2 and the required cutoff stay approximately constant, we expect the allowed window to shrink to zero for large n. Numerically we find that the last allowed multiplet has n = 5.
We refer to Appendix A 2 for a similar argument for scalar WIMPs. These have a slightly different parametric which however results in the same viable EW multiplets of the fermionic case.

D. DM stability
The Lagrangian in Eq. 1 preserves the DM number and as a consequence the DM is automatically stable. Gauge invariant interactions beyond those of Eq. 1 can however induce fast DM decay. Here we discuss whether the effect of these operators can be small enough to allow the DM to be accidentally stable. Throughout this discussion we not only require the DM lifetime to be long enough to circumvent cosmological bounds [37,38] (τ DM 10 19 sec) but also to satisfy the stronger astrophysical bounds on decaying DM [39][40][41] (τ DM 10 28 sec).
For 2 1/2 and 3 1 we can write renormalizable interactions χ c He R and χL c H that break the DM number and lead to a fast DM decay. These EW multiplets require a DM number symmetry, for example a discrete Z 2 -symmetry acting only on the DM field to provide a viable DM candidate.
For even multiplets with Y = 1/2 and n > 2, we can write the following series of higher dimensional operators inducing DM decay where the Higgs bosons are appropriately contracted to make every term an SU (2) singlet. For 4 1/2 , the leading contribution to DM decay comes from the operators with coefficients C 1 and C 4 , while for n ≥ 6 this is given by the operator with coefficient C 5 by closing the χ loop and attaching (n − 2)/2 W µν operators, in a similar fashion to real candidates [1]. In all cases, we need at least Λ UV > 10 10 M DM for O(1) Wilson coefficients to preserve DM stability. This lower bound is incompatible with the upper bound on the UV physics scale required to generate the inelastic splitting δm 0 . This result implies that the DM stability depends upon the properties of the UV physics generating the neutral splitting in Eq. 3. For the 5 1 WIMP the lowest dimensional operators inducing DM decay are (15) These require Λ UV > 10 10 M DM to ensure DM accidental stability which is again incompatible with the upper bound of Λ UV < 20M DM needed to generate the inelastic splitting.
As a result, none of the complex WIMPs with Y = 0 can be accidentally stable in the sense of Minimal Dark Matter [22]. Specific UV completions of the physics generating the inelastic splitting in Eq. 3 might allow for DM accidental stability. This question goes beyond the scope of this study. Conversely, complex WIMPs with Y = 0 can be accidentally stable thanks to the gauging of the unbroken U (1) flavor symmetry in the DM sector [3]. The freeze-out predictions for these millicharged WIMPs are given in Appendix C.

III. MINIMAL SPLITTING
In this section we discuss the DD signal of Complex WIMPs in the minimal splitting scenario, that is when the mass splittings are chosen to be the smallest possible allowed by the requirements of the previous section. For fermions, we set δm 0 = 220 keV for n 1/2 WIMPs with n ≤ 4 as well as for the allowed n 1 WIMPs. For n 1/2 WIMPs with n > 4 the BBN bound in Fig. 1 gives the minimal δm 0 . For scalars instead, we set δm 0 = 4.9 MeV for n 1/2 and δm 0 = 3.7 MeV for n 1 WIMPs, both coming from BBN.
In this minimal setup, 2 1/2 and 3 1 stand out as the only two multiplets where the DM is automatically the lightest state, with a splitting with the Q = 1 state given by the pure EW splitting in Eq. 10 that is 354 MeV for 2 1/2 and 542 MeV for 3 1 . For all the other WIMPs a UV generated splitting is needed to make the DM lighter than the negatively charged states. We fix this splitting to the smallest possible value that gives the DM as the lightest state. As can be seen from Eq. 11 this requires a UV splitting which is of the order of the EW one. We want now to study the phenomenology in direct detection in Sec. III A and at a future muon collider in Sec. III B.  Table I. The light green shaded region is excluded by the present experimental contraints from XENON-1T [36] and PandaX-4T [5], the green dashed lines shows the expected 95% CL reach of LZ/Xenon-nT [8,9] and DARWIN [10,11].

A. Direct Detection prospects
The spin independent scattering cross-section σ SI of DM on nuclei receives two contributions: i) from purely EW loop diagrams ii) from Higgs mediated tree-level diagrams generated by both O 0 and O + . For minimal splitting Higgs mediated scattering is subdominant and σ SI can be computed by considering only EW loop diagrams.
Following [17,42], the Lagrangian describing the spin-independent (SI) DM interactions with quarks and gluons is q is the quark twist-2 operator. The Wilson coefficients are given by [17] f EW where m h is the mass of the Higgs and κ c = 1.32, with c w , s w being the cosine and the sine of the Weinberg angle, respectively. The terms proportional to Y correspond to the exchange of Z bosons inside the EW loops.
After the IR matching of these interactions at the nucleon scale [42], we can express σ SI per nucleon (for where m N is the mass of the nucleon and where the nucleon form factors are defined as (2))/m N , and q(2),q(2) are the second moments of the parton distribution functions for a quark or antiquark inside the nucleon [17]. The values of these form factors are taken from the results of direct computation on the lattice, as reported by the FLAG Col-laboration [43] in the case of N f = 2+1+1 dynamical quarks [44,45]. 3 The direct detection reach on the different complex multiplets for minimal splitting is summarized in Fig. 3. As can be seen by the figure, most of the WIMPs can be probed for minimal splitting by future large exposure direct detection experiments like DARWIN with the notable exception of the 2 1/2 and 5 1 that we discuss below.
Combining Eq. 17 with Eq. 18 the parametric expression for σ SI is where ξ = 16.6 ± 1.3 with the error coming from the lattice determination of the nucleon form factors. This formula makes evident that large cancellations with respect to the natural size of the elastic cross section can take place when Y (n 2 − 1)/ξ. In particular, for n = 2 the exact cancellation takes place at Y 0.44 ± 0.02, which almost matches the exact hypercharge of the doublet leading to the large uncertainty in Fig. 3. Similarly, the cancellation happens at Y 1.2 for n = 5, which explains why in this case the signal entirely lies below the neutrino floor, while it is within DARWIN reach for their real or millicharged counterparts as shown in Ref. [1] and Fig. 7 respectively.
Spin dependent (SD) interactions of DM with the nuclei are also induced by EW loops and can lead to a larger cross section compared to the SI one [17]. Unfortunately, the predicted SD cross section for all the complex WIMPs lies always well below the neutrino floor and it will be impossible to test even at future direct detection experiments.

B. 2 1/2 and 31 at muon colliders
Conversely to the case of real WIMPs studied in Ref. [1], large exposure experiments like DARWIN will not be able to entirely close the window of Complex WIMP candidates. This further motivates us to study the expected reach on these DM candidates at a high energy lepton collider such as a muon collider running well beyond TeV center of mass energy.
In this Section we fix the luminosity as a function of the center of mass energy as [46] and we concentrate on machines running at √ s = 3, 6, 10 TeV [20,21]. In Appendix D we provide results for higher center of mass energies, as well as exhaustive summary tables and plots for all the dark matter candidates at center of mass energies up to 30 TeV, including estimates for scenarios in which the luminosity is a free parameter not fixed by Eq. 20.
We focus on the collider reach for the complex doublet (2 1/2 ) and the complex triplet (3 1 ), that are the lightest WIMPs and have the greater chance to be discoverable at √ s ≤ 10 TeV. Theoretically, these candidates are the most minimal complex WIMPs since they have maximal hypercharge and the neutral component is automatically the lightest one at the renormalizable level. The only required higher dimensional operator is O 0 which generates the inelastic splitting in Eq. 3. The sensitivity of high energy muon colliders to heavier multiplets will be discussed in detail in Appendix D.
The expected signatures at future colliders depend very much on the lifetime of the charged states in the EW multiplet. These might decay back to the neutral DM either promptly or with a macroscopic lifetime on detector scales.
In the prompt case, the decay of the charged state will give rise to soft radiation which would be impossible to disentangle from the SM background and the DM signal will be characterized only by a large missing energy recoiling against EW radiation. This can be further resolved from the SM background at lepton colliders thanks to the knowledge of the center of mass energy by looking for features in the missing mass (MIM) distribution. We systematically studied the different kinds of EW radiation in Ref. [1]. We perform a similar study here (details are given in Appendix D), finding that the MIM reach is dominated by the mono-W and the mono-γ channels. Despite the non-zero hypercharge enhancing the coupling to the Z, the mono-Z channel has a lower reach, similarly to what was found in Ref. [1].
We estimate the reach of MIM searches from cutand-count experiments using the following estimate for the statistical significance where sys parametrizes the expected systematic uncertainties on the measurement. In the long-lived case, depending on the decay length of the charged state it might be beneficial to look at disappearing or "stub" tracks (DT) as studied in Ref. [18], or at long-lived charged tracks that could be distinguished from the SM backgrounds by their energy losses in the material and by the measurement of their time of flight (see for example analogous searches performed by ATLAS and CMS at the LHC [47,48]).
Here we focus on the expected reach from DTs and MIM searches which are relevant for the 2 1/2 and 3 1 , fixing the splitting of the Q = 1 state with respect to the neutral one to be the one generated purely by EW corrections. In this minimal case the rest-frame lifetime of the charged state decaying back to the DM is cτ + ≈ 6.6 mm for the 2 1/2 and cτ + ≈ 0.77 mm for the 3 1 . Our results for the combination of the MIM searches and the DT tracks search are given in Fig. 4 for fixed thermal mass of the 2 1/2 , and similarly for the 3 1 . We provide details on the estimate of the MIM reach, as well as for the DT signature, in App. D, also giving results as a function of the mass M χ of the candidate.
As can be seen from Fig. 4, the MIM search at √ s = 3 TeV with the benchmark luminosity L = 1 ab −1 is not sensitive to the 2 1/2 . A √ s = 6 TeV collider with benchmark luminosity L = 4 ab −1 , instead, is able to probe the doublet WIMP at 2σ C.L. In general the mono-W and mono-γ channels give comparable mass reach for the benchmark luminosity adopted here, full detail is given in Tab. III. The sensitivity of each channel has a specific behavior as a function of the luminosity and for different assumed systematics. The overall combination as function of the total luminosity for fixed thermal mass is shown in Fig. 4. We remark that the reach deviates from a pure rescaling by √ L because the selections, hence the result, have been optimized as a function of L when dealing with sys = 0.
In Fig. 4 we also display results for the DT search as a function of the lifetime of the charge +1 state, that is in a 1-to-1 relation with the mass splitting δm + . We see that a collider of √ s = 6 TeV can probe the 2 1/2 for charge-neutral splitting generated purely by EW interactions. A collider of √ s = 3 TeV can probe a large portion of the allowed lifetimes for the charged tracks corresponding to non-zero UV contributions in Eq. 11.
In Fig. 4 we show similar results for the 3 1 at its thermal mass, which entails interesting results at 6 and 10 TeV center of mass energy machines. For the 10 TeV collider we find that MIM searches are effective probes of this WIMP candidate and can establish a bound at 95% CL with a small luminosity or give a discovery with the nominal luminosity. Mono-W and mono-γ perform similarly well and their combination is worth being done. The DT search cannot probe the EW splitting for the 3 1 , however it can cover a large portion of allowed lifetimes in the non-minimal case in which UV physics is contributing to the chargedneutral splitting.
Finally, it is important to notice that with an O(1) fine tuning of the UV contribution to the charged neutral splitting against the irreducible EW contribution the charged tracks could be extremely long-lived on detector scale for both the 2 1/2 and 3 1 . This observation motivates a detailed detector simulation of a muon collider environment to reliably estimate the expected reach in this channel. We will provide some estimates for the signal yield in Fig. 5 in the next section.

IV. NON-MINIMAL SPLITTINGS
Exploiting the freedom to change the spectrum given by the UV operators in Eq. 1 one can have significant changes in the phenomenology described in Sec. III. We dedicate this section to explore the phenomenological consequences of releasing the minimal splitting assumption.

A. Direct detection
The Higgs portal operators in Eq. 1 generate upon EWSB a linear coupling of the DM to the Higgs boson of the form This coupling mediates tree-level SI scattering processes of DM onto nuclei, therefore it can be constrained by direct detection experiments. As the mass splitting δm 0 has to be sufficiently large to suppress the scattering mediated by the Z boson, we find that the allowed parameter space for the nonrenormalizable couplings is compact. Following Ref. [49] we can integrate out the Higgs boson and write the couplings of the DM to the SM so that the matrix elements f q and f G required to compute σ SI are simply given by where f EW q and f EW G come from EW loops and are given in Eq. 17. We can rewrite the coupling λ D solely in terms of mass splitting as where and ∆M EW Q M is the gauge induced mass splitting in Eq. 10.
Replacing Eq. 24 into Eq. 18 allows us to translate the upper bound on σ SI into an upper bound on δm 0 and δm Q M .

B. Collider searches
The splittings δm 0 , δm Q M determine the lifetime of the charged components of the EW multiplet, which are pivotal to understand the viable collider signatures. In what follows we take as reference the detector geometry proposed in [18] and we classify the pa-rameter space (δm 0 , δm Q M ) in various regions according to the lifetime τ LCP of the Longest Lived Charged Particle (LCP) of a given multiplet.
For cτ LCP > 1 m the LCP gives long charged tracks (CT) with an average length roughly corresponding to the middle layer of the outer tracker. The SM background processes for this "long" track with anomalous properties strongly depends on the properties of the detector, therefore its study is outside the scope of this work. We limit ourselves to estimate the number of expected signal events of charged χ pair production at a muon collider, as shown in Fig. 5. We highlight that this results approximately hold for generic √ s in the domain of tens of TeV, as long as Eq. 20 for the luminosity holds. In Fig. 5 we show the βγ of the produced charged particle, which plays a crucial role to disentangle these special tracks from the SM background.
Our estimates for the LCP signature is based on the counting the number of LCP produced. For the LCP we require p T > 200 GeV, |η χ | < 2, inspired by LHC searches [48]. For the 2 1/2 and 3 1 we include in our counting all charged particles production, assuming that χ 2+ promptly decay to the long lived χ + . For the 4 1/2 and 5 1 , which have more complicated spectra, we stick to the minimal splitting scenario for the estimate of the charged tracks yield, which corresponds to the minimal y + necessary to lift χ − above χ DM . In this spectrum configuration χ − is the most natural and only candidate to make long charged tracks, therefore we only consider this contribution in our result. 4 For 1 m > cτ LCP > 0.33 cm the lightest charged WIMP gives a disappearing track signal. This search has been introduced in Sec. III and it is further described in Appendix D, in particular in Fig. 8 where the reach in the plane mass-lifetime is displayed for judiciously chosen center of mass energies for every WIMP up to the 5 1 .
For cτ LCP < 0.33 cm the lightest charged WIMP decays promptly on collider scales and gives missing energy signatures such as the mono-W and mono-γ discussed in Sec. III. The choice for the threshold value to consider the signal as prompt corresponds to cτ /2 of the charged track for the 2 1/2 WIMP in the minimal splitting scenario, which is a known benchmark where DT reconstruction starts to become challenging. Full details on MIM searches are given in Appendix D and in particular in Fig. 9.

C. Parameter space for complex WIMPs
In the rest of this section, we examine in greater detail the parameter space spanned by δm 0 and δm Q M specifically looking at each of the lightest multiplets of Tab. II up to the 5 1 WIMP. This will allow us to discuss the salient phenomenological features of the complex WIMPs. The 2 1/2 and 3 1 are special WIMPs because they are the only multiplets with maximal hypercharge compatible with our assumptions. In particular, for 2 1/2 requiring δm 0,+ > 0 automatically implies the neutral WIMP candidate is the lightest one. Perturbativity requires δm 0 < 40 MeV for the n = 3 1 WIMP, because of the strong suppression of the 7 dimensional operator generating the neutral splitting. The narrower range of δm 0 for the 3 1 with respect to 2 1/2 was expected from the higher dimensionality of O 0 for Y = 1 as compared to Y = 1/2.
In Fig. 6 we show the constraints on the Dirac 2 1/2 and 3 1 in the parameter space spanned by δm 0 and δm Q M coming from present DD experiments like Xenon1T and PANDAX-4T, as well as prospect from future high exposure Xenon experiments, i.e. LZ, XENONnT, DARWIN and DARWIN/G3. As we can see from Eq. 24, these bounds depend solely on the combination In particular we find . (27) In the region of low mass splitting for both 2 1/2 and 3 1 the direct detection cross section lies below the neutrino floor. For 2 1/2 , Eq. 27 shows that large cancellations between EW loops and tree-level Higgs exchange occur around δµ 300 MeV, while for 3 1 the minimum σ SI is obtained for δµ = 0 and falls below the neutrino floor. To produce a direct detection cross section above the neutrino floor for the 2 1/2 (3 1 ) we need δµ > 1.6 GeV (δµ > 2.0 GeV). PANDAX-4T already excludes mass splitting δµ < 20 GeV (δµ < 30 GeV) for the 2 1/2 (3 1 ).
All in all, DD still leaves a large portion of parameter space unconstrained between the neutrino floor and the PANDAX-4T constraints. Remarkably, this region corresponds to mass splittings that do not require tuned adjustments of the three contributions to δm 0 and δm + . The region of large splittings down to δµ ∼ 1 GeV can be covered by large exposure Xenon experiments while the large portion of the parameter space lying below the neutrino floor should be taken as a major motivation for a future muon collider.
For 2 1/2 WIMP the neutrino floor region can be fully probed only via a combination of charged tracks, DT and MIM searches. The different search strategies become relevant depending on the δm + value, as shown in Fig. 6. For the 3 1 WIMP stub and charged tracks can exclude the entire neutrino floor region, while MIM and DT searches can be complementary to large exposure DD experiments to probe the rest of the parameter space.
While the plot in Fig. 6 show the possible regions in which DT can be reconstructed efficiently, they do not show if the reach is large enough to guarantee exclusion (or discovery) of the WIMP candidate. For these two WIMP candidates a detailed discussion was given in Sec. III B. The general message is that for √ s sufficiently larger than twice the thermal DM mass DT searches will be powerful enough to probe the parameter space.

The Dirac 4 1/2 and 51
Being the first multiplet with non-maximal hypercharge, the constraints coming from DM stability shape the allowed parameter space of the 4 1/2 . In particular, we must require δm − to be positive. This constraint excludes the light blue shaded region in Fig. 6 bottom left.
The region where the direct detection cross section lies below the neutrino floor is reduced to a tiny band for 4 GeV < δµ < 15 GeV. This region can be probed by direct searches at a future high energy muon collider with √ s > 10 TeV because of the thermal mass of the 4 1/2 lies around 4.8 TeV. In particular, stub and charged tracks searches could cover almost the entire neutrino floor region, except for a small portion accessible only to MIM searches. The rest of the viable parameter space for the 4 1/2 can be in principle probed by large exposure direct detection experiments.
We show the parameter space of the 5 1 WIMP in Fig. 6 bottom right. The dominant constraint comes from the perturbativity of the operator generating the neutral splitting, which requires δm 0 < 3 MeV. As a consequence, testing larger splittings for inelastic DM at Xenon1T can already probe large portions of the allowed parameter space of the 5 1 . DM stability requires 1 GeV δm 3+ 2 GeV. The parameter space can be fully excluded by charged and stub track searches if a muon collider of √ s > 20 TeV will be constructed. Besides, all the allowed mass splittings except a small window 1.8 GeV δm 3+ 2.0 GeV can be probed by DARWIN/G3.

V. INDIRECT PROBES OF COMPLEX WIMPS
In this section we briefly comment on indirect probes of complex WIMPs, focusing in particular on the distinctive imprints of their Higgs interactions on precision electroweak observables (Sec. V A) and on the electron dipole moment (Sec. V B).

A. Electroweak precision observables
The addition of new EW multiplets to the SM leads to deviations in EW observables which could be tested with LEP data, at the LHC and at future colliders [50,51]. Since the presence of new EW multiplets affects mainly gauge bosons self energies, their indirect effects can be encoded in the oblique parameters [52,53]. Here we briefly summarize the main features of these contributions following the notation of Ref. [53]. The most relevant oblique parameters can be related to the SM gauge boson vacuum polarizationŝ where Π ij appear in the kinetic terms of the EFT describing the SM vectors interactions at energies smaller than the DM mass In particular, all the EW multiplets (both real and complex) give a universal contribution to theŴ ,Ŷ as previously found in [22] where κ = 1, 1/2, 1/8, 1/16 for Dirac fermions, real fermions, complex scalars and real scalars respectively. In the numerical estimates we normalized the expected contributions to the fermionic 2 1/2 . The rough expectations for heavier WIMPs can be obtained from the equations above by remembering that M DM ∼ n 5/2 . The contribution toŴ scale like as ∼ 1/n 2 while the one toŶ as ∼ 1/n 4 .
W andŶ induce effects in SM observables that grow with energy so that LHC searches have already ameliorated the sensitivity on these operators compared to LEP. In particularŴ has been recently bounded by CMS [54] andŶ is also set to be tested with a similar precision [55]. Even if the current precision is not sufficient to probe the WIMP thermal masses, further improvements are expected in a future high energy muon collider. This could provide interesting indirect tests of EW WIMPs [51,56,57]. The size of the expected WIMP contributions in Eq. 30 and Eq. 31 the results of Ref. [57] indicates that a muon collider of √ s = 30 TeV would be needed to observed these deviations.
As extensively discussed in this paper, complex WIMPs require contact interactions with the SM Higgs to be phenomenologically viable. These in-teractions give an additional contributions to theŜ parameter. Moreover, the mass splittings inside the EW multiplet induced by O 0 and O + break the custodial symmetry in the Higgs sector. As a consequence, complex WIMPs give irreducible contributions also to theT parameter. The explicit formulas for fermionic WIMPs arê where s 2 2W = sin 2 2θ W 0.83 while δµ Q M and δµ 0 are defined in Eq. 25. The corresponding formulas for scalar WIMPs are easily obtained from the fermionic ones by replacingŜ S = 1 4Ŝ F andT S = −T F . Our result forT S is consistent with the result in Ref. [58], first computed in [59].
BothŜ andT are currently constrained at the level of 3 × 10 −3 by LEP data (where the precise value will of course depend on the correlation between these two parameters [53]). High luminosity colliders further exploring the Z pole like FCC-ee will improve this precision typically by a factor of 10 [60] with the ultimate goal of pushing the precision of EW observables to the level of 10 −5 [61].
ConcerningŜ we conclude that the deviations induced by complex WIMPs are unlikely to be visible unless for light multiplets with the ultimate EW precision. Conversely, at fixed large splitting and increasing size of the multiplet, the deviations onT are enhanced to the point that a 12 1/2 multiplet is givinĝ T 3 × 10 −3 , which is within the reach of current LEP data and could help reducing the tension the m W mass extracted from the EW fit and the one measured directly by the CDF collaboration [62]. Of course this extreme scenario requires very large couplings at the boundary of perturbativity for the contact operators inducing the splitting in Eq. 1.

B. Electric dipole moment
As mentioned in Sec. II and detailed in Appendix A 1, the operators generating the splittings for fermionic WIMPs have in general two different chiral structures and different relative phases. This generically induces operators contributing to the electron EDM at 1-loop [63,64]. Following the notation of Ref. [63], the leading operator generated by WIMPs loops is where we neglect for this discussion operators such as |H| 2 BB and H † σ a HWB that would be suppressed by powers of Y /n with respect to the one with two SU (2) field strength insertions. The Wilson coefficient in Eq. 34 is constrained by the recent measurement of the ACME collaboration [65] to be c W W < 4 × 10 −3 g 2 2 Λ UV 10 TeV where the index y I controls the contribution from O 0 or O + and y I,5 the ones from their axial counterpart. If only O 0 is present, like it would be for instance in the minimal splitting case for the 2 1/2 multiplet, the bound in Eq. 35 can be rewritten as an upper bound on the neutral splitting  Table I. Interestingly, future experimental upgrades are expected to be sensitive to electron EDM as small as d e 10 −34 e cm [66] which would provide sensitivity to the WIMP parameter space down to splittings of order ∼ 50 MeV for O(1) couplings and CP-violating phases.
It would be interesting to further explore this direction in models where the O(1) CP-violating phase is motivated by a mechanism producing the observed baryon asymmetry in the Universe such as WIMP baryogenesis [67,68].

VI. CONCLUSIONS
In this study we derived the first full classification of WIMP in complex representations of the EW group. Our only assumption is that below a certain UV scale the SM is extended by the addition of a single EW multiplet containing the DM candidate. Our main results are contained in Table I, where we give all the freeze-out predictions for the WIMP thermal masses, the ranges of phenomenologically viable splittings and the maximal scale separation between the DM mass and the UV physics generating the splittings.
We then inspect the phenomenology of complex WIMPs with the goal of understanding what would be required experimentally to provide the ultimate test on the EW nature of DM. We focus on the interplay between future large exposure direct detection experiments and a future high energy lepton collider, e.g. a muon collider. The reach from the current and forthcoming indirect detection experiments of pure electroweak multiplets are very promising (see e.g. [69,70] for current reaches and [71,72] for prospects). Neverthless a robust assessments of the sensitivities for all the complex WIMPs classified here and the real WIMPs classified in [1] would require a careful evaluation of the expected annihilation cross sections and an assessment of the astrophyisical uncertanties [73]. We leave this work for a future investigation. Similarly, the prospects of discovering WIMPs via kinetic heating of compact astrophysical objects, like neutron stars [74,75] or white dwarfs [76], should be further assessed by the future James Webb Space Telescope infrared surveys [77].
Our main result is that a kiloton exposure Xenon experiment would be able to probe most of the complex and the real WIMPs, with the notable exceptions of the complex doublet with Y = 1/2 whose cross section lies naturally well below the neutrino floor. This result for the complex doublet was of course well known from the many previous studies on the SUSY Higgsino [78], but it gets further substantiated in the context of our WIMP classification. Our findings constitute a major motivation to push forward the research and development of a future high energy lepton collider. Specifically, we find that a future muon collider with √ s = 6 TeV would be able to fully probe the existence of a fermionic complex WIMP doublet at its thermal mass around 1.1 TeV.
We further study the parameter space of the different complex WIMPs and find regions where the direct detection cross section can drop below the neutrino floor because of accidental cancellations of the EW elastic scattering (like for the 5 1 ) or at small (mildly tuned) values of the charged-neutral mass splittings (like for the 3 1 ) or for the destructive interference of EW and Higgs induced scattering (like for the 4 1/2 ). In all these cases a high energy lepton collider might be again the only way of discovering these WIMPs.
Interestingly, the list of WIMP candidates we furnished provides a series of targets that can be probed at successive stages of a future machine. Moreover, our analysis shows that future searches for long-lived charged tracks will be crucial to fully probe the WIMP parameter space. This strongly motivates a detailed collider study assessing the expected sensitivity of these searches in a realistic muon collider environment. We hope to come back to this issue in the near future.

ACKNOWLEDGMENTS
We thank Ranny Budnik and Itay Bloch for useful insights on Xenon1T, Scott Haselschwardt for enlightening discussions about the future of Xenon experiments, Rodolfo Capdevilla, Federico Meloni, Rosa Simoniello and Jose Zurita for comparison on the disappearing track prospects. Simon Knapen, Ennio Salvioni and Andrea Tesi for discussions.

Appendix A: Complex WIMP classification
In this appendix we give further details on our complex WIMP classification. In Sec. A 1 we explicitly classify the possible UV operators, substantiating the result of Sec. II. In Sec. A 2 we give the results for scalar WIMPs.

UV Operators
In Sec. II we showed that in order to make complex WIMPs with Y = 0 viable, new UV sources of splitting between the charged and neutral components and among the neutral components are necessary. In Eq. 1 we showed how the neutral splitting can be generated by O 0 and the charged-neutral splitting by O + . Here we take a step back and investigate the generality of this choice.
The general form of an operator responsible for a Majorana mass term for χ N after EWSB is where the H 4Y is necessary to match the hypercharge of the χ 2 piece, while I is a SU(2) L tensor. We are crucially assuming that the Higgs is the only scalar picking a VEV after EWSB. Under this assumption, since the Higgs is a boson, the only surviving SU(2) L structure is the totally symmetric combination. This can be seen as the symmetric combination of 2Y Higgs pairs in the isotriplet representation (the isosinglet is antisymmetric and vanishes identically), so that we are left with O 0 defined in Eq. 1 and its axial counterpart with the γ 5 insertion. The latter can be shown to give only subleading contributions to the mass splitting between the neutral components of χ. Assuming y to be real, the shift on δm 0 induced by y 0,5 = 0 with respect to its expression in Eq. 3 is which is highly suppressed for Λ UV > M DM > v. The operator O 0 in Eq. 1 is then the dominant contribution to the inelastic splitting among the neutral components.
The isospin structure and field content of O + is already the minimal required to generate additional splitting between the charged components, the only possibility is again to change its chiral structure writing the general operator inducing charged-neutral splitting as χT a (y + + y +,5 γ 5 )χH † σ a H . (A3) Similarly to the neutral case, the γ 5 insertion leads to a subleading shift in the mass splittings with respect to the value of Eq. 11 Finally, we comment on the implications on DD signals. Both operators with the γ 5 insertions give additional contributions to the SI cross-section and match to the effective operators: However, all these operators are strongly momentumsuppressed [79,80], so that their corrections to σ SI are expected to be (q/m N ) 2 ∼ 10 −6 smaller than those from Eq. 16, and thus negligible.

Complex Scalar WIMPs
Following the discussion in Sec. II for the fermions, supported by the previous Appendix A 1, the minimal Lagrangian for a scalar complex WIMP is , Conversely to the fermionic case, no additional operators can be written to generate the fundamental mass splitting. The neutral and charged squared mass splitting, µ 2 Q = M 2 Q − M 2 DM , can be written as The linear mass splittings, δm Q = M Q − M DM , are then given by δm Q = µ 2 Q /(2M DM ). Notice that for Y = 1/2 all the operators are renormalizable and any dependence on the cutoff disappears. The lower bound on δm 0 from DD is identical to that for fermions, since σ SI does not change. Instead, the BBN bound differs since the one-loop decay channel χ 0 → χ DM γ now is now heavily suppressed with respect to the three-body decays. As a consequence, the BBN condition becomes Γ χ0 ≡ Γν ν + Γē e > 6.58 × 10 −25 GeV , which explains why the bound on δm min 0 4 − 5 MeV for scalars is so much stronger than the one for fermions, as shown in Fig. 1. For Y = 1/2, the lower bound on δm 0 sets the upper bound of the allowed window for Λ UV which is the analogous of Eq. 13 for scalar WIMP. Once set y 0 = (4π) 4Y to its NDA maximal value and given the scaling M DM ∼ n 5/2 Eq. A8 can be used to determine the viable EW multiplets. It turns out that the allowed multiplet are the same as for fermions, i.e. all even n 1/2 plus 3 1 and 5 1 .
Finally, we discuss DM stabiity for scalars. All the Y = 1/2 scalars are never accidentally stable, since DM decay can be induced by the renormalizable operator χ 2 χ † H c . Similarly, for 3 1 we can write χ † H 2 , while for 5 1 the lowest dimensional operators are: which ensures stability for Λ UV > 10 21 M DM , way larger than Λ UV < 20M DM from δm 0 perturbativity. Results for the collider reach on the scalars are given in Fig. 10 in Appendix D.

Appendix B: Complex WIMP cosmology
In this Section, we briefly describe the cosmological evolution of the Complex WIMP candidates. In par-ticular, we highlight the main differences with respect to real WIMPs and refer the reader to [1,81,82] for further details.

Sommerfeld Enhancement and Bound State Formation
In the non-relativistic regime, long range interactions in the form of Yukawa or Coulomb potentials significantly affect the WIMP annihilation dynamics. In particular, these interactions can distort the wave function of the incoming particles leading to Sommerfeld Enhancement (SE) or lead to WIMP Bound State Formation (BSF). For WIMPs, the potential induced by EW interactions is attractive for all isospins I √ 2n, I denoting the dimension of the isospin representation of the pair of particles. In these channels, the SE increases the annihilation cross-section by a factor which, for small velocities, can be approximated as Hence, the annihilation cross-section reads where σ I ann is the perturbative, hard annihilation cross-section in the isospin channel I. The SU (2) L symmetric limit approximate better the exact result the more we increase the dimensionality of the EW multiplet. This is because most of the freeze-out dynamics for large n-plets occurs before the EWSB. For smaller multiplets the symmetric approximation generically fails because EWSB become important. In our computation, we always assume the symmetric limit for n ≥ 6, while we compute the SE in the broken phase for smaller multiplets.
Besides, the mass splittings generated by O 0 and O + can alter the pattern of resonances in the SE, so that the thermal mass becomes a function of δm 0 and δm + as well. In Table I we assume the minimal δm 0 , where the effect of the splitting is safely negligible, and account for the dependence on δm + in the theoretical uncertainty. More details about the effect of these UV operators are discussed in Sec. B 3. Another source of uncertainty in the computation of σ SE ann comes from neglecting partial waves higher than the s-wave. We estimate the correction to σ SE ann as In addition to the SE, attractive isospin channels lead to BSF through the emission of one (or more) EW gauge boson in the final state χ i + χ j → BS i j + V a , where V a can be either a weak or hypercharge gauge boson. In the non-relativistic limit, and for single gauge boson emission, the BSF process is induced by an electric dipole moment emission whose dynamics is encoded in the effective Hamiltonian described in Ref. [81,82]. The main difference with respect to real WIMPs are the different selection rules selecting the bound states allowed by symmetries. In the real case, the bound states are of the form χχ and the (anti-)symmetry of the wave function allows only for those BS satisfying P BS ≡ (−1) L+S+ I−1 2 = 1 , while for complex WIMPs the bound states are of the form χχ and the above selection rule no longer applies. The BSF cross-section can be computed in the SU (2) L symmetric limit with good approximation and scales as: where the first and second term account for weak and hypercharge vector boson emmission, respectively, E B I ≈ α 2 eff M DM /4n 2 B -n B being the energy level -is the binding energy and a B = 1/α eff M DM the Bohr radius. Once the BS is formed, it can annihilate into SM, decay to lower lying BS, both with rate Γ ≈ α 2 2,eff α 3 eff M DM or Y 2 α 2 Y α 3 eff M DM , depending on the mediator, or be broken by the interactions with the plasma (ionization). The ionization rate limits the efficiency of DM annihilation through BSF and it is Boltzmann suppressed as ∼ e E B I /T . The Boltzmann suppression gets bigger as we increase the dimensionality of the multiplet, due to the larger binding energies, so that neglecting BS ionizations is a good approximation for large EW multiplets. Under this approximation, the overall DM effective annihilation cross-section reads that is, once a BS is formed it eventually annihilates to SM without being destroyed. In [1] we checked the validity of this approximation for n ≤ 7 and found that it overestimates the thermal mass of about 5 TeV. We assume this holds also for complex WIMPs and include this uncertainty in the estimate of our theoretical uncertainty. In the symmetric limit the effect of the mass splittings generated after EWSB are not taken into account. We expect these to make the heavier components of the multiplet decouple earlier than the lighter ones, thus reducing the cross-section. In order to estimate the error due to this approximation, we compute the thermal mass first setting to zero σ B I for T < max δm Q and then including BSF until the DM abundance saturates. These two effects are the dominant sources of error for n 8. For larger multiplets, the effective cross-section in Eq. B5 suffers from additional theoretical uncertainties. Estimating from the results of Ref. [83] the NLO correction to the potentials controlling the SE we get where β SM 2 = −19/6 is the SM contribution to the SU(2) L beta function. Eq. B6 results in a correction to the Sommerfeld factor S E , which affects both σ SE ann v rel and σ B J v rel as introduced in Eq. B4. NLO corrections start dominating the theoretical uncertainty for n 8.

Perturbative Unitarity Bound
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,84,85]: 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 Bi selects the BS contributions that can be formed by J = 0 initial wave. As we noticed in the previous section, the BSF contribution is larger for complex WIMPs compared to the real case due to the larger multiplicity of bound states that can be formed. While for scalar WIMPs only p-orbitals bound states contribute to the s-wave unitarity, for fermionic WIMPs additional contributions arise from selecting the projection onto the J = 0 partial wave of the BSF cross-section of ns S=1 I states.

Impact of mass splittings
Here we discuss under which conditions the effect of the UV splittings can be neglected in the prediction of the thermal mass. The UV splittings can affect the freeze-out computation in two ways: i) they directly contribute to the DM annihilation cross-section into Higgs bosons, ii) after the EWSB, the interaction Eq. 22 generates a Higgs-mediated Yukawa potential, thus affecting the Sommerfeld enhancement. For definiteness, we focus on fermionic WIMPs and on the O 0 contributions. We checked that similar conclusions hold for scalar WIMPs and for the contrbutions from O + .
For Y = 1/2 WIMPs, the hard annihilation crosssection into higgses induced by O 0 can be estimated as where T R = n(n 2 − 1)/16 is the Dynkin index of the DM SU (2) L representation. This should be compared with the typical EW hard cross-section, which is The condition σ 2H /σ EW < 1 can be translated into an upper bound on the mass splittings. For instance for the inelastic mass splitting we get Since M DM ∼ n 5/2 we get stronger upper bound on the splitting for large multiplets. All in all, we do not expect significant changes in the picture presented in Sec. IV even though for large splittings (above 1 GeV) one can get O(1) effects on the thermal freezeout predictions for small multiplets. For Y = 1 the contribution from O 0 , which now is a 7D operator which controls the 4-body process χχ → 4H whose cross section can be estimated as The condition σ 4H /σ EW < 1 translates into which, for n = 3, 5 is much looser than the upper bound on δm 0 . As a consequence, the UV splittings do not impact the freeze-out prediction for Y = 1 WIMPs. Finally, concerning the contribution to the Sommerfeld enhancement, the potential arising from the Yukawa interaction in Eq. 22 is: Such potential is shorter range with respect to typical EW Yukawa potentials, due to the larger Higgs mass as compared to weak boson masses. We find that V H (r) is negligible compared to the size of the EW potential α 2 /r for the mass splittings of our interest.

Appendix C: Millicharged WIMPs
Complex WIMPs with Y = 0 have an unbroken U (1) flavor symmetry which can be gauged by a new dark photon. Generically the dark photon would mix with the visible one through a kinetic mixing operator F F and the complex WIMP would acquire a EM charge . In this scenario the dark gauge symmetry makes the DM accidentally stable as noticed in Ref. [3] at the price of giving up charge quantization. Here we want to summarize the freeze-out predictions and the basic phenomenology of millicharged WIMPs in the limit of very small (i.e. < 10 −10 ) when their phenomenology resemble the one of the real WIMPs discussed in Ref. [1].
Concerning the freeze-out dynamics, the only difference between real and millicharged WIMPs is in the existence of BS with P BS = (−1) L+S+ I−1 2 = −1 formed byχχ pairs of millicharged WIMPs. These are forbidden by the spin-statistic properties of the χχ wave function for real WIMPs. Since P BS is preserved by dipole interactions for Y = 0, to leading  order no transitions can occur between states with opposite P BS and excited BS with P BS = −1 will dominantly decay to 1s and 2s states with the same P BS . The latter have small decay widths with respect to their P BS = 1 counterparts. ns S=1 1 and ns S=1 5 annihilate into four vectors with a rate annihilates into three vectors with rate: where R n0 ∼ (α eff M DM ) 3/2 is the radial wave function of the BS at the origin. As noticed in Ref. [3] for the millicharged 3-plet a large resonance in the Sommerfeld enhancement leads to two different freeze out predictions. We summarize our freeze out predictions in Table II. For completeness, in Fig. 7 we show the SI scattering cross-section of DM on xenon nuclei for the different millicharged candidates. This cross-section is identical to that computed for real candidates with 10 −10 [3,86], while larger would open up new opportunities for direct detection (see e.g. Fig 1 of [3]). From Fig. 7 we see that even in the worst case scenario of very small millicharge large exposure experiments will be able to fully probe millicharged WIMPs with 200 ton-year exposure. The heavier multiplets could be also firmly discovered at DARWIN with kiloton exposure.

Appendix D: Details on the Collider reach
The kinematic constraint of pair producing the DM particles makes direct searches at colliders of energies √ s ≤ 14 TeV impossible for n ≥ 5. We consider masses M √ s/10, including only DY production channels. Indeed in this range VBF production is subdominant. For the n < 5 candidates, the strategies available are the missing invariant mass (MIM), charged tracks (CT) and disappearing tracks (DT) searches. The first is essentially independent on the splitting values: the collider energies of few TeV needed for DM pair production are much larger than the maximum perturbatively generated splitting. The DT searches in general require, at fixed M DM , a scan in the (δm 0 , δm Q M ) to compute the signal, since the spectrum and the decay width of the particle χ Q is determined uniquely by the splittings δm Q M , δm 0 .
To get a feel of how powerful DT searches can be, we fix the neutral splitting to a representative value δm 0 200KeV for the fermions. Our results are given in Fig. 8 where we display the region of the plane WIMP mass versus cτ where experiments can probe the several WIMPs considered in each panel of the figure. For this result we consider only χ + and its conjugate as candidate long-lived for the 2 1/2 . For the 3 1 , the 4 1/2 and 5 1 all states with charge greater than 1 are assumed to decay promptly to the candidate charged long lived states at the bottom the spectrum. 5 As we have fixed δm 0 much smaller than the mass splitting between the neutral and charged states, we can effectively consider decays into both χ 0 and χ DM as if they were degenerate in mass. Considering the possible decay channels χ − into e − ν e χ 0,DM , µ − ν µ χ 0,DM , π − χ 0,DM and the charge conjugates for χ + we get Γ χ ± → χ 0,DM + SM = Γ e ± + Γ µ ± + Γ π ± , (D1)  Table II. The light green shaded region is excluded by the present experimental contraints from XENON-1T [36] and PandaX-4T [5], the green dashed lines shows the expected 95% CL reach of LZ/Xenon-nT [8,9] and DARWIN [10,11].
where the RHS indicates the decay width into neutral states plus the SM states indicated by the subscript. The relevant widths are given by: where g (n, Y, Q) accounts for the different strengths of the W boson coupling to each of the χ Q : with n the dimensionality of the multiplet and Y its hypercharge, G F is the Fermi constant, f π = 131 MeV is the pion decay constant. In the above formula Φ is the full phase space of the 3-body decay of a massive particle of mass M into a massless lepton (e.g neutrino), a massive lepton with mass m l M (e.g. muon), and a heavy particle with a mass M − δ, for splitting δ M : The massless limit was employed for the electron since in the regions relevant for DTs the splittings allow to neglect m e . In principle, for large splittings, also kaon and tau channels should be considered. However, for such large splittings the decays into the other SM particles are so fast that tracks will never be reconstructed.
If we move away from the δm 0 0 points, we have to consider cases in which χ + is heavier than some of the particles with Q > 1. We checked that for the 3 1 , 4 1/2 , 5 1 there are choices of the couplings y 0 and y + for which M χ ++ < M χ + . In this case, χ ++ cannot decay into χ + , and is forced to decay to the neutral states by emitting twice an off-shell W boson that gives rise to l + ν l or π + . This decay rate has been estimated by summing over all possible final states using NDA. With respect to the decay with just one off-shell W emission, we get an extra factor g(n, Y, 2)G F f π δ/(16π 2 ) and g(n, y, 2)G F δ 2 /(16π 2 ) 2 for any extra π or l + v l pair, respectively. Decays into negatively charged states are even more suppressed, since they need more extra final states. These configurations can give χ ++ DT signals as well as CT signals, which are potentially interesting for their peculiar detector response.
Following [18], we employ two different DT signal regions: i) events with at least one disappearing track in addition to a hard photon with energy E γ > 25 GeV; ii) events with two disappearing tracks, one of which is short, and the hard photon. The former has an estimated background cross section of σ B = 18.8 ab at the 10 TeV collider, leading to a signal-to-noise ratio of order 1. The latter is background free. To compute the number of expected signal events, we do a Monte Carlo simulation with MadGraph5 aMC@NLO [87]. We also got similar results employing the analytic approximation explained in the appendix of [1]. The sensitivity to disappearing tracks is shown in the plots of Fig. 8, for minimal δm 0 imposed by DD and BBN, as a function of the dark matter mass and the lifetime cτ for each of the considered WIMP candidates.
Our results in Fig. 8 show that the thermal mass can be excluded over the entire range of lifetimes that give rise to the DT signature. For the 2 1/2 this corresponds to the full parameter space with splitting larger than the EW one. The contours for the 4 1/2 and 5 1 WIMPs deserve some comment. In this case two singly-charged states are present: at least one state is in the DT region, and DT searches for this long-lived state, shown as solid contours, can probe the thermal mass over the full range of possible lifetimes. We also show the DT reach for the state with shorter lifetime, which could produce a second observable signal in parts of the parameter space. We show with different colors the two different physical solutions where the long-lived mass eigenstate χ + corresponds to the gauge eigenstate χ + or χ c − (corresponding to the regions of low or high δm 2+ in Fig. 2), which have slightly different production cross-sections.
For the MIM searches we have proceeded similarly to our earlier work [1] and we have optimized the selection for the fermion 2 1/2 and 3 1 at √ s = 3, 6 TeV and √ s = 6, 10 TeV, respectively. For the other candidates we used the optimal cuts that we derived from our previous results for real candidates with odd n [1]. For each level of systematics we have interpolated our previous results as functions of n and √ s and we have derived optimal cuts for the new n plets studied in this work. We remark that this procedure is potentially inaccurate as the real odd n-plets contain 2 times fewer degrees of freedom, hence have rates smaller by a factor 2. This can in principle affect the result of the optimization of the selection. We checked that the difference with a dedicated optimization is negligible, thus our event selection should be quite close to the optimal one.
Following the mono-W discussion in [1], we have computed the background rate with Mad-Graph5 aMC@NLO [87] using the full 3-body process for η 6, and matched it to a computation using the Weizsäcker-Williams approximation for larger values of pseudorapidity.
We report in Fig. 9 the sensitivity for the 2 1/2 , 3 1 , 4 1/2 , 5 1 , at colliders with suitable √ s that can provide 95% C.L. exclusion. Full results for all benchmark colliders at the nominal luminosity of Eq. 20 are given in Tables III and IV, for the various search channels  under consideration. Results for generic center-ofmass energy and integrated luminosity are displayed in Fig. 11 for the mono-γ and mono-W channels. The results on this figure have been obtained by rescaling the results at √ s = 3, 6, 10, 14, 30 TeV, under the assumption that all cross-sections and kinematical cuts scale trivially with collider energy, and neglecting systematic errors.
The DT results shown in Fig. 9 are a particular subcase of those shown in Fig. 8, as in Fig. 9 we fixed some choice of the mass splitting governing the spectrum specifically for each WIMP as described in the following. For the 2 1/2 the splitting between the charged and neutral components has been set to its gauge contribution δm + = 354 MeV. For the 3 1 the splitting has been set to its gauge value of δm + = 540 MeV. For the 4 1/2 , the splitting δm + has been taken equal to δm − . This choice implies that the χ + and χ − particles are maximal admixtures of the gauge eigenstates, hence they have same cross-sections and lifetimes. The latter are governed by the mass splitting δm + and δm − , which, from Eq. 11 are constrained to satisfy δm + +δm − = δm 0 +2δ g . This equation implies that for δm 0 δ g in the δm + = δm − case we obtain lifetimes that have relatively challenging detector acceptances. Thus we can consider this choice as a representative but somewhat pessimistic estimate of the mass reach. The χ ++ has been assumed to promptly decay equally into χ + and χ −,c . For the 5 1 WIMP we proceeded similarly to the 4 1/2 case.
A summary of the capabilities of the several stages of a high energy muon collider is given in Fig. 12 under the assumption of luminosity following the scaling of Eq. 20. The upshot of these studies is that, as the center of mass energy of the collider is increased, the higher energy machine gains sensitivity to heavier WIMP candidates. All in all, the list of WIMP candidates that we have described in this work and earlier Ref. [1] provides a series of targets that can be probed at successive stages of a future high energy muon collider.

Charged Tracks
Prompt Thermal mass        FIG. 12. Mass reach in the mono-γ, mono-W and DT channels for fixed luminosity as per Eq. 20 at √ s 3 TeV (yellow), 6 TeV (green), 10 TeV (light blue), 14 TeV (red), and 30 TeV (purple). In the mono-W and mono-γ searches we show an error bar, which covers the range of possible exclusion as the systematic uncertainties are varies from 0 to 1%. The colored bars are for an intermediate choice of systematics at 0.1%. Missing bars denoted by an asterisk * correspond to cases where no exclusion can be set in the mass range Mχ > 0.1 √ s. For such cases it is worth considering VBF production modes at the fixed luminosity Eq. 20 or higher luminosity at potentially smaller √ s as illustrated in Fig. 11