Natural NMSSM after LHC Run I and the Higgsino dominated dark matter scenario

We investigate the impact of the direct searches for SUSY at LHC Run I on the naturalness of the Next-to-Minimal Supersymmetric Standard Model (NMSSM). For this end, we first scan the vast parameter space of the NMSSM to get the region where the fine tuning measures $\Delta_Z$ and $\Delta_h$ at the electroweak scale are less than about 50, then we implement by simulations the constraints of the direct searches on the parameter points in the region. Our results indicate that although the direct search experiments are effective in excluding the points, the parameter intervals for the region and also the minimum reaches of $\Delta_Z$ and $\Delta_h$ are scarcely changed by the constraints, which implies that the fine tuning of the NMSSM does not get worse after LHC Run I. Moreover, based on the results we propose a natural NMSSM scenario where the lightest neutralino $\tilde{\chi}_1^0$ as the dark matter (DM) candidate is Higgsino-dominated. In this scenario, $\Delta_Z$ and $\Delta_h$ may be as low as 2 without conflicting with any experimental constraints, and intriguingly $\tilde{\chi}_1^0$ can easily reach the measured DM relic density due to its significant Singlino component. We exhibit the features of the scenario which distinguish it from the other natural SUSY scenario, including the properties of its neutralino-chargino sector and scalar top quark sector. We emphasize that the scenario can be tested either through searching for $3 l + E_T^{miss}$ signal at 14 TeV LHC or through future DM direct detection experiments.


