Strong constraints of LUX-2016 results on the natural NMSSM

Given the fact that the relatively light Higgsino mass $\mu$ favored in natural supersymmetry usually results in a sizable scattering cross section between the neutralino dark matter and the nucleon, we study the impact of the recently updated direct detection bounds from LUX experiment, including both Spin Independent (SI) and Spin Dependent (SD) measurements, on the parameter space of natural Next-to-Minimal Supersymmetric Standard Model (nNMSSM). Different from the common impression that the SI bound is stronger than the SD one, we find that the SD bound is complementary to the SI bound and in some cases much more powerful than the latter in limiting the nNMSSM scenarios. After considering the LUX results, nNMSSM is severely limited, e.g. for the peculiar scenarios of the NMSSM where the next-to-lightest CP-even Higgs corresponds to the $125 {\rm GeV}$ Higgs boson discovered at the LHC, the samples obtained in our random scan are excluded by more than $85\%$. By contrast, the monojet search at the LHC Run-I can not exclude any sample of nNMSSM. We also investigate the current status of nNMSSM and conclude that, although the parameter points with low fine tuning are still attainable, they are distributed in some isolated parameter islands which are difficult to get. Future dark matter direct search experiments such as XENON-1T will provide a better test of nNMSSM.


Introduction
From various cosmological and astrophysical observations, it has been a well established fact that over 20% of the energy density of the Universe today is composed of Dark Matter (DM) [1]. Among the candidates predicted in new physics models beyond the Standard Model (SM), the Weakly Interacting Massive Particle (WIMP) is a very promising one which has a mass around the electroweak scale and couplings to the SM particles close to the electroweak strength. As a typical example, the Lightest Supersymmetric Particle (LSP) in various Supersymmetric (SUSY) models falls into this category [2]. Motivated by the interactions predicted between DM and SM sector, many direct detection (DD) experiments are on going to search for possible scattering signals of DM particle off the nuclei.
Recently there are updated results from several groups including PICO-2L [3], PandaX-II [4] and LUX [5,6], and the results cover both Spin-Independent (SI) and Spin-Dependent (SD) scattering between DM and nuclei. In many cases the SI scattering is considered to be more promising in detecting DM signals due to its coherent property, of which the scattering cross section is proportional to A 2 of the nucleus and can benefit from the heavy nuclear elements [2]. On the contrary, the SD scattering cross section suffers from the cancelation of the spins of nucleon pairs in the nucleus and thus does not have the A 2 enhancement [7]. However, there are SUSY parameter space inducing cancelations in the SI amplitude and resulting in small and even vanishing SI cross section, the so-called Blind Spot (BS) [8][9][10][11], in which case one has to consider the SD detection. Moreover, nuclei isotopes with un-paired nucleon and high abundance can be good targets to detect SD scattering, e.g. xenon used in XENON and LUX experiments and fluorine used in PICO experiments. Consequently, SI and SD detection methods are complementary to each other and should be considered together if one wants to constrain the parameter space of a certain model.
A particularly interesting SUSY scenario sensitive to DD experiment is natural SUSY (NS) [12][13][14][15][16][17][18][19][20][21]. In this scenario, the Higgsino mass µ tends to be small, and consequently the lightest neutralino χ 0 1 as the DM candidate contains sizable Higgsino components, which enables it to couple rather strongly with CP-even Higgs bosons h i and Z-boson. Given the fact that t-channel exchange of h i (Z boson) is the dominant contribution to SI (SD) cross section for DM-nucleon scattering at tree-level, it is speculated that the continuously improved sensitivity in DD experiments can be promising to test NS if χ 0 1 is fully responsible for the current DM relic density. Since the NS scenario in the Minimal Supersymmetric Standard Model (MSSM) is theoretically unsatisfactory [19], we here investigate this issue in the NS scenario of the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [22]. To be more specific, we study the samples obtained in [19] which predict rather low fine tunings in getting electroweak observables m Z and m h and meanwhile satisfy various experimental constraints, e.g. the DM relic density measured by WMAP and Planck [23,24], the LUX-2015 limit on SI cross section [25] as well as the direct searches for supersymmetric particles at LHC Run-I. Our analyses indicate that the constraints from the LUX experiment in 2016 (LUX-2016), especially those from the upper bounds on SD cross section, are very strong in limiting the NS scenario. Numerically speaking, we find that the LUX-2016 results can exclude more than 90% Type III and Type IV samples obtained in [19]. Although the exact percentage may vary with different scan strategies, this number can exhibit the high sensitivity of the nNMSSM scenarios to the DD experiment. We note that the importance of the SD cross section in limiting SUSY parameter space was not exhibited sufficiently before. This paper is organized as follows. In Section 2, we recapitulate the basics of the NMSSM including the calculation of the SI and SD scattering rates. In Section 3 we illustrate the capability of the LUX-2016 results of limiting the NS scenario of NMSSM (nNMSSM). We also study the status of nNMSSM after the LUX-2016 experiment, which is presented in Section 4. Finally, we draw our conclusions in Section 5.

