Supersymmetry with radiatively-driven naturalness: implications for WIMP and axion searches

By insisting on naturalness in both the electroweak and QCD sectors of the MSSM, the portrait for dark matter production is seriously modified from the usual WIMP miracle picture. In SUSY models with radiatively-driven naturalness (radiative natural SUSY or RNS) which include a DFSZ-like solution to the strong CP and SUSY mu problems, dark matter is expected to be an admixture of both axions and higgsino-like WIMPs. The WIMP/axion abundance calculation requires simultaneous solution of a set of coupled Boltzmann equations which describe quasi-stable axinos and saxions. In most of parameter space, axions make up the dominant contribution of dark matter although regions of WIMP dominance also occur. We show the allowed range of PQ scale f_a and compare to the values expected to be probed by the ADMX axion detector in the near future. We also show WIMP detection rates which are suppressed from usual expectations because now WIMPs comprise only a fraction of the total dark matter. Nonetheless, ton-scale noble liquid detectors should be able to probe the entirety of RNS parameter space. Indirect WIMP detection rates are less propitious since they are reduced by the square of the depleted WIMP abundance.

. Plot of measured Higgs mass squared along with radiative correction and tree-level term 2µ 2 . The latter term is adjusted (fine-tuned) to guarantee that m h = 125 GeV.

Introduction
The discovery of the Higgs boson [1,2] with mass m h = 125.15±0.24 GeV was a great triumph but it brings with it a conundrum: how is it that scalar fields can actually occur in nature? The problem lies in the radiative corrections to their masses: they are quadratically divergent in the energy circulating in the loop diagrams. Since quantum mechanics requires one to sum over a complete set of states, those states with the highest energies bring large quantum corrections which must be compensated by adjusting bare mass terms to maintain the measured value of m h . The situation is depicted in Fig. 1: here we take the SM Higgs potential as V = −µ 2 h |h| 2 + λ h |h| 4 where m 2 h = 2µ 2 h + δm 2 h and m 2 h (tree) = 2µ 2 h . Requiring the quantum corrections not exceed the bare mass (similar to the Gaillard-Lee [3] requirement on ∆m 2 K which predicted the charm quark mass) implies the Standard Model to only be valid at energy scales Q Λ ∼ 1 TeV. These quadratic divergences-which are endemic to scalar quantum fields-led some physicists to ponder whether fundamental scalar fields could really occur in nature [4].
The solution to the above SM naturalness problem was very conservative: expand the fundamental 4 − D spacetime symmetry structure which underlies quantum field theory to its most general structure including graded Lie-algebras [5,6]. The expanded symmetry group-called supersymmetry or SUSY for short-provided once and for all the necessary structure so that scalar field quadratic divergences completely cancelled. Akin to the doubling of particle spectra which occurred when Dirac included Lorentz symmetry into quantum mechanics, SUSY also requires an approximate doubling: under SUSY, for every boson there is a fermion state and vice versa. Since we see e.g. no bosonic electrons with the same mass as electron (similar arguments apply to other SM particle states), SUSY must be a broken symmetry. To stabilize the weak scale, it is expected that SUSY breaking is characterized by soft SUSY breaking terms of weak scale magnitude. In fact, in models bases on local supersymmetry (supergravity or SUGRA), the breakdown of SUSY must occur in a "hidden sector" of the model to maintain phenomenological viability [7]. Taking the limit of M P → ∞ while keeping the gravitino mass m 3/2 fixed, one calculates the soft terms [8] as multiples of m 3/2 where m 3/2 ∼ m 2 hidden /M P . Here M P = 2.4 × 10 18 GeV is the reduced Planck mass. A hidden sector mass scale m hidden ∼ 10 10 GeV gives rise to a weak scale of ∼ 100 GeV. From the above arguments, we arrive at the rough expectation that the new matter particles should inhabit the energy scale Q ∼ 100 − 1000 GeV. Lest one think the above construct is the product of an overly active imagination of theorists, we remark that SUSY is supported by three disparate sets of measurements: • The measured values of the three gauge couplings, when extrapolated to m GUT 2 × 10 16 GeV, very nearly meet at a point [9], as expected in simple unified theories.
• The measured value of the top quark, m t = 173.2 GeV, is in just the right range to drive the up-Higgs soft mass m 2 Hu to negative values, causing the required breakdown of electroweak symmetry [10].
• The measured value of the newly discovered higgs boson, m h 125 GeV, falls squarely within the narrow window m h ∼ 115 − 135 GeV of SUSY requirements which was expected from the pre-LHC era [11]. In contrast, in the SM the Higgs mass could lie anywhere in the 115 − 800 GeV mass range.
In addition, SUSY-as embodied by the MSSM-carries with it several dark matter candidates [12] and several baryogenesis mechanisms [13] whereas the SM contains neither.
In spite of these successes, many authors have proclaimed weak-scale SUSY to be in a state of crisis [14]. While SUSY solves the big hierarchy problem involving quadratic divergences [15], there is a growing Little Hierarchy problem [16] typified by the increasing gap between the W , Z and h masses clustered all around ∼ 100 GeV, and the apparent mass scale of SUSY particles which are seemingly in the multi-TeV range. Presently, LHC8 with 20 fb −1 of data requires mg 1.3 TeV in the case of heavy squark masses and mg 1.8 TeV in the case of comparable squark masses. Furthermore, the value of m h ∼ 125 GeV requires radiative corrections from top-squarks in the tens of TeV range for small top-squark mixing (although few-TeV top squarks are allowed for large mixing induced by trilinear A terms [17]). The lore is that as the mass scale for the soft terms increases, then one must increasingly fine-tune parameters to maintain m W,Z,h ∼ 100 GeV. Since large fine-tuning usually indicates some pathology in any theoretical construct, a number of authors have questioned whether SUSY as we know it is gradually becoming excluded [14]: if so, then new ideas for physics beyond the Standard model are required.
In the following Section 2, we shall refute this point of view. While we shall conclude that many SUSY models are indeed fine-tuned-including the paradigm mSUGRA/CMSSM model-we will find that models characterized by radiatively-driven naturalness [18,19] (radiatively-driven natural SUSY or RNS) are allowed with modest fine-tunings only at the 10% level. Radiatively-driven naturalness occurs in SUSY models with non-universality of Higgs soft terms (as in the NUHM2 model [20]). RNS models are characterized by the presence of light higgsinos with mass µ ∼ 100 − 200 GeV, the closer to m Z the better. The lightest SUSY particle is a candidate for dark matter and is then a higgsino-like WIMP.
We proceed to examine the consequences of RNS for dark matter. In Sec. 3, we require that naturalness occurs also in the QCD sector of the MSSM. This brings to bear the QCD axion albeit as one element of a axion supermultiplet containing also a spin-1/2 R-parity odd axinoã and a spin-0 R-parity even saxion field s. The dark matter then consists of a combination of both axions and higgsino-like WIMPs. In Sec. 4, we present calculations of the expected abundance of axions and WIMPs in RNS SUSY. We display the range in PQ breaking scale f a which is accessible to axion search experiments like ADMX [21]. In Sec. 5, we examine updated prospects for WIMP detection in RNS. While higgsinos may comprise as little as 5-10% of the total dark matter abundance, they should nonetheless be detectable by ton-scale WIMP direct detection experiments owing to their large couplings to the Higgs boson h. Indirect WIMP detection seems less propitious since the detection rate is proportional to the square of the reduced WIMP abundance. We conclude in Sec. 6.