Introduction
In the supersymmetric models such as the Minimal Supersymmetric Standard Model (MSSM) [1,2] and the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [3,4], the Z boson mass is given by [5] where m 2 H d and m 2 Hu represent the weak scale soft SUSY breaking masses of the Higgs fields H d and H u respectively, Σ d and Σ u are their radiative corrections, µ is the Higgsino mass and tan β ≡ v u /v d . As was shown in [6], the corrections Σ d and Σ u can be obtained from the effective Higgs potential at loop level, and in case of a large tan β, their largest contributions arise from the Yukawa interactions of third generation squarks, which are In above formulae, Q denotes the renormalization scale in getting the effective potential, and its optimized value is usually taken as Q = √ mt 1 mt 2 witht 1 andt 2 being the light and heavy top squarks (stop) respectively. Obviously, if the observed value of m Z is obtained without resorting to large cancelations, each term on the right hand side of Eq. (1.1) should be comparable in magnitude with m 2 Z , and this in return can put non-trivial constraints on the magnitudes of µ and mt 1,2 . Numerically speaking, we find that requiring the individual term to be less than 10m 2 Z leads to µ and mt 1,2 upper bounded by about 200 GeV and 1.5 TeV respectively. In history, the scenario satisfying the bounds is dubbed as Natural SUSY (NS) [5].
In the MSSM, the NS scenario is theoretically unsatisfactory due to at least three considerations. First, since the Higgsino mass µ is the only dimensionful parameter in the superpotential of the MSSM, its typical size should be of the order of the SUSY breaking scale. Given that the LHC searches for supersymmetric particles have pushed the masses of gluinos and first generation squarks up to above 1TeV [7,8], µ 200GeV seems rather unnatural. Second, the relic density of the dark matter (DM) predicted in the NS scenario is hardly to coincide with its measured value. Explicitly speaking, it has been shown that if the DMχ 0 1 is Higgsino-dominated 1 , its density is usually about one order smaller than its measured value [5], alternatively if it is Bino-dominated, the correct density can be achieved only in very limited parameter regions of the MSSM [9]. These features make the NS scenario disfavored by DM physics. Third, the NS scenario is further exacerbated by 1 Throughout this work, we denote the mass eigenstates of the neutralinos byχ 0 i with i ranging from 1 to 4 (5) for MSSM (NMSSM), and assume an ascending mass order for theχ 0 i by convention. With such an assumption, the lightest neutralinoχ 0 1 as the lightest supersymmetric particle (LSP) is regarded as the DM candidate.
the uncomfortably large mass of the recently discovered Higgs particle [10,11]: its value m h 125GeV lies well beyond its tree-level upper bound m h ≤ m Z , and consequently stops heavier than about 1TeV must be present to provide a large radiative correction to the mass [12][13][14][15][16][17]. This requirement seems in tension with the naturalness argument of Eq. (1.1). In fact, all these problems point to the direction that the NS scenario should be embedded in a more complex framework. Remarkably, we note that the NMSSM is an ideal model to alleviate these problems.
The NMSSM extends the MSSM by one gauge singlet superfieldŜ, and it is the simplest SUSY extension of the Standard Model (SM) with a scale invariant superpotential (i.e. its superpotential does not contain any dimensionful parameters) [3,4]. In this model, the Higgsino mass µ is dynamically generated by the vacuum expectation value ofŜ, and given that all singlet-dominated scalars are lighter than about v 174GeV, its magnitude can be naturally less than 200GeV. These additional singlet-dominated scalars, on the other hand, can act as the mediator or final states of the DM annihilation [18], and consequently the NS scenario in the NMSSM with a Singlino-dominated DM can not only predict the correct relic density, but also explain the galactic center γ-ray excess [18][19][20][21]. Moreover, in the NMSSM the interaction λŜĤ u ·Ĥ d can lead to a positive contribution to the squared mass of the SM-like Higgs boson, and if the boson corresponds to the next-to-lightest CP-even Higgs state, its mass can be further lifted up by the singlet-doublet Higgs mixing. These enhancements make the large radiative correction of the stops unnecessary in predicting m h 125GeV, and thus stops can be relatively light [22][23][24][25][26][27][28][29].
So far studies on the NS scenario in the NMSSM are concentrated on the assumption thatχ 0 1 is Singlino-dominated [18][19][20][21][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46]. In this case, the branching ratio of the golden channelt 1 → tχ 0 1 in the LHC search for a moderately light stop is highly suppressed. Instead,t 1 mainly decays into the Higgsino-dominatedχ + 1 andχ 0 2,3 in following way [43] where X 0 denotes either Z boson or a neutral Higgs boson. These lengthened decay chains can generate softer final particles in comparison with the golden channel, and consequently weaken the LHC bounds in the stop search. This feature is also applied to other sparticle searches, and it has been viewed as an advantage of the NMSSM in circumventing the tight constraints from the LHC searches for SUSY. In this work, we consider another realization of the NS scenario whereχ 0 1 is Higgsino-dominated. In our scenario, the Higgsinos and the Singlino are degenerated in mass at 50% level, and consequently they mix rather strongly to form mass eigenstatesχ 0 1,2,3 withχ 0 1,2 andχ 0 3 being Higgsino-dominated and Singlinodominated respectively. Since the role of the Singlino component inχ 0 1 is to decrease the DM annihilation rate,χ 0 1 may achieve the relic density measured by Planck and WMAP experiments [47,48] without contradicting the DM direct search experiments such as LUX [49,50]. The phenomenology of our scenario is somewhat similar to that of the popular NS scenario in the MSSM, which was proposed in [5], but our scenario has following advantages • In the parameter regions allowed by current LHC searches for SUSY, it may have a lower fine tuning in getting the Z boson mass. Meanwhile, it has broad parameter regions to predict the right relic density (see our discussions in Sect. III).
• Since all the light particles in our scenario, i.e.χ ± 1 ,χ 0 1 ,χ 0 2 andχ 0 3 , have sizableH u component, the main decay modes oft 1 includẽ and each mode corresponds to different signals. Because the signal oft 1 pair production is shared by rich final states, the LHC bounds on mt 1 are usually weakened. Moreover, we remind that the phenomenology of our scenario is different from that of the NS scenario with a Singlino-dominated DM. This can be seen for example from the decay modes oft 1 , which are presented in Eq.(1.2) and Eq.(1.3) respectively. We also remind that our scenario was scarcely discussed in literatures. In fact, within our knowledge only the work [52] briefly commented thatχ 0 1 may be Higgsino-dominated in the constrained NMSSM.
This work is organized as follows. In Sec. II, we briefly recapitulate the framework of the NMSSM, then we scan its parameter space by considering various constraints to get the NS scenarios in the NMSSM. Especially, we take great pains to implement the constraints from the LHC searches for SUSY by multiple packages and also by detailed Monte Carlo simulations, like what the work [53] did. After these preparation, we exhibit in Sec. III the features of the NS scenario with a Higgsino-dominated DM, including its favored spectrum and the properties of the neutralinos and stops, and subsequently in Sec. IV we take several benchmark points as examples to show the detection of our scenario in future experiments. Finally, we draw our conclusions in Sec.V. The details of our treatment on the LHC searches for SUSY are presented in the Appendix.

