Is natural higgsino-only dark matter excluded?

The requirement of electroweak naturalness in supersymmetric (SUSY) models of particle physics necessitates light higgsinos not too far from the weak scale characterized by mweak∼mW,Z,h∼100GeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{weak}\sim m_{W,Z,h}\sim 100\hbox { GeV}$$\end{document}. On the other hand, LHC Higgs mass measurements and sparticle mass limits point to a SUSY breaking scale in the multi-TeV regime. Under such conditions, the lightest SUSY particle is expected to be a mainly higgsino-like neutralino with non-negligible gaugino components (required by naturalness). The computed thermal WIMP abundance in natural SUSY models is then found to be typically a factor 5–20 below its measured value. To gain concordance with observations, either an additional DM particle (the axion is a well-motivated possibility) must be present or additional non-thermal mechanisms must augment the neutralino abundance. We compare present direct and indirect WIMP detection limits to three natural SUSY models based on gravity-, anomaly- and mirage-mediation. We show that the case of natural higgsino-only dark matter where non-thermal production mechanisms augment its relic density, is essentially excluded by a combination of direct detection constraints from PandaX-II, LUX and Xenon-1t experiments, and by bounds from Fermi-LAT/MAGIC observations of gamma rays from dwarf spheroidal galaxies.

1. The well-tempered neutralino (WTN) [17][18][19][20], wherein the bino and higgsino components were adjusted to comparable values so as to obtain the required relic density, predicted σ S I ( Z 1 p) ∼ 10 −8 pb relatively independently of m Z 1 . The nucleon-WIMP cross section is roughly independent of the WIMP mass because, for heavier WIMPs, the higgsino component of the WIMP needs to be increased to maintain the observed relic density; this increased higgsino component then maintains the direct detection cross-section at roughly a constant value. The light higgsino region is typical of the so-called focus point region/hyperbolic branch [21] of the mSUGRA/CMSSM model [22] and is now solidly excluded [12][13][14][15]23]. 2. The case of wino-like WIMP-only dark matter, which is characteristic of anomaly-mediated SUSY breaking models, predicts rather large rates for WIMP-WIMP annihilation into W W , leading to gamma ray production in areas of the universe where increased WIMP densities are expected (such as galactic cores and dwarf galaxies). Recent limits from Fermi-LAT (at lower m Z 1 ) and HESS (at m Z 1 ∼ TeV-scale) have seemingly excluded this possibility if one includes Sommerfeld enhancement effects in the annihilation cross sections [23][24][25].
Taken all together, the data seem to suggest that weak scale supersymmetry (WSS) [26], if viable, must have at least strongly coupled superpartners with soft SUSY breaking parameters m so f t in the multi-TeV range rather than at the weak scale, m weak ∼ m W,Z ,h ∼ 100 GeV. The confrontation of theory with data then seemingly exacerbates what has become known as the Little Hierarchy problem: why is m weak m so f t ? While the introduction of SUSY can solve the Big Hierarchy problem, avoiding the Higgs mass from blowing up to the Planck scale while avoiding extreme fine-tuning of parameters, now one may expect the Higgs boson mass to inflate to the multi-TeV regime if the heavy superpartners couple directly to the Higgs fields, absent again fine-tuning of SUSY Lagrangian parameters.
The well-known expression for the Z -boson mass obtained from the minimization of the (one-loop) scalar potential of the Higgs fields, serves as a starting point for many discussions of fine-tuning in SUSY models. The last (approximate) equality in Eq. (1) obtains for moderate to large values of tan β required by the measured value of m h = 125 GeV. Here, m 2 H u is the weak scale value of the up-Higgs squared soft mass and μ is the Higgs/higgsino mass parameter occuring in the (SUSY conserving) superpotential. The u u and d d terms contain an assortment of radiative corrections, the largest of which typically arise from the top squarks: x W ) and x W ≡ sin 2 θ W . Also, F(m 2 ) = m 2 log m 2 Q 2 − 1 and g 2 z = (g 2 + g 2 )/8 and f t is the top-quark Yukawa coupling. Expressions for the remaining u u and d d are given in the Appendix of Ref. [27].
Requiring no large unexplained cancellations between the various terms on the right-hand-side of Eq. (1) led us to introduce the electroweak fine-tuning measure EW [27,28] defined as the ratio of the magnitude of the maximal contribution on the right-hand-side (RHS) of Eq. (1) to m 2 Z /2. If the RHS terms in Eq. (1) are individually comparable to m 2 Z /2, then no unnatural fine-tuning is required to generate m Z = 91.2 GeV. We advocate the use of EW for discussions of fine-tuning because it allows for the possibility that model parameters traditionally regarded as independent may turn out to be correlated by the underlying SUSY breaking mechanism, and further, that the most commonly used fine-tuning measure, BG ≡ max i are fundamental parameters of the theory [29,30] reduces to EW [31,32] after appropriate correlations are incorporated. Ignoring the possibility that model parameters (taken to be independent) might turn out to be correlated may lead to prematurely discarding perfectly viable SUSY models.
The most important implications of low electroweak finetuning (which we take to be EW 30) 2 are the following.