Measuring naturalness in SUSY theories
Any serious discussion of naturalness requires the definition of some measure. But first, an important point to be made is that any quantity can look fine-tuned if one splits it into dependent pieces. By re-writing an observable O as O + b − b and allowing b to be large, the quantity might look fine-tuned. In this trivial example, however, combining dependent contributions into independent units (b − b = 0) obviously erases the presumed source of fine-tuning. To avoid such pitfalls, a simple fine-tuning rule has been proposed [22]: When evaluating fine-tuning, it is not permissible to claim fine-tuning of dependent quantities one against another.

Simple electroweak fine-tuning
The simplest relation between the weak scale and the soft SUSY breaking parameters comes from minimizing the scalar potential of the MSSM to determine the vacuum expectation values (VEVs) [6]. The first minimization condition allows one to trade the bilinear soft term B for the more convenient ratio of VEVs tan β ≡ v u /v d . The second condition is given by where m 2 Hu and m 2 H d are the weak scale soft SUSY breaking Higgs masses, µ is the supersymmetric higgsino mass term and Σ u u and Σ d d contain an assortment of loop corrections to the effective potential (for a listing, see Ref. [19]). For naturalness, we require no large unnatural cancellations between independent terms on the right-hand-side of Eq. 2. For instance, if m 2 Hu is driven to multi-TeV negative values at the weak scale, then the completely unrelated value of µ 2 is required to be multi-TeV positive with such high precision as to yield a Z mass of just 91.2 GeV. This fine-tuning occurs on a daily basis by users of SUSY spectrum generator tools [23][24][25][26], but it is hidden in the computer code. While such tuning is logically possible, the overall scenario seems highly implausible, or highly unnatural (in this case, the Z mass would naturally be expected occur in the multi-TeV range).
The quantity ∆ EW measures this implausibility by comparing the largest contribution on the right-hand-side of Eq. 2 to the value of m 2 Z /2. If they are comparable, then no unnatural fine-tunings are required to generate m Z = 91.2 GeV.
The main requirements for EW naturalness can then be read off from Eq. 2. They are the following: • |µ| ∼ 100 − 200 GeV (the closer to m Z the better) [27][28][29]. We note here that the lower bound on µ 100 GeV comes from accommodating LEP2 limits from chargino pair production searches. A low value of ∆ EW yields an upper bound on |µ| depending on how much fine-tuning one is willing to tolerate. A value ∆ EW < 10 (or ∆ −1 EW > 10%) for fine-tuning implies |µ| < 200 GeV.
• The value of m 2 Hu is driven radiatively to small, and not large, negative values [18,19]. In the mSUGRA/CMSSM model, this occurs in the hyperbolic branch/focus point (HB/FP) region [30]. However, the rather large value of m h requires a large trilinear A 0 parameter. Such a large trilinear pushes the HB/FP out to typically m 0 ∼ 10 − 30 TeV [31]. At such high m 0 , then the top squark contributions Σ u u (t 1,2 ) become large and again one is fine-tuned. Alternatively, in models where the Higgs soft terms are non-universal, such as in the two-extra parameter non-universal Higgs model NUHM2 [20], it is possible to have small µ for any m 0 value by simply raising the GUT scale value of m Hu (GUT) ∼ (1.3 − 2)m 0 .
• The top squark contributions to the radiative corrections Σ u u (t 1,2 ) can become large for stops in the multi-TeV region. However, the radiative corrections are minimized for highly mixed (large A 0 ) top squarks [18]. This latter condition also lifts the Higgs mass to m h ∼ 125 GeV.
One advantage of ∆ EW is that-within the context of the MSSM-it is (as discussed in Ref. [19]) 1. model-independent: if a weak scale spectrum is generated within the pMSSM or via some high scale constrained model, one obtains exactly the same value of naturalness. Other virtues of ∆ EW are that it is: 2. the most conservative of the proposed measures, 3. in principle measureable, 4. unambiguous, 5. predictive, 6. falsifiable and 7. simple to calculate.
The principle criticism of ∆ EW is that-since it involves only weak scale parameters-it may not display the sensitivity of the weak scale to variations in high scale parameters. Below we discuss two competing measures, ∆ HS and ∆ BG . Typically, these latter two measures are implemented in violation of the fine-tuning rule. If implemented in accord with the fine-tuning rule, then both essentially reduce to ∆ EW . In this case, ∆ EW portrays the entirety of electroweak naturalness even including high scale physics.

