Current Status of Natural NMSSM in Light of LHC 13TeV Data and XENON-1T Results

In the natural realization of the Next-to-minimal Supersymmetric Standard Model, Higgsinos tend to be lighter than about several hundred GeVs, which can induce detectable leptonic signals at the LHC as well as large DM-nucleon scattering cross section. We explore the constraints from the direct searches for electroweakino and slepton at the LHC Run II and the latest DM direct detection experiments on the scenario with low fine tuning indicator $\Delta_{Z/h} \leq 50$. We find that these experiments are complementary to each other in excluding the scenario, and as far as each kind of experiment is concerned, it is strong enough to exclude a large portion of the parameter space. As a result, the scenario with Bino- or Higgsino-dominated DM is disfavored, and that with Singlino-dominated DM is tightly limited. There are two regions in natural NMSSM parameter space surviving in the current experimental limits. One is featured with a decoupled Singlino-dominated LSP with $\mu \simeq m_{\widetilde{\chi}_1^0}$, which cannot be explored by neither DM detections or collider searches. The other parameter space region is featured by $10^{-47}~{\rm cm^2} \lesssim \sigma^{SI}_{\widetilde{\chi}-p} \lesssim 10^{-46}~{\rm cm^2}$ and the correlation $\mu \simeq m_{\widetilde{\chi}_1^0}$, which will be explored by near future DM detection experiments.