Basics about the NMSSM
In this section we briefly recapitulate the basics of the NMSSM including its Higgs and neutralino sectors, the naturalness argument and the calculation of SI and SD scattering rates. More detailed discussion and complete formulae about the basics can be found in [7,10,19,22] and references therein.

Natural NMSSM
The superpotential of the Z 3 -invariant NMSSM takes the following form [22] where W F is the superpotential of the MSSM without the µ-term, and λ, κ are dimensionless parameters describing the interactions among the Higgs superfields. The Higgs potential of the NMSSM consists of the F-term and D-term of the superfields, as well as the soft breaking terms where H u , H d and S denote the scalar component of the superfieldsĤ u ,Ĥ d andŜ, respectively. In practice, it is convenient to rotate the fields H u and H d by where ε is an antisymmetric tensor with ε 12 = −ε 21 = 1 and ε 11 = ε 22 = 0, and tan β ≡ v u /v d with v u and v d representing the vacuum expectation value of H u and H d fields, respectively. After this rotation, the redefined fields H i (i = 1, 2, 3) have the following form where H 2 corresponds to the SM Higgs doublet with G + , G 0 being the Goldstone bosons eaten by W and Z bosons respectively, while H 1 represents a new SU (2) L doublet scalar field and it has no coupling to W and Z bosons at tree-level. From Eq.(2.4), it is obvious that the Higgs sector of the NMSSM consists of three CP-even mass eigenstates h 1 , h 2 and h 3 , which are the mixtures of S 1 , S 2 and S 3 , two CP-odd mass eigenstates A 1 and A 2 composed by P 1 and P 2 , as well as two charged Higgs H ± . In the following, we assume m h 1 < m h 2 < m h 3 and m A 1 < m A 2 , and call h i the SM-like Higgs boson if its dominant component is composed of the field S 2 . The eigenstates h i are related to the fields S j by with U being the rotation matrix to diagonalize the mass matrix for the S i fields. An interesting feature of NMSSM is that the squared mass term of the filed S 2 in the SM-like Higgs double H 2 is given by where the first term on the right side is the MSSM contribution, and the second term is peculiar to any gauge singlet extension of the MSSM. Moreover, if the relation m 2 S 3 S 3 < m 2 S 2 S 2 holds, the mixing between the fields S 2 and S 3 can further enhance the SM-like Higgs mass. In this case, h 1 is a singlet-dominate scalar and h 2 plays the role of the SM Higgs boson. Benefiting from these features, m h 2 125GeV does not necessarily require a large radiative contribution from stop loops [26,27].
Instead of using the soft parameters m 2 Hu , m 2 H d and m 2 S , one usually trades them for m Z , tan β and µ ≡ λv s by implementing the scalar potential minimization conditions, and uses the following input parameters λ, κ, tan β, µ, M A , A κ , (2.6) where the parameter A λ in Eq.(2.2) is replaced by the squared mass of the CP-odd field P 1 given by Note that M A represents the mass scale of the doublet H 1 and is usually larger than about 300 GeV from the LHC searches for non-standard doublet Higgs bosons. The neutralino sector of the NMSSM consists of the fields BinoB 0 , WinoW 0 , Higgsinos H 0 d,u and SinglinoS 0 , which is the fermion component of the superfieldŜ. Taking the basis , one has the following symmetric neutralino mass matrix where M 1 and M 2 are Bino and Wino soft breaking mass respectively. In the limit of |M 1 |, |M 2 | |µ|, the Bino and Wino components decouple from the mixing, and the remaining three light neutralinosχ 0 i (i = 1, 2, 3) can be decomposed intõ where the elements of the rotation matrix N can be approximated by with mχ0 i denoting the mass ofχ 0 i . In this case, the DM candidate may be either Singlinodominated or Higgsino-dominated [19]. As has been pointed out by numerous studies, DM may achieve its measured relic density in following regions • Higgs boson and Z boson resonance regions, where the Higgs boson may be any of the three CP-even and two CP-odd Higgs bosons.
Naturalness in the NMSSM can be measured by the following two quantities [15] Here h denotes the SM-like Higgs boson, and p i are parameters defined at the weak scale including those in Eq.(2.6) and the top quark Yukawa coupling Y t which is used to estimate the sensitivities of m Z and m h to stop masses [12]. Apparently, ∆ Z (∆ h ) reflects the sensitivity of m Z (m h ) to SUSY parameters at weak scale and a larger value for any of ∆ Z and ∆ h corresponds to more tuning. Formulae of calculating ∆ Z and ∆ h can be found in [12] and [16] respectively.