Large-log measure ∆ HS
The Higgs mass fine-tuning measure, ∆ HS , compares the radiative correction of the m 2 Hu soft term, δm 2 Hu , to the physical Higgs mass m 2 h : If we assume the MSSM is valid up to some high energy scale Λ (which may be as high as m GUT or even M P ), then the value of δm 2 Hu can be found by integrating the renormalization group equation (RGE): By neglecting gauge terms and S (S = 0 in models with scalar soft term universality but can be large in models with non-universality), and also neglecting the m 2 Hu contribution to X t and the fact that f t and the soft terms evolve under Q 2 variation, a simple expression may be obtained by integrating from m SUSY to the cutoff Λ: Here, we take as usual m 2 SUSY mt 1 mt 2 . By requiring [33][34][35][36] ∆ HS 10 (7) then one expects the three third generation squark masses mt GeV. Using the ∆ HS measure of fine-tuning along with m h 125 GeV, one finds some popular SUSY models to be electroweak fine-tuned to 0.1% [37].
Two problems occur within this approach.
1. m 2 Hu (Λ) and δm 2 Hu are not independent: the value of m 2 Hu feeds directly into evaluation of δm 2 Hu via the X t term: the larger the value of m 2 Hu (Λ), then the larger is the cancelling correction δm 2 Hu [38]. It also feeds indirectly into δm 2 Hu by contributing to the evolution of the m 2 Q 3 and m 2 U 3 terms. Thus, the ∆ HS measure as constructed fails the fine-tuning rule [22].
2. In the SM, the SU(2) L ×U(1) Y gauge symmetry can be broken at tree level. However, in the case of SUGRA gauge theories, where SUSY is broken in a hidden sector via the superHiggs mechanism, m 2 Hu ∼ m 2 3/2 > 0. Thus, for SUGRA models, electroweak symmetry is not even broken until one includes radiative corrections. For SUSY models valid up to some high scale Λ m weak , the large log in Eq. 6 is exactly what is required to break EW symmetry in the first place, radiatively driving m 2 Hu to negative values [10].
A simple fix for ∆ HS is to combine the dependent terms into a single quantity. Under such a regrouping [18,19], where now µ 2 and m 2 Hu (Λ) + δm 2 Hu are each independent so each should be comparable to m 2 h in order to avoid fine-tuning. The large log is still present in (m 2 Hu (Λ) + δm 2 Hu ), but now cancellations can occur between the boundary condition and the radiative correction. The regrouping of contributions to m 2 h leads back to the ∆ EW measure since now (m 2 Hu (Λ) + δm 2 Hu ) = m 2 Hu (weak).