Introduction
It is well known that supersymmetric theories provide an elegant solution to the fine tuning problem in the Higgs sector of the Standard Model (SM), where the quadratically divergent contributions to the Higgs mass from the SM fermion loops are canceled exactly by those from corresponding sfermion loops due to supersymmetry, and consequently only relatively mild logarithmic contributions are left in the radiative correction [1,2]. This kind of theories also provide the possibility to unify different forces in nature and feasible dark matter (DM) candidates, which must be present in the universe to explain a large number of cosmological and astrophysical observations. Due to these advantages, supersymmetry has long been regarded as the footstone in building new physics models.
As the most economical realization of supersymmetry in particle physics, the Minimal Supersymmetric Standard Model (MSSM) is theoretically unsatisfactory due to its µ-problem and little hierarchy problem which was firstly discussed in [3,4] and became exacerbated in the last few years by the first run of LHC experiments, especially by the uncomfortable large mass of the discovered Higgs boson m h 125 GeV [5]. Alternatively, its gauge singlet extension called the Next-to-Minimal Supersymmetric Standard Model (NMSSM) has drawn a lot of attention since the first hint of the scalar appeared at the LHC [6][7][8][9][10][11]. In the NMSSM, the µ parameter is dynamically generated by the vacuum expectation value (vev) of the singlet Higgs superfield S, and since the field involves in the electroweak symmetry breaking, the magnitude of µ is naturally at weak scale [12,13]. Moreover, the interaction among Higgs fields λŜĤ u ·Ĥ d can lead to a positive contribution to the squared mass of the discovered Higgs boson, and if the boson corresponds to the next-to-lightest CP-even Higgs state, its mass can be further enhanced by the singlet-doublet Higgs mixing. These effects make the large radiative correction to the mass unnecessary and thus avoid the little hierarchy problem [6,8,10,14,15].
Compared with the MSSM, the introduction of the singlet fieldŜ has profound impacts on the phenomenology of the NMSSM, which is reflected in at least two aspects. One is that the scalar component fields ofŜ will mix with the doublet Higgs fields to form mass eigenstates. Consequently, the properties of the resulting Higgs bosons may deviate significantly from the MSSM predictions [10,16]. In particular, the model predicts singlet-dominated scalars, which can be rather light without conflicting with any experimental constraints [17,18], and they may act as the mediators or final states of DM annihilation [19], and/or as the decay product of heavy sparticles [20][21][22]. The other is that the involvement of the Singlino, the fermionic component field ofŜ, in the electroweakino sector usually extends the decay chain of sparticles [23][24][25]. This case along with the scenario of sparticle decay into the singlet-dominated scalars lead to complicated signal of sparticles at LHC. In the situation that most of the analyses in sparticle search performed by ATLAS and CMS collaborations which are designed for the MSSM, the constraints on the NMSSM can be much weaker [20][21][22][23][24][25]. Besides, due to the presence of the light singlet-dominated scalars and the self interaction of singlet fields κŜ 3 , the Singlino component in the lightest neutralino χ 0 1 makes it a more flexible DM candidate to escape the restriction from DM direct and indirect detection experiments in broad parameter space [19] as well as to explain exotic signals observed by DM experiments in certain scenarios [26][27][28][29]. All these novel features, therefore, necessitate a detailed study of any relevant parameter point in the NMSSM to see whether it is consistent with experimental data.
In the NMSSM, the Z boson mass is given by [30] where m H d and m Hu are the weak scale soft SUSY breaking masses of the Higgs fields H d and H u , d and u denote their radiative corrections, µ is the Higgsino mass and tan β = v u /v d with v u and v d being the vevs of the fields H u and H d . The equation indicates that, in order to get the observed value of Z boson mass m Z without resorting to large cancellations, each term on its right hand side should be comparable in magnitude to m Z . The extent of the comparability can be measured by the quantity [31] with p i denoting any Lagrangian parameter in the NMSSM. Obviously, the smaller value ∆ Z takes, the more natural the theory is in predicting m Z . On the other hand, any upper bound on ∆ Z has non-trivial requirements on the parameter space of the NMSSM, e.g. the Higgsino mass is restricted by µ 300 GeV and the lighter top squark is bounded by m t 1 3 TeV if ∆ Z < 30 [30]. In a similar way, one may define another independent quantitative measure of electroweak naturalness from the expression of the SM-like Higgs boson mass [32] In history, the scenario with ∆ Z O(10 2 ) is dubbed as Natural SUSY (NS) [33] or natural NMSSM so far as the explicit model NMSSM is concerned. In recent years with m h being measured more and more precisely, ∆ h is also considered in defining the NS [30,34]. As for the natural NMSSM scenario, it should be noted that the novel features mentioned above still hold, which make it differ greatly from the natural MSSM scenario. It should also be noted that the scenario prefers relatively light Higgsinos and scalar top quarks, and this preference can be tested at the LHC.
So far the parameter space of the natural NMSSM has been explored relentlessly by considering the constraints from the on-going collider experiments and DM direct and indirect detection experiments [24,25,. These studies indicate that, although the experiments are very effective in excluding the parameter points of the scenario, ∆ Z and ∆ h may still be as low as 2, and the property of the DM candidate is diverse, e.g. it may be either Bino-, Singlino-or Higgsino-dominated [54]. This situation, however, may be changed greatly since experimental search for the production of SUSY particles at LHC and DM-nucleon scattering in DM direct detection experiments has made considerable progress in the last years, which was emphasized in recent works [64,65]. For example, compared with the LHC Run I results, the LHC Run II data have pushed the mass limits on Wino-like χ ± 1 / χ 0 2 from 345 GeV to 650 GeV in simplified model with χ 0 1 = 0 GeV [66], and the recent XENON-1T experiment [67] has improved the sensitivity of the scattering rate by about three times in comparison with the results obtained in 2017 by LUX and PandaX-II experiments [68,69]. So in this work we update previous analyses on the natural NMSSM by including the latest experimental data, and we find that the parameter space of the scenario with ∆ Z/h ≤ 50 shrinks greatly, i.e. some cases become highly disfavored, while some remaining cases will be explored in near future. With the best of our knowledge, these conclusions are not obtained before. This paper is organized as follows: in Sec. 2 we introduce briefly the basics of the NMSSM, and present the results of our exhaustive scans over the scenario with ∆ Z/h ≤ 50 by considering various experiment constraints, including the search for sparticles at the LHC Run I. Then we show the impact of latest LHC and DM direct detection constraints on different cases in natural NMSSM in Sec. 3. The status of the scenario is discussed in Sec. 4. Finally, we draw our conclusion in Sec. 5.

Basics of the NMSSM
As the simplest extension of the MSSM, the NMSSM contains one extra gauge singlet Higgs fieldŜ with the superpotential and soft breaking terms given by [12]: where W F stands the MSSM superpotential without the µ-term,Ĥ u ,Ĥ d andŜ are Higgs superfields with H u , H d and S being their scalar components respectively, the dimensionless coefficients λ and κ parameterize the coupling strength in Higgs sector, and the dimensional quantities m 2 Hu,H d ,S and A λ,κ are soft breaking parameters. The Higgs potential of the model consists of the F -term and D-term of the superfields, as well as the soft breaking terms. After the electroweak symmetry breaking, the fields H u , H d and S acquire the vevs v u , v d and v s , and the soft breaking masses m 2 Hu , m 2 H d and m 2 S can be expressed in terms of v u , v d and v s through the minimization conditions of the scalar potential. In practice, the input parameters of the Higgs sector are usually chosen as instead of the soft masses. Moreover, it is more convenient to consider the field combinations H 1 = cos βH u + ε sin βH * d and H 2 = sin βH u + ε cos βH * d (ε is two-dimensional antisymmetric tensor) in discussion, which take the form [12]: with G + and G 0 corresponding Goldstone bosons and v 2 = v 2 u + v 2 d . In the basis (S 1 , S 2 , S 3 ), the 3 × 3 symmetric CP-even Higgs mass matrix M 2 is given by and consequently the mass eigenstate of CP-even Higgs bosons is h i = j V ij S j with V denoting the rotation matrix to diagonalize the mass matrix M 2 . In a similar way, one can get the CPodd mass eigenstates A 1 and A 2 . In the following, we take 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 the S 2 field. Without the mixing of the S i fields, the squared mass of SM-like Higgs boson gets an additional contribution λ 2 v 2 sin 2 2β in comparison with that of MSSM (see the expression of M 2 S 2 S 2 ), and it can be further enhanced by the mixing effect if M 2 S 3 S 3 < M 2 S 2 S 2 . Consequently, the SM-like Higgs boson does not need large radiative correction to get its measured mass value [10,14,15]. We remind that current experiments have very weak constraints on the S 3 /P 2 dominated scalars, and as a result, these particles may be as light as several GeVs.
At this stage, it is necessary to clarify that the parameter p i in Eq. (2) actually denotes the set of the parameters m 2 Hu , m 2 H d , m 2 S , λ, κ, A λ , A κ and Y t since by the definition of ∆ Z , m Z should be treated as a variable instead of a constant [31]. In this case, m Z , tan β and µ depend on the Lagrangian parameters by the minimization conditions, which enables one to get their derivatives in an analytic formula [31]. Similar treatment is applied to the calculation of ∆ h by noting that m h is related with the parameters by the secular equation det M 2 − m 2 h I 3 = 0 (I 3 denotes a 3 × 3 identity matrix) since m 2 h is one of the eigenvalues of the squared mass matrix M 2 , and the minimization conditions [32].
In the NMSSM, the Singlino S mixes with the gauginos (denoted by B and W respectively) and the Higgsinos H 0 d and H 0 u to form five neutralinos. In the basis of ψ = (−i B, −i W 3 , H 0 d , H 0 u , S), the symmetric mass matrix M 0 is given by [12] where M 1 and M 2 are soft breaking masses of Bino and Wino fields respectively, and g 1 and g 2 are SM gauge couplings. This matrix can be diagonalized by a rotation matrix N so that the mass eigenstates χ 0 i are given by where m χ 0 i is arranged in ascending order of mass, and thus χ 0 1 corresponds to DM candidate. The matrix element N ij measures the component of ψ j field in neutralino χ 0 i , and we call the DM to be ψ j dominated if N 2 1j is larger than the other components. Note that if any two of the five fields are decoupled, one can get the analytic forms of N 1j [26], which are useful to understand intuitively DM physics.
The properties of the other sparticles, such as their masses, are same as those predicted by the MSSM except that they may couple with the singlet fields, which may make their decay product quite complicated and thus increase degree of difficulty in probing them at the LHC [23][24][25]54]. As a result, the exclusion capability of the LHC on the parameter space of the NMSSM is usually weaker than that on the parameter space of the MSSM.