Blind spot in spin independent cross section
In the NMSSM with heavy squarks, the dominant contribution to SI DM-nucleon scattering comes from t-channel exchange of the CP-even Higgs bosons [2,7,28,29]. The SI cross section is then expressed as where n = {p, n} denotes nucleon, µ r is the reduced mass of DM and the nucleon, and 1 representing the coupling of h i with DM (nucleon) [10,11]. The explicit expressions of C h i χχ and C h i nn are given by [22] N n|m q qq|n for q = u, d, s denoting the normalized light quark contribution to nucleon mass, and f (n) related with other heavy quark mass fraction in nucleon [28,29].
The blind spot is defined as the SUSY parameter point for which the SI cross-section vanishes. From the formula in Eq.(2.12), one can get the analytic expression of the BS condition which was studied in detail in [10,11]. This condition correlates the parameters λ, κ, tan β, µ, m A and mχ0 1 in a nontrivial way and in some special cases its expression is quite simple. In the following, we consider two specific situations to illustrate this point.
• Case I: tan β 1, and both the gaugino fields and the singlet field S 3 decouple from the DM-nucleon scattering. In this case, the BS condition takes the simple form (see Eq.(50) in [10]) (2.16) • Case II: h 1 and h 2 correspond to the singlet-dominated and SM-like Higgs boson respectively, and both the gaugino fields and the heavy doublet field S 1 decouple from the DM-nucleon scattering. For this case, the BS condition reads (see Eq.(64) in [10]) with γ, η and A s defined by We remind that the LUX-2016 experiment has imposed an upper bound of the SI cross section at the order of 10 −45 cm 2 . Confronted with such a situation, one can infer that the cancelation among the different h i contributions usually exists, and if no signal is detected in future DD experiments BS will become important for SUSY to coincide with experimental results. . This figure indicates that the upper bound of the LUX experiment on σ SD χ−n is able to exclude the parameter region with µ/(λv) 3, and the tightest limit comes from the region of mχ0