The Structure of the NMSSM
The NMSSM extends the MSSM by adding one gauge singlet superfieldŜ, and since it aims at solving the µ problem of the MSSM, a Z 3 discrete symmetry under which the Higgs superfieldsĤ u,d andŜ are charged is adopted to avoid the appearance of dimensionful parameters in its superpotential. Consequently, the superpotential of the NMSSM can be written as [3] where W F is the superpotential of the MSSM without the µ-term, and the dimensionless parameters λ, κ describe the interactions among the Higgs superfields. The Higgs potential of the NMSSM is given by the usual F-term and D-term of the superfields as well as the soft breaking terms, which are given by where ε is second-order antisymmetric tensor with ε 12 = −ε 21 = 1 and ε 11 = ε 22 = 0, and tan β ≡ v u /v d with v u and v d denoting the vacuum expectation values of H u and H d respectively. In this representation, the redefined fields H i (i = 1, 2, 3) are given by These expressions indicate that the field H 2 corresponds to the SM Higgs doublet with G + and G 0 denoting the Goldston bosons eaten by W and Z bosons respectively, and the field H 1 represents a new SU (2) L doublet scalar field which has no tree-level couplings to the W/Z bosons. These expressions also indicate that the Higgs sector of the NMSSM includes three CP-even mass eigenstates h 1 , h 2 and h 3 which are the mixtures of the fields S 1 , S 2 and S 3 , two CP-odd mass eigenstates A 1 and A 2 which are composed by the fields P 1 and P 2 , as well as one charged Higgs H + . In the following, we assume m h 3 > m h 2 > m h 1 and m A 2 > m A 1 , and call the state h i the SM-like Higgs boson if its dominant component comes from the field S 2 .
In the NMSSM, the squared mass of the filed S 2 is given by where the last term on the right side is peculiar to any singlet extension of the MSSM [3], and its effect is to enhance the mass of the SM-like Higgs boson in comparison with the case of the MSSM. Moreover, if the inequation m 2 S 3 S 3 < m 2 S 2 S 2 holds, the mixing of the field S 2 with the field S 3 in forming the SM-like Higgs boson can further enhance the mass. In this case, h 1 is a singlet-dominate scalar, h 2 acts as the SM-like Higgs boson, and due to the enhancement effects the requirement m h 2 125GeV does not necessarily need the large radiative correction of stops [22,25]. We remind that the singlet-dominated physical scalars (i.e. the mass eigenstates mainly composed by S 3 and P 2 respectively) are experimentally less constrained, and in case that they are lighter than about 200GeV, µ = λv s naturally lies within the range from 100GeV to 200GeV.
In practice, it is convenient to use [3] λ, κ, tan β, µ, M A , A κ , 2) are traded for m Z , tan β and µ by the potential minimization conditions, and A λ is replaced by the squared mass of the CP-odd field P 1 , which is given by Note that M A represents the mass scale of the new doublet H 1 , and it is preferred by current experiments to be larger than about 300GeV. When we discuss the naturalness of the NMSSM, we consider two fine tuning quantities defined by [54] where h represents the SM-like Higgs boson, p i denotes SUSY parameters at the weak scale, and it includes the parameters listed in Eq.(2.5) and top quark Yukawa coupling Y t with the latter used to estimate the sensitivity to stop masses. Obviously, ∆ Z (∆ h ) measures the sensitivity of m Z (m h ) to SUSY parameters at weak scale, and the larger its value becomes, the more tuning is needed to get the corresponding mass 2 . In our calculation, we calculate ∆ Z and ∆ h by the formulae presented in [55] and [56] respectively. The NMSSM predicts five neutralinos, which are the mixtures of the fields BinoB 0 , WinoW 0 , HiggsinosH 0 d,u and SinglinoS 0 . In the basis If |M 1 |, |M 2 | |µ| and 2κ/λ ∼ 1, the Bino and Wino fields are decoupled from the rest fields. In this case, the remaining three light neutralinosχ 0 i (i = 1, 2, 3) can be approximated byχ where the elements of the rotation matrix N roughly satisfy with mχ0 i denoting the mass ofχ 0 i . In the following, we are interested in the parameter region featured by tan β ∼ 10, 2κ/λ ∼ (1 − 1.5) and |µ| ∼ (100 − 200)GeV, which is hereafter dubbed as the NS scenario withχ 0 1 being Higgsino-dominated. From Eq.(2.9) and Eq.(2.10), one can conclude that this scenario has following characters • The lightest two neutralinosχ 0 1 andχ 0 2 are Higgsino-dominated, andχ 0 3 is Singlinodominated. Their masses should satisfy following relations: mχ0 1 < |µ|, mχ0 2 ∼ |µ| and mχ0 3 > | 2κ λ µ| > |µ|.
• As far asχ 0 1 is concerned, its largest component comes fromH 0 u field. If the splitting between mχ0 1 and |µ| is significant, its Singlino component may also be quite large. The importance of the Singlino component is that it can dilute the couplings ofχ 0 1 with W and Z bosons, Higgs scalars and SM fermions, and consequently the density ofχ 0 1 can coincide with the DM density measured by WMAP and Planck experiments.

