Leptonic WIMP Coannihilation and the Current Dark Matter Search Strategy

We discuss the extent to which models of Weakly Interacting Massive Particle (WIMP) Dark Matter (DM) at and above the electroweak scale can be probed conclusively in future high energy and astroparticle physics experiments. We consider simplified models with bino-like dark matter and slepton-like coannihilation partners, and find that perturbative models yield the observed relic abundance up to at least 10 TeV. We emphasise that coannihilation can either increase or decrease the dark matter relic abundance. We compute the sensitivity of direct detection experiments to DM-nucleus scattering, consider indirect detection bounds and estimate the sensitivity of future proton colliders to slepton pair production. We find that current and future experiments will be able to probe the Dirac DM models up to at least 10 TeV. However, current and future searches will not be sensitive to models of Majorana dark matter for masses above 2 or 4 TeV, for one or ten coannihilation partners respectively, leaving around 70 % of the parameter space unconstrained. This demonstrates the need for new experimental ideas to access models of coannihilating Majorana dark matter.

We discuss the extent to which models of Weakly Interacting Massive Particle (WIMP) Dark Matter (DM) at and above the electroweak scale can be probed conclusively in future high energy and astroparticle physics experiments. We consider simplified models with binolike dark matter and slepton-like coannihilation partners, and find that perturbative models yield the observed relic abundance up to at least 10 TeV. We emphasise that coannihilation can either increase or decrease the dark matter relic abundance. We compute the sensitivity of direct detection experiments to DM-nucleus scattering, consider indirect detection bounds and estimate the sensitivity of future proton colliders to slepton pair production. We find that current and future experiments will be able to probe the Dirac DM models up to at least 10 TeV. However, current and future searches will not be sensitive to models of Majorana dark matter for masses above 2 or 4 TeV, for one or ten coannihilation partners respectively, leaving around 70% of the parameter space unconstrained. This demonstrates the need for new experimental ideas to access models of coannihilating Majorana dark matter.

I. INTRODUCTION
Understanding the nature of Dark Matter (DM) is one of the most pressing questions in particle physics. Its existence is well established by a wide range of astrophysical observations and its energy density is measured to 2% accuracy [1]. A thermally produced WIMP (Weakly Interacting Massive Particle) has long been the dominant paradigm. In this picture, dark matter is assumed to have non-negligible interactions with Standard Model (SM) particles. In the early universe, the temperature was very high so the standard model particles and dark matter populated a thermal bath. As the temperature cooled below the mass of the dark matter particle, it self-annihilated more often than it was produced, and so its abundance dropped (it became Boltzmann suppressed). As the universe expanded, the annihilation became inefficient, and the dark matter particles could no longer annihilate -the dark matter froze-out -leaving behind a relic abundance. This picture successfully predicts the observed relic abundance of dark matter if there is a weak-scale interaction cross-section with SM particles. This success, as well as other hints that beyond the standard model physics, such as supersymmetry (SUSY) or new strong dynamics, may be found slightly above the weak scale, have led to a strong theoretical and experimental exploration of the thermal WIMP.
The canonical WIMP is the lightest neutralino of the minimal supersymmetric standard model (MSSM). It is both a well motivated dark matter candidate in its own right and, as an admixture of neutral binos, winos and higgsinos, a powerful parameterisation of a wide range of WIMP models. As such, there has been a large effort to probe its parameter space. Although direct and indirect detection experiments are currently probing this parameter space, no signal has yet been seen. The LHC is also probing the motivated parameter space, but is yet to find signals of a WIMP or any other new physics particles. Although there are many experiments planned for the future, the clear question to answer is: 'will they probe the whole neutralino parameter space?' As such, it is important to identify viable scenarios in the sub-GeV and multi-TeV mass region, and to consider whether the suite of current and proposed future experiments will successfully probe the entire region. It has been shown [2,3] that current and future experiments will be able to probe the neutralino relic surface (where the parameters of the theory are restricted to produce the observed relic abundance via thermal freeze-out) up to masses of 4 TeV, if the sfermions are decoupled. However, once coannihilation with sfermions is taken into account, a larger and more challenging parameter space becomes accessible. It is precisely this scenario we consider in the current work.
Coannihilation [4] has been studied for a long time. It occurs during thermal freeze-out when there are other dark sector particles, φ, similar in mass to the dark matter particle, χ. Freeze-out occurs at temperatures where the abundances of χ and φ are significantly Boltzmann suppressed. In this situation the relic abundance of χ may be reduced if it can effectively annihilate via φ, or its relic abundance may be increased if φ cannot effectively annihilate [5][6][7]. This effect can dramatically change the relic abundance and consequently has an important impact on the relic surface. There has recently been considerable interest in the range of possible coannihilation models [8], their role in producing viable sub-GeV [9,10] and multi-TeV scale [11][12][13][14][15] dark matter candidates and in coannihilating models at the LHC and future colliders [16][17][18]. In this work we use a simplified model framework to explore the impact of coannihilation on multi-TeV dark matter. We consider a bino motivated (gauge-singlet Dirac or Majorana fermion) dark matter candidate accompanied by n dark-sector scalars with unit hypercharge. In the MSSM, a pure bino with no other nearby states cannot efficiently annihilate, resulting in overclosure of the universe. However, when sfermions are included, the observed relic abundance can be recovered for a relatively wide range of masses. We consider the three possible Yukawa couplings with SM electrons, muons and taus individually. This minimal setup lets us study the impact of coannihilation in isolation, and vary the degree of coannihilation by changing the mass difference (between dark matter and the new scalars) and the number of coannihilating partners. We first find, in section II, the relic surface for a range of models, demonstrating that they provide a viable multi-TeV dark matter candidate. We then consider the current and future reach of a range of direct detection experiments (section III), indirect detection experiments (section IV), and proton colliders (section V). We will see that there is a large region of viable parameter space for Majorana dark matter which current and future experiments will be unable to probe, motivating the need for new experimental ideas.
The experimental landscape in dark matter physics consists of colliders, direct detection experiments and indirect detection experiments. The LHC is currently running at 13 TeV and has delivered approximately 100 fb −1 of integrated luminosity to the ATLAS and CMS experiments. This has allowed the experiments to place significant bounds on simplified dark matter models mostly via mono-γ, Z, W, h, t and mono-jet searches. ATLAS and CMS also search directly for the mediators in di-jet or di-lepton plus missing energy searches [19]. For an exhaustive list of possible coannihilating DM searches at the LHC see, e.g,. [8]. A higher energy collider will be required to efficiently produce multi-TeV particles. Currently under discussion are a 27 TeV high energy upgrade to the LHC, dubbed HE-LHC, which would deliver approximately 15 ab −1 of integrated luminosity [20] and a 100 TeV collider, either in Europe or in China, which would deliver approximately 20 ab −1 [21]. Future lepton colliders, such as the ILC [22], CLIC [23] and the FCC-ee [24], are also under consideration. CLIC is designed to reach the highest centre of mass energy among these machines with √ s = 3 TeV. However, this is still not high enough to produce the multi-TeV particles we discuss. It should be noted that although a collider may produce dark sector particles, it cannot determine that any particle is the cosmologically stable dark matter since it can only test particle stability on detector scales.
Direct detection experiments consist of a body of shielded target material. A dark matter particle in our galaxy may interact with the target and deposit energy, which may then be detected as light, heat or ionisation. These experiments have been designed to operate in the dark matter mass range 10 GeV -1 TeV, and the current leaders in sensitivity are LUX [25], PandaX-II [26] and XENON1T [27], which all use xenon as their target. In what follows we take LUX as the illustrative Table I. The new particles we introduce in section II with their respective charges, the number of copies we consider and the number of degrees of freedom.
example. In the future, the most ambitious is the planned DARWIN experiment, which aims to have a sensitivity 100 times better than these experiments with an exposure of 500 ton · years [28]. The final class of experiment is indirect detection of dark matter. Although dark matter stopped annihilating when it froze-out, due to its low number density, gravity interactions have now caused dark matter to cluster, in haloes which encompass galaxies and galaxy clusters. Indirect detection experiments look for annihilation of dark matter where its abundance is expected to be largest. Since thermal dark matter particles could annihilate into any standard model particles, there are a range of strategies looking for photons, neutrinos and a range of anti-matter produced in the galactic centre or in dwarf galaxies. In this work we find the best limits from searches for an excess in continuum photons (as opposed to mono-energetic photons). We find that the strongest constraints are placed on our model by the Fermi-LAT -MAGIC collaboration [29] and HESS [30]. The Fermi-Large Area Telescope (Fermi-LAT) is the principle scientific instrument on the Fermi Gamma Ray Space Telescope spacecraft and is a high-energy gamma-ray telescope covering the energy range from about 20 MeV to more than 300 GeV. MAGIC (Major Atmospheric Gamma Imaging Cherenkov Telescopes) is a system of two ground-based atmospheric Cherenkov telescopes. HESS (High Energy Stereoscopic System) is an array of four ground-based atmospheric Cherenkov telescopes which measure cosmic photons in the energy range from 10s of GeV to 10s of TeV. In the future, the most ambitious planned experiment which improves on these bounds is CTA [31]. CTA (Cherenkov Telescope Array) is the next generation ground-based array, which will operate in a similar energy band with several tens of telescopes.