The EENZ/BG measure
The traditional measure, ∆ BG , was proposed by Ellis, Enquist, Nanopoulos and Zwirner [39] and later investigated more thoroughly by Barbieri and Giudice [40]. The proposal is that the variation in m 2 Z with respect to high scale parameter variation be small: where the p i constitute the fundamental parameters of the model. Thus, ∆ BG measures the fractional change in m 2 Z due to fractional variation in high scale parameters p i . The c i are known as sensitivity coefficients [40].
To evaluate ∆ BG , we first express m 2 Z in terms of weak scale SUSY parameters as in Eq. 2: where the partial equality obtains for moderate-to-large tan β values and where we assume for now the radiative corrections are small. Next, one needs to know the explicit dependence of the weak scale values of m 2 Hu and µ 2 on the more fundamental high scale parameters. These can be obtained from semi-analytic solutions to the renormalization group equations for m 2 Hu and µ 2 and can be found in Ref. [41].
The place where the application of ∆ BG can go wrong is in the identification of the fundamental parameter set p i . Usually, the set p i is taken to be the various soft terms of particular effective theories such as the MSSM, mSUGRA, NUHM2, etc. which arise from integrating out the hidden sector of the underlying SUGRA theory. In these effective theories, variation of the soft SUSY breaking parameters allows for a wide range of possibilities for the (unknown) hidden sector and the dynamics of SUSY breaking. However, recall that in SUGRA gauge theories with SUSY broken in a hidden sector, all soft parameters are actually computed as multiples of the gravitino mass m 3/2 . This means that for any given hidden sector, the soft terms are all correlated: if one increases the value of m 3/2 , then all soft terms increase in magnitude accordingly: i.e. in SUGRA they are not independent. Combining the contributions of the dependent high-scale soft terms to m 2 Z , we arrive at the simple high scale relation Now, to allow for no large unnatural cancellations in Eq. 11, we require µ 2 ∼ m 2 Z (same as ∆ EW ) and also am 2 3/2 ∼ m 2 Z . This latter condition can be fulfilled if m 3/2 ∼ m Z (which now seems highly unlikely in light of LHC8 sparticle search limits and the value of m h ) or if m 3/2 is large but a is small. Since the µ term hardly evolves between m GUT and m weak , we may equate −2m 2 Hu (weak) am 2 3/2 . Since am 2 3/2 ∼ m 2 Z , then also −m 2 Hu (weak) ∼ m 2 Z : i.e. m 2 Hu can start off large with magnitude of order m 3/2 at m GUT , but can be driven radiatively to small values ∼ −m 2 Z at m weak . This is the case of radiatively-driven naturalness.

Naturalness in QCD: the need for axions
If we insist on naturalness in the electroweak sector, then it is only fair to insist as well on naturalness in the QCD sector. In the early days of QCD, it was a mystery why the two-light-quark chiral symmetry U(2) L ×U(2) R gave rise to three and not four light pions [42]. The mystery was resolved by 't Hooft's discovery of the QCD theta vacuum which didn't respect the U(1) A symmetry [43]. As a consequence of the theta vacuum, one expects the presence of a term L θ 32π 2 F AµνF µν A (12) in the QCD Lagrangian (whereθ = θ + arg(det(M)) and M is the quark mass matrix). Measurements of the neutron EDM constrainθ 10 −10 leading to an enormous fine-tuning inθ: the so-called strong CP problem.
The strong CP problem is elegantly solved by Peccei, Quinn, Weinberg and Wilczek (PQWW) [44] via the introduction of PQ symmetry and the concomitant (invisible [45,46]) axion: the offending term can dynamically settle to zero. The axion is a valid dark matter candidate in its own right [47].
Introducing the axion in a SUSY context solves the strong CP problem and renders naturalness to QCD. As a bonus, in the context of the SUSY DFSZ axion model [46] where the Higgs superfields carry PQ charge, one gains an elegant solution to the SUSY µ problem. The most parsimonius implementation of the strong CP solution involves introducing a single MSSM singlet superfield S carrying PQ charge Q P Q = −1 while the Higgs fields both carry Q P Q = +1. The usual µ term is forbidden, but we have a superpotential [48,49] If PQ symmetry is broken and S receives a VEV S ∼ f a , then a weak scale µ term is induced which gives µ ∼ m Z for f a ∼ 10 10 GeV. Although Kim-Nilles sought to relate the PQ breaking scale f a to the hidden sector mass scale m hidden [48], we see now that the Little Hierarchy could emerge due to a mis-match between PQ breaking scale and hidden sector mass scale f a m hidden . For the remainder of this paper, we will assume the SUSY DFSZ axion model holds due to its role in solving the SUSY µ problem.
An elegant model which exhibits this behavior was proposed by Murayama, Suzuki and Yanagida (MSY) [50]. In the MSY model, PQ symmetry is broken radiatively by driving one of the PQ scalars X to negative mass-squared values in much the same way that electroweak symmetry is broken by radiative corrections driving m 2 Hu negative. Starting with multi-TeV scalar masses, the radiatively-broken PQ symmetry induces a SUSY µ term ∼ 100 GeV [51] while at the same time generating intermediate scale Majorana masses for right-hand neutrinos: see Fig. 2. Although we get different solutions for the PQ scale by setting m 3/2 to different masses at Planck scale, the PQ scalar X is driven to negative mass-squared values without m 3/2 dependence at the same Q value. However, the coupling h shifts the position of the Q value where m 2 X becomes negative; increasing h shifts the point to higher energy scales. In models such as MSY, the Little Hierarchy µ m 3/2 is no problem at all but is instead just a reflection of the mis-match between PQ and hidden sector mass scales.