Strategy in Scanning the Parameter Space of the NMSSM
In this part, we perform a comprehensive scan over the parameter space of the NMSSM by considering various experimental constraints. Especially, we take great pains to implement the constraints from the direct searches for SUSY at the LHC. After this procedure, we get the NS scenario with a Higgsino-dominatedχ 0 1 . We begin our study by making following assumptions about some unimportant SUSY parameters: • We fix all soft breaking parameters for the first two generation squarks at 2TeV.
Considering that the third generation squarks can affect significantly the mass of the SM-like Higgs boson, we vary freely all soft parameters in this sector except that we assume m U 3 = m D 3 for right-handed soft breaking masses and A t = A b for soft breaking trilinear coefficients.
• Considering that we require the NMSSM to explain the discrepancy of the measured value of the muon anomalous magnetic moment from its SM prediction, we treat the common value for all soft breaking parameters in the slepton sector (denoted by ml hereafter) as a free parameter.
• We fix gluino mass at 2TeV, and treat the Bino mass M 1 and the Wino mass M 2 as free parameters since they affect the properties of the neutralinos.
Then we use the package NMSSMTools-4.9.0 [57] to scan following parameter space: where all the parameters are defined at the scale of 1TeV. During the scan, we use following constraints to select physical parameter points:  (1), (2) and (3) presented in the text, and that in the second row corresponds to the samples that further satisfy the constraints from the direct search for sparticles at the LHC Run-I, i.e. the constraints (4) and (5). The quantity R in the last row represents the retaining ratio of the samples before and after considering the direct search constraints in our scan. 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. (1) All the constraints implemented in the package NMSSMTools-4.9.0, which include the Z-boson invisible decay, the LEP search for sparticles (i.e. the lower bounds on various sparticle masses and the upper bounds on the chargino/neutralino pair production rates), the B-physics observables such as the branching ratios for B → X s γ and B s → µ + µ − , the discrepancy of the muon anomalous magnetic moment, the dark matter relic density and the LUX limits on the scattering rate of dark matter with nucleon. In getting the constraint from a certain observable which has an experimental central value, we use its latest measured result and require the NMSSM to explain the result at 2σ level.
(2) Constraints from the direct searches for Higgs bosons at LEP, Tevatron and LHC.
Especially we require that one CP-even Higgs boson acts as the SM-like Higgs boson discovered at the LHC. We implement these constraints with the packages HiggsSignal for 125GeV Higgs data fit [58] 3 and HiggsBounds for non-standard Higgs boson search at colliders [61].
(4) Constraints from the preliminary analyses of the ATLAS and CMS groups in their direct searches for sparticles at the LHC Run-I. We implement these constraints by the packages FastLim [62] and SModelS [63]. These two packages provide cut efficiencies or upper bounds on some sparticle production processes in simplified model framework, and thus enable us to impose the direct search bounds in an easy and fast way. In the appendix, we briefly introduce the two packages.
(5) Constraints from the latest searches for electroweakinos and stops by the ATLAS collaboration at the LHC Run-I. We implement these constraints by detailed Monte Carlo simulation. Since we have to treat more than twenty thousand samples at this step, this process is rather time consuming in our calculation by clusters. In the appendix, we provide details of our simulation.
After analyzing the samples that survive the constraints, we find that they 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. In Table 1, we list the favored parameter ranges for each type of samples before and after considering the constraints from the direct search experiments, i.e. the constraints (4) and (5). Note that in the last row of this table, the retaining ratio R is defined by R = N af ter /N bef ore where N bef ore is the number of the samples that satisfy the constraints (1), (2) and (3) in our scan, and N af ter is the number of the samples that further satisfy the constraints (4) and (5). This table indicates that although the direct search experiments are effective in excluding parameter points encountered in the scan, they scarcely change the ranges of the input parameters where ∆ Z and ∆ h take rather low values, or equivalently speaking the NMSSM can naturally predict m Z and m h even after considering the direct search constraints from LHC Run-I. The underlying reason for this phenomenology is that the exclusion capability of the direct search experiments depends not only on sparticle production rate, but also on the decay chain of the sparticle and the mass gap between the sparticle and its decay product. We checked that a large portion of the excluded samples are characterized by Figure 1. ∆ Z and ∆ h predicted by the Type III samples (upper panels) and Type IV samples (lower panels) respectively. Samples in the left panels survive the constraints (1), (2) and (3) presented in the text, and those in the right panels further satisfy the constraints from the direct searches for sparticles at LHC Run I, i.e. the constraints (4) and (5). The panels in each row adopt same color convention for µ, which is presented on the right side of the row. M 2 ≤ 250GeV. In this case, there exist one moderately light neutralino and one moderately light chargino with both of them being Wino-dominated, and their associated production rate at the LHC is quite large so that the 3l + E miss T signal of the production after cuts may exceed its experimental upper bound (see appendix for more information). Among the four types of points, we also find that the lowest fine-tuning comes from type III and type IV samples, for which ∆ Z and ∆ h may be as low as about 2. This character is shown in Fig.1, where we project type III and IV samples on ∆ Z − ∆ h plane. We emphasize that Because the type IV samples were scarcely studied in previous literatures and also because they have similar phenomenology to that of the NS scenario in the MSSM, we in the following focus on this type of samples. In order to make the essential features of the samples clear, we only consider those that satisfy additionally the condition M 1 , M 2 , ml ≥ 300GeV. Hereafter we call such samples collectively as the NS scenario with a Higgsinodominatedχ 0 1 . As we will show below, |µ| in this scenario is upper bounded by about 160GeV, so the condition is equivalent to M 1 , M 2 , ml 2|µ|. In this case, gauginos and sleptons affect little on the properties of the lightest three neutralino, and the lighter chargino.
3 Key features of the NS scenario withχ 0 1 being Higgsino-dominated In this section, we investigate the features of the NS scenario withχ 0 1 being Higgsinodominated. We are particulary interested in neutralino-chargino sector and stop sector since they play an important role in determining the fine tunings of the theory. In Table 2 , and among the ranges, only the lower bound of mt 1 is shifted from 380GeV to 500GeV by the constraints from the direct search experiments. Moreover, we checked that in our scenario, the ratio 2κ/λ is restricted in the range from about 1 to 1.5. In this case, Higgsinos and Singlino are approximately degenerated in mass, and consequently they mix strongly to form mass eigenstates.  In Fig.2, we show the mass spectrum ofχ 0 1 ,χ 0 2 ,χ 0 3 andχ ± 1 in our scenario. This figure indicates that the mass splittings among the particles satisfy 30GeV ∆ ± 70GeV, 50GeV ∆ 0 110GeV and 80GeV (mχ0 3 −mχ0 1 ) 160GeV. We remind that these splittings are induced by the strong mixings between Higgsinos and Singlino, and significantly larger than those amongχ 0 1 ,χ 0 2 andχ ± 1 in the NS scenario of the MSSM [5] . In Fig.3, we show the field components of the statesχ 0 1 ,χ 0 2 andχ 0 3 respectively. As is expected, theH 0 u , S 0 andH 0 d components inχ 0 1 are comparable in magnitude with the largest one coming from theH 0 u component. We emphasize again that the large Singlino component, i.e. N 15 ∼ 0.5, can dilute the interactions of the Higgsino-dominatedχ 0 1 with other fields, and consequentlỹ χ 0 1 can reach its right relic density. In this case, we checked that the main annihilation channels ofχ 0 1 in early universe includeχ 0 1χ 0 1 → W + W − , ZZ, Zh 1 , h 1 h 1 , h 1 h 2 , qq. As forχ 0 2 , its largest component comes from eitherH 0 d field (for most cases) orH 0 u field (in rare cases), and in general the two components are comparable, which can be learned from the figure and also from Eq.(2.10). We checked that due to the spectrum and the mixings, the dominant decay ofχ ± 1 isχ ± 1 →χ 0 1 W * , and that ofχ 0 2 is usuallyχ 0 2 →χ 0 1 Z ( * ) . By contrast, the possible dominant decay modes ofχ 0 3 are rather rich, which includeχ 0 1 Z ( * ) ,χ 0 1 h 1 ,χ ± 1 W ( * ) . Sinceχ 0 3 is Singlino-dominated, its production rate is rather low, and consequently its phenomenology is of less interest.
Next we turn to the properties oft 1 . From the interactions oft 1 presented in [1], one can infer that ift 1 ist R dominated and meanwhile |N 14 | |N 24 |, the relation Br(t 1 → χ + 1 b) : Br(t 1 →χ 0 1 t) : Br(t 1 →χ 0 2 t) 2 : 1 : 1 should hold. On the other hand, ift 1 is t L dominated,t 1 prefers to decay into the Higgsino-dominatedχ 0 1,2 with Br(t 1 →χ 0 1 t) Br(t 1 →χ 0 2 t). These features are exhibited in Fig.4, where we show the correlations between different decay rates oft 1 in our scenario. From Fig.4, one can also learn that the branching ratio oft 1 →χ 0 3 t is less than 10%. This is becauseχ 0 3 is Singlino-dominated and itsH 0 u component is small. Moreover, we note that in our scenario mt 1 is lower bounded by about 500GeV, which is about 100GeV less than that in the NS scenario of the MSSM [64][65][66]. One reason for the difference is thatt 1 in our scenario has richer decay modes.