II. THE MODELS AND THEIR RELIC SURFACES
In this work we focus on several related simplified models. As shown in table I, we introduce dark matter as a Majorana or Dirac fermion, χ, and n copies of an uncoloured scalar coannihilation partner, φ i , with unit hypercharge. The dark matter particle and each coannihilation partner couple to a standard model right-handed charged lepton. In addition to kinetic and mass terms, the Lagrangian only has one new interaction term (ignoring the scalar quartic, which plays no role in our phenomenology) where D µ = ∂ µ −ig Y B µ and the coupling is taken to be universal, i.e., y χ is the same for all φ i . We consider the cases R = e R , µ R and τ R , and assume that all φ i have the same mass, m φ i = m φ , and that m χ < m φ . For illustration, we focus on n ∈ {1, 3, 10}, which will allow us to show the impact of one, several and many coannihilation partners. In a supersymmetric context, the DM particle χ would correspond to a bino and the scalar φ i can be identified with a right-handed slepton. Note that, in SUSY, the number of degrees of freedom of one right-handed slepton corresponds to n = 1, all right-handed sleptons corresponds to n = 3, while all right-and left-handed sleptons correspond to n = 9. Here, we follow a simplified model approach in order to isolate the effect of coannihilation from the added complications due to considering several flavours at once or nontrivial su(2) quantum numbers. Depending on the single lepton flavour involved, we will refer to our models as electron, muon and tau type. We are interested in the slice of parameter space where the models produce the observed relic abundance of χ via thermal freeze-out. In the following we will denote generic standard model bath particles as ψ, whether they are fermions or bosons. In the coannihilation regime, where the Boltzmann equation for the abundance of χ becomes a coupled set of differential equations which also track the abundance of φ i . These can be combined [4] into a single differential equation, the same as the usual Boltzmann equation for a single species, if the χχ →ψψ annihilation cross-section is replaced by where i, j index the DM particle and its coannihilation partners {χ, φ 1 , φ 2 , . . . , φ N }, g i is the number of degrees of freedom of particle i, ∆ i = (m i − m χ )/m χ , the cross-section is σ ij = σ(ij →ψψ) and x = m χ /T . The effective number of degrees of freedom is given by Note that g eff is always greater than g χ . For Majorana and Dirac DM g χ = 2 and 4, respectively. The degrees of freedom of a complex scalar φ i are g φ = 2. Some representative coannihilation diagrams which contribute to σ eff in our model are shown in fig. 1. We can understand some features of these equations on physical grounds. The abundances of both χ and φ are similarly Boltzmann suppressed during freeze-out. This means that the rate of dark sector annihilation χχ →ψψ, χφ →ψψ and φφ * →ψψ, where two rare particles are in the initial state, becomes exponentially smaller than the conversion processes χψ → φψ, which requires only one rare particle and one bath particle. Thus, even if χ annihilation has a very small cross-section, χ can be efficiently depleted by first converting into φ and then annihilating. That is, if σ χφ σ χχ and/or σ φφ σ χχ , and both ∆ and g eff are not too large, then and, since Ω χ h 2 ∼ 1/ σv , coannihilation reduces the relic abundance. This is the usual mechanism used to deplete bino dark matter. As ∆ becomes larger, φ is Boltzmann suppressed to a larger degree, and there is not enough thermal energy to efficiently achieve χψ → φψ, so the mechanism becomes ineffective. We note that although coannihilation is usually thought to increase the effective cross-section (as above), it can also reduce the effective cross-section, increasing the dark matter relic abundance (as noted in, e.g., [5][6][7]). If σ χχ σ χφ , σ φφ but χψ → φψ is still efficient, then the terms in the sum with i = χ and j = χ become negligible and we are left with is a generic standard model particle).
i.e., the cross-section has effectively been reduced by g 2 χ /g 2 eff . We can understand this situation by imagining a temperature above the temperature at which χ freezes out, but below that at which φ would have frozen-out if χ were not present. The energy density which resides in φ cannot go into the thermal bath via φφ →ψψ, since the cross-section is too small, but it can go via φψ → χψ, since this rate is not doubly Boltzmann suppressed. This can be thought to 'top-up' the abundance of χ during freeze-out, resulting in a higher final χ abundance. Again, coannihilation becomes ineffective as ∆ → ∞ as φ can fully annihilate before χ has frozen-out, eliminating the 'top-up'.
To calculate the relic abundances in our models, we use model files written with SARAH v4.12.1 [32] and calculate the relic abundance using micrOMEGAs v4.3.5 [33], which implements coannihilation. We then interpolate the results of a 3-dimensional scan (in m χ , ∆ and y χ ) to determine the value of the coupling y χ which will result in the observed relic abundance of The results for n = 1, 3, 10 are shown in fig. 2. We consider the parameter space 100 GeV < m χ < 10 TeV and 0.01 < ∆ < 0.4 for all 18 cases (Majorana and Dirac χ, each with n ∈ {1, 3, 10}, each for the electron, muon and tau type models). We focus on this region of parameter space since a range of strategies are being pursued for coannihilating dark matter around 100 GeV [34][35][36] and for sub-GeV dark matter [37][38][39], while upcoming proton-proton colliders have little prospect of probing particles heavier than 10 TeV. The ∆ < 0.01 region is extremely fine-tuned, requiring significant theoretical motivation, while the region ∆ > 0.4 will not exhibit significant coannihilation. The relic surfaces for all 18 models can be presented in six plots. Since we are always in the limit m χ m , the impact of the lepton mass on the relic surface is negligible, so the relic surface is independent of the lepton flavour of the model. We see that the observed relic abundance can be reached for perturbative couplings (y χ < 4π/ √ n ≈ {13, 7.3, 4.0} for n = {1, 3, 10}, respectively). Although the coupling remains < 4π/ √ n, it is relatively large in much of our parameter space, which suggests that higher order corrections to our tree-level and one-loop calculations may not be insignificant.
In fig. 2 (top-right) we see the relic surface for the Majorana models with n = 1. As m χ increases, the required coupling increases, as is expected since the annihilation cross-section scales as σ χχ ∼ y 4 χ m −2 χ . Since the relic abundance is roughly inversely proportional to σ χχ , y χ needs to increase as m χ increases to keep σ χχ approximately constant. As ∆ becomes smaller than 0.1, we begin to see the effect of coannihilation. For Majorana χ, σ χχ is velocity suppressed and so is significantly smaller than σ φφ . As coannihilation becomes relevant, σ eff becomes larger, which would reduce the relic abundance, if y χ did not reduce to compensate. On the relic surface, we see the required reduction in y χ . Although coannihilation is not active above ∆ ≈ 0.15, we note that the scalar partners still allow χ to have the correct relic abundance (which would not be the case if they were completely decoupled).
In fig. 2 (top-left) we see the relic surface for the Dirac models with n = 1. Again we see y χ increasing as m χ increases. However, as ∆ reduces below 0.1, y χ now increases. This is both because σ χχ is not velocity suppressed for Dirac dark matter and because in most of the parameter space y χ > 1, which is greater than the electromagnetic coupling of φ. As such, the extra crosssections we add into eq. (3) are small, and the dominant effect is to reduce σ eff due to the increase of g eff . In this situation, coannihilation increases the relic abundance.
In fig. 2 (middle) and (bottom) we see the impact of increasing the number of coannihilation partners. Extra partners change the required y χ both at large and small ∆. We see an effect at large ∆, where the relic abundance is set by σ χχ , since the only contribution to σ χχ is a tchannel diagram with a φ propagator. As such, increasing the number of partners will increase this cross-section by n 2 . To maintain the observed relic abundance, the coupling y χ has to decrease accordingly. As ∆ goes to zero, for Dirac χ, more coannihilating partners simply give a stronger increase in y χ . For Majorana χ, we see a balance between an increase in the effective cross-section due to the extra coannihilation processes and a decrease in σ eff due to an increasing g eff , which is especially pronounced in the n = 10 case. For Majorana dark matter, σ χχ is velocity suppressed while σ φφ grows with n 2 . The χ contribution to g eff , however, is not suppressed in any way, so g 2 eff is not simply proportional to n 2 . At ∆ < 0.05, the dominant effect comes from an increase in σ eff due to σ φφ . At 0.01 < ∆ < 0.2, we see that an increased g eff , which reduces σ eff , is the dominant effect.