Spin dependent cross section
In the heavy squark limit, only t-channel Z exchange diagram contributes to the SD cross section at tree level in the NMSSM. The cross section is then given by [10,11] σ SD χ−(n) with C p ≈ 4.0 and C n ≈ 3.1 for the typical values of f (n) q . With the decoupling of gauginos, we have the following simple expression by using the approximation in Eq.(2.10). From Eqs.(2.18) and (2.19), one can immediately see that σ SD χ−n 0.76 × σ SD χ−p , and σ SD χ−n vanishes with tan β = 1 or pure Higgsino/Singlino DM. Given that the LUX-2016 limit on σ SD χ−n is much stronger than that on σ SD χ−p , we hereafter only consider σ SD χ−n in the following discussion 2 . From the expressions one can 2 Note that although the PICO limit on σ SD χ−p is much stronger than that of the LUX-2016 experiment [3,5], it is still weaker than the LUX-2016 limit on σ SD χ−n in constraining SUSY parameter space after considering the correlation σ SD also get constant contours of σ SD χ−n on µ/(λv) − mχ0 1 /µ plane, which are shown in Fig.1 for tan β = 2, 10 respectively. This figure indicates that for fixed mχ0 1 /µ, the SD cross section decreases monotonously with the increase of µ/(λv), while for fixed µ/(λv) the cross section increases with the increase of mχ0 1 /µ to reach its maximum at mχ0 1 /µ ∼ (0.8 − 0.9). Given that the LUX-2016 limit on σ SD χ−n is at the order of 10 −40 cm 2 [5], one can infer that the experiment can exclude the parameter region with µ/(λv) 3 and the tightest limit comes from the region mχ0 1 /µ ∼ (0.8 − 0.9). Throughout this work, we use the package NMSSMTools [31] to get the particle spectrum of the NMSSM, the package micrOMEGAs [32] to calculate DM relic density and the SI and SD cross sections 3 . We use the default setting of micrOMEGAs, i.e. σ πN = 34 MeV and σ 0 = 42 MeV, to get the values of f (n) q [28,29]. We checked that if we take σ πN = 59 MeV from [33] and σ 0 = 58 MeV from [34][35][36], the SI cross section will be enhanced by a factor from 20% to 40%. Other related discussions can be found in [37][38][39].

