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 m(weak)~ m(W,Z,h)~100 GeV. 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.


Introduction
Supersymmetric models of particle physics have been under assault from both collider search experiments and direct and indirect dark matter detection experiments. From the CERN LHC, the measured value of the Higgs boson mass m h 125 GeV [1] seems to require TeV-scale highly mixed top squarks, at least in the framework of the Minimal Supersymmetric Standard Model, or MSSM [2,3,4,5]. Direct searches for superparticles at LHC have resulted in gluino mass limits mg 2 TeV [6] and top squark limits mt 1 1 TeV [7]. Meanwhile, direct searches for relic WIMP dark matter by LUX [8], PandaX [9] and Xe-1-ton [10] have failed to detect the SUSY WIMP. Indirect WIMP searches from Fermi-LAT/MAGIC [11], expecting to detect WIMP-WIMP annihilation to gamma rays in dwarf spheroidal galaxies, have also placed strong limits on SUSY WIMPs. Taken together, direct and indirect detection limits have eliminated two previously well-regarded candidates for the nature of SUSY WIMP dark matter.
1. The well-tempered neutralino (WTN) [12], wherein the bino and higgsino components were adjusted to comparable values so as to obtain the required relic density, predicted σ SI ( 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 [13] of the mSUGRA/CMSSM model [14] and is now solidly excluded [8,9,10,15].
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 [16,17,15].
Taken all together, the data seem to suggest that weak scale supersymmetry (WSS) [18], if viable, must have at least strongly coupled superpartners with soft SUSY breaking parameters m sof 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 sof 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 Hu 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. [19]. Requiring no large unexplained cancellations between the various terms on the right-handside of Eq. (1) led us to introduce the electroweak fine-tuning measure ∆ EW [20,19] 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 | p i m 2 Z ∂m 2 Z ∂p i | where p i are fundamental parameters of the theory [21] reduces to ∆ EW [22,23,24] 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 fine-tuning (which we take to be ∆ EW 30) 1 are the following. 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 [19,25].
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 [29] (nNUHM2) with parameter space given by 3 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 [31] (nAMSB) with parameter space given by, where in addition to AMSB contributions to soft terms [32], we introduce several bulk induced soft terms to render sleptons non-tachyonic (m 0 (bulk)) and to render the model natural (m Hu (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 Hu (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.
3. 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 well-motivated superstring models [33]. The parameter space is given by [34,35], where α parametrizes the relative gravity-to anomaly-mediation and c m and c m3 coefficients are a continuous generalization of formerly discrete parameters involving modular weights of the relevant fields and a 3 is a 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.
3 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 [30,22,24]. 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.

Dark matter relic density in natural SUSY
For natural SUSY models, we see from Eq. (1) that naturalness requires |µ| and −m 2 Hu 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 thermally-produced relic density for the various SUSY models introduced in Sec. 1. We use the computer code Isajet 7.88 to compute sparticle mass spectra for the nNUHM2, nAMSB and nGMM models [36]. 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 nAMSB model, we scan over 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). 4 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 [37]) 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 [38] 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 higgsinolike 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.
Two well-motivated classes of mechanisms have been proposed to bring 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. [39]). 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-parity-odd axino fieldã. Both saxion and axino are expected to gain masses of order the gravitino mass m 3/2 in supergravity models [40]. In SUSY axion models, the axions can be produced non-thermally 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 [38]. Axinos can be thermally produced in the early universe and then augment the WIMP abundance via decays after thermal WIMP freeze-out. Saxions can be produced both thermally and non-thermally 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 [41]. For low values of the axion decay constant, f a 10 11 GeV, the WIMP abundance is its thermal value since axinos and saxions tend to decay before WIMP freeze-out. If f a 10 11 GeV, then post-freeze out saxion and axino decays may augment the WIMP abundance. 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 and the SUSY particle mass spectrum (which influences the saxion and axino decay branching fractions) [41].
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 dark matter. For spin-independent (SI), spin-dependent (SD) detection rates, and also the neutrino detection rate at IceCube, the target event rates must be scaled by a factor ξ 5 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 . For mixed axion/WIMP dark matter, we will assume ξ = Ω 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 [42] 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 thermallyproduced abundance to gain accord with the measured density of dark matter. This mechanism was originally suggested by Moroi and Randall [43] 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 [44] 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 [45,46,47]. 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.
3 Bounds on natural SUSY WIMPs from direct and indirect WIMP searches

Direct WIMP detection bounds
In Fig. 2, we show the value of ξσ SI ( 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 [48] for our direct and indirect relic scattering calculations. 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. [18]). 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 tonscale 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).
In Fig. 3, we show ξσ SD ( 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 the natural WIMP-only case with ξ = 1. We also show the current SD limits from the PICO-60 experiment [49] and from IceCube [50] (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 [51] (subject to the caveats mentioned above) although some points might still elude SD detection.
In frame b), we show the ξ = 1 case for natural WIMP-only 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. Figure 2: Plot of points in the σ SI ( Z 1 p) vs. m Z 1 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. Figure 3: Plot of points in the σ SD ( Z 1 p) vs. m Z 1 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. Figure 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).

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 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. 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. These m Z 1 > 300 GeV points are excluded by the SI DD bounds from Fig. 2. Likewise, the lone nNUHM2 point with m Z 1 ∼ 250 GeV is excluded by the IDD bounds with ξ = 1.

Concluding Remarks
In this paper we have examined the direct-and indirect-WIMP detection rates 6 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 nonnegligible 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 SD-direct 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 WIMP-only 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 electron-positron collider with √ s ≥ 500 − 600 GeV [52], or a high energy pp collider operating at √ s ∼ 27 − 33 TeV [57] 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-WIMP-only dark matter is indeed excluded. Unnatural higgsino-like WIMP dark matter can still survive as detailed in Ref. [58,59] 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. [60] and in the pMSSM context in the first of Ref. [61], 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. [62].