Relic abundance of axions and WIMPs with implications for axion detection
It is straightforward to calculate the thermal abundance of WIMPs in natural SUSY. To a good approximation, it is given by where s 0 is the current entropy density of the universe, ρ c is the critical density, h is the scaled Hubble constant, x f = m Z 1 /T f is the scaled WIMP freeze-out inverse temperature ∼ 25 and σv is the thermally averaged WIMP annihilation cross section times relative velocity. For a higgsino-like LSP as occurs in RNS, σv is large due to higgsino annihilation into vector boson pairs W W and ZZ. The simple "WIMP miracle" picture seems not to apply to higgsino dark matter where we show in Fig. 3 Ω Z 1 h 2 (from IsaReD [52]) vs. m Z 1 from a scan over NUHM2 parameter space. Here, we keep only solutions with constraints: 1. µ > 100 GeV in accord with LEP2 searches for chargino pair production, 2. 123 < m h < 128 GeV in accord with the CERN Higgs discovery, and allowing for some theoretical error in the RG-improved one loop effective potential computation of m h in Isajet [53] and 3. ∆ EW < 30 (100) as denoted by green stars (blue crosses). The plot shows that for low m Z 1 ∼ 100 GeV (as preferred by naturalness) the predicted thermal abundance of WIMPs is typically a factor 10-30 below the measured value of cold dark matter (CDM) Ω CDM h 2 0.12. If we require ∆ EW < 30, then m Z 1 reaches ∼ 300 GeV maximally with Ω Z 1 h 2 as high as 0.02. At some cost to naturalness, h 2 versus m Z 1 from a scan over NUHM2 parameter space with ∆ EW < 30 (green stars) and ∆ EW < 100 (blue crosses). We also show the central value of Ω CDM h 2 from WMAP9. Ω Z 1 h 2 approaches the measured value for m Z 1 ∼ 600 GeV, where the Z 1 is already frequently a mixed bino-higgsino particle.
Naively, one might expect natural SUSY to be ruled out as being incapable of generating a sufficiently large relic density of WIMPs. However, naturalness in both EW and QCD sectors implies the presence of two dark matter particles: the WIMP and the axion. Axions are expected to be produced dominantly via the non-thermal Bosonic Coherent Motion (BCM) [47,54] yielding where θ i is the initial axion mis-alignment angle, f a is the axion decay constant and N DW is the domain-wall number. Also, f (θ i ) accounts for anharmonicity effects. By proper choice of f a and θ i , BCM-produced axions can always account for the measured CDM abundance. 1 However, as mentioned previously, the axion superfield also contains a spin-1/2 axinoã and a spin-0 saxion s. In SUGRA, one expects m s ∼ m 3/2 while the axino mass is more model-dependent but generally one expects also mã ∼ m 3/2 [57]. The DFSZ axinos can be produced thermally in the early universe at a rate ∝ f −2 a and largely independent of the re-heat temperature T R [58] (in the SUSY KSVZ model, then axino thermal production is proportional to T R ). Once axinos are produced, they undergo 1 Here, we impicitly assume that PQ symmetry is broken before the end of inflation so that topological defects and archioles do not contribute to the ultimate axion relic density [55,56]. (late) decays to sparticle plus particle thereby injecting additional WIMPs into the thermal plasma. If enough WIMPs are produced at the axino decay temeprature, they undergo a process of re-annihilation which still yields an enhanced WIMP abundance [59], but not as much as one-to-one with the population of thermally produced axinos. Of equal importance to WIMP production from axino decays is the axino decay temperature: if axinos decay before WIMP freeze-out, then the injected WIMPs thermalize and one regains the usual thermal WIMP abundance. If axinos decay after WIMP freeze-out, then they always augment the WIMP abundance.
Saxions can also be produced thermally at rates comparable to axino thermal production. In addition, saxions can be produced via BCM which is especially important at large f a . Since saxions are R-parity even, they can decay to pairs of SM particles, thereby injecting extra entropy into the plasma, or they can decay to pairs of SUSY particles, thus also augmenting the WIMP abundance (depending again on the saxion decay temperature). Depending on a combination of PQ charge assignments and VEVs parametrized by ξ s , the saxions may also decay toãã (if kinematically allowed) thus adding to the WIMP abundance, or they may decay to aa thus injecting additional dark radiation into the thermal plasma. Strong limits on dark radiation-parametrized by the effective number of additional neutrinos present in the universe ∆N eff -have been obtained, with a combination of Planck and other data sets finding N eff = 3.15 ± 0.23 [60] (whereas the SM predicts N eff = 3.046). Thus, too much dark radiation from saxion decay can lead to conflict with measured cosmological parameters. In our numerical study, we consider a conservative constraint ∆N eff < 1 (see Fig. 6) at over 3σ with the joint P lanck TT+lowP+BAO result [60]. In addition, if saxions or axinos of sufficient initial abundance decay after the onset of BBN, then they can destroy the successful predictions of light element abundances via BBN, and again the model can be excluded.
The calculation of the mixed axion-neutralino relic abundance can be calculated via semi-analytic techniques [61] or more reliably [62] via the simultaneous solution of eight coupled Boltzmann equations describing the energy densities of 1. radiation, 2. thermally-and decay-produced WIMPs, 3. BCM-produced axion, 4. BCM produced saxions (followed by saxion decay), thermal production and decay of 5. axino, 6. saxions and 7. thermal and decay-induced production of axions and 8. thermal production and decay of gravitinos.
We scan over the following PQ parameters: The result of these calculations were shown in Ref. [62] and in Fig. 4 where the energy densities are tracked as a function of scale factor R from the end of inflation with T = T R to the era of entropy conservation.
In Fig. 5, we show the calculated relic abundance of both WIMPs (blue and red points) and axions (purple points) as a function of f a for a RNS benchmark SUSY model with m 0 = 5000 GeV, m 1/2 = 700 GeV, A 0 = −8300 GeV, tan β = 10, µ = 110 GeV and m A = 1000 GeV. We first take ξ s = 0 so saxion decays to aa andãã are turned off.
At very low f a , axinos are thermally-produced at a large rate but also decay well before neutralino freeze-out so that the WIMP abundance is still given by its expected thermally-produced value. As f a increases, ultimately axinos begin decaying after freeze-out thus augmenting the WIMP abundance. For f a > 10 13 GeV, too many WIMPs are produced and the model parameters are excluded. For very large f a ∼ 10 15 GeV, all points are doubly excluded by producing too much dark matter and violating limits from BBN [63]. We also show the axion abundance. At very low f a , the CDM is axion-dominated [61] although this requires very high values of θ i ∼ π (see frame 5b). which might be considered fine-tuned. For f a ∼ 10 12 GeV, axions can still dominate the CDM abundance but with θ i ∼ 1. For these values of f a , the CDM could also easily be WIMP dominated as well.
In Fig. 6, we show the neutralino and axion relic abundance for the RNS benchmark with ξ s = 1 (saxion decays to axions and axino pairs are turned on). In this case, the additional decay modes allow the saxion to be shorter lived for a given value of m s and f a compared to the ξ s = 0 case. As a consequence, there is a greater range of f a where CDM can be axion-dominated. Ultimately, axinos and saxions decay after freeze-out and the WIMP abundance is enhanced at higher f a ∼ 10 11 − 10 14 GeV. For f a 10 14 GeV, WIMPs are overproduced. Points at very high f a for ξ s = 1 can be triply excluded by producing too many WIMPs and by violating both dark radiation and BBN constraints.
We summarize the results of this section in Fig. 7. We display the range of f a where valid solutions for the relic abundance of mixed axion-higgsino CDM can be found for the RNS benchmark model. The upper bar shows the range of f a for ξ s = 0 while the lower bar shows the range for ξ s = 1. The darker shaded parts of the bars denote θ i values > 3 which might be considered less plausible or fine-tuned. We also show by the bracket the range of f a , assuming the bulk of DM is axion, which is expected to be probed by the ADMX experiment within the next several years [64]. This region probes the most natural region where θ i ∼ 1. We also show a further region of lower f a which might be explored by a new Figure 5. In a) we plot the neutralino relic density from a scan over SUSY DFSZ parameter space for the RNS benchmark case with ξ s = 0. The grey dashed line shows the points where DM consists of 50% axions and 50% neutralinos. In b), we plot the misalignment angle θ i needed to saturate the dark matter relic density Ω Z 1 a h 2 = 0.12. Figure 6. In a) we plot the neutralino relic density from a scan over SUSY DFSZ parameter space for the RNS benchmark case with ξ = 1. The grey dashed line shows the points where DM consists of 50% axions and 50% neutralinos. The red BBN-forbidden points occur at f a 10 14 GeV and are covered over by the brown ∆N eff > 1 coloration. In b), we plot the misalignment angle θ i needed to saturate the dark matter relic density Ω Z 1 a h 2 = 0.12. Figure 7. Range of f a which is allowed in each PQMSSM scenario for the RNS benchmark models. Shaded regions indicate the range of f a where θ i > 3. open resonator technology [65]. About a decade of natural f a ∼ 10 14 GeV seems able to elude ADMX searches for the ξ s = 1 case.