Strong constraints of the LUX-2016 results on nNMSSM
In order to study the constraints of the LUX-2016 results on nNMSSM, we consider the samples discussed in our previous work [19]. These samples were obtained in the following way: • First, we fixed all soft breaking parameters for first two generation squarks and gluino mass at 2 TeV. We also assumed a common value for all soft breaking parameters in slepton sector (denoted by ml hereafter) and m U 3 = m D 3 , A t = A b for soft breaking parameters in third generation squark section.
• Second, we scanned by Markov Chain method following parameter space of the NMSSM 0 < λ ≤ 0.75, 0 < κ ≤ 0.75, 2 ≤ tan β ≤ 60, 100 GeV ≤ ml ≤ 1 TeV, The likelihood function we adopted is In above expressions, Ω obs = 0.1198/h 2 with h being the normalized Hubble constant is the cosmological DM parameter obtained in the latest PLANCK results [24], δΩ is the error including both the observational and theoretical uncertainties as (δΩ) 2 = Figure 2. Type I samples projected on σ SĨ χ−p − mχ0 1 plane (left panel) and σ SD χ−n − mχ0 1 plane (right panel) respectively. The dark green samples are excluded by the LUX-2016 limit on SI cross section, while the red ones survive the limit. Note that the LUX-2016 limit on SD cross section for DM-neutron scattering can not exclude any Type I samples.
• Third, we picked up physical samples from the scan by requiring them to satisfy ∆ Z ≤ 50, ∆ h ≤ 50 and all the constraints contained in the package NMSSMTools-4.9.0 [31], such as various B-physics observables at 2σ level, the DM relic density at 2σ level, and the LUX-2015 limit on SI scattering rate [25]. We also considered the limitations from the direct searches for Higgs bosons at LEP, Tevatron and LHC by using the package HiggsBounds [42] and performed the 125 GeV Higgs data fit with the package HiggsSignal [43]. Moreover, we emphasize here that the constraints from various searches for SUSY at LHC Run-I on the samples were also implemented by detailed Monte Carlo simulations, which is the core of our previous work [19].
According to the analysis in [19], the physical samples can be classified into four types: for Type I samples, h 1 corresponds to the SM-like Higgs boson andχ 0 1 is Bino-dominated, while for Type II, III and IV samples, h 2 acts as the SM-like Higgs boson withχ 0 1 being Bino-, Singlino-and Higgsino-dominated respectively. Among the four types of the samples, the lowest fine-tuning comes from Type III and Type IV samples, for which ∆ Z and ∆ h may be as low as about 2 and therefore they are of particular interest to us. Now we show the impact of the LUX-2016 results on the samples. In Fig.2, we project Type I samples on σ SĨ χ−p − mχ0 1 plane (left panel) and σ SD χ−n − mχ0 1 plane (right panel) respectively. The samples marked by dark green color are excluded by the LUX-2016 limit on SI cross section, while those marked by red color survive the limit. From the figure one can learn that the SI cross sections of the samples are usually larger than 3 × 10 −47 cm 2 , which is within the detection sensitivity of future XENON-1T experiment [44]. One can also learn that among the total 5263 Type I samples, 68.7% of them have been excluded by the LUX-2016 limit on SI cross section, and by contrast the LUX-2016 limit on SD cross section for DM-neutron scattering can not exclude any Type I samples. The underlying reason for these features can be inferred from the formulae that 4 So for moderate light µ which is required to predict a small ∆ Z , the SI cross section is rather large (in comparison with its LUX-2016 limit) given no strong cancelation between h 1 (the SM-like Higgs boson) contribution and the other contributions. On the other hand, because the coupling between DM and Z boson is proportional to m 2 Z /µ 2 cos 2β and thus suppressed in comparison with the h 1χ 0 1χ 0 1 coupling, the SD cross section is not significant with respect to its experimental limit. We remind that in the scan, we do not include any information about the DD experiment in the likelihood function. Although the exact percentage of excluded samples may vary with different scan strategies, the number shown above can still reflect the powerfulness of the LUX-2016 results on Type I samples.
In Figs. 3, 4 and 5, we project Type II, III and IV samples respectively on σ SĨ χ−p − mχ0 1 and σ SD χ−n − mχ0 1 planes. In these figures the black samples are excluded by both the LUX-2016 limit on SI cross section and that on SD cross section for DM-neutron scattering. By contrast the dark green ones and the blue ones are excluded only by either of the limits,   Fig.3, but for Type IV samples. and the red ones survive all the limits. The percentages of each kind of colored samples in the total Type II samples are presented by the disc on the right panel of Fig.3, and so on for the other types of samples. Compared with Type I samples, Type II, III and IV samples exhibit following new features • The SI cross sections for the three types of samples may be much lower than the LUX-2016 limit. The underlying reason for such a behavior is that the singlet dominated h 1 is light, and thus it is able to cancel rather strongly the other h i contributions. However, one should note that for the more general case without fine parameter structure of blind spots, the SI bounds still provide the stronger constraints on the nNMSSM parameter space.
• The LUX-2016 limit on SD cross section can be used to exclude the samples obtained in [19]. Especially for Type III and IV samples, the exclusion capability becomes so strong that more than 80% of the two type samples are physically forbidden. The reason is as follows. For Type II samples, the Zχ 0 1χ 0 1 coupling is similar to that for Type I samples, but the parameter µ for Type II samples may take a significantly lower value (see Table I in [19]) and thus the corresponding SD cross section is enhanced. For Type III and IV samples, the explicit expression of the SD cross section in the heavy gaugino limit is presented in Eq.(2.18), and it is also due to the relatively smallness of µ that the exclusion capability of the SD cross section becomes strong.
About our results presented in above figures, we have the following comments: • In putting the constraints of the LUX-2016 results on the samples, we directly used the upper bounds on the SI and SD rates in [5,6] as inputs. However, from experimentalist point of view, the standard procedure to present the results of a given DM experiment is as follows [30]: First, one of the interaction (SI or SD) is neglected in order to draw conclusions for the other (SD or SI), which is referred to as a pure-SD (resp. -SI) case; then the SI couplings of DM with proton and neutron are assumed to be equal to get the averaged SI limit; finally, for the SD sector the interaction of DM with one type of nucleons is assumed to dominate and thus the SD results can be presented in two independent planes: the pure-proton case (equivalently corresponding to S n = 0 with S n denoting the spin content of neutron in target nuclei) and the pure-neutron one (corresponding to S p = 0).
It is obvious that these assumptions are not pertinent to our model where both σ SI and σ SD χ−p,n are non-zero simultaneously and the LUX experiment where both proton and neutron contribute to the spin of xenon nuclei. In this aspect, we note that an improved way to extract constraints from DD experiments was introduced in [30]. From Fig.3 in [30] and also the relation | S n | | S p | for xenon nuclei, we infer that for most cases, the upper bounds on σ SĨ χ−p and σ SD χ−n provided in [5,6] properly reflect the capability of the experimental data in constraining theory. This testifies the reasonableness of what we did in discussing the constraints of the LUX-2016 results on nNMSSM.
• We note that the LUX-2016 limit on SD cross section is actually based on the analysis of the experimental data collected in 2013. Given that the LUX-2016 limit on SI cross section is about 4 times stronger than the corresponding LUX-2013 result [6, 47], we infer from Eq.(13) of [30] that once the data underlying the LUX-2016 limit are analyzed for SD cross section, the upper bound on σ SD χ−n should be lowered by an approximate factor of 4. If this is true, Type II, III and IV samples in nNMSSM will be further limited.
• From the monojet analyses performed by ATLAS and CMS collaborations [48,49], it seems that the monojet searches at the LHC Run-I are able to put much tighter constraint on DM interactions than the SD limit. Motivated by this observation, we repeat the analyses of the two collaborations by carrying out detailed simulations, and surprisingly we do not find any constraint from the monojet searches on nNMSSM (we note that similar conclusion was obtained in [50]). The difference, we imagine, originates from the assumption on the mediator of the signal. To be more specific, ATLAS considered the process pp → jV → jχχ whereχ denotes a DM particle, and V represents a vector boson with weak couplings to quark and DM pair. As can be seen from Fig.7 in [48], strong constraint can be set only when m V > 2mχ so that the vector boson V is actually produced on-shell and subsequently decays into the DM pair. By contrast, in the NMSSM the monojet signal proceeds from the process pp → jZ ( * ) → jχ 0 1χ 0 1 . If the Z boson in the process is on shell, the upper bound from Z invisible decay set by LEP-I experiments has required the monojet plus E miss T signal sufficiently small, while if the Z boson is produced off-shell, the Z propagator suppression makes the cross section significantly lower than that in the ATLAS analysis. As far as the CMS analysis is concerned, it assumed the contact interaction of DM with quarks (see Fig.5 in [49]) and consequently some kinematic quantities such as E miss T distribute in a way quite different from those of the NMSSM. This will affect greatly the cut efficiencies in the analysis.
• Since the LUX-2016 results have required the SI cross section to be lower than about 10 −45 cm 2 , it is quite general that different h i contributions in Eq.(2.12) cancel each other to get the future measurable cross section, which will introduce another fine tuning problem [51][52][53][54]. To investigate this problem, we define the quantity ∆ DM = max (f p h i /f p ) 2 to measure the extent of the fine tuning. We find that among the samples surviving the LUX-2016 constraints, more than 40% of them predict ∆ DM > 50. This fact again reflects the strong limitation of the latest DD experiment on nNMSSM.