Features of Natural NMSSM
In order to show in detail the features of the natural NMSSM, we repeat the calculation of our previous works [54,57] to get more parameter points than what we obtained in these works. Roughly speaking, we first fix the soft breaking parameters for first two generation squarks and gluino mass at 2 TeV, set a common value M for all soft breaking parameters in slepton sector, and assume M U 3 = M D 3 and A t = A b for the third generation squark section to decrease the number of free parameters. Then we scan by Markov Chain method the rest parameters as follows 0 < λ < 0.75, |κ| < 0.75, 2 < tan β < 60, 100 GeV ≤ M ≤ 1.2 TeV, 100 GeV ≤ µ ≤ 1 TeV, 50 GeV ≤ M A ≤ 2 TeV, |A κ | ≤ 2 TeV, with all the parameters defined at the scale of 1 TeV. In the calculation, the particle spectrum is generated by the package NMSSMTools [70,71], the DM relic density and its spin-independent (SI) and spin-dependent (SD) cross sections are computed with the package micrOMEGAs [72,73], and the likelihood function is taken same as that in [74] except that the limits of LUX-2016 for SI cross section [75] and LUX-2016 for SD cross section [76], instead of the limits of the latest XENON-1T results [67], are adopted since we are going to show the impact of the latest DM detection experiments on the scenario. Note that we take the convention M 2 > 0 in the scan and allow M 1 and κ to be either positive or negative. We keep µM 2 > 0 since this usually leads to a positive contribution from sparticle loops to muon anomalous magnetic moment, which is helpful to alleviate the discrepancy between the measured value of the moment and its SM prediction (we will discuss this issue later) [77]. Due to the differences, the parameter region considered in this work is much broader than that in [54,57].
We further refine the samples obtained in the scan by picking up those which satisfy ∆ Z ≤ 50, ∆ h ≤ 50 and all the constraints implemented in the NMSSMTools, including various Bphysics observables in corresponding experimentally allowed range at 2σ level, DM relic density within ±10% around its measured central value Ωh 2 = 0.1187 [78] 1 , and the upper bounds of LUX-2016 on DM-nucleon scattering cross section at 90% confidence level (C.L.). We also consider the constraints from the direct search for Higgs bosons at LEP, Tevatron and LHC with the package HiggsBounds [79,80] and perform the 125 GeV Higgs data fit with the package HiggsSignals [81][82][83]. Moreover, we implement the constraints from various searches for SUSY at LHC Run I by following procedure: we firstly use the packages FastLim [84] and SModelS [85,86] to obtain preliminary constraints, and then use the package CheckMATE [87][88][89] with all published analyses to limit the rest samples. The Monte Carlo events of relevant SUSY processes are generated by the package MadGraph5 aMC@NLO [90][91][92] with the package PYTHIA [93,94] for parton showering and hadronization.
The scan results before implementing the latest LHC Run II and DM direct detection experimental limits are presented in Fig. 1. In panel (a), we project the samples on the fine tuning indicators ∆ Z − ∆ h plane with colors indicating the value of Higgsino mass µ. One can see that ∆ Z and ∆ h can be as low as about 1.7, and ∆ Z/h ≤ 50 set an upper limit of 547 GeV on the Higgsino mass µ. This conclusion has been obtained in our previous works [54,57], where we aimed to emphasize the importance of the LHC Run I results and DM direct detection results in limiting the scenario. Moreover, in [57] we classified the surviving samples by the dominant component of DM into three types, i.e. Bino-, Singlino-and Higgsino-dominated DM respectively, and found that they show different behaviors to accommodate the constraints from DM detection experiments. In the following, we explore in more detail the features of these types of samples.

