Suppressing the Scattering of WIMP DM with Nucleons in Supersymmetric Theories

The continuously improving sensitivity of dark matter direct detection experiments has limited the interaction between dark matter and nucleons being increasingly feeble, while the dark matter relic density favors it to take part in weak interactions. After taking into account the constraints from the Large Hadron Collider (LHC) search for Higgs bosons and sparticles, it is becoming difficult for the neutralino dark matter in the Minimal Supersymmetric Standard Model and the Next-to-Minimal Supersymmetric Standard Model to possess these two seemingly paradoxical features in their most natural parameter space for electroweak symmetry breaking due to the limited theoretical structure. In contrast, the seesaw extension of the Next-to-Minimal Supersymmetric Standard Model, which was initially proposed to solve the neutrino mass problem, enables the lightest sneutrino to act as a viable dark matter candidate, readily has these features, and thus, it easily satisfies the constraints from dark matter and LHC experiments. Compared with the Type-I seesaw extension, the dark matter physics in the inverse seesaw extension is more flexible, allowing it to be consistent with the experimental results in broader parameter space. We conclude that weakly interacting massive particles (such as the sneutrino in this study) work well in supersymmetric theories as dark matter candidates.


Introduction
Current astronomical observations have confirmed the presence of dark matter (DM) and revealed that it accounts for about 27% of the composition of the universe [1]. Among the various DM candidates, the weakly interacting massive particle (WIMP) is the most promising since it can obtain the measured DM density naturally and straightforwardly. In the popular supersymmetric theories, the lightest supersymmetric particle (LSP) usually has such a property and thus can act as a viable DM candidate [2]. In this case, since the massive supersymmetric particles decay ultimately to the LSP due to R-parity conservation, their signals in colliders will contain missing momentum, a feature widely used when searching for supersymmetry at the Large Hadron Collider (LHC). Although scientists have studied the WIMP as a popular DM candidate for many years, the rapid progress of DM direct detection (DD) experiments in recent years has sharply restricted its interaction with nucleons [3,4], leading more and more people to question its rationality. Since supersymmetric theories usually predicts WIMP DMs, they are also facing unprecedented doubts. It is the primary purpose of this work to discuss whether the WIMP DM predicted by supersymmetry can be naturally consistent with the results of the DM experiments. In particular, although the DM density prefers it to take part in weak interactions, the continuously improved sensitivity of DM DD experiments has limited its interaction with nucleons being increasingly feeble. As shown below, after taking into account the constraints from the LHC search for Higgs bosons and sparticles, it is becoming difficult for the neutralino DM in the Minimal Supersymmetric Standard Model (MSSM) [5][6][7] and the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [8] to possess the two seemingly paradoxical features in their most natural parameter space (i.e., the Natural SUSY scenario in the literature [9]) due to the limited theoretical structure. However, the seesaw extension of the NMSSM, which was initially proposed to solve the neutrino mass problem, enables the lightest sneutrino to act as a viable DM candidate and readily has these features, and thus, it easily satisfies the constraints of both dark matter and LHC experiments.
As a starting point for discussion, we first consider the MSSM [5][6][7]. When the correct DM density is required, the lightest neutralino with a bino field as its dominant component is customarily taken as the DM candidate 1 . Since the bino field transforms non-trivially under the U (1) Y group of the Standard Model (SM), its coupling with the SM-like Higgs boson is approximated as [10][11][12]  if the wino field is very massive, where mχ0 1 denotes the lightest neutralino mass that relates with the bino mass M 1 by mχ0 1 M 1 , µ represents higgsino mass, tan β = v u /v d is the ratio of Higgs vacuum expectation values, and α is the mixing angle of the CP-even Higgs states [7]. In obtaining the final expression, we assume the decoupling limit of the Higgs sector, i.e., m A v, where m A denotes the mass of the CP-odd Higgs boson. Given that the spin independent (SI) cross-section of the DM scattering with nucleons is mainly contributed to by the SM-like Higgs boson and thus approximated by [13] σ SĨ (1. 2) we infer that, if M 1 and µ are of the same sign, µ must be sufficiently large to be consistent with the latest results of the XENON-1T experiment [3]. Numerically speaking, we find |µ| 430 GeV (350 GeV) for tan β = 10 (20) if the DM annihilated by a resonant SMlike Higgs boson to obtain the correct relic density. These formulas also show that, if M 1 and µ are of opposite signs, which can result in the blind spots of the scattering [14][15][16][17], |µ| ∼ 100 GeV seems to be experimentally allowed. However, such a possibility has been limited by the XENON-1T search for spin dependent (SD) DM-nucleon scattering because, regardless of the relative sign between M 1 and µ, small values of |µ| can enhance the scattering cross-section. This can be understood from the following approximation [10][11][12] Cχ0 1χ 0 1 Z e tan θ W cos 2β 2 with C p 1.8 × 10 −40 cm 2 for protons and C n 1.4 × 10 −40 cm 2 for neutrons [18,19]. Furthermore, the LHC search for electroweakinos can play a complementary role in limiting such a possibility. In the Appendix of this work, we study these experimental constraints in the blind spot scenario of the MSSM, and find that the case of |µ| 300 GeV has been excluded for any mχ0 1 by the experimental upper bound of the SD cross-section, and |µ| must be larger than about 390 GeV in the Higgs funnel region due to the LHC experiment. We note that a global fit of the MSSM was recently performed [20], where various experimental constraints, including those from the DM relic density, PandaX-II (2017) results for the SI cross-section [21], PICO results for the SD cross-section [22], and the searches for supersymmetric particles at the 13 TeV LHC with 36 fb −1 data (especially the CMS analysis of the electroweakino production [23]), were considered. The analysis showed that |µ| 350 GeV was favored at a 95% confidence level 2 . Such a µ can induce a tuning of about 3% to predict the Z boson mass [9], and this situation will be further exacerbated if future DM DD experiments fail to detect the sign of the DM and/or forthcoming high luminosity LHC experiments do not find the evidence of the electroweakinos (see, e.g., the discussion in the Appendix about the blind spot of the MSSM).
Next we consider the NMSSM [8]. Since the coupling of the singlino field in this theory with SM particles may be very weak, we discuss the case of the singlino-dominated neutralino as a DM candidate [24][25][26][27], which is feasible if the Higgs self-coupling coefficients λ and κ satisfy λ > 2κ and the gauginos are relatively massive. In the natural NMSSM that essentially requires |µ| 500 GeV [28], this DM may annihilate by the following channels to obtain the correct relic density (see, for example, Ref. [12,13,19,26,27,[29][30][31][32][33][34][35]):χ 0 1χ 0 1 → tt, h s A s , hA s , W W, ZZ, · · · through s-channel exchange of non-resonant Z and Higgs bosons, t/u-channel exchange of electroweakinos (h, h s , and A s denote SM-like, singlet dominated CP-even, and CP-odd Higgs bosons, respectively), the funnels in Z, A s , h and h s , or the co-annihilation with sleptons,χ 0 2 (the next-to-lightest neutralino), and χ ± 1 (the lightest chargino). Among these channels, the cross-sections for the first three are potentially large, so each of them may be fully responsible for the measured DM relic density [13]. However, one can verify with Effective Field Theory (EFT) that such a possibility has been tightly limited by the upper bounds of the SI and SD cross-sections (see the discussion in the Appendix). The basic reason is that all of these channels require a large λ to account for the density, and as a result, the DM couplings with h and Z bosons 2 We emphasize that this conclusion does not mean that the case |µ| < 350 GeV is completely excluded by the experiments. It just implies that the case of |µ| < 350 GeV has been tightly limited and the probability of its occurrence is rather low in frequentist statistics. are sizable, since for λv/|µ| 1 and massive gauginos [13,18], where mχ0 1 denotes the DM mass and v ≡ 174 GeV. The other channels may play an important role in determining the density only with a specific sparticle mass spectrum, and they do not necessarily correspond to a large λ. However, given that the current sensitivities of the DM DD experiments have reached the precision of 10 −47 cm 2 for the SI cross-section [3] and 10 −42 cm 2 for the SD cross-section [4], a large portion of the parameter space for the annihilations in the natural NMSSM can predict SI and/or SD cross-sections that are comparable with or even lager than their bounds. As a result, the DD experiments have tightly constrained the space [33,36]. Furthermore, the LHC experiments can limit the theory, since the annihilations are usually accompanied with light sparticles or Higgs bosons [33]. Recently, researchers performed comprehensive studies of the natural NMSSM by considering various experimental constraints [33,34]. In particular, they included those from the DM and LHC experiments. The conclusion was that only some corners of the parameter space are allowed, which have either of following features: • λ 2κ with λ 0.05. This case corresponds to the decoupling limit of the NMSSM [8], and the DM co-annihilated with the higgsinos to obtain the density [33].
• κ ∼ 0.01, λ 0.2 and there exists at least one light singlet dominated Higgs boson [34]. In this case, the singlino-dominated DM annihilates in certain funnel regions, and the higgsinos decay in a complex way to satisfy the LHC constraints.
We emphasize that this conclusion is valid only for |µ| 500 GeV, and it may be improved by more intensive study. We also emphasize that, as indicated by Eq. (1.5) and (1.6), the increase in |µ| is helpful to relax the constraints of the DM DD experiments 3 .
Based on the discussion above, if the neutralino DM transforms non-trivially under the electroweak gauge group of the SM or it mixes with other fields and obtains the charge of the group, its scattering cross-sections with nucleons tend to be relatively large for the higgsino 3 Very recently we updated a previous study [33]. We no longer required the fine tuning of mZ to be less than 50. Instead, we imposed the condition µ ≤ 1 TeV. We found that for the singlino-dominated DM, in addition to the regions mentioned above, there is a new scenario that is consistent with current experimental constraints, which is characterized by 0.4 λ 0.7, 200 GeV mχ0 1 600 GeV, 400 GeV µ 800 GeV, and the splitting between mχ0 1 and µ is larger than approximately 80 GeV. This conclusion is consistent with a previous report [13]. From this scenario, one can also obtain the bino-singlino well tempered DM scenario reported previously [13] by choosing a negative M1 that satisfies |M1| < µ, and letting |M1| approach the singlino-dominated neutralino mass from below. In either scenario, the SI cross-section can be as low as 10 −50 cm 2 , while the SD cross-section is usually larger than 10 −42 cm 2 . mass µ upper bounded by several hundred GeV. These rates usually contradict relevant experimental bounds after accounting for the rapidly improving sensitivities of the DM DD experiments in recent years. This implies that, if one wants the scattering rate naturally suppressed, the DM should be a gauge singlet field, or its singlet component should at least be naturally far dominant over the other components. Moreover, the situation of the MSSM and NMSSM tells us that the DM is preferably not a neutralino when one extends economically the models. We emphasize that such a DM can still be a WIMP in the sense that it may have weak couplings with the particles beyond the SM, which is necessary to obtain the correct density naturally. These inferences inspire us to augment the NMSSM with a TeV-scale Type-I seesaw mechanism and choose the lightest sneutrinoν 1 as a DM candidate [37]. The resulting theoretical framework not only produces neutrino mass, it also guarantees thatν 1 will be almost purely right-handed, since the tiny neutrino Yukawa couplings suppress the chiral mixing of the sneutrinos significantly [37][38][39] (note that the possibility of a left-handed sneutrino as a feasible DM candidate was ruled out by DM DD experiments one decade ago [40,41]).
Given that the right-handed sneutrino field is a gauge singlet, it can interact directly with the singlet Higgs field by triple or quartic scalar vertexes. This feature allows these fields to form a secluded DM sector, which communicates with the SM sector only by small singlet-doublet Higgs mixing and accounts for the relic density through the annihilation of ν 1 into a pair of singlet dominated Higgs bosons [37]. In addition, the singlet Higgs field could also mediate the transition of the sneutrino pair into a higgsino pair and vice versa in the thermal bath of the early universe. If the sneutrino and the higgsino approximately degenerate in mass, the annihilation of the higgsinos in the freeze-out stage can affect the DM density significantly, which makes the relic density consistent with its experimental measurements [37] (in the literatures, this phenomenon is called co-annihilation [42,43]). Similar to the secluded DM case,ν 1 couples very weakly with the SM particles due to its singlet nature, and this suppresses the DM scattering with nucleons. Last but not least, sinceν 1 is a scalar particle with a definite CP number, the SD cross-section of its scattering with nucleons always vanishes, which is another advantage ofν 1 coinciding with the results of the DM DD experiments. We add that the seesaw extension of the MSSM does not have all of these features (see the introduction in [37]), and the aforementioned theoretical framework extends the field content of the NMSSM only by three generations of the right-handed neutrino superfield. Thus, the NMSSM extension may be the most economical supersymmetric model that can suppress the DM-nucleon scattering naturally.
Similarly, one can embed the inverse seesaw mechanism in the NMSSM, and the resulting theory has similar features to the Type-I seesaw extension in DM physics [44]. Compared with the Type-I extension, this theory has at least two advantages. One is that the neutrino mass is doubly suppressed so that the right-hand neutrino mass can be naturally at the TeV scale even for sizable neutrino Yukawa couplings. The other is that DM physics is very rich due to its more complex structure in the neutrino/sneutrino sector and can be consistent with the experimental results in a more flexible way. The theory is another economic model that naturally suppresses the DM-nucleon scattering. In fact, given that the neutralino DM has been tightly limited by DD experiments in recent years, the sneutrino DM embedded in different extensions of the MSSM has regained broad interest [45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64].
In the NMSSM with any of the seesaw mechanisms, the scattering of the sneutrino DM with nucleons proceeds mainly by the t-channel exchange of CP-even Higgs particles. Although the cross-section of the scattering is usually small, it is still potentially large if the lightest CP-even Higgs boson is significantly lighter than the SM-like Higgs boson discovered at the LHC and it contains sizable doublet Higgs components [37,44]. This case is not only consistent with current Higgs data [65], it also has theoretical advantages of enhancing the SM-like Higgs boson mass and naturally predicting the Z boson mass [66,67]. Thus this case is attractive. Studying the characteristics of the DM-nucleon scattering in this special case for both the seesaw extensions, in particular how it coincides with the latest XENON-1T results, can improve our understanding of the scattering, which is the primary purpose of this paper. In the following, we denote the NMSSM with the Type-I seesaw mechanism as Type-I NMSSM and that with the inverse seesaw mechanism as ISS-NMSSM.
We organize this paper as follows. In Section 2, we introduce the basics of the sneutrino sector in the Type-I NMSSM and ISS-NMSSM, including its mass matrix, its interaction with Higgs bosons and its scattering with nucleons. In Section 3, we consider the particular configuration of the Higgs sector. We vary the parameters in the sneutrino sector and compare the mechanisms of the theories that keep the sneutrino DM compatible with the DD experimental constraints. We also discuss the phenomenology of the extensions. Finally, we draw our conclusions in Section 4.

Theoretical preliminaries
As economic but complete supersymmetric theories to account for neutrino mass, the Type-I NMSSM and ISS-NMSSM adopt the same gauge groups as the NMSSM and extend each lepton generation of the NMSSM only by one and two gauge singlet chiral fields with the lepton number, respectively (see Table 1). These fields can couple directly with the singlet Higgs fieldŝ, and as a result, the singlet-dominated Higgs bosons play extraordinary roles in these theories: in addition to solving the µ problem of the MSSM [8] and enhancing the theoretical prediction of the SM-like Higgs boson mass through the singlet-doublet Higgs mixing [66,68], they are responsible for the massive neutrino mass and make the lightest sneutrino a viable DM candidate [37][38][39]44]. This feature renders the two theories quite different from seesaw extensions of the MSSM in various aspects.
The common feature of the two extensions is that they have same structure in the Higgs and neutralino/chargino sectors as that of the NMSSM, so they predict three CP-even Higgs bosons, two CP-odd Higgs bosons, five neutralinos, and two charginos. Throughout this paper, we label the particles with the same CP and spin quantum numbers in an ascending mass order, e.g., m h 1 < m h 2 < m h 3 for the CP-even Higgs bosons. We have discussed the properties of these particles in our publications [37,44], and here we only emphasize that they play an important role in the annihilation of the sneutrino DM, which includes the following channels [38,39]: (1)ν 1H → XY andHH → X Y , whereH andH denote higgsino-dominated neutralinos or chargino, and X ( ) and Y ( ) represent any possible SM particles (including the massive neutrinos and the extra Higgs bosons if the kinematics are accessible). This annihilation mechanism is called co-annihilation in the literature [42,43], and it is efficient only when the mass splitting betweenH andν 1 is less than about 10%. As pointed out by the Bayesian analysis of the Type-I NMSSM in [37], it is the most important annihilation channel.
(2)ν 1ν1 → ss * , where s denotes a light Higgs boson. This channel proceeds via any relevant quartic scalar coupling, the s-channel exchange of a Higgs boson and the t/u-channel exchange of a sneutrino. It is the second important annihilation channel of the DM by the analysis in [37].
(3)ν 1ν1 → V V * , V s, ff where V and f denote any gauge boson and SM fermion, respectively. This kind of annihilation proceeds mainly by the s-channel exchange of a resonant CP-even Higgs boson.
(4)ν 1ν1 → ν hνh via the s-channel exchange of a Higgs boson and the t/u-channel exchange of a neutralino, where ν h denotes a massive neutrino.
i → XY andν 1ν 1 → X Y whereν 1 denotes a sneutrino with an opposite CP number to that ofν 1 . This annihilation channel is important only whenν 1 and ν 1 are nearly degenerate in mass.
The difference between the extensions and the NMSSM arises only from the neutrino/sneutrino sector, which is evident from their constructions. Since sneutrino DM is the focus of this work, we recapitulate the basics of the sneutrino sector in the following sections.