m 2
H u is driven radiatively from its high scale value to small negative values, comparable to − m 2 Z , at the weak scale [27,28]. 3. The top squark contributions to the radiative corrections u u (t 1,2 ) are minimized for TeV-scale highly mixed top squarks [28]. This latter condition also lifts the Higgs mass to m h ∼ 125 GeV. For EW 30, the lighter top squarks are bounded by mt 1 3 TeV [27,33]. 4. The gluino mass, which feeds into the stop masses at one-loop and hence into the scalar potential at two-loop order, is bounded by mg 6 TeV [27,33].
We will collectively call SUSY models for which EW < 30 natural SUSY models. In the present paper, we examine expectations for SUSY WIMP dark matter from three different, well-motivated classes of natural SUSY models that lead to qualitatively different patterns of gaugino and higgsino masses which in turn determines the nature of the SUSY WIMP.
1. Gravity-mediated SUSY breaking models, as exemplified by the two-extra-parameter natural non-universal Higgs model [37][38][39][40][41][42] (nNUHM2) with parameter space given by 4 For NUHM2, because of the gaugino mass unification assumption, one expects weak scale gaugino masses in the ratio M 1 : M 2 : M 3 ∼ 1 : 2 : 7. 2. A phenomenological generalization of the well-studied anomaly-mediated SUSY breaking model, the natural (generalized) anomaly-mediated SUSY breaking model [44] (nAMSB) with parameter space given by, where in addition to AMSB contributions to soft terms [45][46][47][48], we introduce several bulk induced soft terms to render sleptons non-tachyonic (m 0 (bulk)) and to render the model natural (m H u (bulk), m H d (bulk) and A 0 (bulk)) whilst respecting LHC data. As in NUHM2, we trade the high scale parameter freedom of m H u (bulk) and m H d (bulk) for the more convenient weak scale parameters μ and m A . For nAMSB, one expects weak scale gaugino masses in the ratio M 1 :M 2 :M 3 ∼ 3:1:8 but now with μ < M(gauginos) so that a higgsino-like neutralino (mixed with some wino component) is the lightest SUSY particle (LSP) instead of the neutral wino.

