The last complex WIMPs standing

Abstract

We continue the study of weakly interacting massive particles (WIMP) started in Bottaro et al. (Eur Phys J C 82:31, 2022), 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.

1 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 $\Omega _{\text {DM}}h^2= 0.11933 \pm 0.00091$ [2].

Experimentally, direct detection (DD) experiments robustly excluded the possibility that DM has tree-level 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\ne 0$), where the elastic scattering between DM and the nucleons is forbidden kinematically. These two possibilities 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\ne 0$.¹

Table 1 Thermal masses of complex WIMPs with $Y\ne 0$, obtained including Sommerfeld enhancement and BSF. The upper bound on n for even multiplets comes from the perturbative unitarity bound, as can be seen from the $(\sigma v)^{J=0}_{\text {tot}}/(\sigma v)^{J=0}_{\text {max}}$, where $(\sigma v)^{J=0}_{\text {max}}$ is the maximal allowed annihilation cross section [4]. The loss of perturbativity is also signaled by the Landau pole $\Lambda _\text {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 $\delta m_0$ . For multiplets with $n>5$ the largest UV cutoff $\Lambda _{\text {UV}}^{\text {max}}$ required to generate the minimal viable splitting is smaller than $10M_\text {DM}$. For each candidate we provide the allowed range for the mass splittings. The lower limit on $\delta m_0$ 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 $\mathcal {O}_{0}$ in Eq. (1). Similarly, the lower limit on $\delta m_{ {Q_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 $\mathcal {O}_{+}$ in Eq. (1)

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 1, 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 the-next-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 unspecified UV dynamics).

The inelastic splitting $\delta 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 requirement 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$. Multiplets $Y=1$ and $n>5$ or with higher hypercharges are excluded. At the same time, the inelastic splitting is bounded from above by DD constraints on Higgs-mediated nuclear recoils. This leaves a finite window for the inelastic splitting of every multiplet which we report in Table 1. This window will be further probed by large exposure DD experiments such as LZ [8], XENONnT [9], and ultimately by DARWIN/G3 [10, 11].

The natural value of the charged splitting $\delta 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 $\sqrt{s}=3\text { TeV}$ would only be able to exclude the fermionic EW doublet with a large luminosity from around $\mathcal {L}=5\text { ab}^{-1}$ to $\mathcal {L}=50 \text { 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 $\sqrt{s}=6\text { TeV}$ collider, which can exclude the $2_{1/2}$ with 4 $\text { ab}^{-1}$ and make a discovery with a more demanding 20 $\text {ab}^{-1}$.²

The DD cross section of the $3_1$ can only be probed with large exposure xenon experiment like DARWIN, while a muon collider of $\sqrt{s}=10 \text { TeV}$ would be able to make a discovery at the thermal mass for the fermionic $3_1$ using a luminosity around 10 $\text {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 Sect. 2 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 1 while the more technical aspects of the freeze-out computation are given in Appendix B. In Sect. 3 we show how large exposure DD experiments can probe most of the complex WIMPs with the minimal inelastic and charged splittings. In Sect. 4 we study the full parameter space of the complex WIMPs, in the charged vs neutral splitting plane. We take as examples the doublet and 4-plet with $Y=1/2$ and the triplet and 5-plet with $Y=1$. In Sect. 5 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.

2 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 $Q=T_3+Y$, where $T_3=\text {diag}\left( \frac{n+1}{2}-i\right)$ with $i=1,\dots , 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=\pm \left( \frac{n+1}{2}-i\right)$ for $i=1,\dots , 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\ne 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\ne 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\ne 0$ is:

(1)

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 $\gamma _5$ inside the DM bilinear. We now illustrate the physical consequence of $\mathcal {O}_0$ and $\mathcal {O}_+$ in turn in Sects. 2.1 and  2.2 and derive the implication for complex WIMPs in Sect. 2.3. We comment on DM stability in Sect. 2.4.

Fig. 1

Summary of the lower bounds on the neutral mass splitting. The dark green shaded region is excluded by tree-level Z-exchange in XENON1T [27] for both scalar and fermionic DM. The dashed green line shows what XENON1T could probe by analyzing high recoil energy data. The blue and red lines are the BBN bounds on the splitting for fermionic and scalar DM respectively

2.1 Inelastic splitting

The non-renormalizable operator $\mathcal {O}_{0}$ is required to remove the sizeable coupling to the Z boson of the neutral component $\chi _N$ of the EW multiplet

(2)

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, $\mathcal {O}_{0}$ induces a mixing between $\chi _N$ and $\chi _N^c$. Replacing the Higgs with its VEV, $(H^{c\dagger })\frac{\sigma ^a}{2} H$ is non-zero only if we pick $\sigma ^a=\sigma ^+$, so that the new (pseudo Dirac) mass term in the Lagrangian reads

$$\begin{aligned} \begin{aligned}&\mathscr {L}_m=M_\chi \overline{\chi }_N\chi _N+\frac{\delta m_0}{4}\left[ \overline{\chi }_N\chi _N^c +\overline{\chi ^c}_N\chi _N\right] ,\\&\delta m_0 = 4y_0 c_{nY0}\Lambda _{\text {UV}}\left( \frac{v}{\sqrt{2}\Lambda _{\text {UV}}}\right) ^{4Y}\ . \end{aligned} \end{aligned}$$

(3)

$c_{nYQ}=\frac{1}{2^{Y+1}(4Y)!}\prod _{j=-Y-|Q|}^{Y-1-|Q|}\sqrt{\frac{1}{2}\left( \frac{n+1}{2}+j\right) \left( \frac{n-1}{2}-j\right) }$ contains the normalization of $\mathcal {O}_{0}$ and the matrix elements of the generators. The mass eigenstates are Majorana fermions, $\chi _0$ and $\chi _\text {DM}$, with masses $M_0=M_\chi +\delta m_0/2$ and $M_\text {DM}=M_\chi -\delta m_0/2$, whose coupling to the Z boson is

(4)

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

$$\begin{aligned} \frac{1}{2}\mu v_\text {rel}^2<\delta m_0\ ,\quad \mu =\frac{M_\text {DM}m_N}{M_\text {DM}+m_N}\ , \end{aligned}$$

(5)

where $m_N$ is the mass of the nucleus, $\mu$ is the reduced mass and $v_\text {rel}$ is DM-nucleus relative velocity. In particular, given the upper bound on the relative velocity $v_\text {rel}<v_\text {E}+v_\text {esc}$, where $v_\text {E}=240$ km/s is the Earth’s velocity and $v_\text {esc}=600$ km/s is the assumed escape velocity of DM in the Milky Way, the largest testable mass splitting is $\delta m_0^{\max }=1/2\mu (v_\text {E}+v_\text {esc})^2$ which for xenon nuclei gives $\delta m_0^{\max }\simeq 450$ keV. The splitting for a given recoil energy is

$$\begin{aligned} \delta m_0(E_R)=\sqrt{2m_N E_R}(v_\text {E}+v_\text {esc})-E_R\frac{m_N}{\mu }\ , \end{aligned}$$

(6)

which explain why the maximal constrained splitting experimentally is $\delta m_0^{\max ,\exp }\simeq 240 \text { keV}$ as shown in Fig. 1, given that XENON1T [27] analyzed data only for $E_R<40\text { keV}$. Extending the range of XENON1T to higher recoil energies would be enough to probe splitting up to $\delta m_0^{\max }$ as already noticed in Refs. [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 $\delta m_0$ comes from requiring that the decay $\chi _0\rightarrow \chi _\text {DM}+$SM happens well before BBN. The leading decay channels are $\chi _0\rightarrow \chi _\text {DM}\gamma$, $\chi _0\rightarrow \chi _\text {DM}\bar{\nu }\nu$ and $\chi _0\rightarrow \chi _\text {DM}\bar{e}e$ with decay widths

$$\begin{aligned}&\Gamma _\gamma =\left( 1+\frac{1}{2}\log \left( \frac{m_W^2}{M_\text {DM}^2}\right) \right) ^2\frac{Y^2\alpha _2^2\alpha _{\text {em}}}{\pi ^2}\frac{\delta m_0^3}{M_\text {DM}^2}\ , \end{aligned}$$

(7)

$$\begin{aligned}&\Gamma _{\bar{\nu }\nu }\simeq 6 \Gamma _{\bar{e}e}=\frac{G_F^2\delta m_0^5Y^2}{5\pi ^3}\ . \end{aligned}$$

(8)

Fig. 2

Mass splittings of $3_1$ (left) and $4_{1/2}$ (right) as a function of $\delta m_{2+}$. In the $3_1$ case, no mixing between the charged components of the multiplet can occur and $\delta m_+$ is a monotonic function of $\delta m_{2+}$. In the $4_{1/2}$ case, instead, because of the mixing induced by $\mathcal {O}_{0}$, the positively charged mass eigenstate $\chi ^+$ is always heavier than $\chi ^-$. The splitting between $\chi ^+$ and $\chi ^-$ has a minimum of order $\delta m_0$, which was taken 50 MeV in this plot for display purposes

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 $\chi _0$:

$$\begin{aligned} \Gamma _{\chi _0}\equiv \Gamma _{\bar{\nu }\nu }+\Gamma _{\bar{e}e}+\Gamma _{\gamma }>\tau _\text {BBN}^{-1}\ , \end{aligned}$$

(9)

where $\tau _\text {BBN}^{-1}=6.58 \times 10^{-25}\text { 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 $\chi _0\rightarrow \chi _\text {DM}\gamma$. As a consequence the BBN bounds are stronger for scalar WIMPs.

2.2 Charged splitting

The operator $\mathcal {O}_{+}$ in Eq. (1) is necessary to make the DM the lightest state in the EW multiplet for all the n-plets whose 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\ll M_\chi$ are [12,13,14]

$$\begin{aligned} \Delta M_Q^{\text {EW}}= \delta _g \left( Q^2+\frac{2YQ}{\cos \theta _W}\right) \ , \end{aligned}$$

(10)

where $\delta _g=(167\pm 4)\text { 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 $\mathcal {O}_{+}$, the final splittings $\delta m_{Q}= M_Q-M_\text {DM}$ between the DM and the charged components read

$$\begin{aligned} \delta m_{Q}= & {} \frac{\delta m_{0}}{2}+\delta _g Q^2+\text {sgn}(Q)\nonumber \\{} & {} \times \sqrt{\left( \frac{2Y\delta _g}{\cos \theta _W}-\frac{y_+v^2}{4 \Lambda _{\text {UV}}}\right) ^2Q^2+\frac{\delta m_0^2}{4}\frac{c_{nYQ}^2}{c_{nY0}^2}}\ , \end{aligned}$$

(11)

where $\text {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 $\chi _Q$ and $\chi _{-Q}^c$ induced by $\mathcal {O}_{0}$. This obviously 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 $\delta m_0$ and $\delta m_{ {Q_M} }$, where ${Q_M} \equiv Y+(n-1)/2$ is the largest electric charge in the multiplet. Since $c_{nY {Q_M} }=0$, $\delta 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 $\delta m_{2+}$. In the former case, no mixing occur within the components of the multiplet and $\delta m_+$ is a monotonic function of $\delta m_{2+}$. For the $4_{1/2}$, instead, the mixing induced by $\mathcal {O}_{0}$ between the components with $Q=\pm 1$ makes the positively charged mass eigenstate $\chi ^+$ heavier than $\chi ^-$.

In Sect. 4 we will explore the parameter space spanned by $\delta m_0$ and $\delta m_{ {Q_M} }$, fixing the thermal DM mass of every EW multiplet as shown in Table 1. 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 elastic cross section onto nucleons set by PandaX-4T [5] translates into upper bounds on the neutral and charged splittings as reported in Table 1.

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

$$\begin{aligned} \text {F.T.}\equiv \max _I\left[ \frac{\textrm{d}\log \delta m_Q}{\textrm{d}\log \delta m_I}\right] \ , \end{aligned}$$

(12)

where the index I runs over the three contributions in the definition of $\delta m_Q$ in Eq. (11). Large values of F.T. imply a significant amount of cancellation between two or more parameters.

2.3 Viable complex WIMPs

The EFT approach used to write Eq. (1) is meaningful only if the UV physics generating $\mathcal {O}_{+}$ and $\mathcal {O}_{0}$ is sufficiently decoupled from DM. When $\Lambda _{\text {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\simeq M_{\text {DM}}$ and are likely to modify our freeze-out predictions. To avoid these difficulties we restrict ourselves to $\Lambda _{\text {UV}}\ge 10M_\text {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 $\delta m_0>\delta m_0^{\min }$, we derive the viable window for $\Lambda _{\text {UV}}$:

$$\begin{aligned} 10M_\text {DM}<\Lambda _{\text {UV}}\le \left( \frac{4y_0c_{nY0} v^{4Y}}{2^{2Y}\delta m_0^{\min }}\right) ^{\frac{1}{4Y-1}}. \end{aligned}$$

(13)

We are now interested in estimating for which multiplets the viable window shrinks to zero. Setting $y_0=(4\pi )^{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\le 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 $\sim 10^7\text { TeV}$ as shown in the fifth column of Table 1. 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 $\Lambda _{\text {UV}}$.

For $Y=1$ multiplets the maximal required cutoff is of order $\sim 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_\text {DM}\sim 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.

Fig. 3

Expected SI cross-sections for different complex WIMPs for minimal splitting as defined in Sect. 3. The blue dots correspond to Dirac WIMPs and the red dots to complex scalar WIMPs. The vertical error bands correspond to the propagation of LQCD uncertainties on the elastic cross-section (Eq. 18), while the horizontal error band comes from the uncertainty in the theory determination of the WIMP freeze out mass in Table 1. The light green shaded region is excluded by the present experimental constraints from XENON1T [36] and PandaX-4T [5], the green dashed lines shows the expected 95% CL reach of LZ/XENONnT [8, 9] and DARWIN [10, 11]

2.4 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] ($\tau _{\text {DM}} > rsim 10^{19} \text { s}$) but also to satisfy the stronger astrophysical bounds on decaying DM [39,40,41] ($\tau _{\text {DM}} > rsim 10^{28}\text { s}$).

For $2_{1/2}$ and $3_1$ we can write renormalizable interactions $\overline{\chi ^c}He_R$ and $\overline{\chi }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 $\mathbb {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

$$\begin{aligned} \begin{aligned} \mathscr {L}_{\text {even}}&= \frac{C_1}{\Lambda _{\text {UV}}^{n-3}}\overline{L^c}\chi (H^\dagger H)^{\frac{n-2}{2}}\\&\quad +\frac{C_2}{\Lambda _{\text {UV}}^{n-3}}\overline{L^c}\chi W_{\mu \nu }W^{\mu \nu }(H^\dagger H)^{\frac{n-6}{2}}\\&\quad +\frac{C_3}{\Lambda _{\text {UV}}^{n-3}}\overline{L^c}\chi (W_{\mu \nu }W^{\mu \nu })^{\frac{n-2}{4}}\\&\quad +\frac{C_4}{\Lambda _{\text {UV}}^{n-3}}\overline{L^c}\sigma ^{\mu \nu }\chi W_{\mu \nu }(H^\dagger H)^{\frac{n-4}{2}} \\&\quad + \frac{C_5}{\Lambda _{\text {UV}}^2} \overline{\chi } \chi ^c \overline{L^c}\chi \ , \end{aligned}\nonumber \\ \end{aligned}$$

(14)

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\ge 6$ this is given by the operator with coefficient $C_5$ by closing the $\chi$ loop and attaching $(n-2)/2$ $W^{\mu \nu }$ operators, in a similar fashion to real candidates [1]. In all cases, we need at least $\Lambda _{\text {UV}}>10^{10}M_\text {DM}$ for $\mathcal {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 $\delta 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

$$\begin{aligned} \mathscr {L}_5 = \frac{C_1}{\Lambda _{\text {UV}}^2}(\overline{\chi }L^c\! H)(H^\dagger \! H)+\frac{C_2}{\Lambda _{\text {UV}}^2}(\overline{\chi }\sigma ^{\mu \nu }L^c\! H)W_{\mu \nu }. \end{aligned}$$

(15)

These require $\Lambda _{\text {UV}}>10^{10}M_\text {DM}$ to ensure DM accidental stability which is again incompatible with the upper bound of $\Lambda _{\text {UV}}<20 M_\text {DM}$ needed to generate the inelastic splitting.

As a result, none of the complex WIMPs with $Y\ne 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.

3 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 $\delta m_0= 220\text { keV}$ for $n_{1/2}$ WIMPs with $n\le 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 $\delta m_0$. For scalars instead, we set $\delta m_0=4.9$ MeV for $n_{1/2}$ and $\delta 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\text { MeV}$ for $2_{1/2}$ and $542\text { 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 Sect. 3.1 and at a future muon collider in Sect. 3.2.

3.1 Direct detection prospects

The spin independent scattering cross-section $\sigma _\text {SI}$ of DM on nuclei receives two contributions: (i) from purely EW loop diagrams (ii) from Higgs mediated tree-level diagrams generated by both $\mathcal {O}_{0}$ and $\mathcal {O}_{+}$. For minimal splitting Higgs mediated scattering is subdominant and $\sigma _{\text {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

$$\begin{aligned} \mathscr {L}_{\text {eff}}^{\mathrm{{SI}}} = f_q m_q \bar{\chi } \chi \bar{q} q + \frac{g_q}{M_\text {DM}} \bar{\chi } i \partial ^{\mu } \gamma ^{\nu } \chi \mathcal {O}^q_{\mu \nu } + f_G \bar{\chi } \chi G_{\mu \nu }G^{\mu \nu },\nonumber \\ \end{aligned}$$

(16)

where is the quark twist-2 operator. The Wilson coefficients are given by [17]

$$\begin{aligned} \begin{aligned} f_q^{\mathrm{{EW}}}&\simeq -\frac{\pi \alpha _2^2}{16 m_h^2 m_W} \left[ n^2-1-\left( 1.03+22 (a_q^{V,2}-a_q^{A,2})\right) Y^2 \right] ,\\ g_q^{\mathrm{{EW}}}&\simeq \frac{\pi \alpha _2^2}{24 m_W^3}\left[ n^2-1-\left( 4-18.2 (a_q^{V,2}+a_q^{A,2})\right) Y^2 \right] ,\\ f_G^{\mathrm{{EW}}}&\simeq \frac{\alpha _s\alpha _2^2}{192 m_h^2 m_W} \left[ \left( \sum _{q=c,b,t} \kappa _q +2.6\right) (n^2-1)\right. \\&\quad \left. -\left( 1.03\sum _{q=c,b,t} \kappa _q -7.5\right) Y^2 \right] \, , \end{aligned} \end{aligned}$$

(17)

where $m_h$ is the mass of the Higgs and $\kappa _c = 1.32,\,\kappa _b=1.19,\,\kappa _t=1$. Furthermore we have defined $a_q^V = T_{3q}/2 - Q_q s_w^2$, $a_q^A=-T_{3q}/2$ 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 $\sigma _\text {SI}$ per nucleon (for $M_\text {DM}\gg m_N$) as

$$\begin{aligned} \sigma _{\text {SI}} \simeq \frac{4}{\pi } m_N^4 \vert k_N^{\mathrm{{EW}}} \vert ^2, \end{aligned}$$

(18)

where $m_N$ is the mass of the nucleon and

$$\begin{aligned} k_N^{\mathrm{{EW}}} =\!\!\!\! \sum _{q=u,d,s} f_q^{\mathrm{{EW}}} f_{Tq} + \frac{3}{4} (q(2) + \bar{q}(2)) g_q^{\mathrm{{EW}}} - \frac{8 \pi }{9 \alpha _s} f_{TG} f_G^{\mathrm{{EW}}}\, , \end{aligned}$$

where the nucleon form factors are defined as $f_{Tq} = \langle N| m_q \bar{q} q |N\rangle / m_N$, $f_{TG} = 1-\sum _{q=u,d,s} f_{Tq}$, and $\langle N(p)| \mathcal {O}^q_{\mu \nu }|N(p)\rangle = (p_{\mu } p_{\nu }-\frac{1}{4}m_N^2 g_{\mu \nu })(q(2) + \bar{q}(2))/m_N$, and q(2), $\bar{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 Collaboration [43] in the case of $N_f=2+1+1$ dynamical quarks [44, 45].³

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 the $5_1$ that we discuss below.

Combining Eq. (17) with Eq. (18) the parametric expression for $\sigma _{\text {SI}}$ is

$$\begin{aligned} \sigma _{\text {SI}} \approx 10^{-49} \text { cm}^2 (n^2-1-\xi ~ Y^2)^2, \end{aligned}$$

(19)

where $\xi =16.6\pm 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\simeq \sqrt{(n^2-1)/\xi }$. In particular, for $n=2$ the exact cancellation takes place at $Y\simeq 0.44\pm 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\simeq 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.

3.2 $2_{1/2}$ and $3_{1}$ at muon colliders

As opposed 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]

$$\begin{aligned} \mathcal {L}=10 \text { ab}^{-1}\cdot \left( \frac{\sqrt{s}}{10\text { TeV}}\right) ^2\,, \end{aligned}$$

(20)

and we concentrate on machines running at $\sqrt{s}=3,6,10\text { 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 $\sqrt{s}\le 10~\textrm{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 $\mathcal {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-$\gamma$ 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].

Fig. 4

Top two rows: Collider reach for the $2_{1/2}$ of mass $M_\chi =1.1$ TeV and $\sqrt{s}=3$ TeV (left) and $\sqrt{s}=6$ TeV (right). Upper row: Combined reach for MIM (mono-$\gamma$+mono-W) as a function of the luminosity $\mathcal {L}$ for different value of the expected systematic uncertainties parametrized by $\epsilon _{\text {sys}}$: $1\%$ (purple), $0.1\%$ (cyan), 0 (green). The red lines show the benchmark luminosity following Eq. (20). Lower row: Reach of DT for varying lifetime of the charged track. The red vertical line shows the benchmark lifetime for EW splitting. The blue and orange curves correspond to the reach with 1 and 2 disappearing tracks. For the single track search the solid (dashed) line assumes 0 (10%) systematics on the background. Fine-Tuning values computed on $\delta m_+$ following Eq. (12) are displayed as color code from green (no fine-tuning) to red (higher fine-tuning). Bottom two rows: same as for the top rows, but for the case of the $3_{1}$ of mass $M_\chi =2.85$ TeV and $\sqrt{s}=6$ TeV (left) and $\sqrt{s}=10$ TeV (right)

We estimate the reach of MIM searches from cut-and-count experiments using the following estimate for the statistical significance

$$\begin{aligned} \text {significance}=\frac{S}{\sqrt{B+S+\epsilon _{\text {sys}}^2 (B^2+S^2)}}\ , \end{aligned}$$

(21)

where $\epsilon _{\text {sys}}$ parametrizes the expected systematic uncertainties on the measurement which we take to be zero, $0.1\%$ or $1\%$.

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]). For the DT we compute the reach using Eq. (21) with different assumption on the background systematics. For the single DT we rescale the background from the study of Ref. [18] and show the reach taking the systematic uncertainties up to $10\%$. The double DT is assumed to be background free and it is hence insensitive to systematics on the background. The $2\sigma$ reach correspond to four signal events in this channel. For the charged tracks we will limit ourselves to give the signal yield, given the lack of dedicated background studies at the moment.

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\tau _+\approx 6.6\text { mm}$ for the $2_{1/2}$ and $c\tau _+\approx 0.75\text { 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_\chi$ of the candidate.

As can be seen from Fig. 4, the MIM search at $\sqrt{s}=3$ TeV with the benchmark luminosity $\mathcal {L}=1\; \textrm{ab}^{-1}$ is not sensitive to the $2_{{1/2}}$. A $\sqrt{s}=6$ TeV collider with benchmark luminosity $\mathcal {L}=4\; \textrm{ab}^{-1}$, instead, is able to probe the doublet WIMP at $2\sigma$ CL. In general the mono-W and mono-$\gamma$ channels give comparable mass reach for the benchmark luminosity adopted here, full detail is given in Table 3. 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 $\sqrt{\mathcal {L}}$ because the selections, hence the result, have been optimized as a function of $\mathcal {L}$ when dealing with $\epsilon _{\text {sys}}\ne 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 $\delta m_+$. We see that a collider of $\sqrt{s}=6\text { TeV}$ can probe the $2_{{1/2}}$ for charge-neutral splitting generated purely by EW interactions. A collider of $\sqrt{s}=3\text { TeV}$ can probe a large portion of the allowed lifetimes for the charged tracks corresponding to non-zero UV contributions in Eq. (11).

Fig. 5

Signal yield of the charged tracks for the $2_{1/2}$, $3_1$, $4_{1/2}$, $5_1$. The plot holds for all $\sqrt{s}$ in the tens of TeV range, assuming the luminosity scales as Eq. (20). The cuts applied on the produced charged $\chi$’s are $p_T>200$ GeV, $|\eta _\chi |<2$

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-$\gamma$ 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 charged-neutral splitting.

Finally, it is important to notice that with an $\mathcal {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.

4 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 Sect. 3. We dedicate this section to explore the phenomenological consequences of releasing the minimal splitting assumption.

Fig. 6

Direct Detection and collider signatures in the plane $\delta m_0$ vs. $\delta m_{ {Q_M} }$ from the $2_{1/2}$ (upper left), $3_1$ (upper right) $4_{1/2}$ (lower left) and $5_1$ (lower right). The green shaded regions are excluded by current DD constraints on inelastic DM at small $\delta m_0$ and on Higgs-exchange elastic scattering at large splittings. Gray hatched regions are excluded by BBN constraints on the longest lived unstable particle in the multiplet. Red dashed lines delimit the range of perturbative mass splitting at fixed $\Lambda _{\text {UV}}/M_\text {DM}$ ratio. Light blue shaded patches are not viable because the lightest WIMP in the n-plet is not the neutral one. Dashed green lines show prospects from future high exposure xenon experiments: LZ, DARWIN and DARWIN/G3 (arrows pointing to the direction of the expected probed region). In the gray shaded region the DD signal falls below the neutrino floor of xenon experiments [11, 16]. The vertical XENON1T-high energy line show the ultimate reach of xenon experiments on inelastic DM (see Fig. 1). Blue dot-dashed lines correspond to different expected lengths of charged tracks. Accordingly, different hues of brown distinguish regions where different signatures at future colliders are expected, as defined in Sect. 4.2. In every panel the DM thermal mass gives an idea of the center-of-mass energy of an hypothetical collider which would be needed to produce a DM pair. The precise center-of-mass energy for exclusion/discovery requires a more detailed study (see Sect. 3.2 and App. D) The dashed gray line shows the EW value of $\delta m_{ {Q_M} }$. Purple dashed lines show the contours of the fine-tuning among the different mass splittings as defined in Eq. (12). Mass splittings above the F.T.=1 line are not fine tuned

4.1 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

$$\begin{aligned} \mathscr {L}_{D,h}=-\frac{\lambda _Dv}{2\Lambda _{\text {UV}}}\chi _\text {DM}^2h\, . \end{aligned}$$

(22)

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 $\delta 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 non-renormalizable couplings is compact.

Following Ref. [49] we can integrate out the Higgs boson and write the couplings of the DM to the SM

$$\begin{aligned} \mathscr {L}_{\text {eff,h}}^{\text {SI}}=\frac{\lambda _D}{2m_h^2\Lambda _{\text {UV}}}\chi ^2\left( m_q \bar{q}q-\frac{\alpha _s}{4\pi }G^{a}_{\mu \nu }G^{a\mu \nu } \right) , \end{aligned}$$

(23)

so that the matrix elements $f_q$ and $f_G$ required to compute $\sigma _{\text {SI}}$ are simply given by

$$\begin{aligned} f_q = f_q^{\text {EW}}+\frac{\lambda _D}{2m_h^2\Lambda _{\text {UV}}},\, f_G = f_G^{\text {EW}}-\frac{\alpha _s}{8\pi }\frac{\lambda _D}{2m_h^2\Lambda _{\text {UV}}}, \end{aligned}$$

where $f_q^{\text {EW}}$ and $f_G^{\text {EW}}$ come from EW loops and are given in Eq. (17). We can rewrite the coupling $\lambda _D$ solely in terms of mass splitting as

$$\begin{aligned} \frac{\lambda _D}{\Lambda _{\text {UV}}}=-\frac{2Y}{v^2}(\delta \mu _{ {Q_M} }+n\delta \mu _0)\, , \end{aligned}$$

(24)

where

$$\begin{aligned} \!\!\!\! \delta \mu _{ {Q_M} }\!\equiv \! \frac{2\delta m_{ {Q_M} }-\delta m_0-2\Delta M_{ {Q_M} }^{\text {EW}}}{2 {Q_M} }\, ,\ \delta \mu _0\equiv \frac{\delta m_0}{n} \end{aligned}$$

(25)

and $\Delta M_{ {Q_M} }^{\text {EW}}$ is the gauge induced mass splitting in Eq. (10).

Replacing Eq. (24) into Eq. (18) allows us to translate the upper bound on $\sigma _{\text {SI}}$ into an upper bound on $\delta m_0$ and $\delta m_{ {Q_M} }$.

4.2 Collider searches

The splittings $\delta m_0, \delta 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 parameter space $(\delta m_0, \delta m_{ {Q_M} })$ in various regions according to the lifetime $\tau _{\textrm{LCP}}$ of the Longest Lived Charged Particle (LCP) of a given multiplet.

For $c\tau _{\textrm{LCP}}>1\;\textrm{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 $\chi$ pair production at a muon collider, as shown in Fig. 5. We highlight that this results approximately hold for generic $\sqrt{s}$ in the domain of tens of TeV, as long as Eq. (20) for the luminosity holds. In Fig. 5 we show the $\beta \gamma$ 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, $|\eta _\chi |<2$, inspired by LHC searches [48]. For the $2_{1/2}$ and $3_1$ we include in our counting all charged particles production, using the fact that $\chi ^{2+}$ always promptly decays to the long lived $\chi ^+$. 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 $\chi ^{-}$ above $\chi _{\text {DM}}$. In this spectrum configuration $\chi ^{-}$ is the most natural and only candidate to make long charged tracks, therefore we only consider this contribution in our result.⁴

For $1\; \textrm{m}>c\tau _{\textrm{LCP}}>0.33\;\textrm{cm}$ the lightest charged WIMP gives a disappearing track signal. This search has been introduced in Sect. 3 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\tau _{\textrm{LCP}}<0.33\;\textrm{cm}$ the lightest charged WIMP decays promptly on collider scales and gives missing energy signatures such as the mono-W and mono-$\gamma$ discussed in Sect. 3. The choice for the threshold value to consider the signal as prompt corresponds to $c \tau /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.

4.3 Parameter space for complex WIMPs

In the rest of this section, we examine in greater detail the parameter space spanned by $\delta m_0$ and $\delta m_{ {Q_M} }$ specifically looking at each of the lightest multiplets of Table 2 up to the $5_1$ WIMP. This will allow us to discuss the salient phenomenological features of the complex WIMPs.

4.3.1 The Dirac $2_{1/2}$ and $3_1$

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 $\delta m_{0,+}>0$ automatically implies the neutral WIMP candidate is the lightest one. Perturbativity requires $\delta m_0<40\text { 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 $\delta m_0$ for the $3_{1}$ with respect to $2_{1/2}$ was expected from the higher dimensionality of $\mathcal {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 $\delta m_0$ and $\delta 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

$$\begin{aligned} \delta \mu \equiv 2\delta m_{ {Q_M} }+(2 {Q_M} -1)\delta m_0\ . \end{aligned}$$

(26)

In particular we find

$$\begin{aligned} \sigma _\text {SI}\approx 10^{-48}\text {cm}^2\left\{ \begin{aligned}&\left( 0.3-\frac{\delta \mu }{1 \text { GeV}}\right) ^2,\ [2_{1/2}]\\&\left( 0.2+\frac{\delta \mu }{1 \text { GeV}}\right) ^2, \ [3_{1}] \end{aligned}\right. \, . \end{aligned}$$

(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 $\delta \mu \simeq 300\text { MeV}$, while for $3_1$ the minimum $\sigma _\text {SI}$ is obtained for $\delta \mu =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 $\delta \mu >1.6$ GeV ($\delta \mu > 2.0$ GeV). PANDAX-4T already excludes mass splitting $\delta \mu <20$ GeV ($\delta \mu <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 $\delta m_0$ and $\delta m_+$. The region of large splittings down to $\delta \mu \sim 1\text { 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 $\delta 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 Sect. 3.2. The general message is that for $\sqrt{s}$ sufficiently larger than twice the thermal DM mass DT searches will be powerful enough to probe the parameter space.

4.3.2 The dirac $4_{1/2}$ and $5_1$

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 $\delta 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\text { GeV}<\delta \mu < 15\text { GeV}$. This region can be probed by direct searches at a future high energy muon collider with $\sqrt{s}>10 \text { 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 $\delta m_0<3\text { 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\text { GeV}\lesssim \delta m_{3+}\lesssim 2\text { GeV}$. The parameter space can be fully excluded by charged and stub track searches if a muon collider of $\sqrt{s}>20$ TeV will be constructed. Besides, all the allowed mass splittings except a small window $1.8\text { GeV}\lesssim \delta m_{3+}\lesssim 2.0\text { GeV}$ can be probed by DARWIN/G3.

5 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 (Sect. 5.1) and on the electron dipole moment (Sect. 5.2).

5.1 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 polarizations

$$\begin{aligned} \begin{aligned}&\hat{Y}= \frac{g_Y^2m_W^2}{2}\Pi ''_{BB}(0), \quad \hat{W}= \frac{g_2^2 m_W^2}{2}\Pi ''_{33}(0)\ ,\\&\hat{S}= \Pi _{3B}'(0)\quad ,\quad \hat{T}= \frac{\Pi _{33}(0)-\Pi _{WW}(0)}{m_W^2}\ , \end{aligned} \end{aligned}$$

(28)

where $\Pi _{ij}$ appear in the kinetic terms of the EFT describing the SM vectors interactions at energies smaller than the DM mass

$$\begin{aligned} \mathscr {L}_\text {oblique}=-\frac{1}{2}V_i^\mu \Pi _{ij}V_{j,\mu }\ . \end{aligned}$$

(29)

In particular, all the EW multiplets (both real and complex) give a universal contribution to the $\hat{W}$, $\hat{Y}$ as previously found in [22]

$$\begin{aligned}&\hat{W}=\frac{\alpha _\text {em} \cot ^2\theta _W}{180\pi }\frac{m_Z^2}{M_\text {DM}^2}\kappa n(n^2-1)\nonumber \\&\quad \simeq 3.8\kappa \times 10^{-7}\left( \frac{1\text { TeV}}{M_\text {DM}}\right) ^2\!\!n(n^2-1)\ , \end{aligned}$$

(30)

$$\begin{aligned}&\hat{Y}=\frac{\alpha _\text {em}}{15\pi }\frac{m_Z^2}{M_\text {DM}^2}\kappa nY^2\!\simeq 3.4\kappa \times 10^{-7} \left( \frac{Y}{1/2}\right) ^2\nonumber \\&\quad \quad \times \left( \frac{1\text { TeV}}{M_\text {DM}}\right) ^2n \ , \end{aligned}$$

(31)

where $\kappa = 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_{\text {DM}}\sim n^{5/2}$. The contribution to $\hat{W}$ scale like as $\sim 1/n^2$ while the one to $\hat{Y}$ as $\sim 1/n^4$.

$\hat{W}$ and $\hat{Y}$ 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 $\hat{W}$ has been recently bounded by CMS [54] and $\hat{Y}$ 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]. Given the size of the expected WIMP contributions in Eqs. (30) and (31) the results of Ref. [57] indicates that a muon collider of $\sqrt{s}=30\text { 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 interactions give an additional contributions to the $\hat{S}$ parameter. Moreover, the mass splittings inside the EW multiplet induced by $\mathcal {O}_{0}$ and $\mathcal {O}_{+}$ break the custodial symmetry in the Higgs sector. As a consequence, complex WIMPs give irreducible contributions also to the $\hat{T}$ parameter. The explicit formulas for fermionic WIMPs are

$$\begin{aligned} \hat{S}_F&=-\frac{\alpha _{\text {em}}Yn(n^2-1)}{9\pi s_{2W}M_\text {DM}}\delta \mu _{ {Q_M} }\simeq -1.6 \times 10^{-6}\left( \frac{Y}{1/2}\right) \nonumber \\&\quad \times \left( \frac{1 \text { TeV}}{M_\text {DM}}\right) \left( \frac{\delta \mu _{ {Q_M} }}{10\text { GeV}}\right) n(n^2-1)\, , \end{aligned}$$

(32)

$$\begin{aligned} \hat{T}_F&=-\frac{\alpha _{\text {em}}}{ 18\pi s_{2W}^2}n(n^2-1)\frac{\delta \mu _{Q_M}^2+\delta \mu _{0}^2}{m_Z^2}\nonumber \\&\simeq -2.3 \times 10^{-6}\frac{\delta \mu _{ {Q_M} }^2+\delta \mu _0^2}{100 \text { GeV}^2}n(n^2-1)\, , \end{aligned}$$

(33)

where $s^2_{2W}=\sin ^22\theta _W\simeq 0.83$ while $\delta \mu _{ {Q_M} }$ and $\delta \mu _0$ are defined in Eq. (25). The corresponding formulas for scalar WIMPs are easily obtained from the fermionic ones by replacing $\hat{S}_S=\frac{1}{4}\hat{S}_F$ and $\hat{T}_S=-\hat{T}_F$. Our result for $\hat{T}_S$ is consistent with the result in Ref. [58], first computed in [59].

Both $\hat{S}$ and $\hat{T}$ are currently constrained at the level of $3\times 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 $\hat{S}$ 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 on $\hat{T}$ are enhanced to the point that a $12_{1/2}$ multiplet is giving $\hat{T}\simeq 3\times 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).

5.2 Electric dipole moment

As mentioned in Sect. 2 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

$$\begin{aligned} \mathscr {L}_{\text {EDM}}\supset \frac{c_{WW}}{\Lambda ^2_{\text {UV}}} \vert H\vert ^2 W\tilde{W}\ , \end{aligned}$$

(34)

where we neglect for this discussion operators such as $\vert H\vert ^2 B\tilde{B}$ and $H^{\dagger }\sigma ^a H W\tilde{B}$ 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

$$\begin{aligned} c_{WW}< 4\times 10^{-3}g_2^2\left( \frac{ \Lambda _{\text {UV}}}{10\text { TeV}}\right) ^2\frac{d_e}{1.1\times 10^{-29}\text {e\,cm}}\ . \end{aligned}$$

(35)

The effective operator in Eq. (34), is generated at 4Y-loops from the $\mathcal {O}_0$ operator and its axial partner and at 2-loop from the $\mathcal {O}_+$ operator and its axial partner. For example the 2-loops Wilson coefficients can be estimated as

$$\begin{aligned} c_{WW}\vert _{2-\text {loop}}\simeq \frac{g_2^2 n (n^2-1)^2\mathcal {I}[y_I y_{I,5}^{*}]}{32 (16\pi ^2)^2\Lambda _{ \text {UV}}^2}\ , \end{aligned}$$

(36)

where the index $y_I$ controls the contribution from $\mathcal {O}_0$ or $\mathcal {O}_+$ and $y_{I,5}$ the ones from their axial counterpart. If only $\mathcal {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

$$\begin{aligned} \delta m_0<16\text { GeV}\left( \frac{2}{n}\right) ^{\frac{3}{2}}\left| \theta _{\text {CP}}\right| ^{\frac{1}{2}}\sqrt{\frac{|d_e|}{1.1 \times 10^{-29}\text {e}\,\text {cm}}}\ , \end{aligned}$$

(37)

where we defined $\theta _{\text {CP}}=y_I^2/\mathcal {I}[y_I y_{I,5}^{*}]$. A similar bound is found for $\delta m_{ {Q_M} }$ in the minimal splitting case, where $\mathcal {O}_{0}$ can be neglected. Hence, for $\mathcal {O}(1)$ couplings and CP-violating phases, the EDM already gives competitive if not stronger bounds with respect to those shown in Table 1. Interestingly, future experimental upgrades are expected to be sensitive to electron EDM as small as $d_e\simeq 10^{-34}\text {e}\, \text {cm}$ [66] which would provide sensitivity to the WIMP parameter space down to splittings of order $\sim 50\text { MeV}$ for $\mathcal {O}(1)$ couplings and CP-violating phases.

It would be interesting to further explore this direction in models where the $\mathcal {O}(1)$ CP-violating phase is motivated by a mechanism producing the observed baryon asymmetry in the Universe such as WIMP baryogenesis [67, 68].

6 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 1, 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 uncertainties [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$ and the fiveplet with $Y=1$, whose cross sections lie 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 $\sqrt{s}=6 \text { 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.

Appendices

Appendix A: Complex WIMP classification

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

1.1 1. UV Operators

In Sect. 2 we showed that in order to make complex WIMPs with $Y\ne 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 $\mathcal {O}_{0}$ and the charged-neutral splitting by $\mathcal {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 $\chi _N$ after EWSB is

$$\begin{aligned} \overline{\chi }(y_0+y_{0,5}\gamma _5)\mathcal {I}\chi ^c H^{4Y}, \end{aligned}$$

(A1)

where the $H^{4Y}$ is necessary to match the hypercharge of the $\chi ^2$ piece, while $\mathcal {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 $\mathcal {O}_{0}$ defined in Eq. (1) and its axial counterpart with the $\gamma _5$ insertion. The latter can be shown to give only subleading contributions to the mass splitting between the neutral components of $\chi$. Assuming y to be real, the shift on $\delta m_0$ induced by $y_{0,5}\ne 0$ with respect to its expression in Eq. (3) is

$$\begin{aligned} \delta m_0 \rightarrow \delta m_0\sqrt{ 1 + 4\Re [y_{0,5}]^2\frac{\Lambda _{\text {UV}}^2}{M_\text {DM}^2}\left( \frac{v}{\sqrt{2}\Lambda _{\text {UV}}}\right) ^{8Y}}\, , \end{aligned}$$

(A2)

which is highly suppressed for $\Lambda _{\text {UV}}>M_\text {DM}>v$. The operator $\mathcal {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 $\mathcal {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

$$\begin{aligned} \overline{\chi }T^a(y_+ + y_{+,5}\gamma _5)\chi H^\dagger \sigma ^a H. \end{aligned}$$

(A3)

Similarly to the neutral case, the $\gamma _5$ insertion leads to a subleading shift in the mass splittings with respect to the value of Eq. (11)

$$\begin{aligned} M_Q^2\rightarrow M_Q^2 +y_{+,5}^2\frac{v^4}{16M_\text {DM}^2\Lambda ^2}(Q-Y)^2. \end{aligned}$$

(A4)

Finally, we comment on the implications on DD signals. Both operators with the $\gamma _5$ insertions give additional contributions to the SI cross-section and match to the effective operators:

$$\begin{aligned} \begin{aligned} \mathscr {L}_{\text {eff}}^{\mathrm{{SI}}} =&\tilde{f}_q m_q \bar{\chi } i\gamma _5\chi \bar{q} q + \frac{\tilde{g}_q}{M_\text {DM}} \bar{\chi } i \partial ^{\mu } \gamma ^{\nu } i\gamma _5\chi \mathcal {O}^q_{\mu \nu }\\&+ \tilde{f}_G \bar{\chi } i\gamma _5\chi G_{\mu \nu }G^{\mu \nu }. \end{aligned} \end{aligned}$$

However, all these operators are strongly momentum-suppressed [79, 80], so that their corrections to $\sigma _{\text {SI}}$ are expected to be $(q/m_N)^2 \sim 10^{-6}$ smaller than those from Eq. (16), and thus negligible.

1.2 2. Complex scalar WIMPs

Following the discussion in Sect. 2 for the fermions, supported by the previous Appendix A 1, the minimal Lagrangian for a scalar complex WIMP is

$$\begin{aligned}&\mathscr {L}_{\text {S}} =|D_\mu \chi |^2-M_\chi |\chi |^2 +\frac{y_0}{\Lambda _{\text {UV}}^{4Y-2}}\mathcal {O}_{0}^S + y_+\mathcal {O}_{+}^S+\text {h.c.}\, ,\nonumber \\&\mathcal {O}_{0}^S=\frac{1}{2(4Y)!}\left( \chi ^\dagger (T^a)^{2Y}\chi ^c\right) \left[ H^{c\dagger }\frac{\sigma ^a}{2} H\right] ^{2Y}\, , \\&\mathcal {O}_{+}^S=-\chi ^\dagger T^a\chi H^\dagger \frac{\sigma ^a}{2} H\, .\nonumber \end{aligned}$$

(A5)

Conversely to the fermionic case, no additional operators can be written to generate the fundamental mass splitting. The neutral and charged squared mass splitting, $\mu _Q^2=M_Q^2-M_\text {DM}^2$, can be written as

$$\begin{aligned} \begin{aligned} \mu _0^2&=4y_0 c_{nY0}\Lambda _{\text {UV}}^2\left( \frac{v}{\Lambda _{\text {UV}}\sqrt{2}}\right) ^{4Y}\, ,\\ \mu _{Q}^2&=\frac{\mu _0^2}{2}+2M_\text {DM}\delta _gQ^2+\text {sgn}(Q)\\&\quad \times \sqrt{\left( \frac{4Y\delta _gM_\text {DM}}{\cos \theta _W}-\frac{y_+v^2}{4}\right) ^2Q^2+\frac{\mu _0^2}{4}\frac{c_{nYQ}^2}{c_{nY0}^2}}. \end{aligned}\nonumber \\ \end{aligned}$$

(A6)

The linear mass splittings, $\delta m_Q=M_Q-M_\text {DM}$, are then given by $\delta m_Q=\mu ^2_Q/(2M_\text {DM})$. Notice that for $Y=1/2$ all the operators are renormalizable and any dependence on the cutoff disappears.

The lower bound on $\delta m_0$ from DD is identical to that for fermions, since $\sigma _\text {SI}$ does not change. Instead, the BBN bound differs since the one-loop decay channel $\chi _0\rightarrow \chi _\text {DM}\gamma$ now is now heavily suppressed with respect to the three-body decays. As a consequence, the BBN condition becomes

$$\begin{aligned} \Gamma _{\chi _0}\equiv \Gamma _{\bar{\nu }\nu }+\Gamma _{\bar{e}e}>6.58 \times 10^{-25}\text { GeV}\, , \end{aligned}$$

(A7)

which explains why the bound on $\delta m_0^{\text {min}}\simeq 4--5\text { MeV}$ for scalars is so much stronger than the one for fermions, as shown in Fig. 1. For $Y\ne 1/2$, the lower bound on $\delta m_0$ sets the upper bound of the allowed window for $\Lambda _{\text {UV}}$

$$\begin{aligned} 10\, M_{\text {DM}}<\Lambda _{\text {UV}}\le \left( \frac{2y_0c_{nY} v^{4Y}}{2^{2Y}\delta m_0^{\min }M_\text {DM}}\right) ^{\frac{1}{4Y-2}}\, , \end{aligned}$$

(A8)

which is the analogous of Eq. (13) for scalar WIMP. Once set $y_0=(4\pi )^{4Y}$ to its NDA maximal value and given the scaling $M_\text {DM}\sim 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 stability for scalars. All the $Y=1/2$ scalars are never accidentally stable, since DM decay can be induced by the renormalizable operator $\chi ^2\chi ^\dagger H^c$. Similarly, for $3_1$ we can write $\chi ^\dagger H^2$, while for $5_1$ the lowest dimensional operators are:

$$\begin{aligned} \mathscr {L}=\frac{C_1}{\Lambda _{\text {UV}}}\chi ^\dagger H^2 (H^\dagger H) + \frac{C_2}{\Lambda _{\text {UV}}}(\chi ^\dagger \chi )(\chi ^\dagger H^2)\, , \end{aligned}$$

(A9)

which ensures stability for $\Lambda _{\text {UV}}>10^{21}M_\text {DM}$, way larger than $\Lambda _{\text {UV}}<20M_\text {DM}$ from $\delta 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 particular, we highlight the main differences with respect to real WIMPs and refer the reader to [1, 81, 82] for further details.

1.1 1. 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) [81,82,83] . For WIMPs, the potential induced by EW interactions is attractive for all isospins $I\lesssim \sqrt{2}n$, 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

$$\begin{aligned} S_E^I \approx \frac{2\pi \alpha _\text {eff}}{v_\text {rel}}\ , \end{aligned}$$

(B1)

where $\alpha _\text {eff}= \alpha _2(I^2+1-2n^2)/8 + Y^2 \alpha _Y\equiv \alpha _{2,\text {eff}} + Y^2 \alpha _Y$. Hence, the annihilation cross-section reads

$$\begin{aligned} \langle \sigma _{\text {ann}}^\text {SE}v_\text {rel}\rangle =\sum _I \langle S_E^I\sigma _{\text {ann}}^Iv_\text {rel}\rangle \ , \end{aligned}$$

(B2)

where $\sigma _{\text {ann}}^I$ 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\ge 6$, while we compute the SE in the broken phase for smaller multiplets.

Besides, the mass splittings generated by $\mathcal {O}_{0}$ and $\mathcal {O}_{+}$ can alter the pattern of resonances in the SE, so that the thermal mass becomes a function of $\delta m_0$ and $\delta m_+$ as well. In Table 1 we assume the minimal $\delta m_0$, where the effect of the splitting is safely negligible, and account for the dependence on $\delta m_+$ in the theoretical uncertainty. More details about the effect of these UV operators are discussed in Sect. B 3. Another source of uncertainty in the computation of $\sigma _{\text {ann}}^\text {SE}$ comes from neglecting partial waves higher than the s-wave. We estimate the correction to $\sigma _{\text {ann}}^\text {SE}$ as

$$\begin{aligned} \frac{\sigma _{\text {ann},\text {p-wave}}^\text {SE}}{\sigma _{\text {ann},\text {s-wave}}^\text {SE}}\simeq v_\text {rel}^2+\alpha _\text {eff}^2\ . \end{aligned}$$

(B3)

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 $\chi _i + \overline{\chi }_j \rightarrow \text {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 $\chi \chi$ and the (anti-)symmetry of the wave function allows only for those BS satisfying $P_\text {BS}\equiv (-1)^{L+S+\frac{I-1}{2}}=1$ , while for complex WIMPs the bound states are of the form $\overline{\chi }\chi$ 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:

$$\begin{aligned} {\sigma _{B_I}v_\text {rel}\simeq \frac{E_{B_I}a_B}{M_\text {DM}^2 n^2}(\alpha _{2,\text {eff}}S_E^{I\pm 2}+Y^2\alpha _Y S_E^{I})}\ , \end{aligned}$$

(B4)

where the first and second term account for weak and hypercharge vector boson emission, respectively, $E_{B_I}\approx \alpha _\text {eff}^2M_\text {DM}/4n_B^2$ - $n_B$ being the energy level - is the binding energy and $a_B=1/\alpha _\text {eff}M_\text {DM}$ the Bohr radius.

Once the BS is formed, it can annihilate into SM, decay to lower lying BS, both with rate $\Gamma \approx \alpha _{2,\text {eff}}^2\alpha _\text {eff}^3M_\text {DM}$ or $Y^2\alpha _{Y}^2\alpha _\text {eff}^3M_\text {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 $\sim 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

$$\begin{aligned} \langle \sigma _{\text {eff}}v_\text {rel}\rangle \equiv \langle \sigma _{\text {ann}}^\text {SE}v_\text {rel}\rangle +\sum _{B_J}\langle \sigma _{B_J}v_\text {rel}\rangle , \end{aligned}$$

(B5)

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\le 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, see also Ref. [84]. 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 $\sigma _{B_I}$ for $T<\max \delta m_Q$ and then including BSF until the DM abundance saturates. These two effects are the dominant sources of error for $n\lesssim 8$.

For larger multiplets, the effective cross-section in Eq. (B5) suffers from additional theoretical uncertainties. Estimating from the results of Ref. [85] the NLO correction to the potentials controlling the SE we get

$$\begin{aligned} \frac{\Delta V_{\text {NLO}}}{V_{\text {LO}}}\sim \beta _2^{\text {SM}}\frac{\alpha _2}{4\pi }\log \left( \frac{m_W}{M_\text {DM}}\right) \, , \end{aligned}$$

(B6)

where $\beta _2^{\text {SM}}=-19/6$ is the SM contribution to the SU(2)$_L$ beta function. Equation (B6) results in a correction to the Sommerfeld factor $S_E$, which affects both $\sigma _{\text {ann}}^\text {SE}v_\text {rel}$ and $\sigma _{B_J}v_\text {rel}$ as introduced in Eq. (B4). NLO corrections start dominating the theoretical uncertainty for $n > rsim 8$.

1.2 2. 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, 86, 87]:

$$\begin{aligned} (\sigma _\text {eff} v_\text {rel})^{J}\le \frac{4\pi (2J+1)}{M_\text {DM}^2 v_\text {rel}}\ , \end{aligned}$$

(B7)

where $\vec {J} = \vec {L} + \vec {S}$ is the total angular momentum. The stronger inequality comes from the s-wave channel (i.e. $J=0$) which can be written as

$$\begin{aligned} (\sigma _\text {ann} v_\text {rel})+\sum _{B_J} f_{B_J}^0 (\sigma _{B_J} v_\text {rel})\le \frac{4\pi }{M_\text {DM}^2 v_\text {rel}}\ , \end{aligned}$$

(B8)

where $f_{B_i}^0$ 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_I^{S=1}$ states.

1.3 3. 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 $\mathcal {O}_0$ contributions. We checked that similar conclusions hold for scalar WIMPs and for the contributions from $\mathcal {O}_+$.

For $Y=1/2$ WIMPs, the hard annihilation cross-section into Higgses induced by $\mathcal {O}_{0}$ can be estimated as

$$\begin{aligned} \sigma _{2H} v_\text {rel}\simeq \frac{3y_{0,+}^2}{32n^2\pi \Lambda _{\text {UV}}^2}T_R\, , \end{aligned}$$

(B9)

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

$$\begin{aligned} \sigma _{\text {EW}} v_\text {rel}\simeq \frac{\pi \alpha _2^2(n^2-1)^2}{32nM_\text {DM}^2}\, . \end{aligned}$$

(B10)

The condition $\sigma _{2H}/\sigma _{\text {EW}}<1$ can be translated into an upper bound on the mass splittings. For instance for the inelastic mass splitting we get

$$\begin{aligned} \delta m_0< 1.6\text { GeV}\left( \frac{1\text { TeV}}{M_\text {DM}}\right) n^2\ . \end{aligned}$$

(B11)

Table 2 Freeze-out mass predictions for millicharged WIMP DM. The annihilation cross-section includes both the contribution of SE and BSF. For the triplets, a second prediction, denoted with a *, for the thermal mass is present due to the emergence of a large resonance in the Sommerfeld enhancement for $M_\text {DM}\approx 2.4$ TeV. We provide a measure of how close the DM annihilation cross-section is to the unitarity bound for s-wave annihilation $(\sigma v)^{J=0}_{\text {max}}=4\pi /M_{\text {DM}}^2v$. We derive the scale where EW gauge coupling will develop a Landau pole by integrating-in the WIMP multiplet at its freeze-out mass

Since $M_\text {DM}\sim 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 Sect. 4 even though for large splittings (above 1 GeV) one can get $\mathcal {O}(1)$ effects on the thermal freeze-out predictions for small multiplets.

For $Y=1$ the contribution from $\mathcal {O}_{0}$, which now is a 7D operator which controls the 4-body process $\overline{\chi }\chi \rightarrow 4H$ whose cross section can be estimated as

$$\begin{aligned} \sigma _{4H} v_\text {rel}\approx \frac{\pi y_0^2M_\text {DM}^4}{(16\pi ^2)^3\Lambda _{\text {UV}}^6}\frac{(n^2-1)^2}{64n}\, . \end{aligned}$$

(B12)

The condition $\sigma _{4H}/\sigma _{\text {EW}}<1$ translates into

$$\begin{aligned} \delta m_0<0.45\text { GeV}\left( \frac{1\text { TeV}}{M_\text {DM}}\right) ^3 n^2, \end{aligned}$$

(B13)

which, for $n=3,5$ is much looser than the upper bound on $\delta 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:

$$\begin{aligned} V_H(r)=-\left( \frac{\lambda _D v}{8\pi \Lambda _{\text {UV}}}\right) ^2\frac{e^{-m_h r}}{r}. \end{aligned}$$

(B14)

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 $\alpha _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 $\epsilon F F'$ and the complex WIMP would acquire a EM charge $\epsilon$. 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 $\epsilon$ (i.e. $\epsilon <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_\text {BS}= (-1)^{L+S+\frac{I-1}{2}}=-1$ formed by $\bar{\chi }\chi$ pairs of millicharged WIMPs. These are forbidden by the spin-statistic properties of the $\chi \chi$ wave function for real WIMPs. Since $P_\text {BS}$ is preserved by dipole interactions for $Y=0$, to leading order no transitions can occur between states with opposite $P_\text {BS}$ and excited BS with $P_\text {BS}=-1$ will dominantly decay to 1s and 2s states with the same $P_\text {BS}$. The latter have small decay widths with respect to their $P_\text {BS}=1$ counterparts. $ns_1^{S=1}$ and $ns_5^{S=1}$ annihilate into four vectors with a rate

$$\begin{aligned} \Gamma (ns_{1,5}^{S=1}\rightarrow VVVV)\simeq \frac{\alpha _\text {eff}^4}{16\pi ^2 M_\text {DM}^2}\frac{|R_{n0}(0)|^2}{M_\text {DM}^2}\, , \end{aligned}$$

(C1)

while $ns_3^{S=0}$ annihilates into three vectors with rate:

$$\begin{aligned} \Gamma (ns_{3}^{S=0}\rightarrow VVV)\simeq \frac{\alpha _\text {eff}^3}{4\pi M_\text {DM}^2}\frac{|R_{n0}(0)|^2}{M_\text {DM}^2}\, . \end{aligned}$$

(C2)

where $R_{n0}\sim (\alpha _\text {eff}M_\text {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 2.

Fig. 7

Expected SI cross-sections for different millicharged complex WIMPs with $\epsilon \lesssim 10^{-10}$. The blue dots correspond to Dirac WIMPs and the red dots to complex scalar WIMPs. The vertical error bands correspond to the propagation of LQCD uncertainties on the elastic cross-section (Eq. 18), while the horizontal error band comes from the uncertainty in the theory determination of the WIMP freeze out mass in Table 1. The light green shaded region is excluded by the present experimental constraints from XENON1T [36] and PandaX-4T [5], the green dashed lines shows the expected 95% CL reach of LZ/XENONnT [8, 9] and DARWIN [10, 11]

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 $\epsilon \lesssim 10^{-10}$ [3, 88], while larger $\epsilon$ 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 $\sqrt{s}\le 14$ TeV impossible for $n\ge 5$. We consider masses $M > rsim \sqrt{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_\text {DM}$, a scan in the $(\delta m_0, \delta m_ {Q_M} )$ to compute the signal, since the spectrum and the decay width of the particle $\chi ^Q$ is determined uniquely by the splittings $\delta m_ {Q_M}$, $\delta m_0$.

To get a feel of how powerful DT searches can be, we fix the neutral splitting to a representative value $\delta m_0 \simeq 200 \textrm{KeV }$ for the fermions. Our results are given in Fig. 8 where we display the region of the plane WIMP mass versus $c\tau$ where experiments can probe the several WIMPs considered in each panel of the figure. For this result we consider only $\chi ^+$ 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.⁵ As we have fixed $\delta m_0$ much smaller than the mass splitting between the neutral and charged states, we can effectively consider decays into both $\chi ^0$ and $\chi ^{\textrm{DM}}$ as if they were degenerate in mass. Considering the possible decay channels $\chi ^-$ into $e^- \nu _e \chi ^{0,\textrm{DM}}\,$, $\mu ^- \nu _\mu \chi ^{0,\textrm{DM}}\,$, $\pi ^- \chi ^{0,\textrm{DM}}$ and the charge conjugates for $\chi ^+$ we get

$$\begin{aligned} \Gamma \left( \chi ^\pm \rightarrow \chi ^{0,\textrm{DM}}+\textrm{SM} \right) = \Gamma _{e^\pm }+\Gamma _{\mu ^\pm }+\Gamma _{\pi ^\pm } \;, \end{aligned}$$

(D1)

Fig. 8

Sensitivity at 95% CL from DT searches as a function of $M_\chi$ and the proper lifetime of the charged particle $\chi ^+$. The mass splitting $\delta m_+$ is also shown in the extra vertical axis in two upper panels. In the lower panels, the two different contours for each search channel correspond to the different physical situations where the heavy mass eigenstate corresponds to one of the two gauge eigenstates $\chi _+$ or $\chi _{-}^{c}$. Top Left: Dirac $2_{1/2}$ at $\sqrt{s}=6$ TeV, $\mathcal {L}= 4\; \textrm{ab}^{-1}$. Top Right: Dirac $3_1$ at $\sqrt{s}=10$ TeV, $\mathcal {L}= 10\; \textrm{ab}^{-1}$. Bottom Left: Dirac $4_{1/2}$ at $\sqrt{s}=14$ TeV, $\mathcal {L}= 20\; \textrm{ab}^{-1}$. Bottom Right: Dirac $5_1$ at $\sqrt{s}=30$ TeV, $\mathcal {L}= 90\; \textrm{ab}^{-1}$

Fig. 9

Different bars show the 2 $\sigma$ (solid wide) and 5$\sigma$ (hatched thin) reach on the fermionic WIMP mass at a muon collider for different search channels. The first three bars show the channels discussed in Sect. 4 where DM would appear as missing invariant mass (MIM) recoiling against one or more SM objects: mono-$\gamma$, inclusive mono-W, and the combination of all these MIM channels (blue). The mono-Z channel is not reported since it gives results below the minimum mass shown. The last two bars show the reach of disappearing tracks, requiring at least 1 disappearing track (red), or at least 2 tracks (orange). All the results are shown assuming systematic uncertainties to be 0 (light), 1‰(medium), or 1% (dark) for MIM searches, while they are taken 0 (light), 1% (medium), 10% (dark) for the single DT search. The double DT search is not affected by systematics.. The vertical red bands show the freeze-out prediction. Tracks have been computed assuming minimally split neutral states, and a charged-neutral splitting $\delta m_+$ equal to the gauge contribution for the $2_{1/2}$ and $3_1$. For the $4_{1/2}$ and $5_1$ the splitting-inducing couplings $y_0$ and $y_+$ have been tuned so that $\delta m_+ = \delta m_-$ (up to minimal $\delta m_0$ corrections), which represents a least favorable condition for DT searches. Top left: Dirac $2_{1/2}$ for $\sqrt{s}=6$ TeV and $\mathcal {L} = 4 \; \textrm{ab}^{-1}$. Top right: Dirac $3_{1}$ for $\sqrt{s}=10$ TeV and $\mathcal {L} = 10 \; \textrm{ab}^{-1}$. Bottom left: Dirac $4_{1/2}$ for $\sqrt{s}=14$ TeV and $\mathcal {L} = 20 \; \textrm{ab}^{-1}$. Bottom right: Dirac $5_{1}$ for $\sqrt{s}=30$ TeV and $\mathcal {L} = 90 \; \textrm{ab}^{-1}$

Table 3 95% CL reach on the mass of the Dirac $2_{1/2}$ and $3_{1}$ 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 $\eta _X$ and $p_{T,X}$,where X is the single vector boson.The numbers are shown for different collider energies $\sqrt{s}$ and systematic uncertainties $\epsilon _{\textrm{sys}}$

Table 4 95% CL reach on the mass of the Dirac $4_{1/2}$ and $5_{1}$ 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 $\eta _X$ and $p_{T,X}$,where X is the single vector boson.The numbers are shown for different collider energies $\sqrt{s}$ and systematic uncertainties $\epsilon _{\textrm{sys}}$

where the RHS indicates the decay width into neutral states plus the SM states indicated by the subscript. The relevant widths are given by:

$$\begin{aligned} \begin{aligned}&\Gamma _{e^\pm }= g\left( n,Y, \pm 1 \right) \frac{G_F^2 \delta m_\pm ^5}{60 \pi ^3}\;,\\&\Gamma _{\mu ^\pm }=g\left( n,Y,\pm 1 \right) \frac{G_F^2 \delta m_\pm ^5}{2\pi ^3} \Phi (\delta m_\pm ,m_\mu ) \;,\\&\Gamma _{\pi ^\pm }=g\left( n,Y, \pm 1 \right) \frac{G_F^2 f_\pi ^2 |V_{ud}|^2 \delta m_\pm ^3}{4\pi }\sqrt{1-\frac{m_\pi ^2}{\delta m_\pm ^2}} \; , \end{aligned} \end{aligned}$$

(D2)

where $g\left( n,Y,Q \right)$ accounts for the different strengths of the W boson coupling to each of the $\chi ^{Q}$:

$$\begin{aligned} g\left( n,Y,Q \right) =n^2-1-4(Q-Y)(Q-Y-{\text {sgn}}(Q) )\,,\nonumber \\ \end{aligned}$$

(D3)

with n the dimensionality of the multiplet and Y its hypercharge, $G_F$ is the Fermi constant, $f_\pi =131$ MeV is the pion decay constant. In the above formula $\Phi$ 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 \ll M$ (e.g. muon), and a heavy particle with a mass $M-\delta$, for splitting $\delta \ll M$:

$$\begin{aligned} \Phi (\delta ,m_l)= & {} \frac{1}{60}\sqrt{\delta ^2-m_l^2}\left( 2 \delta ^4-9\delta ^2m_l^2-8m_l^4 \right) \nonumber \\{} & {} +\frac{1}{4}m_l^4\delta {\text {ArcCoth}}\left( \frac{\delta }{\sqrt{\delta ^2-m_l^2}}\right) \;. \end{aligned}$$

(D4)

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.

Fig. 10

Results for scalar WIMPs. Same notation as in Fig. 9

If we move away from the $\delta m_0 \simeq 0$ points, we have to consider cases in which $\chi ^+$ 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_{\chi ^{++}}<M_{\chi ^+}$. In this case, $\chi ^{++}$ cannot decay into $\chi ^+$, and is forced to decay to the neutral states by emitting twice an off-shell W boson that gives rise to $l^+ \nu _l$ or $\pi ^+$. 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_\pi \delta /(16\pi ^2)$ and $g(n,y,2)G_F \delta ^2/(16\pi ^2)^2$ for any extra $\pi$ 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 $\chi ^{++}$ 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_\gamma >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 $\sigma _B=18.8$ ab at the 10 TeV collider, leading to a signal-to-noise ratio of order 1. The latter has negligible background. We use these cross sections in order to compute the number of background events needed to get the significance in Eq. (21). To compute the number of expected signal events, we do a Monte Carlo simulation with MadGraph5_aMC@NLO [89]. For a more detailed discussion on how to compute DT signal events and reach, we refer to the Appendix of [1]. Results are consistent with the analytic approximation explained in the same work.

The sensitivity to disappearing tracks is shown in the plots of Fig. 8, for minimal $\delta m_0$ imposed by DD and BBN, as a function of the dark matter mass and the lifetime $c \tau$ 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 $\chi ^+$ corresponds to the gauge eigenstate $\chi _+$ or $\chi _-^c$ (corresponding to the regions of low or high $\delta 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 $\sqrt{s}=3, 6$ TeV and $\sqrt{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 $\epsilon$ we have interpolated our previous results as functions of n and $\sqrt{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 MadGraph5_aMC@NLO [89] using the full 3-body process for $\eta \lesssim 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 $\sqrt{s}$ that can provide 95% CL exclusion. Full results for all benchmark colliders at the nominal luminosity of Eq. (20) are given in Tables 3 and 4, for the various search channels under consideration. The collider reach for scalar candidates are instead given in Fig. 10. Results for generic center-of-mass energy and integrated luminosity are displayed in Fig. 11 for the mono-$\gamma$ and mono-W channels. The results on this figure have been obtained by rescaling the results at $\sqrt{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.

Fig. 11

Mass reach (in TeV) in the mono-$\gamma$ and mono-W channels as a function of collider center-of-mass energy and luminosity, for Dirac fermion doublets, triplets, and 4-plets (blue lines). The thermal freeze-out mass is shown in red. Blue shades show the expected values of the ratio of the signal rate over background. Systematic uncertainties are set to zero. The muon collider luminosity Eq. (20) is shown as a white line, with the benchmark values of $\sqrt{s}$ highlighted by the colored squares

Fig. 12

Mass reach in the mono-$\gamma$, mono-W and DT channels for fixed luminosity as per Eq. (20) at $\sqrt{s}$ 3 TeV (yellow), 6 TeV (green), 10 TeV (light blue), 14 TeV (red), and 30 TeV (purple). In the mono-W and mono-$\gamma$ searches we show an error bar, which covers the range of possible exclusion as the systematic uncertainties are varied from 0 to 1%. The colored bars are for an intermediate choice of systematics at 0.1%. In the 1 DT search the error bar corresponds to uncertainties from 0 to 10%, and the colored bars to 1%. The 2 DT search is not affected by systematics. Missing bars denoted by an asterisk * correspond to cases where no exclusion can be set in the mass range $M_\chi >0.1\sqrt{s}$. For such cases it is worth considering VBF production modes at the fixed luminosity Eq. (20) or higher luminosity at potentially smaller $\sqrt{s}$ as illustrated in Fig. 11

The DT results shown in Fig. 9 are a particular sub-case 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 asymptotic gauge contribution $\delta m_+=354$ MeV. Its 1-loop value computed at the thermal mass is 343 and 335 MeV for the fermion and scalar respectively. Using the latter values increases the 2$\sigma$ mass reach by 10% (20% for the two DT search for the scalar). This does not change appreciably the minimum collider energy required for exclusion. For the $3_1$ the splitting has been set to its gauge value of $\delta m_+\simeq 540$ MeV. For the $4_{1/2}$, the splitting $\delta m_+$ has been taken equal to its value at the maximal mixing point, for which it is equal to $\delta m_-$, up to corrections of order $\delta m_0$. This choice implies that the $\chi ^+$ and $\chi ^-$ particles are maximal admixtures of the gauge eigenstates, hence they have same cross-sections and lifetimes. The latter are governed by the mass splitting $\delta m_+$ and $\delta m_-$, which, from Eq. (11) are constrained to satisfy $\delta m_+ + \delta m_- = \delta m_0 +2\delta _g$. This equation implies that for $\delta m_0 \ll \delta _g$ in the $\delta m_+ =\delta 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 $\chi ^{++}$ has been assumed to promptly decay equally into $\chi ^+$ and $\chi ^{-,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.