Detection at 14 TeV LHC
From the analysis in last section, one can learn that the NS scenario withχ 0 1 being Higgsinodominated is characterized by predicting µ ≤ 160GeV and sizable mass splittings among the Higgsino-dominated neutralinos and chargino, i.e. 30 GeV ≤ ∆ ± ≤ 70 GeV and 50 GeV ≤ ∆ 0 ≤ 110 GeV. Although this kind of spectrum is allowed by the direct searches for the electroweakinos at LHC Run I, it is expected to be tightly constrained at the upgraded LHC. Table 3. Mass and decay information of the benchmark points P1, P2, P3 and P4 in our study. We investigate this issue by considering the neutralino and chargino associated production processes at 14 TeV LHC. For simplicity we adopt 4 benchmark points listed in Table  3, which are discriminated by the values of µ (or equivalently mχ± 1 ), ∆ 0 and ∆ ± . Sinceχ ± 1 for these points decays intoχ 0 1 plus an off-shell W boson, andχ 0 2 decays mainly intoχ 0 1 plus a Z boson (on-shell or off-shell), the signal region SR0τ a in the ATLAS direct searches for electroweakinos by trileptons and large E miss T signal [67], which was proposed in the analy- Table 5. The best two signal bins and corresponding significances for the benchmark points with 30 fb −1 and 300 fb −1 integrated luminosity data respectively at 14 TeV LHC. sis [68] and also briefly introduced in the appendix of this work, is most pertinent to explore those points. In our analysis, we simulate the processes pp →χ ±  Table 4. We also present in the table the backgrounds of the bins at 14 TeV LHC, which were obtained by detailed simulations done in [9]. With these results, we evaluate the significance S = s/ b + ( b) 2 for each bin, where s and b correspond to the number of signal and background events and = 10% is the assumed systematical uncertainty of the backgrounds. Assuming 30 fb −1 and 300 fb −1 integrated luminosity data at 14 TeV LHC, we present the best two signal bins and corresponding expected significances for P1, P2, P3 and P4 in Table 5. This table reveals following information: • With 30 fb −1 integrated luminosity data, P1 and P2 can be excluded at 2σ confidence level, and with 300 fb −1 data P3 can also be excluded. In any case, the point P4 is hard to be excluded.
• For each point, which signal bin is best for exclusion depends on the mass splittings among the neutralinos and chargino. For example, since ∆ ± < m W for all the four points, the most effective bins usually require m T < 80 or 110 GeV. For points P1, P2 and P3, the bins satisfying m SF OS < m Z are preferred for exclusion since ∆ 0 < m Z , and by contrast the bins with |m SF OS − m Z | < 10 GeV (such as bins 14 and 15) are favored by point P4 since in this caseχ 0 2 can decay into an on-shell Z boson. Note that for bins 14 and 15, the backgrounds are relatively large, and that is why the point P4 can not be excluded at 14 TeV LHC after including the systematic uncertainties.
• With 300 fb −1 data, the point P2 can be discovered at 14 TeV LHC. This is partially becauseχ 0 2 andχ ± 1 are relatively light so that the rate of their associated production is large, and also because they have sizable mass splittings fromχ 0 1 to result in moderately energetic decay products.
• Since the bins in the SR0τ a are disjoint, in principle their results can be statistically combined to maximize the significance. We did this, but we found that the improvement is not significant.
Since the four benchmark points stand for the typical situation in our scenario, we conclude that the future LHC experiments can exclude most part of the parameter space for the scenario, and consequently the fine-tuning of the NMSSM can be pushed to higher level. We will discuss such an issue extensively in our forthcoming work.