Natural generalized mirage mediation model (nGMM)
where gravity-and anomaly-mediated contributions to soft SUSY breaking terms are comparable. The nGMM mass pattern is expected to emerge from several wellmotivated superstring models [49,50]. The parameter space is given by [51,52], where α parametrizes the relative gravity-to anomalymediation and c m and c m3 co-efficients are a continuous generalization of formerly discrete parameters involving modular weights of the relevant fields and a 3 is a 4 Models such as mSUGRA/CMSSM where all soft scalar masses are set to m 0 are no longer natural for m h = 125 GeV while respecting LHC sparticle mass limits [31,43]. Historically, the common mass m 0 originated in the assumption of a flat Kähler potential in minimal supergravity models. The CMSSM requirement that m Hu = m H d = m 0 appears ad-hoc and artificial given that the Higgs multiplets belong to different GUT representations from the matter scalars. The NUHM2 model rectifies the artificial degeneracy requirement and then allows for naturalness whilst respecting LHC Higgs mass and sparticle search constraints. Indeed the NUHM3 generalization where third generation scalars are treated differently from scalars of the first two generations has also been examined in the literature.
continuous generalization of the formerly discrete trilinear gravity-mediated A term. For nGMM models, one expects the weak scale gaugino masses with M 1 < M 2 < M 3 but with compressed spectra (depending on the value of α since the scale of mirage unification μ mir = m GU T e −8π 2 /α ) since they appear to unify at some intermediate scale rather than m GU T 2 × 10 16 GeV.

Dark matter relic density in natural SUSY
For natural SUSY models, we see from Eq. (1) that naturalness requires |μ| and − m 2 H u to be not too far above m 2 Z , but imposes only loop-suppressed restrictions on other soft SUSY breaking parameters. Hence, one expects the LSP to be dominantly higgsino-like, but with a non-negligible gaugino component (lest u u (m W 2 ) becomes large for too large wino masses). The first question then is: do the natural SUSY models produce the measured relic abundance of dark matter in the universe given by DM h 2 ≡ ρ DM ρ c h 2 where ρ c is the critical closure density of dark matter and h is the scaled Hubble parameter. Of course, since higgsinos annihilate with full gauge strength in the early Universe, we do not expect that the relic density of thermally produced, light higgsinos to saturate the observed relic density, but it is nonetheless instructive to examine the expectations for the thermal relic density in well-motivated natural SUSY models.
To answer this question, we next compute the thermallyproduced relic density for the various SUSY models introduced in Sect. 1. We use the computer code Isajet 7.88 to compute sparticle mass spectra for the nNUHM2, nAMSB and nGMM models [53,54]. For nNUHM2 and nAMSB models, we have performed a broad random scan as well as an additional focused scan (over the parameter ranges shown in parenthesis below) in an attempt to further zero in on the natural SUSY region of the parameter space, while for the nGMM model, our scan is already quite focussed. For the NUHM2 model we scan over the parameter range: For the nGMM model, we scan over For each solution, we require the light Higgs boson m h : 122-128 GeV (allowing for ± 3 GeV error in the Isajet m h calculation). 5 To enforce naturalness, we require of each solution EW < 30. We also require mg > 2 TeV and mt 1 > 1 TeV in accord with LHC sparticle search limits.
The results of our calculations of the thermal LSP relic density T P Z 1 h 2 (using the Isajet subcode IsaReD [55]) are shown versus m Z 1 in Fig. 1 for the three natural SUSY models. We plot points from our scan that yield EW ≤ 30 and also satisfy the Higgs boson mass and LHC sparticle mass constraints as blue pluses (nGMM model), green stars (nAMSB model) and yellow crosses (nNUHM2 model). We see first that m Z 1 is bounded from below by m Z 1 100 GeV due to LEP2 limits on m W 1 100 GeV (which we set as the lower limit on the μ parameter scan). Also, m Z 1 is bounded from above by m Z 1 350 GeV from the naturalness constraint, EW < 30. For the lower range of m Z 1 values, then T P Z 1 h 2 is typically a factor ∼ 20 below the measured value CDM = 0.1199 ± 0.0022 [56] while for the high range of m Z 1 then the calculated relic abundance is about a factor ∼ 4 below the measured result. The range of under-abundance just mentioned applies to all three models with the possible exception of nAMSB where some of the green stars lie at even lower T P Z 1 h 2 values. The reason for this is that in nAMSB models, for a lower range of m 3/2 values then the wino can range down to M 2 : 200-300 GeV so that for this model the Z 1 can be mixed higgsino-wino variety: then the neutralino annihilation rate in the universe is enhanced even beyond the higgsino-like case leading to even lower relic density. Thus, natural SUSY models typically predict an under-abundance of thermally produced neutralinos in standard Big Bang cosmology by a factor ∼ 5-25. Other mechanisms are required to bring the expected DM abundance into accord with data. 5 In our previous studies of naturalness we had used a ±2 GeV window on m h . We have checked that our conclusions are quite insensitive to this wider and more conservative window. Two well-motivated classes of mechanisms have been proposed to bring the thermally-produced under-abundance of neutralinos into accord with the measured dark matter abundance. In the first class, the dark matter is multi-component with thermal higgsinos comprising only a fraction of the observed dark matter, with the remainder consisting of other particle(s). The axion is perhaps the best-motivated candidates for the remainder of the dark matter (for a review, see e.g. Ref. [57]). In the second class of models, the dark matter is all neutralinos, with a non-thermal component from late decays (to neutralinos) of heavy particles making up the balance of the observed relic density. We will see below that if the neutralino is dominantly the higgsino of natural SUSY, the second class of models is essentially ruled out by the data.