III. DIRECT DETECTION
Current and future direct detection experiments will place important constraints on our models. In these experiments, dark matter in the neighbourhood of the earth may pass through the detector and interact with the nucleus of one of the target atoms. The energy deposited causes the emission of light, charge (electrons) and heat. Direct detection experiments are typically sensitive to two of these three signals, and use them to place a limit on the rate of interactions seen in the target material. They then translate these into a bound on the DM-nucleon interaction cross-section, assuming a contact interaction. However, as we will see, the models we are discussing do not have a contact interaction but instead dipole and anapole interactions. As such, we will consider the expected interaction rate between our DM models and the target nuclei and compare them to the rates that can be probed in experiments. The DM-nucleus scattering rate per unit target mass is Figure 4. The anapole moment (left) for electrons (orange), muons (green) and taus (blue) running in the loop. For the electron case, we take a representative E r = 50 keV. The anapole moment is larger for smaller lepton masses. The differential rate (right) as a function of the recoil energy for a dipole moment (black) and anapole moment for electrons (orange), muons (green) and taus (blue), including the nuclear form factors and LUX efficiency (see text for details). We assume m χ = 1 TeV and ∆ = 0.1 and restrict y χ to lie on the relic surface.
given by where ρ 0 ∼ 0.3 GeV/cm 3 is the local DM density, m N is the target nucleus mass, v min = E r /2(m χ + m N )/(m χ m N ) is the minimal DM velocity required to give a recoil energy E r and f M W (v + v e ) describes the DM velocity distribution in the rest frame of the detector. The particle physics interactions are contained within σ χN , which we now discuss.
Although the dark matter in our models is uncharged, it can interact with a nucleon, and hence the nucleus, at one-loop level. In our model, the dominant contribution to direct detection comes from the one-loop diagrams shown in fig. 3. Note that for Majorana DM there are two additional diagrams with crossed χ legs. The loop diagrams can be mapped onto effective DMphoton interactions, where the most general effective Lagrangian for our interaction is given by where d M and d E are the magnetic and electric dipole moments and A denotes the anapole moment. We see that the dipole operator appears at dimension five whereas the anapole operator is dimension six. For Majorana DM the magnetic and electric dipole moments are identically zero (which ultimately leads to a dramatically reduced rate for the Majorana models). The loop diagrams and their contribution to the dark matter-nucleon scattering cross-section for one coannihilation partner in the n = 1 case have been computed in [40][41][42][43]. Here we briefly summarise the results relevant for our work.
For Majorana dark matter, the one-loop contribution to the anapole moment for the tau and muon models is given by for momentum transfers q much smaller than the lepton mass, |q 2 | m 2 , and where µ = (1 + ∆) 2 and = m 2 /m 2 χ . The factor n accounts for the sum of diagrams when n coannihilation partners are present. Since the momentum transfer is typically larger than the electron mass, we take the limit |q 2 | m 2 for the anapole moment of the electron models, In fig. 4 (left) we show the anapole moment as a function of ∆ for electrons, muons and taus. For the electron case the anapole moment depends on the momentum transfer, which is given by q = √ 2E r m N , and we show the moment for the exemplary value E r = 50 keV. Assuming m ≈ 0, the anapole moment for all models has a log divergence as ∆ → 0, reflecting the fact that all particles in the loop can be on-shell simultaneously (when m φ ≈ m χ ). The anapole moment has a further divergence which is regularised by either the lepton mass or the momentum transfer, which explains why it is larger for smaller lepton masses. Finally, the anapole moment tends to zero as ∆ → ∞, since the coannihilation partners decouple in this limit.
For Dirac DM, the anapole moment is half as large as the moment for Majorana DM, eqs. (10) and (11), for ∈ {e, µ, τ }. Additionally, it generates a dipole moment given by Since we are only interested in the limit where m m χ , the above expression simplifies to which is independent of the lepton mass.
The anapole and dipole moments contribute to the differential DM-nucleus cross-section, dσ Dip where the first term in both expressions corresponds to the spin-independent part (where the DM scatters on the nuclear charge Z), while the second term parameterises the spin-dependent interaction (where the DM scatters on the nuclear magnetic moment, d A ). Here, α is the finestructure constant, J and m N are the spin and the mass of the nucleus, respectively, v is the velocity of the incoming DM particle and E r the recoil energy. Note that eq. (16) and eq. (17) are written for one nucleon isotope only and the spin-dependent and spin-independent parts have to be summed separately over the relevant isotopes. Each summand is weighted by the isotope abundance. In the present paper we focus on Xenon as the target material, with Z = 54 and sum over the isotopes given in [44,45]. For the nuclear charge and spin form factors we use [44] F Z (E r ) = 3 exp −q 2 s 2 /2 sin(qr) − qr cos(qr) (qr) 3 , where q = √ 2E r m N , s = 1 fm, r = √ R 2 − 5s 2 , R = 1.2A 1/3 fm, R s = A 1/3 fm and A is the nuclear mass number.
Turning now to the astrophysical quantities in eq. (8), we assume the standard halo model (an isotropic and isothermal sphere) for the DM distribution, which leads to a Maxwell-Boltzmann velocity distribution in the galactic frame smoothly truncated at the galaxy escape velocity, v esc = 550 km [46], given by where N is a normalisation constant and v 0 = 220 km/s [46] describes the velocity of the sun about the centre of the galaxy. Since the velocity distribution is given in the rest frame of the Milky Way, we use a Galilean transformation to move to the rest frame of the detector, where v e is the velocity of the Earth relative to the galactic centre. For simplicity, we neglect the velocity of the Earth with respect to the sun and take v e = 220 km/s. Note that all terms in eq. (16) and eq. (17) are either independent of the velocity or proportional to 1/v 2 [47]. The integral for velocity dependent terms can be solved analytically (see for instance [46,48]), while the integral with constant terms can be solved analytically, once numerical values for v esc , v 0 and v e have been provided.
We now have all the ingredients necessary to calculate dR/dE r , eq. (8), and we show an example spectrum in fig. 4 (right). The differential rate for the dipole moment dominates over the rate for the anapole moment by ∼ 5 orders of magnitude, reflecting the fact that the dipole operator is We are now ready to compute bounds from existing LUX data and derive projections for future experiments, such as DARWIN, using the statistical procedure outlined in [49]. The current LUX exposure is 3.35 × 10 4 kg · days [25], while DARWIN aims at exposures of 2 and 500 ton · years [28]. We take into account the efficiency of nuclear recoil event detection in LUX given in [25] and assume the same efficiency profile for future experiments. Figure 5 shows current LUX (solid lines) and projected DARWIN (dashed and dotted lines) bounds for Dirac (top-left) and Majorana (top-right and bottom) DM as a function of ∆ and m χ . The line colours correspond to one, three and ten coannihilation partners. As shown in fig. 4 (right), the dipole contribution dominates over the anapole contribution which can thus be neglected for Dirac DM. Since the dependence of the dipole on the lepton mass is negligible, the bounds for Dirac DM are the same for all leptons. Current LUX results (solid lines) exclude Dirac DM up to masses of 0.5 TeV for large values of ∆, but up to 10 TeV in the strong coannihilation region of ∆ 0.02. LUX bounds are slightly more stringent for models with more coannihilation partners. An exposure of 2 ton · years at DARWIN (dashed line) can exclude masses up to 6 TeV for large ∆ and 10 TeV even for ∆ ∼ 0.3. With the nominal exposure of 500 ton · years, this region of parameter space for Dirac DM will be probed completely.
Majorana DM contributes only to the anapole moment, which depends on the lepton mass. The top-right, bottom-left and bottom-right panels show the bounds for Majorana DM coannihilating with an e, µ and τ -type partner, respectively. The bounds are significantly weaker than for Dirac DM, and LUX does not currently constrain any of the parameter space. An exposure of 2 ton·years can access masses up to 1 TeV and ∆ < 0.1 for the model with ten electron-type coannihilation partners, but is barely sensitive to the models with fewer partners or models that are muon or tau type. The bounds are strongest for electrons, since the anapole moment is larger for smaller lepton masses. An exposure of 500 ton · years can exclude Majorana DM masses ∼ 0.5 TeV for all ∆, and up to a maximum of 8 TeV in the coannihilation region ∆ < 0.1.