Dark Matter Direct Search
Sinceχ 0 1 as the DM candidate has quite largeH 0 u ,H 0 d andS 0 components in our scenario, its interactions with the Higgs bosons and Z boson are unsuppressed. As a result, the spin-independent (SI) and spin-dependent (SD) cross sections of theχ 0 1 -nucleon scattering are sizable, and may reach the sensitivities of future DM direct detection experiments such as XENON-1T and LZ-7.2T experiments [69]. In this section, we investigate such an issue. In Fig.5, we project the samples of our scenario on mχ0 1 − σ SĨ χ−p and mχ0 1 − σ SD χ−p planes with σ SĨ χ−p and σ SD χ−p denoting the SI and SD cross sections respectively. The blue lines, red lines and green lines are the sensitivities of LUX, XENON-1T and LZ experiments respectively. For the SI cross section, one can learn from Fig.5 that the future XENON-1T experiment is able to probe a large portion of the samples, and the LZ experiment can test even more. Anyhow, there still exist some samples remaining untouched by these future experiments. We numerically checked that for the untouched samples, there exists rather strong cancelation among the contributions induced by different CP-even Higgs bosons. On the other hand, the story for σ SD χ−p is quite different. From the right panel of Fig.5 one can see that the future XENON-1T experiment can test nearly all of the samples, let alone the more sensitive LZ experiment. The underlying reason is that in the NMSSM with heavy sfermions, the SD cross section gets contribution mainly from the t-channel Z-mediated diagram. As a result, the size of the cross section is determined by the Zχ 0 1χ 0 1 coupling, which is given by In the last step of the equation, we have used the information of N 13 and N 14 presented in Fig.3. Moreover, it is interesting to see that although the benchmark point P4 in Table  3 is hard to be excluded by sparticle direct search experiments at the LHC, its SD cross section is quite large and so the point will be tested by future dark matter direct search experiments. In getting Fig.5, we use the package micrOMEGAs [70] to calculate the cross sections. We choose its default setting σ πN = 34MeV and σ 0 = 42MeV as input. We checked that if we take σ πN = 59MeV from [71] and σ 0 = 58MeV from [72][73][74], the SI cross section will be enhanced by a factor from 20% to 40%, and this does not affect the conclusions presented in this work.