Mixed axion/WIMP dark matter
As mentioned, one possibility is that the total WIMP abundance does not saturate the measured relic density but that, like visible matter, the dark matter is comprised of several particles. A very natural choice for a second dark matter particle is the QCD axion which also seems to be required to solve the strong CP problem in QCD. In a supersymmetric context, then the axion should occur as but one element of an axion superfield which would also necessarily contain a spin-0 R-parity even saxion field s and a spin-1 2 R-parityodd axino fieldã. Both saxion and axino are expected to gain masses of order the gravitino mass m 3/2 in supergravity models [58,59].
In SUSY axion models, the axions can be produced nonthermally via (1) vacuum misalignment, (2) thermally, and also (3) non-thermally via (late time) saxion decay s → aa. The latter two may lead to relativistic axions whose population is limited by strict bounds on the effective number of relativistic degrees of freedom N eff = 3.15 ± 0.23 derived from fits to CMB and other cosmological data [56]. Axinos can be thermally produced in the early universe and then aug-ment the WIMP abundance via decays after thermal WIMP freeze-out. Saxions can be produced both thermally and nonthermally and then decay to SM particles (resulting in entropy dilution of all relics from their value at the time of decay), SUSY particles (which augment the WIMP abundance) or to axions as mentioned above. WIMPs can be produced thermally or non-thermally via axino, saxion or gravitino decay. The resultant mixed axion-WIMP abundance has been evaluated by solving eight-coupled Boltzmann equations [60][61][62][63]. The Boltzmann equations track the interrelated abundances of • thermally and non-thermally produced WIMPs, • thermally and decay-produced axions, • axions from vacuum mis-alignment/bosonic coherent motion (BCM), • thermal production and decay of axinos, • thermal production and decay of saxions, • saxion BCM production and decay, • thermal gravitino production and decay and • production of radiation at re-heat and from saxion/axino decay.
The exact rates also depend on the underlying SUSY axion model assumed (KSVZ or DFSZ), as well as on other parameters such as mã, m s , θ s , m 3/2 , the re-heat temperature T R , the initial axion mis-alignment andgle θ i and the SUSY particle mass spectrum (which influences the saxion, axino and gravitino decay branching fractions) [60][61][62][63]. Points in parameter space may become excluded via overproduction of WIMPs or axions, or by increasing N e f f via relativistic axion production from saxion decays or by violation of BBN constraints. For low values of the axion decay constant, 6 f a 10 11 GeV, the WIMP abundance is its thermal value since axinos and saxions tend to decay before WIMP freezeout so that ξ = ξ T P ≡ T P Z 1 h 2 /0.12. If f a 10 11 GeV, then post-freeze out saxion and axino decays may augment the WIMP abundance so that ξ T P < ξ < 1. For very large f a 10 14 GeV, then almost always WIMPs are overproduced via saxion and axino decays (ξ > 1), N e f f becomes too large and BBN constraints on late-decaying neutral relics are violated.
For the non-excluded points ξ ≤ 1, the upshot is that the expected rates for direct and indirect WIMP detection now depend on the fractional WIMP abundance denoted by ξ = Z 1 h 2 /0.12 < 1 since now there are fewer target WIMPs compared to the WIMP-only hypothesis for 6 The axion decay constant f a is defined via its coupling to two gluons: Aμν where f a can range from ∼ 10 9 GeV (from SN1987A energy loss rate) up to possibly beyond the Planck scale (for tiny initial mis-alignment angle θ i ). Here, N DW = 6 is the domain wall number for the DFSZ axion. dark matter. For spin-independent (SI), spin-dependent (SD) detection rates, and also the neutrino detection rate at Ice-Cube, the target event rates must be scaled by a factor 7 ξ while for indirect WIMP detection (IDD) via WIMP-WIMP annihilation into gamma-rays or particle-antiparticle pairs, the event rates must be scaled by a factor ξ 2 . To be conservative, for mixed axion/WIMP dark matter, we will assume ξ = ξ T P = T P Z 1 h 2 /0.12 which is usually the lower bound on ξ . For special cases at high f a , bosonic collective motion (BCM) produces a large saxion abundance in the early universe. If parameters are adjusted properly (the saa coupling is tiny or zero to avoid relativistic axion production and m s < 2m Z 1 so s decays only to SM particles) then it is possible to have large entropy dilution of all relics [64] and even lower ξ values; this seems rather contrived, and we will ignore this possibility in this paper.