Sneutrino sector of Type-I NMSSM
With the field content in Table 1, the superpotential and the soft breaking terms of the Type-I NMSSM are given by [38,39] where W MSSM and L MSSM represent the corresponding terms of the MSSM without including those for the Higgs sector, the coefficients λ and κ parameterize the interactions between the Higgs fields, Y ν andλ ν are neutrino Yukawa couplings with the flavor index omitted, m i (i = H u , H d , · · · ) denotes the soft breaking masses, and A i (i = λ, κ, · · · ) are soft breaking parameters for the trilinear terms. Noting that the soft masses m 2 Hu , m 2 H d , and m 2 S are related to the vacuum expectation values (vev) of the fields H u , H d and S, by the minimization conditions of the Higgs potential after electroweak symmetry breaking [8], we take λ, κ, 2 as theoretical inputs in the following discussion. In the Type-I NMSSM, the active neutrino mass matrix is formulated by representing the heavy neutrino mass matrix [69]. Since the active neutrino masses are at the 0.1 eV order, the magnitude of Y ν should be around 10 −6 if the massive neutrino masses are taken at the TeV order. Moreover, to reproduce neutrino oscillation data, m ν must be flavor non-diagonal. This can be realized by choosing a flavor non-diagonal Y ν and a diagonalλ ν [70]. If one further assumes that the soft breaking parameters in the sneutrino sector, such as ml,mν, and A λν , are flavor diagonal, the flavor mixing in the sneutrino mass matrix is induced only by the off-diagonal elements of Y ν , which is greatly suppressed. Given these facts, it is sufficient to consider only one generation sneutrino in studying the sneutrino DM [37]. In the following, we concentrate on the third generation case and use the symbols λ ν , A λν , and mν to denote the 33 element ofλ ν ,Ā λν , andmν, respectively.
After rephrasing the sneutrino fields by CP-even and CP-odd parts, the sneutrino mass matrix in the bases (ν L1 ,ν R1 ,ν L2 ,ν R2 ) is given by If all the parameters in the matrix are real, i.e., there is no CP violation, the real and imaginary parts of the sneutrino fields will not mix. In this case, the 4 × 4 mass matrix splits into two 2 × 2 matrices: where i = 1 and 2 denote CP-even and CP-odd states, respectively, and the minus signs in the matrix elements are relevant to the CP-odd states. These formulas indicate that the chiral mixings of the sneutrinos are proportional to Y ν , and hence, they are negligible. Thus, the sneutrino mass eigenstate and the chiral state coincide. In our study, we selected the lightest right-handed state as the DM candidate. The formulas also indicate that the mass splitting between CP-even and CP-odd right-handed states is given by ∆m 2 ≡ m 2 even − m 2 odd = 4m 2 RR , which implies that the sneutrino DM has a definite CP number, i.e. it is CP-even if m 2 RR < 0 and CP-odd for the other case. In our study of the sneutrino DM, we consider both CP possibilities.
The coupling strengths of the CP-even sneutrino DM with Higgs bosons are given by where Z ij (i, j = 1, 2, 3) and Z mn (m, n = 1, 2) are the elements of the rotations to diagonalize the CP-even Higgs mass matrix in the bases (Re

Re[S]) and the CP-odd Higgs mass matrix in the bases (
, respectively. Those for the CP-odd DM case are obtained from the formulas for the CP-even state by the substitution λ ν → −λ ν . These expressions indicate that Cν 1ν1 h i is suppressed by a factor λλ ν cos β if h i is the SM Higgs boson, which corresponds to setting Z i1 = cos β, Z i2 = sin β, and Z i3 = 0, and all three couplings may be moderately large only if the Higgs bosons are singlet dominant. Moreover, among the parameters in the sneutrino sector, λ ν and A λν affect both the couplings and masses of the sneutrinos, while m 2 ν only affects the masses. These features are helpful for understanding the behavior of the DM-nucleon scattering discussed below.

Sneutrino sector of ISS-NMSSM
Compared with the Type-I NMSSM, the ISS-NMSSM is much more complex in its neutrino/sneutrino sector. With the assignment of the quantum number for the fields in Table  1, its renormalizable superpotential and soft breaking terms take the following form [44] where the terms in the first bracket on the right side make up the Lagrangian of the NMSSM, and those in the second bracket are needed to implement the inverse seesaw mechanism. The coefficients λ ν and Y ν in the superpotential are neutrino Yukawa couplings, A λν and A ν in L sof t are coefficients for the soft breaking trilinear terms, and all of them are 3 × 3 (diagonal or non-diagonal) matrices in the flavor space. Moreover, among the parameters in the superpotential, only the matrix in the flavor space µ X is dimensional. This matrix parameterizes the effect of the lepton number violation (LNV), which may arise from the integration of heavy particles in an ultraviolet high-energy theory with LNV interactions (see, for example, [71][72][73]), so the magnitude of its elements should be suppressed. Similarly, the soft breaking parameter B µ X tends to be small.
for an arbitrary matrix M , one can approximate the mass matrix of the light active neutrinos by under the condition µ X M D M R . In this approximation, F ≡ M T D M T −1 R , and the magnitudes of its elements are of the order M D / M R . Thus, in the inverse seesaw mechanism, the active neutrino masses are suppressed in a double way, i.e., by the smallness of the LNV matrix µ X and also by the suppression factor if the mass scale of the massive neutrinos M R is taken at the TeV order.
Given the expression in Eq. (2.6), one can solve the matrix µ X with neutrino masses m ν i and the unitary Pontecorvo-Maki-Nakagawa-Sakata matrix U PMNS extracted from low energy experiments [74] to obtain [75,76]: . This formula indicates that if the Yukawa couplings Y ν and λ ν are flavor diagonal, the neutrino oscillation observed in the low energy experiments is attributed only to the non-diagonality of µ X . In this case, being unitary in the neutrino sector requires [77] [  22 . In addition to the inputs λ, κ, µ, and tan β in the Higgs sector, the sneutrino sector in the ISS-NMSSM involves the parameters Y ν , λ ν , A ν , A λν , µ X , and B µ X and the soft breaking masses ml, mν, and mx. As a result, the squared mass of the sneutrino fields is given by a 9 × 9 matrix in three-generation (ν L ,ν * R ,x) bases, whose form is quite complicated. However, we note the fact that among these parameters, only µ X must be flavor non-diagonal to predict the neutrino oscillations, but since its elements are usually less than 10 KeV [75], it can be safely neglected in calculating the sneutrino mass. Thus, if there are no flavor mixings for the other parameters, the matrix is flavor diagonal, and one generation (ν L ,ν * R ,x) bases can be used to study the properties of the sneutrino DM. In this work, we take the third generation sneutrinos as the DM sector, which is motivated by the fact that both the unitary bound and the constraints of the LHC search for sparticles in this sector are significantly weaker than those in the other generations [44]. In the discussion below, when we refer to the parameters Y ν , λ ν , A ν , A λν , mν, mx, and ml, we are actually referring their 33 elements, which is the same as what we did for the Type-I NMSSM.
If the sneutrino field are decomposed into CP-even and CP-odd parts, the squared mass matrix of the CP-even fields is given by This matrix shows that the mixing of the φ 1 field with the other fields is determined by the parameters Y ν and A ν . As Y ν approaches zero, m 12 and m 13 vanish, and consequently, φ 1 no longer mixes with the fields φ 2 and φ 3 . This situation is quite similar to that of the Type-I NMSSM. Moreover, if the first terms in m 22 and m 33 are far dominant over the other terms in their respective expressions, then m 22 m 33 . This results in maximal mixing between the φ 2 and φ 3 fields. In this case,ν 1 is approximated byν 1 if the left-handed field is decoupled [44]. Such a situation is frequently encountered in the ISS-NMSSM.
The mass matrix of the CP-odd sneutrino fields can be obtained from that of the CP-even fields by the substitution µ X → −µ X and B µ X → −B µ X . Since µ X and B µ X represent the degree of the LNV, their effect on m 33 should be much smaller than the other contributions. In the extreme case of µ X = 0 and B µ X = 0, any CP-odd sneutrino state is accompanied by a mass-degenerate CP-even sneutrino state. Consequently, the sneutrino particle as a mass eigenstate corresponds to a complex field, and it has an antiparticle [78]. Alternatively, if B µ X takes a naturally suppressed value, the mass splitting between the CP-even and CP-odd states is usually tiny, e.g., less than 0.2 GeV for B µ X = 100 GeV 2 and mν 1 ∼ 100 GeV. Such a sneutrino is called a pseudo-complex particle in the literature [73,[79][80][81]. This feature of the sneutrino DM is quite different from that in the Type-I NMSSM.
In the ISS-NMSSM, theν * 1ν 1 h i coupling strength is given by ) denotes the coupling ofν 1 with the scalar field s. For one generation sneutrino case, it is given by where V denotes the rotation matrix to diagonalize the squared mass matrix in Eq. (2.9). This formula indicates that, among the parameters in the sneutrino sector, Y ν , λ ν , A ν , and A λν affect not only the interactions of the sneutrinos but also the mass spectrum and the mixing of the sneutrinos. In contrast, the soft breaking masses m 2 ν and m 2 x affect only the latter property. Given the typical value of the quantities in Eq. (2.11), i.e., tan β 1, 100 GeV. This estimation reflects the fact that |Cν * 1ν 1 Re[S] | is usually much larger than the other two couplings. The basic reason for this is thatν 1 is a singlet-dominated scalar, so it can couple directly with the field S, while in the case of V 11 = 0, the other couplings emerge only after electro-weak symmetry breaking.

DM-nucleon scattering
Since the sneutrino DM in the NMSSM extensions is a singlet-dominated scalar with definite CP and lepton numbers, its interaction with nucleon N (N = p, n) is mediated mainly by the CP-even Higgs bosons h i (i = 1, 2, 3), yielding the effective operator where the coefficient f N is given by [82] and C h i N N denotes the Yukawa coupling of the Higgs boson h i with the nucleon N that relies on the nucleon form factors . This operator does not contribute to the SD cross-section for theν 1 − N scattering, while it predicts the SI crosssection as follows [82]: where µ red = m N /(1 + m N /mν 1 ) is the reduced mass of the nucleon with mν 1 , and the quantities a ui and a di are defined by The expressions of a u i and a d i reveal the following important features of the scattering: 4 In the case that the DM candidate is a Majorana fermion, e.g., the lightest neutralino in the MSSM and NMSSM, the scattering cross-section takes the same form as Eq. (2.13) except that aq i is obtained from Eq. (2.14) by the replacement Cν * 1ν 1 h i /mν 1 → Cχ0 1χ 0 1 h i [44]. This similarity was used to compare the scattering rate in different models [44].
• Although the SD cross-section vanishes, the SI cross-section depends not only on the parameters in the Higgs sector but also on those in the sneutrino sector, which include λ ν , A λν , and mν for the Type-I NMSSM and λ ν , Y ν , A λν , A ν , mν, and mx for the ISS-NMSSM. This feature provides the theories (especially the ISS-NMSSM) with a great deal of freedom to be consistent with the experimental results. In particular, a small |µ| no longer enhance the cross-section, and it is quite often that, after fixing the parameters in the Higgs sector, an experimentally allowed DM candidate can be predicted by only adjusting the inputs of the sneutrino sector. In contrast, in the MSSM or the NMSSM with the lightest neutralino being a DM candidate, the SI and SD cross-sections rely only on the DM mass and the parameters in the Higgs sector. Due to different dependencies of the cross-sections on the parameters and also due to the constraints of the LHC experiments on the parameters, it is not easy to suppress the two cross-sections simultaneously in the light higgsino case [14,15,[17][18][19].
• Each of the h i contributions to the SI cross-section is naturally suppressed. Explicitly speaking, the mass of the Re[H 0 d ]-dominated Higgs particle is usually at the TeV order, so its contribution is suppressed by the large mass. The Re[H 0 u ]-dominated scalar corresponds to the SM-like Higgs boson and its coupling withν 1 is suppressed by the factor λλ ν cos β and the small Y ν . Moreover, as far as the ISS-NMSSM is concerned, the accidental cancellation between the different terms in Cν * 1ν 1 Re[H 0 u ] can further suppress the coupling. In most cases, the contribution from the singletdominated scalar is the most important, but it vanishes if there is no singlet-doublet mixing in the CP-even Higgs sector.
Due to these features, the extensions can easily satisfy the constraints of the DM DD experiments, and this was proven by the Bayesian analysis of the Type-I NMSSM [37].
To further illustrate the behavior of the scattering rate, we consider a special case where Cν * 1ν 1 H d = Cν * 1ν 1 Hu = 0 and m H ± 1TeV. We first integrate out the heavy doublet Higgs field so that the CP-even Higgs sector contains only the SM Higgs field sin βRe[H 0 u ]+ cos βRe[H 0 d ] and the singlet field Re[S]. We then calculate the scattering amplitude by the mass insertion method. The result takes the following form where δ denotes the splitting between m 2 h 1 and m 2 h 2 normalized by the squared mass of the singlet dominated Higgs boson, and θ is the mixing angle of the SM Higgs field and the singlet field. This formula indicates that, besides reducing Cν * 1ν 1 S in the specific parameter space of the sneutrino sector, a small mixing angle obtained by adjusting the parameters in the Higgs sector can suppress the scattering. This small mixing, however, is favored by the Higgs data at the LHC.
At first glance, it appears that the formula in Eq. (2.15) may be applied to the NMSSM with the singlino-dominated neutralino being the DM candidate by the replace-ment Cν * 1ν 1 S /mν 1 → Cχ0 1χ 0 1 S κ. This speculation is incorrect because, in order to obtain the correct density, the neutralino must contain sizable higgsino components, and this will induce the direct coupling of the neutralino with the SM Higgs field. As a result, the scattering cross-section is quite large, even for the case of sin θ = 0 (this is evident from the fact that λ > 2κ, and it is demonstrated in the example of Case II in the Appendix). In the seesaw extensions of the NMSSM, however, the sneutrino DM may naturally correspond to an almost pure singlet field. In this case, its coupling with the SM Higgs field emerges only after electroweak symmetry breaking and is suppressed by the factors λλ ν v d and Y 2 ν v u . This is a significant difference between these the extensions and the NMSSM.
Before concluding the introduction of the theories, we point out that DM physics in the B µ X = 0 case of the ISS-NMSSM is slightly different from the previous description in two aspects. One is that the sneutrino DM corresponds to a complex field, and its anti-particle also acts as a DM candidate with an equal contribution to the relic density. As such, this case is actually a two-component DM theory. The other is that the Z-boson can mediate the elastic scattering of the DM with nucleons, and consequently, it contributes to the SI cross-section. Since the total SI cross-section in such a theory is obtained by averaging over theν 1 N andν * 1 N scatterings and the interferences between the Z and the Higgs exchange diagrams for the two scatterings have opposite signs [91], the SI cross-section can be written as [41] where σ h N is the same as before, and the Z-mediated contributions are given by where G F denotes the Fermi constant and θ W is the weak angle. In this case, σ SI n may differ significantly from σ SI p . Consequently, the effective cross-section of the coherent scattering between the DMs and xenon nucleus (defined as the averaged cross-section σ SĨ ν 1 −Xe /A 2 , where A denotes the mass number of the xenon nucleus) is given by the following general form: where the three coefficients on the right side of the equation are obtained by considering the abundance of different xenon isotopes in nature. This effective cross-section has the property σ SI eff = σ SI N if σ SI p and σ SI n are equal, and it can be compared directly with the bound of the PandaX-II and XENON-1T experiments if the SD scattering cross-section is negligible [91].
Using the micrOMEGAs code [83][84][85], we confirmed that for B µ X 200 GeV 2 , the DM observables, such as its relic density and its current annihilation rate σv 0 , are insensitive to the value of B µ X [44]. However, as far as the SI cross-section is concerned, the predictions for the B µ X = 0 and B µ X = 0 cases differ significantly due to the reason discussed above. ). All of the masses and Z ij in this table are obtained by the setting Y ν = λ ν = 0. Non-zero Y ν and λ ν may slightly alter their values due to sneutrino loop effects. In addition, the properties of h 2 are consistent with the Higgs results from the recent ALTAS analysis, which are based on 80 fb −1 data at the 13 TeV LHC [96].
Numerically speaking, we found that the DM DD experiments require |V 11 | 0.01 for the B µ X = 0 case, since the Z-mediated contribution is usually much larger than the Higgsmediated contribution for sizable V 11 [92]. In contrast, |V 11 | may be as large as 0.1 for the B µ X = 0 case (we will present the results in our forthcoming work). We also confirmed that, when building the model file for the micrOMEGAs with the package SARAH [93][94][95], takingν 1 as a complex field (corresponding to B µ X = 0) or as a real field (corresponding to B µ X = 0) significantly affects the time cost for calculating the relic density: in the B µ X = 0 case, more Feynman diagrams must be calculated, and consequently, the computation is quite time expensive.
Throughout this work, we neglect the Z-mediated contribution by setting B µ X = 100 GeV 2 .

Numerical results
In this section, we compare how the sneutrino DM in the two extensions stays consistent with the bound of the XENON-1T experiment on the scattering cross-section under the premise of predicting the correct density and the photon spectrum of the DM annihilation in dwarf galaxies compatible with the Fermi-LAT observation. We choose the benchmarksetting of Region I in [65] for the parameters in the Higgs sector with detailed information presented in Table 2. This setting predicts m h 1 96 GeV, a TeV magnitude v s (v s ≡ √ 2µ/λ = 1273.5 GeV), and large singlet-doublet Higgs mixing. Thus, the DM-nucleon scattering may be quite large by the expression of the Cν 1ν1 S and the formula in Eq. (2.15). We emphasize that, although the setting is consistent with the latest data of the LHC related to the discovered Higgs boson and the search for extra bosons at the LEP and LHC [65], it is actually a rare case since the mixing angle is usually small in the broad parameter space of the NMSSM after considering the data. However, studying this extreme case is very helpful for improving our understanding of the scattering. In particular, it may reveal the mechanisms that keep the theories to be consistent with the tight XENON-1T bound.
The procedure in our study is as follows. We first construct the likelihood function for the DM physics and perform sophisticated scans over the parameters of the sneutrino sector for either theory by requiring the lightest sneutrino as a DM candidate. We then plot the map of the profile likelihood (PL) [37,104] 5 in different two-dimensional planes to illustrate their features and determine the underlying physics. The likelihood function we adopt is expressed as follows: where L Ων 1 , L DD , and L ID account for the relic density, the XENON-1T experiment and the Fermi-LAT observation of dwarf galaxies, respectively. Their expressions are as follows: • L Ων 1 is Gaussian distributed, i.e., where Ω th denotes the theoretical prediction of the density Ων 1 h 2 , Ω obs = 0.120 represents its experimental central value [1], and σ = 0.1×Ω obs is the total (including both theoretical and experimental) uncertainty of the density.
• L DD takes a Gaussian form with a mean value of zero [97]: where σν 1 −p denotes the theoretical prediction of the DM-proton scattering rate, δ σ is evaluated by δ 2 σ = U L 2 σ /1.64 2 +(0.2σν 1 −p ) 2 , U L σ denotes the upper limit of the latest XENON-1T results on the scattering cross-section at a 90% C.L. [3], and 0.2σν 1 −p parameterizes the theoretical uncertainty of σν 1 −p .
• L ID is calculated by the likelihood function proposed in [98,99] with the data of the Fermi-LAT collaboration presented in [100,101]. 5 In frequentist statistics, the PL for a likelihood function L represents the most significant likelihood value in a specific parameter space. With the two dimensional (2D) PL as an example, it is defined by where Θ = (ΘA, ΘB, · · · ) is a set of parameters that L depends on, and one obtains the maximization by changing parameters other than ΘA and ΘB. In our study, we utilized the package SARAH-4.11.0 [93][94][95] to build the models, the SPheno-4.0.3 code [102] to generate the particle spectrum, and the package MicrOMEGAs 4.3.4 [83,85,103] to calculate the DM observables. We set the soft breaking masses for the first two generation sneutrino fields at 2 TeV throughout this work.

Features of sneutrino DM in Type-I NMSSM
We performed two independent scans over the following parameter space of the Type-I NMSSM 6 0 < mν < 200 GeV, 0.025 < λ ν < 0.5, |A λν | < 1 TeV, (3.4) 6 We note that the SM-like Higgs boson may decay into a massive neutrino pair if it is lighter than about 60 GeV, and the branching ratio may be significantly large since the boson in the benchmark setting contains a sizable singlet Higgs component. To avoid such a possibility, we require λν 0.025, so the decay is kinematically forbidden.
with the MultiNest algorithm [105,106] by taking the prior probability density function (PDF) as uniformly distributed andν 1 as a DM candidate to be CP-even and CP-odd, respectively. With the samples obtained in the scan, we show the map of the PL for the function L DM on different planes in Fig. 1. From the results of the left panels for the CP-even case, following facts are obtained • mν 1 is concentrated on the range from 125 to 135 GeV, which is close to mχ0 1 . Since λ ν is small, the annihilationsν 1ν1 → h 1 h 1 , h 1 h 2 , h 2 h 2 are unimportant (discussed below), and the sneutrino DM mainly co-annihilated with the higgsino-dominated electroweakinos to achieve its measured density. In this case, the relic density is much more sensitive to the splitting between mν 1 and µ than to λ ν .
• λ ν in the 2σ confidence interval (CI) is upper bounded, e.g., λ ν 0.09 for mν 1 128 GeV and λ ν 0.12 for mν 1 134 GeV. This feature is mainly due to the constraints from the XENON-1T experiment, which can be understood as follows. First, for the parameters in Table 2 (in particular the remarkable feature that both the ratio µ/λ and the singlet-doublet Higgs mixing are quite large), the Cν 1ν1 h 1 and Cν 1ν1 h 2 couplings in Eq. (2.5) are contributed to mainly by the terms in the second brackets. Thus, one can conclude from Eq. (2.15) that This approximation reflects that the SI cross-section is roughly proportional to 1/m 2 ν 1 , and it increases monotonically with λ ν . Second, the soft-breaking term of the Type-I NMSSM in Eq. (2.1) and the sneutrino mass matrix in Eq. (2.3) show that the parameter mν does not introduce any interactions, and it is related to mν 1 by As a result, one may substitute mν with mν 1 as an input parameter. In this way, mν 1 no longer relies on λ ν and A λν except that, to make the CP-even sneutrino state as a DM candidate lighter than its CP-odd partner, A λν should be negative and satisfy |A λν | κ/λµ with κ/λµ 100 GeV in our parameters. However, a negative A λν ensures that the second term in the curl brackets of Eq. (3.5) always interferes constructively with the first term to strengthen the constraint of the DD experiment. Finally, we emphasize that a relaxed experimental bound on the scattering crosssection with the increase in mν 1 is also an important factor for determining the upper bound.
• The annihilationν 1ν1 → h 1 h 1 for m h 1 96 GeV is unable to account for the measured DM density. This is because the density requires λ ν 0.17 for the parameters in Table 2 if the channel is fully responsible for the density through the quartic scalar interactions (this conclusion was obtained by the formula of the relic density in [107,108]). Such a large λ ν , however, is strongly disfavored by the XENON-1T experiment.
As for the upper left panel of Fig. 1, it may appear that the dependence of the SI crosssection on the parameter λ ν becomes a step function at λ ν 0.09. A proper understanding of the panel involves the concept of the two dimensional (2D) PL L(λ ν , mν 1 ), which is defined by (footnote 5) In plotting the panel, the maximization over A λν was obtained through the following procedure: we split the λ ν −mν 1 plane into 80×80 equal boxes (i.e., we divided each dimension of the plane by 80 equal bins), fitted the samples obtained in the scan into each box so that the samples in each box corresponded roughly equal λ ν and mν 1 (but A λν might differ greatly), and finally selected the maximum likelihood value from the samples in each box as the PL value. It is obvious that L(λ ν , mν 1 ) reflects the preference of the theory on the parameters λ ν and mν 1 , and for a given point in the λ ν − mν 1 plane, its value represents the capability of the point to account for experimental data. Sequentially one can define the 2D CI as the region on the plane satisfying where χ 2 ≡ −2 ln L(λ ν , mν 1 ), and χ 2 min is the χ 2 value for the best sample (best point) obtained in the scan. In our study, χ 2 min 0 because the DM experimental data are independent and consistent with each other, and the Type-I NMSSM can explain the data well. With this knowledge, one can infer that the CIs are not necessarily contiguous [37,104]. In fact, the panel actually reveals the following conclusions: the 2D CI is mainly located in two isolated parameter islands in the λ ν − mν 1 plane, which are featured by λ ν 0.09, 125 GeV mν 1 135 GeV, and λ ν 0.08, mν 1 ∼ 134 GeV, respectively, and the two islands are connected at mν 1 ∼ 132 GeV with 0.08 ≤ λ ν ≤ 0.09. We confirmed that at the bridge, A λν may vary from −260 to −410 GeV to predict χ 2 ≤ 6.18.
Next, we concentrate on the CP-odd sneutrino DM case, where a q i is given by (discussed below Eq. (2.5)) (3.9) Compared with the CP-even case, the results on the right panels of Fig. 1 show similar features except for three aspects. The first is that λ ν in the co-annihilation region can take a much larger value than the CP-even case. The reason is that A λν in the CP-odd case may be either positive or negative, and consequently, the second term in the curly brackets of Eq. (3.9) can interfere destructively with the first term to weaken the constraint of the DM DD experiment. The second aspect is that the 1σ CI in the λ ν − mν 1 plane includes some separated islands with mν 1 125 GeV and λ ν 0.17. For samples located in these islands, the annihilationsν 1ν1 → h 1 h 1 , h 1 h 2 can account for the measured density and the scattering cross-section can be consistent with the bound of the XENON-1T experiment due to the cancellation. The last aspect is that there is a broad 2σ CI characterized by mν 1 125 GeV and λ ν 0.11. In this region, the correct DM density cannot be achieved by the co-annihilation and the annihilation into Higgs bosons, although the constraint from the DM DD experiment can be satisfied. This region is a compromise of the two constraints that maximizes the PL in Eq. (3.7). The information of the best point for the two cases is as follows: • CP-even case: whereχ 2 ≡ −2 ln L DM , and the fine-tuning quantities ∆ Ω and ∆ σ are defined by 7 where p i denotes any of the input parameters of the theory. These quantities reflect the adaptation of the parameter point to be consistent with the experimental measurements. The larger they are, the more finely tuned the theory will be to coincide with the measurements. From our numerical results, we found that Ω th is most sensitive to the mass splitting betweenν 1 andχ 0 1 , i.e., to mν, λ ν , and A λν . In contrast, σν 1 −p is sensitive only to λ ν and A λν . Given the values of ∆ Ω and ∆ σ for the best points, we concluded that the Type-I NMSSM is natural in DM physics for either case. Moreover, we also calculated the Bayesian evidence Z and found that ln Z = −8.98 for the CP-even case and ln Z = −8.27 for the CP-odd case. Since the Jeffrey's scale was only 0.71, we concluded that the theory does not show a strong preference of the CP-odd case over the CP-even case [109,110].
Before concluding this section, we emphasize that the 1σ CIs of our results for either case also include the regions characterized by mν 1 m h 1 /2 or by mν 1 m h 2 /2. In these regions, the DM was annihilated by a resonant Higgs boson to obtain its measured density. Quite similar to the co-annihilation case, λ ν may be very small, and consequently, σν 1 −p can be as small as 10 −50 cm 2 . These regions, however, have been excluded by the LHC search for the 2τ + E Miss T signal, which is induced by the process pp →χ ± 1χ ∓ 1 , due to the large mass splitting betweenχ ± 1 andν 1 [37]. As such, we do not show them on the panels. 7 In the definition of ∆σ, the unit of σν 1 −p, i.e. 10 −47 cm 2 , represents its current experimental sensitivity. Since the relic density has been precisely measured while the scattering cross-section is only upper bounded, we chose different definitions of ∆Ω and ∆σ to parameterize the theory's ability to account for the experimental results.  Table 3. Information of the benchmark points in the ISS-NMSSM.

Features of sneutrino DM in ISS-NMSSM
Similar to our approach for the Type-I NMSSM, we scanned the following parameter space of the ISS-NMSSM: 0 < Y ν , λ ν < 0.7, 0 < mν, mx < 250 GeV, |A ν |, |A λν | < 1TeV, (3.11) by taking the lightest sneutrino as a DM candidate and requiring that the Yukawa couplings Y ν and λ ν satisfy the unitary constraint in Eq. (2.7). The results of the CIs are projected onto the λ ν − mν 1 and λ ν − σν 1 −p planes in Fig. 2. From this figure, one can learn following facts: • For mν 1 ranging from about 127 to 133 GeV, the DM mainly co-annihilated with the higgsino-dominatedχ 0 1 to achieve the measured density. In this case, σν 1 −p can be suppressed to 10 −51 cm 2 . In contrast, the cross-section is usually larger than 10 −50 cm 2 in the Type-I NMSSM for the same mass range. Moreover, the XENON-1T experiment limits λ ν 0.44, which is significantly weaker than the limitation on λ ν in the Type-I NMSSM 8 .
• Similar to the Type-I NMSSM, 1σ CIs is allowed to be in the mass range 110 GeV mν 1 125 GeV, where the DM obtains its density mainly through the process ν 1ν1 → h 1 h 1 , h 1 h 2 . The difference is that the CI of the ISS-NMSSM is much broader than that of the Type-I NMSSM, and the cross-section can be significantly lower than the prediction of the Type-I NMSSM.
Correspondingly, the DM density requires λ ν ∼ 0.35 (see the formula for the relic density in [107,108]). This phenomenon is absent in the Type-I NMSSM.
The differences reflect the fact that the ISS-NMSSM has greater degrees of freedom for adjusting its parameters to be consistent with relevant experimental constraints than the Type-I NMSSM, and consequently, its DM physics are more flexible. In addition, the left-handed slepton soft-breaking mass ml also plays a role in determining the DM observables through the matrix element V 11 , which is different from the situation for the Type-I NMSSM. In presenting the contours in Fig. 2, we fixed ml = 2 TeV. We verified that the conclusions were not affected by this specific choice. Instead, adopting a smaller ml might slightly improve the fit of the ISS-NMSSM to the DM observables. For example, in the region with mν 1 125 GeV, we obtained χ 2 min 0.45 for the best point with ml = 2 TeV and χ 2 min 0.02 for ml = 400 GeV. In Table 3, we list the detailed information of the best point in Fig. 2 (labelled as Point I with corresponding χ 2 0) and a point with mν 1 96 GeV and χ 2 1.4 (labelled as Point II). From this table, the followings are evident: • For Point I, the DM particle is dominated by thex field in its component, while for Point II, it is an equal mix of theν R andx fields 9 . In contrast, the DM in the Type-I NMSSM is mainly composed of theν R field.
• The values ∆ Ω = 25 and ∆ σ = 0.08 for Point I reflect the insensitivity of these quantities to the input parameters. The implication is that there is a large parameter space around the point that satisfies the experimental constraints. We checked that the mass splitting ofν 1 andχ 0 1 was mainly determined ∆ Ω , which was similar to the best points of the Type-I NMSSM.
• Given ∆ Ω = 167 and ∆ σ = 99, Point II requires considerable tuning to coincide with the experimental results, and thus, it is difficult to obtain in the scan. Our results indicate that, as far as the point is concerned, Ω th is most sensitive to λ ν and A λν , while σν 1 −p is most sensitive to λ ν and Y ν .

Effective natural NMSSM scenario
In the seesaw extensions of the NMSSM, the sparticle's signals may be distinct from those in the NMSSM, and so is the strategy to look for them at the LHC. This feature is relfected in two aspects. One is that, since the sneutrino DM carries a lepton number and has very weak interactions with particles other than the singlet-dominated Higgs boson and massive neutrinos, the decay chain of the sparticles is usually long. Moreover, the decay branching ratio depends not only on the particle mass spectrum but also on new Higgs couplings, such as Y ν and λ ν . As a result, the phenomenology of the sparticles is quite complicated. The other aspect arises because both the DM search experiments and the collider experiments relax significantly their constraints on the extensions [65]. Consequently, broad parameter spaces in the NMSSM, which have been excluded by the experiments, are resurrected as experimentally allowed. In particular, the higgsino mass maybe around 100 GeV to predict the Z-boson mass naturally [37,44]. This fact makes the phenomenology of the extensions quite rich. Despite these differences, we will show in the following that the phenomenology of the extensions may still mimic that of the NMSSM in a particular case, which we dub the effective natural NMSSM scenario (ENNS).
The ENNS contains only the fields of the NMSSM, with its parameter space automatically satisfying the constraints from DM physics and its potentially significant collider signals self-contained in the framework. The following observations motivate this scenario: • From the discussion in previous sections and our Bayesian analysis of the Type-I NMSSM in [37], the co-annihilation of the sneutrino DM with the lightest neutralino is the most important mechanism to obtain the measured DM density 10 .
• The co-annihilation is insensitive to the parameter λ ν , and it is consistent with the constraints from the DM DD and ID experiments given that λ ν is not too large.
• The co-annihilation has the distinct kinematic feature mν 1 mχ0 1 , which implies that the DM mass is roughly determined by mχ0 1 in the NMSSM. In this case,χ 0 1 corresponds to missing momentum at the LHC.
• As indicated by the best points in previous discussions, λ ν and Y ν are preferred to be less than 0.1. In this case, theν R -orx-dominated sneutrino couples very weakly with the other particles. For such a sneutrino DM, it intervenes in the phenomenology of the theories mainly by appearing in the decay chain of heavy sparticles.
The key features of the ENNS are as follows • The lightest neutralino is higgsino-dominated, so the sneutrino DM can co-annihilate with it to obtain the right density. In this case, the higgsino-dominated particlesχ 0 1,2 andχ ± 1 act as the lightest supersymmetric particles of the scenario and are shown as missing momentum at the LHC if the splitting of their masses with the sneutrino 10 Very recently, we performed an analysis of the ISS-NMSSM that was siimlar to our analysis of the Type-I NMSSM. After studying the posterior PDFs of the samples, we determined that this conclusion also applies to the ISS-NMSSM.
DM mass is less than several tens of GeV. This situation is the same as the natural MSSM [9].
The decay of the singlino-dominated neutralino is somewhat complicated. If its higgsino component is small, it may decay dominantly into DM and a massive neutrino, which will subsequently decay into W τ , Zν, and hν states. This case, however, is of less theoretical interest since this neutralino couples weakly with the other heavier sparticles and has little effect on their decay. Alternatively, if the neutralino contains sizable higgsino components, it prefers to decay into the higgsino-dominated neutralino/chargino plus a vector boson or a Higgs boson.
• The decays of the gaugino-dominated particles and the colored sparticles are the same as the prediction of the NMSSM with the lightest neutralino being a DM candidate because they do not interact directly with the sneutrinos.
• The decay modes of charged sleptons are scarcely changed. This can be understood as follows. The extensions predict singlet-dominated particles ν h (massive neutrinos) andν (sneutrinos), so the left-handed slepton has additional decay channelsl ± L → νH ± ,νW ± ,χ ± 1 ν h . Since H ± is preferred to be heavy by the LHC search for charged Higgs bosons and also by B-physics measurements, the decayl ± L →νH ± is usually kinematically forbidden for a moderately lightl L . For the processl ± L →νW ± , it proceeds by the smallν L component inν, so its width is suppressed. The channel l ± L →χ ± 1 ν h is induced by the Yukawa coupling Y ν and can be negligible if |Y ν | is much smaller than the magnitude of the gaugino component inχ 0 1 , which enables the decaỹ l ± L → l ±χ0 1 to be dominant [11]. The right-handed sleptons may also decay into these final states. Compared with the left-handed sleptons, the decays must be proceeded by an additional chiral flippingl R →l L , so their widths are further suppressed.
One may also discuss the decay of the left-handed sneutrino, and the conclusion is that its decay pattern changes little for a small Y ν .
• The singlet dominated Higgs bosons may be light, which is one of the interesting features in the NMSSM. In this case, one can adjust the parameter λ ν to enhance m ν h so that the decay of the Higgs intoν h ν h is kinematically forbidden. In this case, the decay of the Higgs boson is the same as the NMSSM predictions.
We emphasize that the ENNS can be realized by only requiring 2κ λ, and it corresponds to the case of the extensions that is most favored by DM experiments. This situation encourages the study of the phenomenology of the NMSSM without considering the constraints from DM experiments.

Conclusions
With the rapid progress in the DM DD experiments, the sensitivity to DM-nucleon scattering has reached unprecedented precision, which is at the order of 10 −47 cm 2 for the SI cross-section and 10 −42 cm 2 for the SD cross-section. These experiments, together with the fruitful LHC experiments, have tightly limited the light higgsino scenario of popular supersymmetric theories such as the MSSM and NMSSM, which customarily take the lightest neutrino as the DM candidate, and thus deteriorate their naturalness in predicting electroweak symmetry breaking. Therefore, it is necessary to seek new mechanisms to suppress the scattering naturally. It is notable that such a mechanism usually alters the DM properties and consequently changes the phenomenology of the traditional theories. This property may be helpful for the underlying theory to survive the LHC experiments.
After analyzing the fundamental origin of the tight constraints on the MSSM and NMSSM, we realized that, if the DM corresponds to a singlet field under the SM gauge group or at least its singlet component is naturally highly dominant over the other components, the DM-nucleon scattering can be spontaneously suppressed. Based on such observations and our previous work, we pointed out that, in the NMSSM augmented with the Type-I seesaw mechanism or the inverse seesaw mechanism, the lightest sneutrino as a DM candidate can automatically possess this property. Furthermore, due to the simplicity of the theoretical structure, these extensions are economical supersymmetric theories that naturally suppress the scattering. In addition, we showed that the singlet-dominated Higgs boson in this framework plays a vital role in multiple ways: besides generating the higgsino and massive neutrino masses as well as significantly affecting the SM-like Higgs boson mass, it also mediates the annihilation of the DM or acts as the final state of the annihilation. Consequently, the sneutrino in the extensions is a viable WIMP DM.
Although it is known that the scattering is usually suppressed in the framework, we considered a particular configuration in the Higgs sector that was able to enhance the scattering cross-section to improve our understanding of the DM-nucleon scattering. Our study of this dangerous case revealed that the sneutrino DM could co-annihilate with the lightest neutralino to obtain the measured density. In this case, the constraint from the DD experiment is readily satisfied, and the induced tunings in the DM physics are not serious. Our study of the dangerous case also revealed that the DD experiments were able to exclude the annihilation of the sneutrino DM into a pair of singlet-dominated Higgs bosons as the dominant channel for the DM relic density in the Type-I NMSSM. However, they still allowed the process to account for the density in the ISS-NMSSM. This conclusion shows that the ISS-NMSSM is more flexible in DM physics than the Type-I NMSSM, since the ISS-NMSSM has more theoretical parameters and a more complex structure.
We also discussed the phenomenology of the seesaw extensions of the NMSSM and showed that it might be quite similar to that of the NMSSM in the co-annihilation case, which readily satisfies the constraints from DM experiments. This fact implies that one may ignore the constraints of DM physics on the NMSSM when studying its phenomenology.
In summary, the WIMP DM (such as the sneutrino discussed herein) in supersymmetric theories is still a good DM candidate that can be naturally consistent with the DM DD and ID experimental results and the LHC search for sparticles even when the higgsino mass is around 100 GeV. Consequently, relevant supersymmetric theories deserve an intensive study. We comment that, although we assume that the sneutrino DM accounts for the total DM density, the conclusion that the seesaw extensions can naturally suppress the DM-nucleon scattering is valid in the multi-component DM case. In this case, this kind of  [112]. In the left panel, constraints from the CMS search for electroweakinos with 36 fb −1 data at the 13 TeV LHC as well as the estimated exclusion capabilities with 300 fb −1 data at the 14 TeV LHC in the same search are also shown, which correspond to the green region and the orange region, respectively. In the right panel, parameter points on the curve labeled by Case I predict the right DM relic density through the annihilation channelχ 0 1χ 0 1 → tt, while those on the curve labeled by Case II obtain the correct density through the processχ 0 1χ 0 theory is valuable if one of the DM components corresponds to a WIMP, and its scattering with nucleons is found to be tiny.
conditional [10][11][12] and cannot be applied to all the cases that arise, we obtain the values of the elements by numerically diagonalizing the neutralino mass matrix in this appendix. In the theories under consideration, the SI scattering is induced by the SM-like Higgs boson and the other CP-even Higgs bosons. If the latter contribution to the cross-section is negligible, the blind spot condition of the MSSM with a bino-dominated DM can be simplified as sin 2β = −mχ0 1 /µ [14]. In the heavy wino limit, mχ0 1 = M 1 by the formula for mχ0 1 in [12]. Likewise, the blind spot of the NMSSM with a singlino-dominated DM is characterized by sin 2β = mχ0 1 /µ and mχ0 1 = 2κµ/λ in the heavy gaugino limit [12,18]. We will use these formulas in our discussion.

Constraints on blind spot of MSSM
Given sin 2β = −M 1 /µ at the blind spot of the MSSM, the neutralino mass matrix depends only on µ and the gaugino masses M 1,2 . In the left panel of Fig. 3, we plot the SD constraint of the XENON-1T experiment [4] in the mχ0 1 − µ plane of the MSSM (the dark blue region) as well as the constraint of the future LZ experiment [111,112] in the plane (the purple region) by setting M 2 = 5 TeV. These results show that the current bound of the SD cross-section requires |µ| at the blind spot to be larger than about 300 GeV regardless of the DM mass, and the future experiments can further exclude |µ| up to 800 GeV.
As a useful complement to the DD experiments, we also show collider constraints at the 95% confidence level on the blind spot of the MSSM, which arise from the CMS search for electroweakinos through multi-lepton signals of the process pp →χ 0 2,3χ ± 1 at the 13 TeV LHC with 36 fb −1 data [23]. We mark the excluded region with green color in the left panel of Fig. 3. In obtaining this region, we performed detailed simulations of the electroweakino production process in the same way as previously reported [33]. We compared our results with the upper bounds of the cross-section for the electroweakino production, which were provided by the CMS collaboration, and found that they roughly coincided. Exclusion capabilities with 300 fb −1 data at the 14 TeV LHC in the same search were estimated in Fig. 16 of [16], and they are plotted as the orange region in the left panel of Fig. 3. These results show that |µ| has been excluded by the LHC experiment up to about 390 GeV for mχ0 1 = 0 and will be further excluded up to 700 GeV in the future. We emphasize that, in both calculations of the exclusion capability, the branching ratios Br(χ 0 2,3 →χ 0 1 Z) and Br(χ 0 2,3 →χ 0 1 h) were computed in the decoupling limit m A v, and their effects on the signal event number were taken into account in the simulations.

Constraints of SD cross-section on blind spots of NMSSM
We discuss three cases of the natural NMSSM where the channelsχ 0 1χ 0 1 → tt, h s A s , hA s are responsible, respectively, for the measured DM relic density. The common feature of the cases is that the density requires a moderately large |λv/µ| through an explicit or indirect dependence on the combination, and under the circumstance that the h s induced contribution to σ SĨ χ 0 1 −N is negligible, rather strong cancellation between the two terms on the right side of Eq. (1.5) is necessary to render the case consistent with the bound of the SI cross-section (discussed below). This fact is the basic reason that we considered the blind spot of the SI cross-section and studied the SD constraint for the cases when discussing their compatibility with the DM DD experiments.
The neutralino sector of the NMSSM involves the parameters λ, κ, tan β, µ, and the gaugino masses M 1,2 [8]. After taking M 1,2 = 5 TeV and tan β = mχ0 1 /µ = 2κ/λ, the SD cross-section depends only on λ, κ, and µ. Considering that the approximation of theχ 0 1χ 0 1 Z coupling in Eq. (1.6) relies on the combinations λv/µ and mχ0 1 /µ and that a large µ can relax the constraint, we project the SD constraint of the XENON-1T experiment [4] onto the λv/µ−mχ0 1 /µ plane by setting µ = 500GeV, which is shown on the right panel of Fig. 3. This panel indicates that the region of λv/µ 0.15, which corresponds to λ 0.45 and is marked with dark blue color, has been excluded. We also show the exclusion capability of the future LZ experiment in the same plane, which corresponds to the purple region of the figure. The results indicate that the region of λv/µ 0.06 (corresponding to λ 0.18) will be excluded. Alternatively, if we take a lower value of µ, e.g., µ = 300 GeV, the region of λv/µ 0.14 (corresponding to λ 0.23) has been excluded, and a broader region bounded by λv/µ 0.05 (corresponding to λ 0.09) will be excluded.
Next, we consider the three cases mentioned above.
• Case I: the annihilationχ 0 1χ 0 1 → tt is responsible for the DM density. This case usually occurs when mχ0 1 > m t [13,19]. Barring the s-channel exchange of a resonant Higgs boson, the annihilation is dominated by the Z-mediated contribution, and the density requires the coupling of the DM pair with Goldstone boson to satisfy [13]  In the right panel of Fig. 3, we project this correlation onto the λv/µ − mχ0 1 /µ plane to obtain the curve labeled by Case I. This curve shows that, although the parameter points on it are consistent with the SI constraints, they do not satisfy the density and the SD constraints simultaneously. This conclusion has been pointed out explicitly in [19].
As a supplement to the discussion above, we comment on the SI cross-section without the blind spot condition. Eq. (5.2) implies that the first term on the right side of Eq. (1.5) should be larger than 0.1, while the XENON-1T bound on the SI crosssection requires |Cχ0 1χ 0 1 h | 0.02 for mχ0 1 = 300 GeV. Thus, strong cancellation in Eq. (1.5) is necessary and that is why we considered the blind spot for Case I.
• Case II: the annihilationχ 0 1χ 0 1 → h s A s is responsible for the DM density. Without possible resonant contributions, this annihilation proceeds mainly through the singlino-dominatedχ 0 1 in the t/u channel, and it is usually the dominant channel when mχ0 1 < m t [13,32]. Similar to Case I, the relic density requires [13] |Cχ0 1χ 0 1 hs | = |Cχ0 which translates into the correlation λκ(1 − 0.24λ 2 ) 0.067 (5.5) at the blind spot for µ = 500 GeV. In the λv/µ − mχ0 1 /µ plane of Fig. 3, we show the correlation with the curve labelled by Case II. Again, the case cannot satisfy all the constraints of DM physics.
In analogy to the study of Case I, we discuss the situation of the SI cross-section. Eq. (5.4) shows that κ 0.1 for mχ0 1 = 150 GeV is able to predict the correct density. Given mχ0 1 2κµ/λ in case of λv/µ 1 [12], the first term in Eq. (1.5) is approximately 4 √ 2κ 2 v/mχ0 1 0.066. This fact together with the XENON-1T bound on the SI cross-section (i.e., |Cχ0 1χ 0 1 h | 0.016 for mχ0 1 = 150 GeV) again inspired us to consider the blind spot.
• Case III: the annihilationχ 0 1χ 0 1 → hA s is responsible for the DM density. This channel proceeds mainly by the t/u-channel exchange of higgsino-dominated neutralinos [13,32]. The relic density for this case requires [13] λ 3 sin 2β µ 700 GeV 2 , (5.6) which corresponds to λ 2 κ 0.25 at the blind spot for µ = 500 GeV. This correlation needs an even larger λ than Case II, i.e., λv/µ 0.4 for mχ0 1 < m t , and is not shown on the panel of Fig. 3. In addition, regardless of the blind spot condition, Eq. (5.6) implies that the magnitude of the second term in Eq. (1.5) is approximated by 0.25/λ for µ = 500 GeV. Since the SI constraint of the XENON-1T experiment is equivalent to |Cχ0 At the end of this section, we emphasize that the conclusions are based on several hypotheses, i.e., the negligible smallness of the non SM-like Higgs contribution to the SI cross-section, one single annihilation channel responsible for the relic density, and |µ| 500 GeV. For a general case of the NMSSM, these hypotheses are usually violated. In particular, µ may be chosen at 1 TeV if one does not consider the fine tuning. In such a situation, these annihilation channels may still be most important to obtain the density [13,19]. We also emphasize that, although we take µ = 500 GeV as an example to show the tension between the relic density and the SI/SD cross-section, it holds for µ < 500 GeV.