Results for variable µ
We may convert the RNS benchmark point into a model line by allowing for variable µ. In this case, we have variable higgsino mass with the lower bound given by the LEP2 limit µ ∼ m W 1 > 103.5 GeV while the upper bound is determined by how much fine-tuning one is willing to tolerate with In Fig. 8, we show the thermally-produced relic density of neutralinos along the variable µ RNS model line. Here, Ω Z 1 h 2 ∼ 0.007 for low µ but increases as µ increases since the Z 1 becomes increasingly bino-like. At µ ∼ 300 GeV, the Z 1 becomes more bino-like than higgsino-like, and at µ ∼ 340 GeV, naively too much neutralino dark matter is produced. As we have seen, it is easy to increase the neutralino abundance from its thermal expectation by allowing for axino and saxion production with decay taking place after neutralino freeze-out. It is much harder to reduce the neutralino abundance from its thermal value: the three most common ways include 1. entropy dilution from saxion decay to SM particles only at very high f a ∼ 10 15 GeV, 2. allowing for R-parity violation (in which case one must somehow stabilize the proton) or 3. allowing for a lighter LSP than the neutralino (e.g. a light axino or gravitino into which the neutralino may decay).
In Fig. 9, we plot the contours of allowed regions (allowed below the contours) in the f a vs. µ plane by varying µ along the RNS model line for the ξ s = 1 case. We show the boundaries for three different assumptions on the axino/saxion masses: mã ,s = 5, 10 and 20 TeV. The lower bound is always f a 10 9 GeV from supernovae/red giant astrophysical cooling limits [66] (although the lower range Figure 9.
Contours of allowed f a region as µ is varies along the RNS benchmark model-line.
requires some tuning on θ i ∼ π). For our canonical case where we expect m 0 mã ,s ∼ m 3/2 = 5 TeV, f a can range up to 10 13 GeV beyond which too much neutralino mass density is produced. As µ increases, the upper bound tends to decrease because the neutralino thermal abundance is increasing and there is less "room" for additional neutralino production from axino/saxion decay. As mã ,s increase, the upper bound on f a increases. This is because as mã ,s become more massive, their widths increase and their lifetimes decrease: for a given f a value, they are more likely to decay at earlier times and so re-annihilation from decay-produced neutralinos occurs at higher axino/saxion decay temperature (and the re-annihilation yield is inversely proportional to decay temperature [59]). For the ξ s = 0 case, we have more constrained upper f a boundaries since saxion decays into axions and axinos are turned off and hence the saxion is longer lived.