Conclusions
With the great discovery of the 125 GeV Higgs boson and the increased bounds on sparticle masses at LHC Run I, the NS scenario in MSSM has become theoretically unsatisfactory. By contrast, the situation may be improved greatly in the NMSSM. This motivates us to scrutinize the impact of the direct searches for SUSY at LHC Run I on the naturalness of the NMSSM.
We start our study by scanning the vast parameter space of the NMSSM to get the region where the fine tuning measures ∆ Z and ∆ h at electroweak scale are less than about 50. In this process, we considered various experimental constraints such as DM relic density, LUX limits on the scattering rate of DM with nucleon and the 125 GeV Higgs data on the model. We classify the surviving samples 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. After these preparations, we specially study the influence of the direct searches for SUSY on the samples. We implement the direct search constraints by the packages FastLim and SModelS and also by simulating the electroweakino and stop production processes. Our results indicate that although the direct search experiments are effective in excluding the samples, the parameter intervals for the region and also the minimum reaches of ∆ Z and ∆ h are scarcely changed by the constraints, which implies that, contrary to general belief, the fine tuning of the NMSSM does not get worse after LHC Run I. Our results also indicate that the lowest fine-tuning comes from type III and type IV samples, for which ∆ Z and ∆ h may be as low as about 2 without conflicting with any experimental constraints.
Considering that the type IV samples were scarcely studied in previous literatures and that they have similar phenomenology to that of the NS scenario in the MSSM, we investigate the essential features of this kind of samples. We find that they are characterized by strong mixings between Higgsinos and Singlino in forming mass eigenstates called neutralinos. As a result, the lightest neutralinoχ 0 1 as the DM candidate has significant Singlino component so that it can easily reach the measured DM relic density, and meanwhile the mass splittings among the Higgsino-dominated particlesχ 0 1 ,χ 0 2 andχ ± 1 are usually larger than 30GeV. These features make the samples rather special. For example, we show that due to the rich decay products of the lighter scalar top quark, its lower mass bound is decreased by about 100GeV in comparison with that in the NS scenario of the MSSM, and that the neutralino-chargino sector of the samples can be readily tested either through searching for 3l + E miss T signal at 14 TeV LHC or through future dark matter direct detection experiments.
In summary, we conclude that so far the fine tuning of the NMSSM is scarcely affected by the direct searches for SUSY at LHC Run I, and it can still predict Z boson mass and the SM-like Higgs boson mass in a natural way. This conclusion, however, may be altered by the upcoming 14 TeV LHC experiments and DM matter direct search experiments.   In this section, we briefly introduce the packages Fastlim [62] and SModelS [63], which can be used to implement the constraints from the direct searches for sparticles at LHC Run I in an easy and fast way.

A.1 Fastlim
Fastlim is a package that limits SUSY parameter space using the LHC-8TeV data. It incorporates in its database cut efficiency tables for different sparticle production processes in simplified model framework. Now it supports the processes of gluino pair production and third generation squark pair productions, and it involves many decay modes of gluino and the squarks, such asg →t * 1,2 t →χ 0 1 tt org → qqχ 0 1 andb 1,2 → bχ 0 1 where q stands for the first two generation quarks. Any processes that contain the same final states at the detector level are combined by Fastlim to improve the signal significance.
In this work, we only use the experiments of searching for gluino and squarks, which are listed in Table 6. This is because that Fastlim gives better limits for strong SUSY productions in comparison with the package SModelS [63].

A.2 SModelS
SModelS have the same function as that of Fastlim, but it could give better limits for slepton productions and electroweakino productions. SModelS has implemented the information about following experiments: • Direct slepton searches (ATLAS): ATLAS-CONF-2013-049 [83].
• Electroweakino searches (CMS): SUS-12-022 [84], SUS-13-006 [85], SUS-13-017 [90], and it contains the cross section limits at 95% confidence level for above analyses. If one parameter point predicts cross sections larger than those presented in the analyses, it will be excluded. Otherwise, it is allowed. Note that the efficiencies or the upper bounds on SUSY signals in the database of the two packages are usually based on certain assumptions, which may not be applied to some parameter points encountered in our scan. In this case, the encountered point is considered to be experimentally allowed. Also note that the packages are based on the preliminary analyses of the ATLAS and CMS groups, which were done in 2013. Given that most of these analyses have been updated in past two years, a more powerful exclusion capability may be obtained if one repeats the updated analyses by detailed Monte Carlo simulation. Anyhow, these two packages can serve as useful tools to exclude some SUSY parameter points.