Non-thermally produced WIMP-only dark matter
Another option is to assume WIMP-only dark matter where the additional WIMP abundance is assumed to arise from non-thermal processes. The prototypical non-thermal WIMP production process occurs from light modulus field φ production in the early universe via the BCM (which also occurs for saxion and cold axion production). If the modulus field (of mass m φ ) then decays after WIMP freeze-out but before the onset of BBN, then it may augment the thermally-produced abundance to gain accord with the measured density of dark matter. This mechanism was originally suggested by Moroi and Randall [65] to account for how wino-like LSPs from AMSB models could account for the observed dark matter. It was later emphasized by Gondolo and Gelmini [66] that the measured relic density could be achieved for any value of T P Z 1 h 2 > 10 −5 (100 GeV/m Z 1 ) by adjusting just two parameters: b/m φ and T R2 where b is the number of neutralinos produced per φ decay and T R2 is the (second) reheat temperature arising from φ decay. This reheating temperature is related to the φ field energy density as T R2 ∼ ρ −1/4 φ . Nonthermal WIMP production has also been recently invoked to reconcile an underproduced WIMP relic density with measured value in string-motivated models with a wino-like LSP [67][68][69][70]. For the case of natural WIMP-only dark matter, we will assume the thermal and non-thermal relic density contributions sum to the measured dark matter density so that ξ = 1 for this case.