Direct and indirect detection of WIMPs
In this Section, we update our previous projections [67] for direct and indirect detection of higgsino-like WIMPs from radiatively-driven natural SUSY. Our current results contain several improvements: 1. Our previous scan over NUHM2 parameter space was restricted to a range of m A : 0.15 − 1.5 TeV.
However, low ∆ EW solutions can be found for much higher m A values [68] and so here we expand the m A range to as far as 20 TeV so that the bounds on our scanned parameter space is dictated by the value of ∆ EW rather than an arbitrary parameter cutoff.
2. We have updated the nucleon mass fraction parameters which enter the quark and gluon matrix elements in IsaReS [69] to values given in Table 1 of Ref. [70]. These mainly lessen the contribution from strange quarks from older estimates of the spin-dependent scattering cross section. In our case, the computed values of σ SI ( Z 1 p) decrease by typically a factor of two.
3. We have increased our sampling statistics in NUHM2 parameter space. Figure 10. Plot of rescaled higgsino-like WIMP spin-independent direct detection rate ξσ SI ( Z 1 p) versus m Z 1 from a scan over NUHM2 parameter space with ∆ EW < 30 (green) and 30 < ∆ EW < 100 (blue). We also show the current reach from the LUX experiment and projected reaches of Xe-1-ton, LZ(10) and Darwin.
In Fig. 10, we show the spin-independent neutralino-proton scattering rate in cm 2 as calculated using the updated IsaReS [69]. The result is rescaled by a factor ξ = Ω std Z 1 h 2 /0.12 to account for the fact that the local relic abundance might be less than the usually assumed value ρ local 0.3 GeV/cm 3 , as suggested long ago by Bottino et al. [71] (the remainder would be composed of axions). Green stars denote points with ∆ EW < 30 while blue crosses denote points with 30 < ∆ EW < 100.
The higgsino-like WIMP in our case scatters from quarks and gluons mainly via h exchange. The Z 1 − Z 1 − h coupling involves a product of both higgsino and gaugino components. In the case of RNS models, the Z 1 is mainly higgsino-like, but since m 1/2 is bounded from above by naturalness, the Z 1 contains enough gaugino component that the coupling is never small: in the notation of Ref. [6] L −X h where and where v 4 are the bino and wino components of the lightest neutralino, α is the Higgs mixing angle and g and g are SU(2) L and U(1) Y gauge couplings. Thus, for SUSY models with low ∆ EW 30 − 100, the SI direct detection cross section is also bounded from below, even including the rescaling factor ξ.
From Fig. 10, we see that the current reach from the LUX experiment (solid contour) has begun sampling the upper limits of predicted ξσ SI ( Z 1 p) values. The projected reach of Xe-1-ton, a ton scale liquid Xenon detector, is also shown. It is seen to cover nearly all the predicted parameter space points. We also show the projected reach of LZ (10), an upgrade to LUX. The projected LZ reach is seen to cover the entire set of points generated. Thus, the projected ton scale noble liquid detectors (or other comparable WIMP detectors) can make a complete exploration of the RNS parameter space. Since deployment of these ton-scale detectors is ongoing, it seems that direct WIMP search experiments may either verify or exclude RNS models in the near future. These searches should either verify or rule out a very essential aspect of natural SUSY models.
In Fig. 11, we show the rescaled spin-dependent neutralino-proton scattering cross section ξσ SD ( Z 1 p). Here we show recent limits from the COUPP [72] detector. Current limits are still about an order of magnitude away from reaching the predicted rates from RNS models. We also show limits from the IceCube experiment. IceCube searches for high energy neutrinos which could be produced from WIMP annihilations in the solar core. The IceCube expected rates depend on the Sun's ability to capture WIMPs which in turn depends on a product of spin-dependent neutralino-proton scattering cross section times the local WIMP abundance. 2 The IceCube limits have barely entered the RNS parameter space and excluded just the largest values of ξσ SD ( Z 1 p).
In Fig. 12, we show the rescaled thermally-averaged neutralino annihilation cross section times relative velocity in the limit as v → 0: ξ 2 σv | v→0 . This quantity enters into the rate expected from WIMP halo annihilations into γ, e + ,p ord. The rescaling appears as ξ 2 since limits depend on the 2 In a previous work [67], it was mistakenly suggested that the IceCube detection rate was independent of local abundance due to equilibration between solar capture rate and WIMP annihilation rate. Figure 11. Plot of rescaled spin-dependent higgsino-like WIMP detection rate ξσ SD ( Z 1 p) versus m Z 1 from a scan over NUHM2 parameter space with ∆ EW < 30 (green stars) and 30 < ∆ EW < 100 (blue crosses). We also show current reach from the COUPP and IceCube detectors. Figure 12. Plot of rescaled ξ 2 σv | v→0 versus m Z 1 from a scan over NUHM2 parameter space with ∆ EW < 30 (green stars) and 30 < ∆ EW < 100 (blue crosses). We also show current reach from Fermi LAT, Ref. [78].
square of the local WIMP abundance [76]. Anomalies in the positron and γ spectra have been reported, although the former may be attributed to pulsars [77], while the latter 130 GeV gamma line may be instrumental. On the plot, we show the limit derived from the Fermi LAT gamma ray observatory [78] for WIMP annihilations into W W . These limits have not yet reached the RNS parameter space due in part to suppression from the squared rescaling factor.