IV. INDIRECT DETECTION
We now consider the bounds from indirect detection. In regions of high dark matter density such as the galactic centre or dwarf galaxies, χ andχ can annihilate and form pairs of highenergy opposite-sign leptons at tree-level, which may decay. In the process, photons, positrons and anti-protons may be produced, along with other SM particles. At one-loop level, pairs of mono-energetic photons may also be produced. Current and future experiments search for excesses of these particles above the expected astrophysical backgrounds, which may be interpreted as a signal of dark matter. We first focus on continuum photon searches and consider the constraints placed on our model by Fermi-LAT, HESS and the region of parameter space that will be probed by CTA. These will turn out to be the most important constraints on our models.
Indirect search strategies aim to maximise the potential signal and minimise the background. The two main targets commonly considered are the Galactic Centre (GC) and dwarf spheroidal galaxies (dSphs). Dark matter density is highest at the GC, and the signal from the GC is expected to be several orders of magnitude larger than that from dwarf galaxies. However, the GC suffers from both large background sources of gamma rays and significant uncertainty in the local DM density. For this reason, the GC is considered a likely target to provide the first measurement of a DM signal, but subsequent measurement of a signal in dwarf galaxies will usually be needed to make the claim that the signal is unambiguously due to DM. The DM density in the Milky Way is well measured away from the centre, but is poorly known in the inner ∼ 2 kpc. DM models are either cusped, e.g., Einasto profile and Navarro-Frenk-White (NFW) profile, or cored, e.g., Burkett profile. Although N -body simulations suggest a cuspy profile, interactions with baryonic matter could lead to a cored profile (where the dark matter density is constant below ∼ 2 kpc). HESS have developed strategies for both situations [30,50,51], but their sensitivity is much higher if the distribution is cusped. In our GC analysis we assume the profile is cusped. If the distribution is in fact cored, neither CTA nor HESS can place any limits on the models we consider. We also determine weaker but more robust bounds obtained from observing dwarf galaxies.
Annihilating dark matter can produce photons both via direct emission (primary) and by secondary production via Inverse Compton Scattering (ICS) of ± on the ambient photon background. The differential flux of photons from direct emission, in a solid angle ∆Ω, is given by where σv is the thermally averaged dark matter annihilation cross-section, N γ is the photon yield per annihilation and the J-factor integrates the square of the dark matter density along the line of sight over the solid angle ∆Ω. The spectrum of secondary photons produced via ICS can be calculated by convolving the e ± injection spectrum with a halo function for the inverse Compton process [52]. For electrons and muons, a non-negligible photon flux is generated from ICS. However, the precise contribution depends strongly on the assumed halo function and on the spatial region considered in the analysis. For the HESS limits, we use the results in [50] which ignore the contribution from ICS. For the prospective CTA limits we use the results in [52] which include this contribution, strengthening these limits. All of the dependence on the particle physics model is contained within the σv term, the rest being dependent on astrophysical quantities. Since DM in the Milky Way is travelling relatively slowly (v ≈ 10 −3 c), we can take the non-relativistic limit for the annihilation cross-section and assume that all dark matter particles are travelling with the same speed. For the Dirac case of our model, taking the limit m m χ , the annihilation cross-section is which agrees with [53] for n = 1. For the Majorana case we have where v is the relative velocity between the two annihilating dark matter particles, which agrees with [54,55] for n = 1. We see that in the Majorana case, the cross-section is velocity suppressed, as is well known. This suppression means that indirect detection is a poor probe of this case. The strongest indirect detection constraints come from measurements of continuum photons from the galactic centre, assuming a cusped halo profile. To calculate the current limits from HESS, we use the 95 % C.L. upper limits for χχ → τ + τ − and χχ → µ + µ − presented in [50]. 1 To produce their limit, a difference measurement is preformed between a region near the GC and another region further away, and a 2D binned Poisson maximum likelihood analysis is used to distinguish signal from background using both spatial and spectral information. As mentioned above, photons produced by ICS are ignored. To estimate the sensitivities to the electron model, we rescale the muon limit using the integrated flux of prompt photons between 160 GeV (the threshold for HESS) and m χ , using the results presented in [57]. We use the muon limit since the prompt spectrum from muons is closest to that from electrons. To calculate the prospective limit from CTA, we use the 95 % C.L. sensitivity limits on χχ → τ + τ − , χχ → µ + µ − and χχ → e + e − presented in [52]. These limits assume a cuspy Einasto profile and 500 hours of observation, and use a likelihood ratio statistical test to derive the 95% C.L. sensitivity limits. These limits include photons produced via ICS and ignore systematic uncertainties in the datasets and the galactic diffuse emission.
As mentioned above, limits derived from observations of dwarf galaxies are more robust as they do not depend on assumptions about the dark matter halo in the centre of the Milky Way. We consider 95 % C.L. upper limits set by Fermi-LAT on χχ → τ + τ − and χχ → µ + µ − [29], and rescale the muon limit to electrons using the integrated flux of prompt photons. Neither the HESS [58] nor CTA [31] constraints from dwarf galaxies are large enough to place any constraints on our models.
In fig. 6 we show current and prospective limits from HESS and CTA on our models. Only the models with Dirac dark matter are shown as even in the most optimistic scenarios, CTA cannot probe any of the parameter space for Majorana dark matter. This is due to the velocity suppression in the annihilation cross-section, seen in eq. (23). As mentioned above, if a cored dark matter profile is assumed, HESS and CTA place no limits on any of the parameter space of any of the models discussed in this paper.
In fig. 6 (left) we see the constraints on the electron type models for n = 1, 3, 10. HESS excludes the region of low dark matter mass and small ∆. On dimensional grounds, the strongest constraints are always expected at low dark matter masses, since this is the characteristic scale of dark matter annihilation and σv ∼ m −2 χ . For Dirac dark matter, coannihilation leads to a relic surface with large Yukawa couplings at small mass splittings, so stronger constraints are also seen at low ∆. Since the diagram responsible for the indirect detection signal is the same as that responsible for dark matter freeze-out, the significant variation in y χ between the models is mostly cancelled by the n 2 in eqs. (22) and (23). However, it is not completely cancelled due to the interplay between numerator and the denominator in the coannihilation formula, eq. (3). This means that slightly stronger constraints are seen with models with more coannihilation partners. The CTA limits are significantly stronger and probe the electron type models up to m χ ∼ 7 TeV for large ∆, and even higher for ∆ 0.1. In fig. 6 (middle) we show the limits on the muon type model. We see that the results are broadly similar to the electron type models, but the region of parameter space which can be probed is slightly smaller. This is because muons produce fewer primary and ICS photons than electrons [52]. In fig. 6 (right) we show the limits on the tau type model. Here the limits are significantly stronger. HESS can probe the parameter space up to m χ ∼ 4 TeV for any ∆ while CTA probes the whole parameter space. This is due to the large number of primary photons produced by taus. Figure 6 also shows the limits from Fermi-LAT -MAGIC observations of dwarf galaxies on our Dirac dark matter models. The bounds are weak and only constrain the n = 10 electron model (m χ < 0.5 TeV, ∆ < 0.02), the n = 3 tau model (m χ < 0.75 TeV, ∆ < 0.02) and the n = 10 tau model (m χ < 1.5 TeV, ∆ < 0.06). However, these constraints do not depend on any assumptions about the DM halo profile or ICS and are therefore more robust than those from the galactic centre. Future CTA bounds for dwarf galaxies [59] are too weak to place any constraints on any of our models. Furthermore, all the searches we consider are too weak to constrain any of the Majorana DM models. If CTA were strengthened by a factor of ≈ 10, it would begin to probe this parameter space [55].
The lepton pairs produced in the DM annihilation can also lead to primary or secondary positrons and anti-protons. Experimental limits on these final states [60][61][62] extend only to DM masses of ∼ 0.5 TeV. Since we are primarily interested in heavier DM, we do not show these bounds in detail. Although AMS-02 may see an excess in the positron fraction [60], there are many uncertainties in determining the background and a detailed analysis is beyond the scope of this paper. A preliminary study shows, however, that these positron bounds are very weak and only constrain the Dirac models with ten coannihilation partners at m χ < 0.5 TeV and ∆ < 0.01. Furthermore, this model can only produce a signal larger than the astrophysical background, at any positron energy, for m χ < 0.5 TeV or ∆ < 0.01.
Finally, we mention mono-energetic photons. These may be produced via a one-loop diagram, with φ and running in the loop. We find that the loop suppression is too large for Fermi-LAT, HESS or CTA observations [63] to provide any constraints. Internal bremsstrahlung also produces a sharp feature in the gamma ray spectrum. Although this process is not velocity suppressed, even for Majorana dark matter, the cross-section is too small for the experiments we consider to constrain the models, due to phase-space suppression. For a detailed discussion of these gamma-ray constraints, see [55].