Direct WIMP detection bounds
In Fig. 2, we show the value of ξσ S I ( Z 1 p) vs. m Z 1 for a) the case with ξ = T P Z 1 h 2 /0.12 < 1 (corresponding to mixed axion/WIMP DM with no non-thermal WIMP production or dilution) while in frame b) we show the case with natural WIMP-only DM and ξ = 1. We use the Isajet subcode IsaReS [22] for our direct and indirect relic scattering calculations. 8 In both frames, we also plot the current SI DD bounds from LUX, PandaX and Xe-1ton (solid curves), along with a future projected bound from Xe-1ton (dashed). From frame a), we see that present bounds already exclude many natural SUSY model points even with ξ < 1, if we assume that the neutralino relic density is given by its thermal value. Especially, a large fraction of nAMSB model points are excluded. This is because in nAMSB the winos can be relatively light compared to mg and the h Z 1 Z 1 coupling occurs as a product of gaugino times higgsino components (see Eq. (8.117) of Ref. [26]). The enhanced Z 1 p scattering rate for nAMSB more than compensates for the somewhat diminished relic abundance. For the nNUHM2 and nGMM models, the major portion of model points survive the current SI DD bounds. But future ton-scale noble liquid search experiments will cover the remainder of parameter space, assuming that the neutralino relic density is not diluted from its thermal value by entropy injection in the early Universe.
In frame (b), for WIMP-only DM with ξ = 1, then we see that current bounds exclude almost every point of all three models. A single point from the scan with m Z 1 ∼ 250 GeV has survived. The surviving point lies within the future reach of ton-scale noble liquid detectors. Thus, it appears from this plot alone that natural WIMP-only DM appears to be essentially excluded (but for one nNUHM2 point which, we have checked, has gaugino masses close to their naturalness upper limit, and hence a reduced gaugino content and correspondingly reduced neutralino coupling to h). 9 We also show in frame (b) the latest Xe-1ton bound with an added factor of two uncertainty in the experimental bound. In this case, one additional point with m Z 1 ∼ 205 GeV could barely be allowed as it is just inside the limit band.
In Fig. 3, we show ξσ S D ( Z 1 p) vs. m Z 1 . Again, in frame (a) we take ξ = T P Z 1 h 2 /0.12 < 1 while in (b) we show 8 The IsaReS SI direct detection cross sections depend sensitively on the strange quark content of the proton [71,72]. For IsaReS, we use the central values of updated quark mass fractions and moments as tabulated by Hisano et al. [73]. 9 The outlier point with m Z1 250 GeV was generated with mg = 6.2 TeV and EW = 29 and m h = 122.5 GeV. Thus, it inhabits the outermost extremity of the naturalness and Higgs mass allowed regime. ( Z 1 p) vs. m Z1 plane from a scan over the natural NUHM2, nGMM and nAMSB model parameter space for a ξ < 1, assuming the neutralino relic density is given by its thermal value, and b ξ = 1 the natural WIMP-only case with ξ = 1. We also show the current SD limits from the PICO-60 experiment [74] and from IceCube [75] (the latter assuming dominant WIMP annihilation within the solar core into W W final states). From frame (a), we see that, save for a few points around m Z 1 100 GeV, all points avoid the present SD DD bounds. We also see that the bulk of natural SUSY points will be probed by PICO-500 [76] (subject to the caveats mentioned above) although some points might still elude SD detection.

Fig. 2 Plot of points in the σ SI
In frame (b), we show the ξ = 1 case for natural WIMPonly DM. In this case, we see that a combination of PICO-60 and IceCube have already ruled out a significant fraction of natural SUSY model points. The projected reach of PICO-500 should probe the remaining possibilities.

Indirect WIMP detection bounds
In Fig. 4, we show the quantity ξ 2 σ v , the thermally averaged WIMP-WIMP annihilation cross section times velocity, evaluated as v → 0, scaled by the square of the depleted relic abundance, vs. m Z 1 . In this figure, the mixed axion/WIMP dark matter points with ξ 1 (lower set of points), again assuming the thermal neutralino relic density is close to its real value, are neatly separated from the ξ = 1 points  ( Z 1 p) vs. m Z1 plane from scans over the parameter space of the the natural NUHM2, nGMM and nAMSB models for a ξ < 1, assuming the neutralino relic density is given by its thermal value, and b ξ = 1 for WIMP-only dark matter (upper set of points). We also show the present bounds from the combined Fermi-LAT and MAGIC collaborations derived from observations of gamma rays from dwarf spheroidal galaxies. 10 We have added a possible "factor of two" uncertainty in the experimental limit so that it can also be interpreted as a limit band. Corresponding limits from HESS are relevant only for higher, unnatural values of m Z 1 , and not shown in the figure. We see that all of the mixed axion/WIMP dark matter points fall well below the experimental bounds. However, we also see that all the natural WIMP-only points with ξ = 1 points are excluded by present bounds save for a few points with m Z 1 > 300 GeV. If we instead use the limit band, then points with m Z 1 250 GeV are still allowed. We have checked that the m Z 1 > 250 GeV points are excluded by the SI DD band from Fig. 2b). Likewise, we have checked that the two nNUHM2 point with m Z 1 ∼ 200, 250 GeV are excluded by the IDD limit band with ξ = 1. 10 For our case of mainly higgsino-like WIMPs, we have checked that the WIMP-WIMP annihilation takes place almost entirely into the W W and Z Z channels along with a smaller component into Zh. Thus, the Fermi-LAT/MAGIC channel to be compared against is their result for annihilation ito W W since the gammas come primarily from V → qq → π 0 s → γ γ and these configurations are similar for V = W or V = Z . Fig. 4 The scaled values ξ 2 σ v from scans over the parameter space of the the natural NUHM2, nGMM and nAMSB models for ξ < 1, assuming the neutralino relic density is given by its thermal value (lower set), and ξ = 1 (upper set). We plot the Fermi-LAT + MAGIC bound including the central value along with a possible factor of two uncertainty