Conclusions
We have found in this paper, and in previous works, that if one insists on naturalness-in both the electroweak and the QCD sectors-then the simple picture of SUSY WIMP dark matter changes radically. Naturalness in the electroweak sector implies a low value of the superpotential µ parameter: the closer to m Z the better. In models with gaugino mass unification, as favored in simple GUTs, this implies the LSP is a higgsino-like neutralino with a predicted thermal abundance a factor of 10-15 below the measured dark matter density. This seeming disaster is in fact an attribute if one also insists on naturalness in the QCD sector, i.e solving the strong CP problem. In this case, the most compelling solution invokes a PQ symmetry with its concommitant axion. In this situation, the axion makes up the remaining abundance, and in fact over most of parameter space the axion is the dominant CDM component while WIMPs are subdominant.
Invoking the axion in a SUSY context brings along both the axino and the saxion. The dark matter abundance calculation becomes more intertwined since axions can be produced via BCM, via thermal production and via saxion decay. WIMPs can be produced thermally but also via axino, saxion and gravitino decays. If WIMPs are produced via decays at sufficient rates, then WIMP re-annihilation occurs. Additional entropy can be produced at late times by the decays of heavy unstable states, thus diluting all relics which are present. The ensuing abundance calculation is more complicated than the simple WIMP miracle picture, but in many ways it is more elegant and compelling. Our abundance calculations here have used the SUSY DFSZ axion model which provides an elegant solution to the SUSY µ problem. We have outlined the range of f a values which are allowed in RNS, and shown the regions which ADMX and other experiments hope to probe in the near future.
With regard to WIMP detection, higgsino-like LSPs which contain significant gaugino components (as is required in natural SUSY) generally have large rates for both direct and indirect detection, at least compared to binos. However, the propitious detection rates are off-set by the fact that now the WIMPs might comprise only a small fraction of the local abundance instead of the entirety of CDM. To compensate, one must temper detection rates by the ξ = Ω Z 1 h 2 /0.12 factor. For instance, direct detection via SI or SD scattering are both reduced by a factor ξ. Nonetheless, ton-scale noble liquid WIMP detectors are projected to probe the entirety of RNS parameter space: if a WIMP signal is not ultimately seen, then the RNS picture will have to be seriously modified or abandoned. Detection rates for indirect WIMP searches via halo WIMP annihilation into gammas or antimatter are suppressed by a factor of ξ 2 . This suppression will make detection of WIMPs in these channels more difficult, except in the cases where WIMPs still comprise the bulk of dark matter.