V. COLLIDER CONSTRAINTS
It is challenging to search for our dark matter models directly at a hadron collider, since the dark matter is a gauge singlet which only couples to leptons. The coannihilation partner, φ ± i , however is a charged scalar of similar mass. It will be pair produced in the process pp → φ + i φ − i with a subsequent decay of φ ± i to a lepton, , and χ, depicted in fig. 7 (left), where BR(φ ± i → χ ± ) = 1. We focus on final states containing two opposite-sign same-flavour leptons and missing energy. As τ reconstruction at future colliders is particularly challenging to model, we do not provide collider limits for the τ models. However, we can assume that the collider reach on τ models will be somewhat worse than the limits on the models involving electrons and muons. For the tables and figures in this section we focus on the muon type model, and provide the limits for the electron type models in appendix A.
Since we are interested in multi-TeV dark matter, the LHC at 13 TeV only provides weak constraints. E.g., [19] excludes our n = 1 models only for m χ < 0.3 TeV. We therefore present sensitivity projections for the HE-LHC with √ s = 27 TeV assuming an integrated luminosity of 15 ab −1 [20] and for the FCC-hh with √ s = 100 TeV and 20 ab −1 [21]. We estimate the sensitivity Figure 7. The leading order partonic process contributing to pp → φ + φ − →χχ + − (left) and its crosssection at 13 TeV, 27 TeV and 100 TeV including a K-factor of 2 (right).
of future colliders to our models by adapting the analysis used in [64] to search for slepton pair production with subsequent decay to neutralinos and leptons.
The signal pp → φ + φ − is simulated using a custom SARAH v4.12.1 [32] model, we generate the signal and background parton level events using MadGraph5 v2.6.2 [65], simulate the showering using Pythia6 v6.4.28 [66] and perform the detector simulation with Delphes v3.3.3 [67]. For our 27 TeV simulations, we use the default Delphes card. For the simulations at 100 TeV we use the FCC Delphes card implementing the configurations proposed by the FCC working group [68]. For the signal simulation, we adapt the card to treat the DM particle as missing energy. We use the LO partonic production cross-sections and multiply by a generous K-factor of 2, as we want to find the exclusion limits in the optimistic case, fig. 7 (right). To validate our analysis, we reproduce the relevant backgrounds in [64] and find good agreement when using σ W W = 72 pb [64,[69][70][71], σ W Z = 26 pb [64,72], σ ZZ = 9.0 pb [64,73], σ tt = 230 pb [74] and σ W t = 23 pb [75].
The main SM backgrounds to our signal are W W , V V , W V , tt, W t and V +jets, where V = Z, γ. While only W W and V V are irreducible backgrounds, W V , tt and W t contribute if a lepton or one or two b-jets are missed. The V +jets background is important at low values of m T2 , but is negligible above m T2 ≈ 100 GeV. In order to isolate the signal, we impose the following cuts. Two opposite-sign same-flavour light leptons are required with p T > 35 GeV and p T > 20 GeV for leading and subleading leptons, respectively. We veto events with any other leptons, which reduces the W V background. Removing events with m µµ < 20 GeV and |m µµ − m Z | < 10 GeV significantly reduces backgrounds with a Z-boson in the final state. Finally we cut on the transverse mass [76,77], m T2 > 200 GeV , where we use   where p µ + T and p µ − T are the transverse momenta of the leptons, q T is an arbitrary two-vector which represents the unknown transverse momenta of the dark matter particle associated with µ − , and For a process where two particles each decay to a lepton and missing energy, the m T2 distribution will have an end point at the mass of the heavier particle [78]. Although in [64] a cut of m T2 > 90 GeV is used, we increase this to m T2 > 200 GeV. This has a small effect on our signal efficiency, as we are mostly interested in dark matter candidates with mass larger than 200 GeV, while strongly reducing the background from tt, W t. However, even with this large cut, we find a significant background from W W , W V and V V , where at least one of the vector bosons is extremely off-shell.
To include this effect in MadGraph we simulate pp → + − ν allν all and pp → + − all ν, where ν all is ν e , ν µ or ν τ and all is any charged lepton. We do not find a similar large contribution from off-shell particles in the tt and W t channels. Even though the cross-section of these gluon initiated channels grows faster than the di-boson processes as the collider energy is increased, they remain a subdominant background as the t is narrower and as this background only passes the cuts if a jet is missed, reducing the m T2 endpoint. Finally, we checked that the contribution from jets faking muons is negligible. In table II we show the cross-sections at each stage in the analysis for the background and for an example signal, m χ = 0.6 TeV, ∆ = 0.34 (27 TeV) and m χ = 0.8 TeV, ∆ = 0.2 (100 TeV), both for the n = 10 muon type model.
In fig. 8 we show the differential distribution in m T2 for the events passing all cuts, for the background and example signal. We see that µ + µ − ν allν all is the dominant background, and µ + µ − all ν is around an order of magnitude smaller. This is due to both the smaller initial cross-section and the smaller efficiency. We see that both the background and the example signal falls sharply from m T2 = 200 GeV to m T2 ≈ 500 GeV. However, the signal will continue to higher values of m T2 for other points in our parameter space. We also see that at 27 TeV, the µ + µ − ν allν all continues out to higher values of m T2 , while at 100 TeV the situation is reversed.
To estimate the expected exclusion limit, we use a Poisson counting procedure for the signal and background events which pass all the cuts, based on a frequentist framework [1,79]. We use the likelihood ratio λ(0) = L(0)/L(μ), where  Figure 8. The m T2 distribution for background events passing all cuts for the muon model, and an example signal for n = 1, 3, 10, at 27 TeV (left) and 100 TeV (right). We do not use this information in determining the reach, but simply perform a cut-and-count analysis based on these events.
n is the observed number of events, s is the number of signal events, b is the number of background events andμ = (n−b)/s. In the large sample limit, the significance Z 0 is given by Z 0 = −2 ln λ(0). Replacing n with the expected value, s + b, we find that the significance is given by In the limit s b, this reduces simply to s/ √ b. The significance is related to the p-value by Z 0 = Φ −1 (1 − p), where Φ is the standard Gaussian cumulative distribution. The 90% confidence limit is given by Z 0 ≈ 1.64.
In fig. 9 we present the 90% C.L. sensitivity for the muon type models at a 27 TeV and a 100 TeV proton-proton collider. The reach on electron type models is shown in Appendix A. The parameter space probed is where m χ is small and ∆ is relatively large. The reach is independent of whether dark matter is Majorana or Dirac, since it depends on the φ-pair production cross-section and the fact that BR(φ ± i → χ ± ) = 1. The large m χ region is not probed as m φ increases with m χ , and the φ-pair production cross-section decreases rapidly as m φ increases, fig. 7 (right). We see that in both cases the limits are strongest when there are more coannihilation partners. This is because the pp → χχ + − cross-section scales as n 2 . For n = 1, the 27 TeV (100 TeV) machine can probe m χ < 0.75 TeV (1.2 TeV), for n = 3 it can probe m χ < 1.3 TeV (2.3 TeV) while for n = 10 the limits are m χ < 2.0 TeV (4.0 TeV). The small ∆ region is not probed as in this region the momentum of the leptons is small and they are not efficiently reconstructed. This is a well known problem in the coannihilation region. The gap for lower ∆ can be closed, e.g., by looking for ISR [34,35] or for disappearing charged tracks [36]. A thorough study of the reach in the challenging small ∆ region is not pursued here as in these models there is a nice complementarity with direct detection experiments, which can be seen to cover this region.
We also overlay the direct and indirect detection bounds from sections III and IV, to give a summary of all the relevant current and future experimental constraints. We see that the situation is dramatically different for Dirac and Majorana χ. For Dirac χ, small masses and mass splittings Figure 9. The reach of future colliders, at 90% C.L., current and future direct detection experiments, at 90% C.L., and current and future indirect detection experiments, at 95% C.L., for Dirac (left) and Majorana (right) DM interacting with a muon and one (top), three (middle) and ten (bottom) coannihilation partners in the ∆-m χ plane. The lightly, moderately and strongly shaded regions correspond to the direct detection limits by the future DARWIN experiment with 500 ton · years, 2 ton · years and the LUX limits, respectively, which are discussed in detail in section III. The circle and cross signify our example signals. have already been excluded by LUX. In the future, DARWIN will probe the full parameter space, while colliders and indirect detection will be sensitive for relatively low masses and large or small ∆, respectively. We see that the challenging small ∆ region at colliders is excluded by the existing bound from LUX.
For Majorana χ, on the other hand, DARWIN, with the maximum exposure, is limited to small masses and small ∆, while indirect constraints do not feature. This is due to the velocity suppression of both the DM-nucleus and the annihilation cross-sections. The collider bounds are the same as in the Dirac case, since the mass term of χ does not enter into the production and decay of φ-pairs. In this case, future colliders are essential for probing the large ∆ region of the parameter space.

VI. CONCLUSIONS
As nature has not yet provided clues to the mass scale of dark matter, we have to explore many orders of magnitude in mass and coupling strengths. One of the first and best studied regions in parameter space is 10 GeV -1 TeV dark matter, with a large fraction of the experimental progress of the past decade targeting this mass range. Although many of the viable WIMP models in this mass range can be exhaustively probed at the next generation of direct and indirect detection and collider experiments, we point out that large regions of viable parameter space are inaccessible when coannihilation is taken into account.
Assuming a minimal and versatile dark matter model, we restricted to the relic surface and studied the reach of current and future experiments on the dark matter parameter space. We take the dark matter particle to be a gauge-singlet Majorana or Dirac fermion and introduce n charged scalars as coannihilation partners, which together couple to SM right-handed leptons. The relic surfaces of these models demonstrate a viable, perturbative multi-TeV dark matter candidate, whose relic abundance either decreases (in the Majorana case) or increases (in the Dirac case) with coannihilation.
Direct detection experiments are sensitive to our models via loop-induced dipole and anapole interactions. We compute the interaction cross-section between dark matter and target nuclei, taking into account nuclear form factors. By explicitly integrating over the DM velocity distribution and including experimental efficiencies we calculate the expected event rate, which we compare to bounds from the current LUX experiment and use to derive projected limits from DARWIN. We find starkly different result for Dirac and Majorana dark matter. Over the parameter space we consider, LUX excludes the Dirac models at m χ < 0.3 TeV (for any ∆) and ∆ < 0.02 (for any m χ ). With 2 ton · years of exposure these constraints strengthen to cover around 95% of the parameter space. DARWIN, with an exposure of 500 ton · years, can probe the entire region for Dirac dark matter. The models with Majorana dark matter are, however, currently unconstrained, due to velocity suppression in the direct detection cross-section. Only DARWIN with an exposure of 500 ton · years can make progress, with the ability to exclude Majorana DM up to 0.5 TeV for large ∆ and between 2 TeV and 8 TeV in the region ∆ < 0.1. The remaining ∼ 90% of the parameter space, however, remains unconstrained. In all cases, increasing coannihilation (both by reducing ∆ and increasing the number of partners) tends to make direct detection a better probe of the models.
In indirect detection experiments, we only find significant limits via continuum photons if the dark matter is Dirac, and if the Milky Way halo profile is cusped. The models coupled to taus have the best limits, since many continuum photons are produced when the tau decays. If the Milky Way halo profile is cusped, HESS has excluded the tau models up to 4 TeV for all values of ∆, and at least up to 10 TeV for n = 10 and ∆ < 0.02. Again if it is cusped, CTA will be able to probe our whole parameter space for Dirac dark matter which couples to taus via n ≤ 10 scalars. For the muon (electron) case, even if the Milky Way profile is cusped, CTA can only probe the Dirac models to m χ < 5 TeV (m χ < 7 TeV), unless ∆ < 0.15 in which case it can probe up to at least 10 TeV. If the Milky Way halo profile is not cusped then HESS and CTA provide no bounds and the best limits come from Fermi-LAT -MAGIC observations of dwarf spheroidals, although they are relatively weak, only covering a region with m χ < 2 TeV and ∆ < 0.05. Finally, due to velocity suppression, there are no indirect detection constraints on any of the Majorana models.
The collider bounds we find are insensitive to the difference between Majorana and Dirac DM, as we consider pair production of the coannihilation partner, which subsequently decays into two same-flavour opposite-sign leptons and missing energy with a branching ratio of 1. We only provide bounds for the models where dark matter couples to electrons and muons since tau reconstruction is difficult to model at a future collider. In any case, the exclusion reach for the tau models will be worse than for the electron or muon models. We simulate signal and background, apply a set of cuts to isolate the signal and derive the expected reach of the HE-LHC with √ s = 27 TeV and an integrated luminosity of 15 ab −1 and a future collider with √ s = 100 TeV and 20 ab −1 . We find that, as the production cross-section is proportional to n 2 , the bounds strengthen as the number of coannihilation partners increases. For the muon model, a 27 TeV (100 TeV) proton-proton collider can exclude the n = 1 model up to 0.74 TeV (1.2 TeV), the n = 3 model up to 1.3 TeV (2.2 TeV) and the n = 10 model up to 2 TeV (4 TeV). Although the analysis does not target the more challenging low ∆ region, small ∆ searches are unlikely to have a larger mass reach and this region is covered by current and future direct detection experiments.
We thus see that while viable, perturbative models of Dirac dark matter will be well probed by DARWIN, future colliders, and, if the Milky Way Halo is cuspy, by CTA, the future suite of dark matter searches will leave around 70% of the parameter space of our viable models of Majorana dark matter untouched. We are limited in one direction by velocity suppressed interactions, and in another direction by the large mass scales involved. For direct detection, overcoming these limitations could involve the optimisation of current experiments for larger DM masses or a greater emphasis on electron recoils (for probing the electron type models at tree level), by developing new target materials, e.g, [80], the application of novel techniques, e.g., [81], or scenarios where the velocity suppression is lifted, e.g., due to infall into neutron stars [82]. The situation in indirect detection would be considerably improved with a better determination of the dark matter density in the galactic centre, and by analyses which do not depend on a cuspy profile, e.g., [30]. However, key to probing the velocity/loop/phase-space suppressed Majorana models is substantially improving bounds from gamma-ray observations [55]. We see that future colliders will be able to produce a handful of 10 TeV particles in optimistic scenarios, but the analysis considered here retains thousands of background events. Improved search strategies, which choose varying cuts depending on the new physics parameter point under consideration, or which optimise the signal against the background using a wider spectrum of information (for instance, using neural networks) could push the sensitivity to larger masses, even though it is intrinsically limited by the steeply falling cross-section with mass. Exploring displaced vertex analyses and using ISR will be useful in the coannihilation region. We see that, to conclusively test the WIMP paradigm, new experimental techniques will need to be developed to surmount these challenges.

VII. ACKNOWLEDGMENTS
We would like to thank Tilman Plehn and Joachim Kopp for collaboration in the early stages of the project, Tim Tait for discussions on coannihilation, Joachim Kopp for advice on direct detection, Alaettin Serhan Mete from ATLAS and Rakhi Mahbubani for useful discussions on the collider analysis, and Paolo Panci and Lucia Rinchiuso for helpful advice on indirect detection. MJB was supported by the German Research Foundation (DFG) under Grant Nos. KO 4820/1-1 and FOR 2239, by the European Research Council (ERC) under the European Union's Horizon Figure 10. The reach of future colliders, at 90% C.L., current and future direct detection experiments, at 90% C.L., and current and future indirect detection experiments, at 95% C.L., for Dirac (left) and Majorana (right) DM interacting with an electron and one (top), three (middle) and ten (bottom) coannihilation partners in the ∆-m χ plane. The lightly, moderately and strongly shaded regions correspond to the direct detection limits by the future DARWIN experiment with 500 ton · years, 2 ton · years and the LUX limits, respectively, which are discussed in detail in section III.