Status of nNMSSM after the LUX-2016 results
In this section, we investigate the status of nNMSSM after the LUX-2016 results. In order to get more physical samples for study, we repeat the scan introduced in Section 3 by including the LUX limits in the likelihood function as what was done in [40]. In the renewed Markov Chain scan, we find that it becomes rather difficult to get Type II, III and IV samples, and that the surviving samples are distributed in some isolated parameter islands. These facts reflect the strong limitation of the LUX-2016 results on nNMSSM.
In Fig.6, we show the impacts of the LUX-2016 limitations on ∆ Z and ∆ h for Type I samples. Samples on the left panel are obtained without considering the LUX-2016 constraints, while those on the right panel satisfy the constraints. From this figures, one can learn that if the LUX-2016 limits are not considered, ∆ Z and ∆ h can be as low as about 5, while the LUX experiment has pushed them up to be larger than 10. One can also infer that the experiment has excluded samples with µ 230GeV, and by contrast the electroweakino searches at LHC Run-I can not do this. Moreover, we remind from   Fig.2 that if the future XENON-1T experiment does not detect any DM signal, ∆ Z and ∆ h for Type I scenario must be larger than 50, and correspondingly the lower bound of the parameter µ will be larger than 400GeV.
In Figs. 7, 8 and 9, we show similar graphs to Fig. 6 for Type II, III and IV samples respectively. These figures indicate that even after considering the LUX-2016 limits, ∆ Z and ∆ h can still be as low as about 2. This implies that until now the NMSSM is still a viable theory in naturally predicting the quantities.
Since Singlino-dominated or Higgsino-dominated DM is peculiar to the NMSSM, in the following we concentrate on the Type III and IV samples and study the annihilation of DM in early universe. For this purpose, we project the surviving samples on the λ − κ plane  with different color representing the dominant annihilation mechanism for each sample. The corresponding results are shown in Fig.10 with the left panel for Type III samples and the right panel for Type IV samples. This figure indicates that for the Type III samples, DM may annihilate via the co-annihilation introduced in section II, large Higgsino and Singlino mixing or resonant funnel to get its measured relic density, while the Type IV samples achieve the relic density only through the sizable mixing. This figure also indicates that for the Type III samples, there are points with λ < 0.1. We examined the property of these samples and found that they are characterized by v s > 1 TeV while basically all singletdominated particles are lighter than about 250 GeV. Obviously, such a configuration of the Higgs potential is somewhat unnatural since the vacuum expectation value v s is much larger than the related scalar masses. Interestingly, we checked the vacuum stability of the potential by the package Vevacious [55,56] and found that the vacuums are usually long lived.