Bino-dominated DM
For this type of samples, the DM annihilated mainly through three channels to get its measured relic density, which are In this case, the annihilation cross section is given by [12] σ( χ 0 where X and X denote SM particles, is the coupling between χ 0 1 s and h 1 (Z) given by and f s (g s ) is the generic functions for h 1 (Z) funnel depending on the s-channel momentum and the involved masses [95]. If one further assumes that the Wino and Singlino fields decouple from the rest of the neutralino sector, N 12 , N 15 ∼ 0, the other component coefficients of the DM roughly satisfy [26] N 11 : N 13 : N 14 (m 2 This relation implies that |N 11 | ∼ O(1), given that tan β 1 and µ m χ 0 1 , and consequently C h i χ 0 1 χ 0 1 and C Z χ 0 1 χ 0 1 are suppressed by µ −1 and µ −2 respectively. Since the annihilation cross section in Eq. (9) must be moderately large to get right DM relic density, µ should be upper bounded by certain values for the two annihilation channels. In Fig. 1 (b), we show the surviving samples with Bino-dominated DM on m χ 0 1 − m χ ± 1 plane with colors indicating slepton mass m . From this figure, one can see clearly that µ 480 GeV and µ 440 GeV for the Higgs funnel and Z funnel region respectively (we checked that the lighter chargino is Higgsino dominated for m χ ± 1 400 GeV). This situation is similar to the case of MSSM, which, according to the recent study of [65], is strictly limited by the latest LHC Run II result for electroweakino searches.
corresponding to the SM-like Higgs boson. This annihilation cross section can be written as where C h i χ 0 j χ 0 1 is the coupling of h i with χ 0 1 χ 0 j state, and h s denotes an auxiliary function encoding the complex mass dependence [95]. We checked that only few samples in Fig. 1  • Co-annihilation with sleptons.
With the assumption of a common slepton mass scale m , this annihilation cross section depends only on m χ 0 1 and m [96]. As indicated by Fig. 1 (b), such coannihilation channel occurs over a broad range of m χ 0 1 from 40 GeV to 220 GeV, and numerical results show the difference of the two masses varying from 60 GeV to 5 GeV with the increase of m χ 0 1 . Moreover, we note that either h 1 (in most case) or h 2 may act as the SM-like Higgs boson in this case.

Singlino-dominated DM
For this type of samples, only h 2 (h 1 ) can act as the SM-like Higgs boson for m χ 0 1 150GeV (m χ 0 1 220GeV). The properties of the DM differ from those of the Bino-like DM in following aspects: • Besides the three channels for the Bino-dominated DM, the Singlino-dominated DM may also annihilate by the process χ 0 1 χ 0 1 → W + W − through t-channel exchange of a chargino χ ± 1 . This requires the mass splitting between χ ± 1 and χ 0 1 to be about 10 GeV for Higgsino-dominated χ ± 1 and about 45 GeV for Wino-dominated χ ± 1 , which is shown on Fig. 1 (c) for Singlino-dominated DM case. We note that for the Higgsino case, the co-annihilation of the Higgsinos with χ 0 1 is also important since the mass splitting is less than 10% [97,98].
• For the Singlino-dominated DM, the elements of the matrix N have following relationship [26,54]: in the limit of |µ| |M 1 |, |M 2 |. This implies that the DM has less Higgsino components than the Bino-dominated DM, and consequently the h χ 0 1 χ 0 1 and Z χ 0 1 χ 0 1 coupling strengths may be significantly smaller than those for the Bino-dominated DM case if m χ 0 1 , µ, λ and tan β are taken same values (see the expressions in Eq. (10)). That is why the Higgsino masses in Fig. 1 (c) are visibly smaller than that in Fig. 1 (b) for the funnel regions.
• Compared with the Bino-dominated DM case, we found more samples that the DM annihilate significantly by the channel χ 0 The underlying reason is that the Singlino-dominated DM prefers certain parameter space of the NMSSM, such as 2|κ| < λ and moderately light µ, so that h 2 prefers to be the SM-like Higgs boson for m χ 0 1 ∼ 100 GeV. By contrast, for Bino-like DM case h 1 prefers to be the SM-like Higgs boson.
• The DM may annihilate by a resonant singlet-dominated A 1 , which has long been considered as a viable annihilation mechanism in literature [99], but after considering the constraints such a case becomes rare. This fact can be understood from the sum rule for the masses of the singlet fields [61,63] M 2 0,55 and the approximation M 2 , which is valid for most cases. Then the resonant annihilation condition m A 1 2m χ 0 1 implies that Since m 2 h 1 > 0, the equation holds only when M 2 0,55 is significantly larger than m χ 0 1 , which can be achieved by a large λ to induce sizable mixing between Higgsinos and Singlino in the neutralino mass matrix. Such a parameter space then predicts a light h 1 as well as a singlino-dominated DM whose Higgsino component is also sizable. Obviously, this situation tends to predict a large DM-nucleon scattering rate, and is therefore limited by DM direct detection experiments [61,63,100,101]. In fact, we find that only when the DM mass lies within a range from 88 GeV to 122 GeV can the situation survive the constraint. Moreover, we note that some samples with m χ 0 1 < 10 GeV and the A 1 funnel as DM dominant annihilation channel are presented in [100]. We checked the properties of these samples and found that they have been excluded by the 3 + E miss T search at the LHC [102].

Higgsino-dominated DM
This scenarios is characterized by approximately degenerated Higgsinos and Singlino in mass, and consequently H 0 u , S and H 0 d components of the DM are comparable in magnitude with the largest one coming from the H 0 u [54]. The non-negligible Singlino component N 2 15 , which is around 30%, can dilute the interactions of the DM with other fields so that DM density can coincide with the measured value of WMAP and Planck experiments [54]. We checked that the main annihilation channels of the DM in early universe include χ 0 where h 2 always denotes the SM-like Higgs for this type of samples, as well as the co-annihilation with sleptons. As was shown in [57], this scenario is tightly restricted by LUX-2016 on SD cross section for DM-nucleon scattering, and only samples with m χ 0 1 80 GeV and tan β ∼ O(1) are experimentally allowed. Part of these features are illustrated in Fig. 1 (d) of this work as well as in Fig. 2 with n OSSF < 2 and 2 τ h yes ----- Table 1: The summarisation of signal region categories defined in the CMS search for electroweakinos in final state with two light leptons of the same charge or with three or more leptons [102]. "OSSF" "OS" and "SS" stand "opposite sigh same flavor", "opposite sign" and "same sign" leptons, respectively. τ h denotes tau-tagged jet. "yes" means the corresponding variable is used to category bins. All quantities with mass dimension are given in units of GeV.

Impact of LHC Run II and DM direct detection results
Due to the requirement on naturalness, the DM is bounded by m χ 0 1 < 440 GeV with either moderately light chargino or light slepton. This feature motivates us to study the direct searches for sleptons and neutralinos/charginos pair production at LHC Run II and DM direct detections.

Sparticle searches at LHC Run II
To implement the LHC Run II limits on the slepton and electroweakino of samples, we add the following LHC Run II experimental analyses to CheckMATE: • The CMS search for electroweakinos in final state with either two or more leptons of the same charge, or with three or more leptons [102]. In simple terms, the target processes of this analysis are pp → χ ± i χ 0 j with different decay models into 2/3/4 + E miss T final state. The decay models can be classified into light slepton scenario and heavy slepton scenario. In light slepton scenario, the dominated decay chain of neutralino is χ 0 i → → + − χ 0 1 with i > 1, and main decay chain of chargino is χ ± i → ν ± / ν ± → ± ν χ 0 1 . The mass of slepton m and the flavor of the slepton in the decay chain both directly affect the property of final state. In the heavy slepton scenario, decay models jj/ ± ν and h → will lead to two/threelepton final states. Here h refers to the 125 GeV SM Higgs boson. Our natural NMSSM samples cover both scenarios.
After passing the basic selections, the signal events are categorized into 158 bins which are summarized into 12 signal regions (SRs) categories shown in Tab. 1. The first SR category, SS, is designed to the compressed scenarios in which one of the leptons from the decay chain of neutralino can be very soft, and therefore requires 2 same-sign (SS) leptons. The SR categories requiring three reconstructed leptons can be further classified by the number of τ h candidate. For three-leptons final state without τ h , signal events with (without) an opposite-sign same flavor (OSSF) lepton pair are categorized into SR category SRA (SRB). For three-leptons final state with one τ h candidate, SRs are defined as SRC, SRD and SRE by the signal events with OSSF lepton pair, opposite-sign (OS) lepton pair and (SS) lepton pair respectively. The SRF requires two τ h candidates of three reconstructed leptons. The events with final state of four or more than four leptons are classified into SRG to SRK by the number of OSSF pair n OSSF and the number of τ h . They aim for the production of a Z boson or h Higgs boson in the decay chain, which finally decays into two light flavor leptons or two taus.
• The CMS searches for electroweakinos with compressed mass spectra using events including two soft OS leptons and missing transverse energy [103]. The analysis is conceived to provide sensitivity to the process pp → χ 0 2 χ ± 1 → ( χ 0 1 W * )( χ 0 1 Z * ) for mass differences between χ 0 2 and χ 0 1 (∆m) of less than 20 GeV, where Z * and W * stand virtual Z and W bosons. The analysis requires an OS pair of light leptons, moderate E miss T and at least one jet. No significant excess was reported in the 12 SRs defined based on dilepton invariant mass and E miss T . In simplified model, Wino-like χ 0 2 / χ ± 1 masses masses up to 230 GeV are excluded for ∆m of 20 GeV. This analysis should be sensitive to the Singlino-dominated DM annihilating through t-channel chargino.
• The CMS search for electroweakinos in events with a lepton, two b-tagged jets and significant imbalance in the transverse momentum [104]. This search targets the neutralino and chargino pair production pp → χ 0 2 χ ± 1 → ( χ 0 1 h)( χ 0 1 W ± ) with decay models h → bb and W → ν . The kinematic variables used in this analysis including E miss T , the invariant mass of the two b jets M bb , the transverse mass of the lepton-E miss T system M T and the contransverse mass variable where p b1 T and p b2 T are the transverse momenta of the tow b jets, and ∆φ bb is the azimuthal angle between the b jets pair. After requiring 90 GeV < M bb < 150 GeV, M T > 150 GeV and M CT > 170 GeV, two exclusive SRs of 125 GeV < E miss T < 200 GeV and E miss T > 200 GeV are defined to enhance sensitivity to signal models with different mass spectra. The results show no significant excess in the two SRs, and exclude Winolike χ 0 2 / χ ± 1 between 220 and 490 GeV when the χ 0 1 is massless in simplified model. This analysis should be sensitive to the Bino-dominated DM scenario in which Higgsino-like χ 0 2,3 can decay to χ 0 1 h with large branch ratios.
• The CMS search for electroweakinos in final states with two leptons consistent with a Z boson and E miss T [105]. This search is designed for both strong and electroweak SUSY production leading to the on-Z signature, by selecting events with exactly one OSSF lepton pair consistent with the Z boson mass, two non b-tagged jets consistent with the W boson mass and large E miss T . Two Electroweak-production on-Z SRs, HZ and VZ, were defined with the invariant mass of two jets M jj , the variable M T2 ( ) [106,107] using the two selected leptons and M T2 ( b b) using two combinations of one lepton and one b-tagged jet as the visible object. The SRs are then divided into bins in E miss T . The analysis excludes Wino-like χ 0 2 / χ ± 1 masses between approximately 160 and 610 GeV for massless χ 0 1 with decay branch ratios Br( χ ± 1 → W ± χ 0 1 ) = Br( χ 0 2 → Z χ 0 1 ) = 100%. Thus it is sensitive to both the Bino-dominated DM scenario and the Singlino-dominated DM scenario.
• The CMS search for sleptons in final states with one OSSF lepton pair, no jet and large missing transverse momentum [108]. This search is optimized on the production of selectron pair and smuon pair in simplified model that Br( → χ 0 1 ) = 100%. In order to suppress tt and W W backgrounds, the SR selects events with 20 GeV • The CMS search for stau in the semi-leptonic and all-leptonic final state [109]. This search is targeting for direct τ pair production process in final state with two different flavor leptons formed one OS pair, which could be divided into eµ, eτ and µτ channels. The kinematic variable used in this search to bin SRs include E miss T , M T2 and Dζ, where Dζ is defined as: here ζ is the bisector between the direction of the two leptons, p T ( 1 ) and p T ( 2 ) are the transverse momenta of two leptons. In this search, signal events are binned into 144 SRs. Since the data from collider are consistent the SM expectations, no mass point in direct τ production can be excluded. For a τ mass of 90 GeV and a χ 0 1 of 1 GeV with decay mode Br( τ → τ χ 0 1 ) = 100%, the 95% C.L. upper limit for direct τ pair production cross section is up to 0.66 pb.
• The CMS search for stau pair production in the all-hadronic final state [110]. This search examines events with two hadronically decaying τ leptons and large E miss T . In this search, the angle between two τ h candidates ∆φ(τ 1 , τ 2 ), M T2 (τ 1 , τ 2 ), E miss T and ΣM T are used in the signal selection criteria to reduce the SM background, where ΣM T = M T (τ 1 , p miss T ) + M T (τ 2 , p miss T ). Three exclusive SRs are used to improve the sensitivity towards signal models with different stau masses. This analysis is most sensitive to a scenario with left-handed stau of around 125 GeV and a massless χ 0 1 .
We have submitted the implementations of above analyses to the CheckMATE database. The validations of cut-follows can be found in the website and Appendix A, which shows that our simulations agree with the corresponding experimental results within a 20% uncertainty.
For the surviving samples described in Sec. 2.2, we generate MC events of following processes at 13 TeV LHC, using MadGraph5 aMC@NLO [90][91][92] with the package PYTHIA [93,94] for parton showering and hadronization. Although the cross section of slepton pair production is much smaller than the cross section of electroweakino pair production, two high p T leptons from sleptons decay are always appeared in the final state. Process pp → ν ν, for example, can provide a SS lepton pair with a large E miss T if sneutrino pair decay through ν ν → ( χ ± 1 ∓ )( χ ± 1 ∓ ), which sensitive to the SS category in analysis [102]. And then the events are passed into CheckMATE which includes Delphes-3.2.0 [111] for detector fast simulation. The cross section are normalized to NLO using PROSPINO2 [112].
We firstly use the R values obtained from CheckMATE to apply the constraints from above searches. Here R ≡ max{(S i − 1.96∆S i )/S 95 i,obs } for individual analysis, where S i is the number of simulated signal events in i th SR or bin of the analysis, ∆S i stands the uncertainty of S i and S 95 i,obs represents the 95% C.L. upper limit of the event number in the SR. The samples that the R value of any above analysis is larger than 1 are deemed to be excluded by searches at LHC Run II at 95% C.L. in the following text. Then we combine the first four CMS electroweakino searches [102][103][104][105] though CLs method [113] with RooStats [114], because the SRs of them are mutually exclusive [66]. We use the likelihood function described in [65] for the combination, in which relative uncertainties of signal event is assumed to equal 5% and covariance matrices are not included.

DM direct detection
Complementary to the LHC experiments, DM direct detection experiments can also limit tightly the natural NMSSM scenario by measuring the SI and SD cross section for DM-nucleon scattering. In the NMSSM with heavy squark limit, the dominant contribution to the SI scattering comes from t-channel exchange of CP-even Higgs bosons [115][116][117], and the cross section is expressed as [118] where (n) denotes nucleon, µ r is the reduced mass of DM and nucleon, and C h i χ 0 1 χ 0 1 (C h i nn ) represents the coupling of h i with DM(nucleon). Note that light Higgs boson mass appearing in the denominator of Eq. ( (17)) can enhance the SI cross section, while on the other hand the cancellation between the contributions of different Higgs boson can suppress greatly the cross section [118]. The SD cross section is induced by the exchange of Z boson, which is given by [57,119] where n and p denote neutron and proton respectively, C p ≈ 4.0 and C n ≈ 3.1 for the typical values of form factor f (n) q . So far the tightest bound on the SI and SD cross sections comes from the XENON-1T experiment in 2018 [67] and the LUX measurement of DM-neutron scattering in 2017 [120] respectively. Both experiments improve the limits adopted in [57] by more than six times, so we think it mandatory to update the constraints on the scenario discussed in [57] with the latest limits.

Numerical results
Now we study the impact of the LHC experiments and the DM detection experiments on the three types of samples in natural NMSSM scenario.
In Fig. 2 Since there is no sample surviving both the constraints, it is fair to say that, at least for the assumptions made in this work, the natural NMSSM scenario with Bino-dominated DM and ∆ Z/h < 50 is strongly disfavored by current experiments.
In order to show more details about the results, we also divide the samples into two cases by different symbols: those marked by dot denote the case of m < µ, and the others marked by triangle denote the case of m > µ. The difference of the cases is that for the former case, Higgsinos prefer to decay into slepton first, which can enhance the branching ratio for leptonic final state. With the division, one can infer from Fig. 2 following facts:  the funnel regions, which can suppress the leptonic signal of the dominant electroweakino production process pp → χ ± 1 χ 0 2 . The net result of these facts is that the CLs values of the samples are slightly larger than 0.05, which means that they are at the edge of being excluded by the LHC analyses at 95% C.L.. On the other hand, since the annihilation mechanisms set an upper bound on µ so that the coupling C h χ 0 1 χ 0 1 is not suppressed too much, the SI cross section is moderately large, and consequently these samples are excluded by the XENON-1T experiment.
• For most green color samples, χ ± 1 is Higgsino-dominated with m χ ± 1 250 GeV, and h 2 acts as the SM-like Higgs boson. In this case, the h 2 contribution to the SI cross section can be cancelled by the h 1 contribution to a great extent [118] so that the SI cross section may be as low as 10 −48 cm 2 . Moreover, the SD cross section may also be suppressed by the cancellation between |N 13 | 2 and |N 14 | 2 , which is reflected by Fig. 2 (d) and Eq. (18). We remind that it is actually a common case in the NMSSM with a moderately light µ that only one of the cross section is suppressed [57], and the rare situation that both the cross sections are suppressed simultaneously was recently discussed in [59].  similar to those of the green color samples except that the cancellation between |N 13 | 2 and |N 14 | 2 is not strong to result in a sizable SD cross section.
In Fig. 3, we illustrate the features of the Singlino-dominated DM case in a similar way to that of Fig. 2 with additional blue points standing for those which survive all the experimental constraints. From this figure, one can learn following facts: • All the samples with h 2 acting as the SM-like Higgs boson satisfy µ 300 GeV, and some of them also satisfy M 2 180 GeV or M 400 GeV. While for the samples with h 1 corresponding to the SM-like Higgs boson, they satisfy µ 450GeV with µ m χ 0 1 or m m χ 0 1 . These features entail following conditions for the samples to be consistent with the experimental constraints: moderately strong cancellation between the h 1 and h 2 contributions to the SI cross section, |N 13 | 2 |N 14 | 2 as well as the suppressed spectrum of the sparticles with χ 0 1 [54,57] 3 .
• Similar to the Bino-dominated DM case, samples featured by m χ 0 1 m Z/h /2 or µ > m are completely excluded by the current experimental limits. Constrains from the LHC electroweakino searches play critical roles.
• Nearly all the samples with an approximate degeneracy of Wino-like χ ± 1 with χ 0 1 in mass are excluded. Some of them may survive the constraints from the LHC experiments and  Table 2: Detailed information about two benchmark points P1 and P2 for Singlino-dominated DM case. All quantities with mass dimension are given in units of GeV.
the XENON-1T experiment, but are excluded by the measurement on the SD cross section. These samples correspond to the yellow color samples in the four panels of 200 GeV in Fig. 3 (a) and m < 400 GeV in Fig. 3 (b), the LHC experiments have no exclusion capability. The SI cross section is sizable in comparison its detection limit, which varies from 2 × 10 −47 cm 2 to 3 × 10 −46 cm 2 , while the SD cross section is suppressed too much to be less than 2 × 10 −42 cm 2 . For the green samples, they satisfy µ < m , and the constraints of the LHC Run II mainly come from the associated production of Wino-like chargino and neutralino.
• Most important, there exist samples that survive all the constraints, which correspond to the co-annihilation region of the DM with Higgsinos to get the right relic density, and are marked by blue color in Fig. 3. In addition, some of these samples may also co-annihilate with slepton, and consequently the mass splitting between the DM and the Higgsinos can be slightly larger in getting the right DM relic density. Compared with the green samples discussed above, sleptons and Wino-like neutralino/chargino are heavier to escape the constraints from LHC Run II. Moreover, we note that there are surviving samples with high singlet purity (N 2 15 > 0.99). In this case, the DM decouples with SM particles so that both SI and SD cross sections are lower than the future LZ detection limits [121]. This case was recently emphasized in [59]. For the other samples without such high singlet purity, the SI cross section may be at the order of 10 −47 cm 2 , which will be explored by near future DM direction detection experiments, and the SD cross section is usually less than 5 × 10 −43 cm 2 which is far below its current detection limits.
In order to emphasize the property of the samples with Wino dominated χ ± 1 and Higgsino dominated χ ± 1 in the co-annihilation region, we choose two benchmark points P1 and P2 with their detailed information presented in Tab. 2. Both the points pass the LHC constraints, but their behaviors confronted the DM detection limits are quite different: for the Wino dominated χ ± 1 case (point P1), the SI cross section is far below its detection limit while the SD cross section is around its detection limit, and the situation is reversed for the Higgsino dominated χ ± 1 case (point P2).
Finally we consider the Higgsino-dominated DM case. This kind of samples predict a light CP-even Higgs boson with m h 1 < 125 GeV, 70 GeV m χ 0 1 100 GeV, µ 160 GeV and moderately large mixing between Higgsino and Singlino in forming neutralino mass eigenstates [54]. In Fig. 4, we project the samples on different planes like what we did in Fig. 2. From this figure, one can learn following facts: • Although the mass splittings between χ 0 2/3 / χ ± 1 and χ 0 1 are relatively small, the quite large cross section of neutralino/chargino pair production leads to the exclusion of all the samples by the electroweakino searches described in Sec. 3.1. Some of the samples can also be excluded by the slepton searches at LHC.

Status of the natural NMSSM
The results in previous sections reveals that in the natural NMSSM scenario with ∆ Z/h ≤ 50, only the Singlino-dominated DM case can survive the tight experimental constraints if the correlation µ m χ 0 1 holds. This has non-trivial impact on the parameter space of the NMSSM and also on the fine tunings of the theory. In Fig. 5, we project the samples obtained in the scan on λ − κ plane and µ − tan β plane with the grey (blue) color samples being experimentally excluded (allowed). This figure indicates that, after considering the constraints, the scenario is restricted in certain narrow corners of the NMSSM parameter space, which is featured by λ/κ 2.5 with λ ≈ 0.02, 100 GeV µ 200 GeV and 8 tan β 32 for κ > 0 and by λ/κ −2.5 with λ 0.05, µ 460 GeV and 36 tan β < 60 for κ < 0. In Fig. 6, we show the fine tuning indicators of the scenario before and after considering the LHC Run II and DM detection results with different colors representing the values of µ (see the color bar on the right side of the figure). This figure shows again that the experimental constraints are very powerful in limiting the scenario and have reduced significantly the range of ∆ Z/h . Given the status of the natural NMSSM, it is interesting to ask following questions: 1. Since the DM relic density is another precisely measured quantity, what is the tuning needed to get its measured value?
2. Is the natural NMSSM scenario able to explain the discrepancy of muon anomalous magnetic moment?
3. What are the effects on the conclusions given above if one relaxes the requirement on the fine tuning measurements by ∆ Z , ∆ h ≤ 100?
4. What will happen if one takes the value of the DM relic density measured by Planck just as an upper bound?
In order to answer the first question, we define the fine tuning measurement of the density as where p i denotes the variables in Eq. (8), and present ∆ Ωh 2 for the surviving samples on m χ 0 1 − ∆ Ωh 2 plane in Fig. 7 with the color bar denoting the mass splitting between m χ ± 1 and m χ 0 1 . This panel indicates that for the samples in the co-annihilation region of the Singlino-dominated DM with Higgsinos, ∆ Ωh 2 35 which is insensitive to the DM mass, while for those in the co-annihilation region with sleptons, ∆ Ωh 2 can be as large as 95. We note that our results about ∆ Ωh 2 coincide with those in [122]. As for the second question, we categorize the surviving samples by whether they can explain the muon g − 2 anomaly at 2σ level or not, and present them on the ∆ Z − ∆ h plane of Fig. 7. The samples marked by blue color are able to explain the anomaly, while those marked by grey color fail to do so. This panel indicates that the explanation of the anomaly places additional restrictions on the scenario, and consequently, due to the shrink of the allowed parameter space from µ 100 GeV to µ 150 GeV, the lower bound on ∆ Z (∆ h ) is shifted from 2(2) to 7 (5). In getting the results, we use the default setting of the NMSSMTools to take into account the theoretical and experimental uncertainties of the anomaly. With respect to the third question, we note that relaxing the constraint on ∆ Z and ∆ h will allow the parameter µ to vary over a broader range since both fine tuning indicators are sensitive to µ. A larger µ can suppress the rate of electroweakino pair production at the LHC as well as the DM-nucleon scattering rate, which is helpful for the theory to escape the experimental constraints. Our results from an additional scan of the parameter space in Eq. (8) indicate that allowing ∆ Z/h ≤ 100 can increase the samples in the slepton co-annihilation region for the Binodominated DM case and the Higgsino co-annihilation region for the Singlino-dominated DM case without violating the constraints. The results also show that the Higgsino-dominated DM case is scarcely affected by relaxing the fine tuning measurements. Finally, we point out that taking a lower value of the density Ω h 2 is equivalent to relax the upper bound on the cross section of the DM-nucleon scattering by a factor (Ω h 2 )/0.1187, and consequently the constraints from the DM detection experiments are weakened. As far as the Bino-dominated DM case and the Singlinodominated DM case are concerned, a lower relic density can be achieved by narrowing the mass gap between the DM and its co-annihilating particles. This will not affect the constraints from the LHC experiments. Moreover, without the right relic density, the DM may be a pure Higgsino particle. In this case, its relic density is less than 0.01 [30,34], and its scattering with nucleon is suppressed greatly since there is no triple doublet-Higgs interaction in the superpotential of the NMSSM.
Before we end this section, we have following comments about our results: • In our discussion, we do not consider the constraints from the direct search for top squarks at the LHC Run II [123]. We checked that the surviving samples in Fig. 5 may predict the lighter stop mass as low as about 600 GeV, and part of those samples are sure to be tested by the search. This will further shrink the parameter space of the scenario.
• We note that the discovery potential for the electroweakino production process pp → χ ± 1 χ 0 2 at future high luminosity LHC has been estimated by ATLAS collaboration [124,125] and CMS collaboration [126] in trilepton and WH channels. Using the relevant analysis codes for ATLAS collaboration at 14 TeV LHC [125], which was provided by the package CheckMATE, we find the analysis has no exclusion capability for the surviving samples even for the luminosity as high as 3000 fb −1 .
• As we mentioned before, in order to satisfy the strong constraint of XENON-1T experiment on the SI cross section, the h 2 contribution must be cancelled greatly by the h 1 contribution. This induces another kind of fine tuning in DM physics which is different from the tuning in the electroweak symmetry breaking and was discussed in [36]. The In the left panel, the color indicates the mass difference between χ ± 1 and χ 0 1 , and in the right panel, the blue (gray) points stand for the samples which are able to (unable to) explain the discrepancy of the muon anomalous magnetic moment at 2σ level.
origin of the tuning comes from two aspects. One is that the Higgsino mass µ should be of O(10 2 GeV) to predict m Z in a natural way. Such a light µ can enhance the SI cross section greatly. The other is that the parameters in the Higgs sector have been tightly limited by the LHC search for Higgs bosons, and this determines the relative size of each h i contribution to the cross section [118]. Take the heavy doublet dominated Higgs boson as an example, its contribution to the SI cross section can be neglected safely in most cases since the search for extra Higgs bosons at the LHC has required its mass at TeV scale, which can suppress the contribution greatly.

Summary
In this work, we explore the constraints from the direct searches for electroweakino and slepton at the LHC Run II and the latest DM direct detection experiments on the natural NMSSM scenario for three types of samples, namely those with Bino, Singlino and Higgsino as dominant DM component respectively. We have following observations: • Moderately light Higgsinos are favored by this scenario, which usually results in detectable leptonic signals at the LHC as well as large DM-nucleon scattering rate. Moreover, in some cases Wino and sleptons with mass around several hundred GeVs are also predicted. This situation makes the scenario to be testable readily by the experiments, and surviving these experiments necessitates great cancellation among different Higgs contributions to the SI cross section of DM-nucleon scattering, |N 13 | 2 |N 14 | 2 and suppressed sparticle spectrum. This, on the other hand, induces a kind of tuning of the theory which is other than the fine tuning in electroweak symmetry breaking sector.
• The signal of the electroweakino/slepton pair productions at the LHC Run II and the SI and SD cross section for DM-nucleon scattering are sensitive to different parameter space of the NMSSM, and their constraints are complementary to each other in excluding the samples of the natural NMSSM scenario. As far as each kind of the experiments is concerned, its individual constraint is strong enough to exclude most samples of the scenario.
• With the assumptions made in this work, the samples with Bino-or Higgsino-dominated DM are completely excluded by the experiments, and most samples for Singlino-dominated DM case are also excluded. As a result, some input parameters of the natural NMSSM scenario are restricted in certain narrow corners of the NMSSM parameter space.
• Although future LHC experiments and DM detection experiments can further limit the parameter space of the natural NMSSM scenario, there exist special parameter regions where the Singlino-dominated DM decouples from the SM sector. In this case, neither LHC experiments nor DM direct detection experiments can probe the scenario.
In summary, given the tight experimental constraints on the natural NMSSM scenario, its charm is fading, and one may either accept the current situation of the theory or insist on the fine tuning criteria as a guidance of new physics to construct more elaborated theories. For the latter choice, the seesaw extensions of the NMSSM, which is motivated by neutrino mass, provide an economical solution to the problem of the strong constraints by choosing the lightest sneutrino as the DM candidate [74,122,127]. As was shown in [74,122], a moderately light µ in this framework is favored not only by predicting naturally Z boson mass, but also by predicting right DM physics. The signals of sparticles at the LHC may be quite different from those in the MSSM or NMSSM, which is helpful to evade collider constraints [74,127].  Table 3: Cut-flow validation for signal region categories SRA and SRB in analysis [102]. The yields in "3 tight e, µ or τ h " of "CheckMATE" are normalized to "3 tight e, µ or τ h " of "CMS". "efficiency" is defined as the ratio of the event number passing though the Cut-flow to the event number of the previous one.