B Details of our simulation
In our study, any point that passed Fastlim and SModelS is further tested by simulations to see whether it survives the constraints from the direct search experiments. In detail, we first use MadGraph/MadEvent [91] to generate parton level events for certain sparticle production processes, and feed them into Pythia [92] for parton showering and hadronization. Then we use the package CheckMATE [94] where a well-tuned Delphes [93] is provided for the detector simulation and analyses. We define R ≡ max S i /S 95 i,obs to decide whether the point survives the analysis, where S i stands for the simulated signal events in the ith signal region of the analysis, and S 95 i,obs represents the 95% C.L. upper limit of the event number in the signal region. If R > 1, the parameter point is excluded by the analysis and otherwise it is allowed.
Since the fine tuning of the NMSSM at the electroweak scale is mainly affected by its chargino-neutralino sector and stop sector, we repeat by simulation the ATLAS analyses in [95], [67], [97] and [99]. All these analyses are based on 20.3fb −1 data at the LHC-8TeV with the former two presenting so far the strictest limits on electroweakino productions, and the latter two providing the tightest constraints on the pair productions of third generation squarks. In the following, we briefly introduce these analyses.

B.1 Search for Electroweakino at the LHC
The analysis [95] targets final states with two leptons and large E miss T . In our simulation of the analysis, we mainly focus on the signal region named "SR-Zjets", which is specifically designed for the process pp →χ 0 2χ ± 1 → Zχ 0 1 Wχ 0 1 → χ 0 1 qqχ 0 1 . This signal region requires that the two leading leptons should be same flavor but opposite sign (SFOS), and that their invariant mass locates at Z-peak. As was pointed out in [95], it provides the strongest constraints on the chargino-neutralino sector among the anlayses with dilepton final state.
The analysis in [67] also searches for electroweakino productions but with final states of three leptons and large E miss T . Here we concentrate on the signal region named "SR0τ a" which is optimized for processes pp →χ 0 2χ ± 1 → Z ( * )χ0 1 W ( * )χ0 1 → χ 0 1 νχ 0 1 . This region requires a pair of SFOS leptons in its signal, and utilizes the transverse mass m T = 2 p T E miss T − 2 p T · E miss T (here p l T is the transverse momentum of the lepton not forming the SFOS lepton pair) to suppress the SM background. It considers twenty bins, which are categorized by the SFOS leptons' invariant mass, m T and E miss T , to maximize its sensitivity to different mass spectrums ofχ 2 ,χ ± 1 andχ 0 1 . We remind that in the NS scenarios of the NMSSM, allχ 0 iχ ± 1 associated production processes with i = 2, 3, 4, 5 may contribute sizably to the trilepton signal, so in our simulation we include all these contributions. By contrast, SModelS does not combine processes that have the same final states at the detector level. Moreover, in our analysis we combine the signal region "SR-Zjets" with the bins in "SR0τ a" to maximize the discovery significance by the CL s method in RooStats as we did in [9].

B.2 Search for stops at the LHC
In the NS scenario of the NMSSM,t 1 usually decays liket 1 →χ 0 2,3 t →χ 0 1 Z ( * ) t ort 1 → χ + 1 b →χ 0 1 W ( * ) b. Considering that these two-step topologies haven't been included in the database of Fastlim, we repeat recent ATLAS analyses on stop pair productions, which were presented in [97] and [99].
The analysis in [97] searches for stops in final states containing exactly one isolated lepton, at least two jets and a large E miss expected. These facts motivates physicists to define five signal regions, which are called "SR2(A,B,C)" and "SR3(A,B)" respectively, by the number of leptons in the signal, the jet multiplicity and whether the Z boson is boosted. For example, the signal regions "SR2(A,B,C)" require that the signal events contain exactly two signal leptons and the Z boson is boosted. The regions "SR2A" and "SR2B" are optimized for low jet multiplicity, while the region "SR2C" is designed for high jet multiplicity case.
[76] "Search for strong production of supersymmetric particles in final states with missing transverse momentum and at least three b-jets using 20.1fb −1 of pp collisions at √ s = 8TeV with the ATLAS Detector", ATLAS-CONF-2013-061.
[78] "Search for direct third generation squark pair production in final states with missing transverse momentum and two b-jets in √ s = 8TeV pp collisions with the ATLAS detector", ATLAS-CONF-2013-053.
[79] The ATLAS collaboration, "Search for direct top squark pair production in final states with two leptons in √ s = 8TeV pp collisions using 20fb −1 of ATLAS data", ATLAS-CONF-2013-048.
[80] "Search for squarks and gluinos with the ATLAS detector in final states with jets and missing transverse momentum and 20.3fb −1 of √ s =8TeV proton-proton collision data", ATLAS-CONF-2013-047.
[81] "Search for direct top squark pair production in final states with one isolated lepton, jets, and missing transverse momentum in √ s = 8TeV pp collisions using 21fb−1 of ATLAS data", ATLAS-CONF-2013-037.