Concluding remarks
In this paper we have examined the direct-and indirect-WIMP detection rates 11 for three different natural SUSY models with very different gaugino spectra: nNUHM2, nAMSB and nGMM. The three models all have higgsino-like LSPs but qualitatively different and non-negligible gaugino components. They have suppressed values of thermally produced neutralino relic abundances -lower than the measured abundance of CDM by factors ranging from 5-25. For the three models, we have examined their WIMP SI-and SDdirect detection rates and also their indirect detection rates for two different possibilities: (1) mixed axion-WIMP dark matter where only a fraction ξ , determined by the thermal neutralino relic abundance, is assumed to be due to WIMPs, while the remainder is axions, and (2) the case of WIMPonly dark matter where the thermal relic abundance is supplemented by non-thermal production from processes like modulus field decay in the early universe. In this second case, then we take the fractional WIMP abundance ξ = 1.
From our scans of the parameter space of natural SUSY models, we find that models where the WIMP relic density (taken to be its thermal value) forms just ∼ 5-20% of the measured CDM density comfortably survive constraints from LHC as well as those from direct and indirect searches. Direct searches at ton-sized detectors (Xenon-nT or LZ) will probe the entire natural SUSY parameter space, assuming that the relic abundance is given by its thermal expectation. In this case, future experiments such as PICO-500 -designed to measure the spin-dependent neutralino-nucleon scattering -will also probe a large part (but not all) of the parameter space. Otherwise, future colliders such as an electronpositron collider with √ s ≥ 500-600 GeV [77], or a high energy pp collider operating at √ s ∼ 27-33 TeV [83] will be necessary for a definitive probe of the natural SUSY scenario with multi-component dark matter.
The situation for natural SUSY models where the neutral higgsino-like WIMP saturates the observed relic density is qualitatively different. These scenarios are essentially excluded both by bounds from direct detection experiments as well as by independent bounds from Fermi-Lat + Magic observations of high energy gamma rays from dwarf galaxies. More correctly, while a few points from our scans survive the indirect searches, these are excluded by direct detection, and vice-versa. Such models would also be decisively probed by spin-dependent direct-detection at PICO-500.
Thus, the answer to the question posed in the title is: yes, it appears the case of natural higgsino-like-WIMPonly dark matter is indeed excluded. Unnatural higgsinolike WIMP dark matter can still survive as detailed in Refs. [84,85] although these models would have a difficult time explaining why it is that the weak scale is a mere 100 GeV instead of lying in the multi-TeV range. Another possibility is to have models with non-universal gaugino masses where M 3 > 2 TeV to satisfy LHC gluino mass bounds but where M 1 ∼ 50-150 GeV with |M 1 | < |μ|. This case, explored with running non-universal gaugino masses in Ref. [86] and in the pMSSM context in the first of Ref. [87] (for other studies with non-universal gaugino masses, see also [88][89][90][91][92]), has a mainly bino-like LSP while still satisfying naturalness bounds. It is unclear as to the origin of the rather large mass gap between bino and gluino.
As a whole, our results seem to bolster the case for a second dark matter particle such as the axion. While the remainder of the dark matter could be in the hidden sector, the axion is a very well motivated candidate which may well constitute the bulk of dark matter in our Universe. Prospects for the complementary axion searches in SUSY axion models have been examined in Ref. [93].