Conclusions
In this work we studied the constraints from the latest LUX results including both SI and SD measurements on natural NMSSM (nNMSSM). Since this theoretical framework usually favors light Higgsino mass µ which can induce a sizable rate for the neutralino DM-nucleon scattering, the constraints are expected to be rather tight. Our main observations include: • The SD bound is complementary to the SI bound in limiting nNMSSM since they have different dependence on SUSY parameters. Especially for the blind spots where the SI cross section vanishes due to strong cancelations among different contributions, only the SD DM-nucleon scattering contributes to the DM signal in the DD experiments and is therefore vital for the detection. We note that since the LUX-2016 experiment has set an upper bound on the SI cross section at about 10 −45 cm 2 , most samples that survive the bound should lie near the blind spots.
• For the peculiar scenarios of the NMSSM where the next-to-lightest CP-even Higgs corresponds to the 125 GeV Higgs boson discovered at the LHC, more than 85% of the samples obtained in our random scan are excluded by the latest LUX results. Although the exact percentage may vary with different scan strategies, this number can exhibit the high sensitivity of the nNMSSM scenarios to the DD experiment. By contrast, the constraint from the monojet searches at LHC Run-I on nNMSSM is rather weak due to the small production cross section, which is constrained by the invisible Z boson width or suppressed by the off-shell Z propagator.
• Quite distinctively, the SD bound is much more powerful than the SI bound in excluding SUSY parameter for Type III and Type IV samples in our study. This is opposite to the common impression that the SI bound is stronger than the SD one. However, one should also note that for the more general case without fine parameter structure of blind spots, the SI bounds still provide the stronger constraints on the nNMSSM parameter space.
• Low fine tuning samples are strongly limited within some isolated regions of the NMSSM parameter space and difficult to obtain. Future dark matter direct search experiments such as XENON-1T will provide a better test of nNMSSM.
Finally, we remind that the LUX-2016 limit on SD cross section is actually based on the analysis of the experimental data collected in 2013. Given that the LUX-2016 limit on SI cross section is about 4 times smaller than the corresponding LUX-2013 result [6,47], it is expected that, once the data underlying the LUX-2016 limit are analyzed for SD cross section, the upper bound on σ SD χ−n might also be lowered by an approximate factor of 4. In this case, the constraint from the SD bound will become stronger than what we obtained in this work.