Sneutrino DM in the NMSSM with inverse seesaw mechanism

In supersymmetric theories like the Next-to-Minimal Supersymmetric Standard Model (NMSSM), the lightest neutralino with bino or singlino as its dominant component is customarily taken as dark matter (DM) candidate. Since light Higgsinos favored by naturalness can strength the couplings of the DM and thus enhance the DM-nucleon scattering rate, the tension between naturalness and DM direct detection results becomes more and more acute with the improved experimental sensitivity. In this work, we extend the NMSSM by inverse seesaw mechanism to generate neutrino mass, and show that in certain parameter space the lightest sneutrino may act as a viable DM candidate, i.e. it can annihilate by multi-channels to get correct relic density and meanwhile satisfy all experimental constraints. The most striking feature of the extension is that the DM-nucleon scattering rate can be naturally below its current experimental bounds regardless of the higgsino mass, and hence it alleviates the tension between naturalness and DM experiments. Other interesting features include that the Higgs phenomenology becomes much richer than that of the original NMSSM due to the relaxed constraints from DM physics and also due to the presence of extra neutrinos, and that the signatures of sparticles at colliders are quite different from those with neutralino as DM candidate.


Introduction
From recent cosmological and astrophysical measurements with unprecedented precision, it has been a robust fact that over 20% of the energy density of the Universe today is composed of Dark Matter (DM) [1]. Among various kinds of DM candidates, the massive neutral stable particle with weak couplings to quarks/leptons is a promising one, and has been widely discussed in different new physics models for past decades. In the popular supersymmetric models such as the Minimal Supersymmetric Standard Model (MSSM) [2,3], the lightest neutralino with bino field as its dominant component has such properties [4][5][6], and is customarily treated as DM candidate in phenomenological study. In this setup, the interactions of the DM with Higgs bosons are inversely proportional to higgsino mass µ [7], and the lighter the higgsino is, the stronger the couplings become. This in return results in an increased DM-nucleon scattering rate 1 . On the other hand, the higgsino mass determines at tree level the Z boson mass, and naturalness favors light higgsinos up to several hundred GeV [8]. Obviously, with the rapidly improved sensitivity of DM direct detection (DD) experiments such as PandaX-II [12,13], LUX [14] and XENON-1T [15] to DM-nucleon scattering rate in recent years, there emerges increasing tension between naturalness and the DD experiments [16,17]. Confronted with such a situation, some authors recently emphasized the role of blind spots in escaping the strong constraints from the DD experiments [18][19][20][21]. These parameter points, however, require subtle cancelation among different contributions to the scattering rate, and hence lead to a certain degree of fine tuning. Another long-standing problem that the MSSM fails to account for comes from the observation of neutrino oscillation, which can be explained only if neutrinos have tiny masses [22][23][24]. Given the fact that the MSSM with R-parity conservation has no built-in mechanism to generate the masses, neutrino oscillation indicates unambiguously the existence of extra physics.
In this work, we intend to seek for the theory that can address the origin of neutrino mass and the nature of DM simultaneously. To be more specific, we require it to have following properties: • predicting in a natural way the masses of active neutrinos and also the recently discovered Higgs boson with its field content as economical as possible; • providing a testable mechanism to generate sterile neutrino masses; • easily satisfying the experimental data such as the neutrino oscillation data, the electroweak precision measurements, and the lepton-flavor violation; 1 We emphasize here that we only consider the case of one-component DM with its mass at electroweak scale. In this case, the lightest neutralino in the MSSM is the admixture of gaugino and higgsino with bino as its largest component field in order to predict the right relic DM density. Alternatively if the lightest neutralino is an almost pure higgsino which can be realized in natural SUSY [8] or mirage mediation scenarios [9], its current density will fall far short to account for the measured value of DM density [8,9], and meanwhile the corresponding DM-nucleon scattering rate is usually suppressed too [10]. Note that the tendency of a light µ to enhance DM-nucleon scattering rate is also applied to the NMSSM where the lightest neutralino is usually bino-dominated or singlino-dominated [11].
• easily coinciding with the observations in DM physics even for light higgsinos, especially that DM-nucleon scattering rate should be naturally suppressed to satisfy the very tight constraints from the recent XENON-1T experiment.
In constructing such a theory, we note that among the ideas to generate the tiny neutrino masses, the inverse seesaw mechanism [25] is rather attractive since in its configuration, the smallness of neutrino masses is attributed to lepton number violation (LNV) and a doubly suppressed ratio, all involved dimensional parameters are at weak scale and the Yukawa couplings of the neutrinos may be moderately large, all of which indicate that the mechanism can provide a natural, simple and testable way to realize the small neutrino masses at low energy [26]. We also note that the gauge singlet extensions of the MSSM like the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [27] have great theoretical advantages, e.g. their capability of generating dynamically the higgsino mass µ and enhancing the SM-like Higgs boson mass by the singlet-doublet interaction among the Higgs fields in the theory and/or by the singlet-doublet Higgs mixing effect [28,29]. These features motivate us to incorporate the inverse seesaw mechanism in the NMSSM as an attempt at weak scale to solve the problems mentioned above. Interestingly, we find that the resulting theory not only inherits all the merits of the NMSSM and the seesaw mechanism, but also exhibits following new features: • Except for the tiny Majorana masses for extra family of sterile neutrino fields, which violates lepton number by two units and is naturally small according to 't Hooft's naturalness criterion [30], there is no dimensional parameters in its superpotential. As a result, the mass for any new particle beyond the Standard Model, such as sterile neutrinos and supersymmetric particles, is determined by the vacuum expectation values (vev) of Higgs fields and/or by soft supersymmetry breaking coefficients .
• The lightest sneutrinoν 1 may act as a viable DM candidate. In more detail, unlike some pioneer studies in this direction [31][32][33][34][35][36][37][38][39][40], the sneutrino DM in our framework has two attractive characters. One is thatν 1 can mainly annihilate into a pair of singlet dominated Higgs bosons to get the right relic density and meanwhile satisfy all experimental constraints. This process is determined by the interactions ofν 1 with the singlet Higgs fields for a given Higgs boson spectrum, and consequently DM observables are sensitive only to the parameters in sneutrino sector. The other is owe to the fact that the singlet field can mediate the transition betweenν 1 pair and the higgsino pair, which implies thatν 1 and the higgsinos can be in thermal equilibrium in early Universe before their freeze-out. If their mass splitting is less than about 10%, the number density of the higgsinos can track that ofν 1 during freeze-out, and consequently the higgsinos played an important role in determining DM relic density [41] (in literature such a phenomenon was called coannihilation [42]). As a result, even for very weak couplings ofν 1 with SM particles,ν 1 may still reach the correct relic density by coannihilating with the higgsino-dominated particles. Again, this translates to the constraints only on the parameters in sneutrino sector if the higgsino mass is less than the other neutralino masses.
Due to the mentioned properties ofν 1 , the DM-nucleon scattering rate in our model can be naturally suppressed by the small mixing between singlet-doublet Higgs fields (corresponding to the former case) or by the highly sterile nature ofν 1 (the latter case). This suppression is independent of the parameter µ, and hence there is no tension any more between the weak scale naturalness and DM physics.
• Due to potentially relaxed DM constraints on the theory and also due to the presence of possible light sterile neutrinos, Higgs physics is enriched greatly compared with that of the unextended NMSSM. Moreover, the signature of sparticles at colliders is greatly changed for sneutrino DM instead of the customary neutralino DM.
With respect to these features, we have more explanations. One is that in the original MSSM and NMSSM, only left-handed sneutrinos are predicted, and consequently the sneutrinoν 1 as DM candidate is excluded by DD experiments due to its sizable coupling with Z boson [43]. In Type-I seesaw extended models, however, sneutrino may be a viable DM because the inclusion of right-handed (RH) neutrino superfields in the theory enables the DM to be RH sneutrino dominated, which can reduce the coupling strength greatly [44]. In the inverse seesaw extension, beside the RH fields an extra family of sterile neutrino fields are also introduced, which is able to further suppress the left-handed sneutrino component ofν 1 to get a smaller DM-nucleon scattering rate. In fact, this is one of the reasons that we are interested in the incorporation of the inverse seesaw mechanism within supersymmetric theories. The other is that the features mentioned above are not unique to the inverse seesaw extension of the NMSSM. In fact, the Type-I seesaw extension of the NMSSM also possesses these properties, and in particular it has an advantage over our framework in that it corresponds to a more economical field assignment [44]. However, as we will discuss at the end of this work, our framework provides more flexibility to accommodate low energy data (such as the neutrino oscillation data, the electroweak precision measurements and the lepton-flavor violation) and richer phenomenology than the Type-I seesaw extension, which make it worthy of an intensive study.
The main purpose of this work is to illustrate the properties of the sneutrino DM in the NMSSM with inverse seesaw mechanism (ISS-NMSSM). For this end, we vary the parameters in sneutrino sector to obtain physical parameter points, and show howν 1 annihilated to get the right relic density and meanwhile avoids the constraints from Fermi-LAT search for DM annihilation in dwarf galaxies. In particular, we pay great attention to study DM-nucleon scattering, and exhibit suppression mechanisms of the theory on the rate. We note that the inverse seesaw mechanism has been intensively studied in the framework of the MSSM [45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64] and the supersymmetric B-L models [65][66][67][68][69][70][71][72][73][74][75][76], and that most seesaw extensions of the NMSSM focused on the augmentation of simple Type-I mechanism to study the spectral characters of gamma-ray from DM annihilations [44,[77][78][79][80][81][82][83][84][85]. These studies usually concentrated on the parameter region which predicts a large DM-nucleon scattering rate and hence has been excluded by current DD experiments. By contrast, only several works have been done to study the theory and phenomenology of the ISS-NMSSM [86][87][88][89][90][91]. In particular, we note that only the work [87] adopted same symmetries as our model, and it studied the effect ofν 1 on the properties of a O(10GeV) CP-odd Higgs boson. This situation necessitates our study as a helpful attempt to explore the nature of DM. Obviously, our result on DM physics may be distinct from the previous ones since they are based on different theoretical assumptions and also for different purposes. This paper is organized as follows. In section 2, we introduce the basics of the ISS-NMSSM, including the annihilation mechanisms of sneutrino DM and the features of the spin-independent (SI) cross section for DM-nucleon scattering. In Section 3 we scan the parameter space of the model by considering relevant experimental constraints to get viable parameter points, and analyze numerically the key features of sneutrino DM. In section 4, we study the constraints of the LHC experiment on our choice of the NMSSM parameters. For this purpose, we simulate the neutralino/chargino production processes, and point out that current experimental analyses on sparticle search can not exclude the light higgsinodominated particles due to their unconventional signatures. Section 5 is devoted to a brief exploration of the phenomenology of the ISS-NMSSM, and we will show that the phenomenology is quite rich and distinct. Finally, we draw our conclusions in section 6.

NMSSM with Inverse Seesaw Mechanism
In this section we first introduce the basics of the ISS-NMSSM, including its Lagrangian and neutrino physics, then we concentrate on sneutrino DM case. We analyze the features of sneutrino mass matrix, and present useful formula to calculate the cross sections for DM annihilations and also that for DM-nucleon scattering.

Model Lagrangian
Depending on field assignment and the symmetry adopted in model construction, there are various ways to implement the inverse seesaw mechanism in the NMSSM [86][87][88][89][90][91]. Here we consider the minimal framework which extends the NMSSM by only two gauge singlet chiral fields ν and X for each generation with lepton numbers L = −1 and L = +1 respectively. We assume that the lepton number L and Z 3 symmetry are broken slightly, while the R-parity and (−1) L parity are still good symmetries. With these assumptions, we write down the theory of the ISS-NMSSM with its field content presented in Table 1, and its superpotential and corresponding soft breaking terms given by [87] In above formulae, the coefficients λ and κ parameterize the interactions among the Higgs fields, Y f (f = u, d, e, ν) and λ N are Yukawa couplings for quarks and leptons, m i (i = u, d, · · · ) denote soft breaking masses, and A i are soft breaking coefficients for trilinear terms. About the Lagrangian in Eq.(2.1), five points should be noted. First, the terms in the first bracket of the superpotential W correspond to that of the NMSSM with Z 3 symmetry [27], and those in the second bracket are for the newly added neutrino superfields. The expression of L sof t has same structure. Second, we have neglected flavor indices in writing down the expressions for the sake of simplicity. So all the parameters except for those in the Higgs and gaugino sectors are actually 3 × 3 (diagonal or non-diagonal) matrices in flavor space. Third, among the parameters in the superpotential only µ ν and µ X are dimensional. These coefficients parameterize the effect of LNV, which may arise from the integration of heavy particles in an ultraviolet high energy theory with LNV interactions (see for example [86,88] and also discussions in [87]), so the magnitude of their elements should be suppressed. Similarly, the coefficients B µν and B µ X tend to be small. Fourth, the fields H 0 u,d and S acquire their vevs after electroweak symmetry breaking, i.e. H 0 These vevs are related with the soft breaking squared masses m 2 Hu , m 2 H d and m 2 S by the minimization conditions of the Higgs potential [27], and in practice one may use m Z , tan β ≡ v u /v d and µ = λ √ 2 v s instead of the squared masses as input parameters of the ISS-NMSSM. Finally, we emphasize that the last two terms in the W can induce three/four scalar interactions involving sneutrinos and Higgs bosons, and their corresponding soft breaking terms induce only three scalar interactions. These interactions, as we mentioned before, play an important role in DM physics. We also emphasize that the Yukawa coupling Y ν can introduce extra interactions for the superfieldŝ l andĤ u , and consequently the signature of left-handed sleptons and higgsinos at the LHC may be altered greatly.
Obviously the Higgs sector of the ISS-NMSSM is same as that of the NMSSM. In this work, we adopt the convention of the NMSSM that h i with i = 1, 2, 3 (A j with j = 1, 2) denote mass eigenstates of CP-even Higgs bosons (CP-odd Higgs bosons) with their mass satisfying m h 1 < m h 2 < m h 3 (m A 1 < m A 2 ). Since this sector has been introduced in detail in [27], we in the following only consider the neutrino and sneutrino sectors. As we will show below, the singlet Higgs fields can play an important role in these sections.

Neutrino Sector
In the ISS-NMSSM, the neutrino Yukawa interactions take the following form and they generate the neutrino masses after the electroweak symmetry breaking. In the interaction basis (ν L , ν * R , x), the 9 × 9 neutrino mass matrix reads with the 3 × 3 Dirac mass matrices given by Since this mass matrix is complex and symmetric, it can be diagonalized by a 9 × 9 unitary matrix U ν according to This gives three light neutrino mass eigenstates and six heavy neutrino mass eigenstates, which are related with the interaction state ν by ν m = U ν ν. Without loss of generality, the matrix U † ν can be decomposed into the blocks where the 3 × 3 matrix U is responsible for the oscillations of active neutrinos, and the value of its elements can be extracted from relevant experimental data. With the definition M ≡ Tr(M † M ) for an arbitrary matrix M and in the limit µ X , µ ν M D M R , one can extract the mass matrix of the light active neutrinos from the expression in Eq.(2.3), which is given by and the magnitude of its elements is of the order M D / M R . So in inverse seesaw mechanism, the smallness of the active neutrino masses is not only due to the small elements of the lepton-number violating matrix µ X , but also due to the suppression factor This usually leads to observable lepton flavor violation (LFV) signals as discussed in literatures [45,47,48,50,58,60,62,92]. Note that although both µ X and µ ν are naturally small, µ X controls the size of the light neutrino masses, while µ ν is irrelevant. In view of this, for the sake of simplicity we set the matrix µ ν (and also its soft breaking parameter B µν ) to be zero and do not discuss its effect any more. Also note that the mass scale of the heavy neutrinos is determined by the magnitude of M R .
The symmetric effective light neutrino mass matrix M ν can be diagonalized by the unitary Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix where m ν 1 , m ν 2 and m ν 3 are the masses of the three lightest neutrinos. Generally speaking, due to the mixings among the states (ν L , ν * R , x), the matrix U in Eq.(2.5) does not coincide with U PMNS , instead in the limit µ X M D M R , they are related by In this sense, η = 1 2 F F † is a measure of the non-unitarity of the matrix U , which is obtained from neutrino experiments. On the other hand, since current experiments have tightly limited the violation of the unitarity [93], one can assume U PMNS U and use the data of neutrino experiments to limit the parameters in M ν . So far two parameterizations schemes are adopted in literature (see for example [94]) in doing this. One is to express the Yukawa coupling matrix Y ν in terms of U PMNS by using a modified Casas-Ibarra parameterization [96], which is given by and R is a complex orthogonal matrix given by where c i ≡ cos θ i , s i ≡ sin θ i and θ 1 , θ 2 , and θ 3 are arbitrary angles. In this scheme, the neutrino Yukawa coupling Y ν is usually flavor non-diagonal. The other scheme utilizes the fact that once the matrix Y ν and M R are given, µ X alone can be responsible for neutrino experimental data. In this case, µ X is given by [94,95] Note that for this scheme, one may set Y ν and λ N to be flavor diagonal, and this choice can simplify greatly our study on the properties of sneutrino DM (see following discussion about sneutrino mass matrix).

Sneutrino Dark Matter
In the ISS-NMSSM, the lightest sneutrinoν 1 may be a better DM candidate than the customary lightest neutralino after considering the negative result in recent DM DD experiments, which is the main standpoint of this work. In the following, we will present in detail the features ofν 1 , including its mass, its annihilation channels as well as its scattering with nucleon.

Sneutrino mass matrices
After decomposing sneutrino fields into CP-even and CP-odd parts with i = 1, 2, 3 representing flavor index, one can write down the mass matrix for the CP-odd sneutrinos in the basis σ j (j = 1, · · · 9) as follows (2.14) where and all the m ij are 3 × 3 matrices in flavor space. From the expression of m 2 ν I , one can get following conclusions • In the case of no flavor mixing in the matrix m ij , which can be obtained by neglecting the small flavor non-diagonal matrix µ X presented in Eq.(2.12) (and also the coefficient of the bilinear term B µ X ) and is the situation considered in this work, one can rearrange the basis σ j by the order (σ 1 , σ 4 , σ 7 , σ 2 , σ 5 , σ 8 , σ 3 , σ 6 , σ 9 ) so that m 2 ν I is flavor diagonal. In this case, there are only the mixings between (ν L ,ν R ,x) for same generation sneutrinos. If the lightest sneutrino comes from a certain generation, e.g. the third generation, and at same time it is significantly lighter than the other generation sneutrinos, one only needs to consider the mass matrix for this generation sneutrinos in discussing the properties of the DM 2 . This will greatly simply our analysis. In the following, we only consider one generation of sneutrinos in studying the property ofν 1 .
• Among the parameters in sneutrino sector, Y ν , A ν , λ N and A λ N affect not only the interactions of the sneutrinos, but also the mass spectrum of the sneutrinos. By contrast, the soft breaking masses m 2 ν and m 2 x and the small bilinear coefficient B µν only affect the spectrum. Considering that the former four parameters are tightly limited by various experiments (see below), one can conclude that the spectrum is mainly determined by the soft breaking masses for heavy sneutrino case; on the other hand, since v s is usually much larger than v u , the spectrum is more sensitive to λ N and A λ N than to the other parameters for the case of light sneutrinos, mν i ∼ v u .
• The mixing ofν L with the other fields is determined by the parameters Y ν and A ν .
In the limit Y ν = 0, m 12 and m 13 vanish, and consequentlyν L does not mix withν R andx any more. Furthermore, if the first term in m 22 is far dominant over the rest terms in m 22 and so is m 33 , m 22 m 33 and this results in a maximal mixing betweeñ ν R andx. In this case,ν 1 is approximated byν 1 . This situation is frequently encountered in our results.
In a similar way one may discuss the mass spectrum of the CP-even sneutrinos. We find that their mass matrix m 2 Since the Majorana mass µ X and the bilinear coefficient B µ X reflect the effect of LNV, they should be suppressed greatly. In the limit µ X = 0 and B µ X = 0, any CP-even sneutrino particle must be accompanied with a mass-degenerate CP-odd sneutrino. In this case, one may say that the sneutrino as an mass eigenstate corresponds to a complex field, and it has its anti-particle [54]. If alternatively B µ X takes a moderately small value and consequently the mass splitting between the CP-even sneutrino particle and its corresponding CP-odd particle is at eV order, one may call such a sneutrino pseudo-complex particle. This case has interesting implication in DM physics [59,88].
In this work, we only consider the case that B µ X is moderately large, B µ X = 20GeV 2 , so that the CP-odd state is lighter than its corresponding CP-even state by ∼ 0.1GeV, and sneutrinos as mass eigenstates have definite CP number. We note that the lightest CP-even sneutrinoν R 1 can decay intoν 1 γ with a width around the order of 10 −8 GeV, and it usually coannihilated in early Universe withν 1 to get the right DM relic density. We numerically checked by the code micrOMEGAs [97][98][99] that the observables such as the relic density and DM-nucleon scattering rate discussed in this work are insensitive to the choice of B µ X .

Relic density of sneutrino DM
In the cosmological standard model, the abundance of a thermal DM Y (T ) is defined as the number density divided by entropy density s(T ), and its Boltzmann equation is [100] where g * is an effective number of degrees of freedom (dof) derived from thermodynamics describing state of the Universe, M p is Plank mass, Y eq is thermal equilibrium abundance, and σv is the relativistic thermally averaged annihilation cross section with v denoting the relative velocity between the annihilating particles. With the aid of present day abundance Y (T 0 ), the DM density today can be written as [100] where s(T 0 ) is the entropy density at present time and h is the normalized Hubble constant.
These formulae indicate that in order to get the right relic density, one has to solve the evolution equation of Y (T ), which is usually a complicated work and has to be done numerically.
As far as the ISS-NMSSM is concerned, its influence on the relic density ofν 1 enters through the cross section σv , which includes all annihilation and coannihilation channels predicted by the model, and is given by [ where g i is the number of dof, σ ij;kl is the cross section for the annihilation of a pair of supersymmetric particles with masses m i , m j into SM particles k and l, p ij is the momentum of incoming particles in their center of mass frame with total energy √ s, and K 1 and K 2 are modified Bessel functions. In practice, the potentially important contributions to σv include following annihilation channels (1)ν 1ν1 → ss with s denoting either a CP-even or CP-odd singlet dominant Higgs boson.
This annihilation proceeds via a four-point scalar coupling, s-channel mediation of a Higgs boson and t/u-exchange of a sneutrino, which are depicted in Fig.1 for the case that s is the lightest CP-odd Higgs A 1 .
(2)ν 1ν1 → ηη * with η denoting a SM particle or any of the heavy neutrinos. This annihilation is mediated by any of the CP-even Higgs bosons, and since the involved interactions are usually weak in getting the right relic density, one of the bosons must be at resonance. (1) and (2). Note thatν R 1 plays an important role in determining the relic density sinceν R 1 is always nearly degenerate withν 1 in mass.
withχ denoting a higgsino-like or wino-like electroweakino. These annihilations are called coannihilation in literature [41,42], and to make the effect significant, the mass splitting betweenχ andν 1 should be less than about 10%.
In the following, we consider for illustration purpose the cross section of the annihilation channel shown in Fig.1 with collision energy √ s, which is given by (2.21) In above formulae, a i = s 2 + m 2 , Cν 1ν1 A 1 A 1 denotes the coupling of twoν 1 s with two A 1 s which is mainly determined by the parameter λ N , and the other coefficients C XY Z are the triple scalar couplings involving the particles X, Y and Z. In getting the approximation in Eq.(2.19), we use the relation s = 16m 2 ν 1 /(4−v 2 ) with v denoting the relative velocity of the twoν 1 s, and assume following conditions: (1) v ∼ O(0.1), which means that the collision is nonrelativistic; (2)ν 1 is much heavier than m A 1 ; (3) mν 1 /|2mν 1 − m h i | is at most a O(1) quantity, which excludes the possibility that the mediating Higgs boson is resonant. With these conditions, the coefficient b is usually smaller than the coefficient a.
The thermal averaged cross section of the annihilation at freeze-out temperature T f and that at present time are then given by [42] σv T f a + 6b This implies that, if the annihilationν 1ν1 → A 1 A 1 is fully responsible for current relic density so that σv Obviously, such a large σv 0 is tightly limited by the Fermi-LAT search for DM annihilation from dwarf spheroidal galaxy (dSph). In order to avoid the constraint, one may consider following cases as pointed out by the classical paper [42] • Coannihilation, or more general mixed annihilations. In this case, the annihilatioñ ν 1ν1 → A 1 A 1 plays a minor role in contributing to the relic density, and consequently σv 0 can be lowered significantly.
•ν 1 is slightly lighter than A 1 , which is called forbidden annihilation in [42]. In this case, since the freeze-out occurs at a temperature T f mν 1 /25, and also since the velocity ofν 1 is Boltzmann distributed with the temperature, the annihilation may proceed in early Universe, but can not occur at present time. So σv 0 is suppressed greatly.
• Resonant annihilation mediated by h i as the main contribution to the relic density.
In this case, σv 0 can be significantly lower than σv T f if 2mν 1 < m h i [42,101].
As we will show below, these cases are frequently encountered in our scan over the parameter space of the ISS-NMSSM to escape the constraints from the dSph. Throughout this work, we use the package micrOMEGAs [97][98][99] to evaluate observables in DM physics, including the relic density, photon spectrum from DM annihilation in the dSph which is used for DM indirect detections, and also the cross sections for DM-nucleon scattering. The package solves the equation for the abundance in Eq.(2.16) numerically without any approximation. In addition, it also estimates the relative contribution of each individual annihilation or coannihilation channel to the relic density at the freeze-out temperature.

Direct detection
Since the DMν 1 in the ISS-NMSSM is a scalar with certain lepton and CP numbers, its interaction with nucleon N (N = p, n) is mediated only by CP-even Higgs bosons to result in the effective operator In this formula, C h i N N is the Yukawa coupling of the Higgs boson h i with nucleon N , S ij is the (i, j) element of the matrix S which is used to diagonalize the CP-even Higgs mass matrix in the basis (H d , H u , s), F . This operator indicates that the spin-dependent cross section forν 1 scattering with proton vanishes, whereas the SI cross section is given by [102] σ SĨ where µ red = m p /(1 + m 2 p /m 2 ν 1 ) is the reduced mass of proton with mν 1 , and the quantities a ui and a di are defined by to facilitate our analysis. By contrast, we note that a ui and a di in the MSSM with the lightest neutralinoχ 0 1 acting as DM candidate take following form [85] where Cχ0 (2.26) We will return to this issue later. In our numerical calculation of σ SĨ ν 1 −p , we use the default setting of the package mi-crOMEGAs [97][98][99] for the nucleon form factors, σ πN = 34MeV and σ 0 = 42MeV, and obtain F (p) u 0.15 and F (p) d 0.14 3 . In this case, Eq.(2.23) can be approximated by 3 We remind that different choices of σπN and σ0 can induce an uncertainty of O(10%) on F d . For example, if we take σπN = 59MeV and σ0 = 57MeV, which are determined from [103] and [104] respectively, we obtain F where Cν 1ν1 S (S = H d , H u , s) denotes the coupling ofν 1 with the scalar field S, and for one generation sneutrino case it is given by Z denoting the rotation matrix to diagonalize the CP-odd sneutrino mass matrix and consequentlyν i = Z i1νL + Z i2νR + Z i3x .
In the following, we analyze the features of σ SĨ ν 1 −p . From Eq.(2.27), we learn that the dependence of σ SĨ ν 1 −p on the parameters of the ISS-NMSSM comes from the expression in the bracket, which is quite complicated. To simplify the analysis, we assume that the left-handed sneutrino component inν 1 is suppressed greatly, e.g. |Z 11 | 0.05, and that Then the couplings Cν 1ν1 S can be approximated by For this case, S 12 sin β ∼ 1, S 11 cos β, a u1 a d1 , and (2.30) This formula indicates that the cross section is determined by Y ν and T ν ≡ Y ν A ν , and may be suppressed ifν 1 isx dominated. We remind that a small Y ν is not only favored by the recent XENON-1T constraints on σ SĨ ν 1 −p , but also consistent with the limitation on the non-unitarity of the U matrix in neutrino sector.
As a comparison, one may also discuss the DM-nucleon scattering rate in the MSSM, which can be obtained from σ SĨ  Taking mν 1 = mχ0 1 m 1 , we conclude that the ratio is about 10×(Y ν Z 12 µ/m 1 ) 4 after neglecting unimportant terms. This fact indicates that σ SĨ ν 1 −p in the ISS-NMSSM can be easily much lower than that in the MSSM.
II. h 1 acts as the SM-like Higgs boson, h 2 is singlet dominated with m h 2 v, and h 3 is decoupled.
In this case, S 12 sin β ∼ 1, S 11 cos β and S 23 ∼ 1. At same time, |S 13 | and |S 22 | are usually moderately larger than |S 21 |, but all of them should be less than about 0.1 to coincide with the 125GeV Higgs data. Consequently, a d1 a u1 , a d2 is much larger than a u2 since tan β 1, and a d3 , a u3 0 since they are suppressed by 1/m 2 h 3 . σ SĨ ν 1 −p is then given by which indicates that the magnitude of σ SĨ ν 1 −p is partially decided by the mixing S 13 . The implication of this special case is that the interaction ofν 1 with the singlet fields alone can be responsible for DM physics in the ISS-NMSSM, namely predicting correct relic density and also possibly sizable DM-nucleon scattering cross section.
We remind that Eq.(2.33) also holds if the element S 13 is not suppressed too much so that S 13 Cν 1ν1 s Cν 1ν1 Hu . We will study in detail this situation later.
-In Eq.(2.32), the first term comes from the interchange of h 1 , and the second term denotes the contribution of h 2 . These two contributions are usually comparable in size since a d1 /a d2 ∼ O(1), and in some cases the latter may be more important. We will show that the two contributions may interfere destructively or constructively in contributing to the cross section. Note that since we have set the higgsino mass at 120GeV, which is motivated by naturalness argument, all the higgsino-dominated particles such asχ 0 1 ,χ 0 2 andχ ± 1 are light with mass around 120GeV. Also note that the masses for the Higgs bosons are slightly altered by the parameters in sneutrino sector through loop effects, so their values in this table are actually obtained for the case of Y ν = λ N = 0.

Numerical Results
In this section, we study the property of the sneutrino DMν 1 by presenting some numerical results. In order to illustrate the underlying physics as clearly as possible, we first fix the parameters in the NMSSM sector, and give in Table 2 the values of some quantities which are relevant to our study. Then we adopt the Metropolis-Hastings algorithm 4 implemented in the code EasyScan HEP [107] to scan following parameter space in sneutrino sector 5 In the calculation, we utilize the package SARAH-4.11.0 [108][109][110] to build the model and the code SPheno-4.0.3 [111] to generate the particle spectrum, and we consider following constraints 4 To be more explicit, we adopt the likelihood function L = Lm h 1 × L Ωh 2 × L Br(B→Xsγ) × L Br(Bs→µ + µ − ) for the Markov Chain Monte Carlo scan where Lm h 1 , L Ωh 2 , L Br(B→Xsγ) and L Br(Bs→µ + µ − ) are likelihood functions for experimentally measured SM-like Higgs boson mass, DM relic density, Br(B → Xsγ) and Br(Bs → µ + µ − ) respectively, which are taken to be Gaussian distributed [105,106]. 5 Since we concentrate on the property ofν1 instead of on neutrino oscillations, we set µX = 0 for simplicity, and only consider the effects of the third generation sneutrinos by setting Yν = 0 and the diagonal elements of mν and mx at 1TeV for the other two generation sneutrinos. With such a treatment, λN in Eq.  • 123GeV ≤ m h 1 ≤ 127GeV, which is the most favored range of the SM-like Higgs boson mass by current LHC results [112]. This constraint arises from the fact that the parameters in sneutrino sector can alter the Higgs boson mass spectrum through loop effects [54,57,67].
• consistence of the Higgs properties with the data from LEP, Tevatron and LHC experiments. This is due to the consideration that the non-standard neutrinos may serve as the decay products of the SM-like Higgs boson, and thus change the branching ratios of its decay into SM particles, and also the consideration that the moderately light h 2 may induce sizable signals at the colliders. We implement the requirement by the packages HiggsBounds-5.0.0 [113] and HiggsSignal-2.0.0 [114].
• low energy flavor observables, such as B → X s γ, B s → µ + µ − and ∆M Bs , and muon anomalous magnetic momentum within 2σ range around its experimental central value. These observables can be calculated automatically by the code SPheno-4.0.3 under the instruction of the package SARAH-4.11.0.
• mν 1 < mχ0 1 , and 0.107 < Ωh 2 < 0.131 in order to account for the Planck measurement of DM relic density at 2σ level [1].  As we introduced in last section, the sneutrino sector of the ISS-NMSSM provides great flexibility to account for DM physics. In our study, we consider about three thousand samples obtained from the scan with the constraints considered. We find |Z 11 | < 0.1 for all the samples, and |Z 12 | |Z 13 | 1/ √ 2 for a sizable portion of the samples. In Fig.2, we project the samples on σv 0 − mν 1 plane with the dark blue color to represent those that satisfy both the dSph constraint and the XENON-1T constraint, and the lime color, the cyan color and the golden yellow color to denote those which are excluded by either the XENON-1T constraint or the dSph constraint, or the both respectively. In implementing the dSph constraint, we use the data provided by Fermi-LAT collaboration [115], and adopt the likelihood function proposed in [116,117], while in imposing the XENON-1T constraint, we use directly the 90% exclusion limits on the SI DM-nucleon scattering cross-section of the recent XENON-1T experiment [15].
From Fig.2, one can infer that for the samples around the green, red and yellow vertical lines,ν 1 annihilated in early Universe mainly by the resonant Higgs mediated processes ν 1ν 1 , ff , W + W − respectively, and for the samples near the gray line, it achieves right relic density mainly by the annihilation of the higgsinos. Moreover, the annihilation channelν 1ν1 → A 1 A 1 opens up in the early Universe for the samples close from left to the blue line, and it soon becomes the dominant one with the increase of mν 1 up to about 100GeV. Since this annihilation is a s-wave dominant process, the dSph constraint is rather strong to exclude a large portion of the samples. While on the other hand, there still exist various ways to escape the constraint as we introduced in last section. Different from the left panel in Fig.2 where only the constraints listed in the text are considered, the right panel further considers the constraint Y ν v u /(λ N v s ) < 0.1 on the samples. This constraint is motivated by the limitations on the non-unitary of neutrino mixing matrix [54] and the electroweak precision data [57]. We remind that the masses of A 1 , h 1 and h 2 are slightly altered by the radiative correction from the sneutrino sector, and the positions of the vertical lines only act as a rough indication of the masses.
Next we consider the SI DM-nucleon scattering rate, which is the focus of this work. In Fig.3, we project the samples of Fig. 2 on σ SĨ ν 1 −p − mν 1 plane with the same color convention as that of Fig. 2. As expected in last section, the constraint from the XENON-1T experiment is rather weak on the sneutrino DM in the ISS-NMSSM, and only a small portion of the samples are excluded. Especially if we further require Y ν v u /(λ N v s ) < 0.1, only few samples are excluded. We emphasize that in the ISS-NMSSM, the SI cross section can be lower than the neutrino background even for light higgsinos, and consequently the DM may never be probed in DD experiments.
In the following, we present more information about the SI cross section. In Fig.4, we only consider the samples in the left panel of Fig. 3, and project them on Cν 1ν1 h 1 − Cν 1ν1 h 2 plane (left panel) and a d1 − a d2 plane (right panel) respectively, where the green samples are excluded by the XENON-1T experiment, and the dark blue ones are not. The left panel indicates that |Cν 1ν1 h 1 | 4GeV, and Cν 1ν1 h 2 varies in a much wider range from −20GeV to 50GeV for the surviving samples. It also indicates that the couplings Cν 1ν1 h 1 and Cν 1ν1 h 2 seem to be roughly linear dependent for most samples. The underlying reason for the correlation is that Cν 1ν1 h 2 Cν 1ν1 s and Cν 1ν1 h 1 Cν 1ν1 Hu + S 13 Cν 1ν1 s S 13 Cν 1ν1 s where we used the fact |Cν 1ν1 s | |Cν 1ν1 Hu | (for similar discussions, see Eq. 2.33). The right panel shows that −2GeV −3 a d1 1GeV −3 and −3GeV −3 a d2 1.5GeV −3 for the surviving samples, and a similar correlation between a d1 and a d2 exists for most samples. About Fig.4 three points should be noted. First, for the typical setting of the NMSSM parameters in Table 2, the coefficients a i obey the relations: a u1 a d1 , |a d2 | is several times larger than |a u2 | due to the large tan β, and |a d1 | |a u3 |, |a d3 |. As for a d1 and a d2 , their magnitudes may be comparable, and they can interfere constructively or destructively in contributing to the cross section. Second, since a d2 ∝ Cν 1ν1 h 2 /m 2 h 2 , the  range of Cν 1ν1 h 2 must be narrowed correspondingly to survive the XENON-1T constraints if we choose a lighter h 2 . In this case, more parameter space of the ISS-NMSSM will be limited by the DD experiments. Third, in case of Cν 1ν1 Hu 0 where the correlation holds, only the interactions ofν 1 with the singlet Higgs field are significant. These interactions alone can be responsible for the right relic density, and meanwhile contribute to the cross section. This cross section, however, is usually lower than the bound of the XENON-1T experiment, and is thus experimentally favored.
Finally, we consider the dependence of the cross section on the parameters in sneutrino sector. In Fig.5, we project the samples in Fig.4 on T ν − Y ν plane, where the colors correspond to the values of a d1 and the circled samples are excluded by the XENON-1T experiment. This figure indicates that the sample with a large Y ν and/or a large T ν tends to predict large |a d1 | and σν 1 −p , and the parameter region preferred by the experiment is Y ν 0.15 and |T ν | 100GeV. This fact can be understood from Eq.(2.30) by noting that the h 2 contribution is insensitive to the two parameters. We note that for the samples along the black line direction, their predictions on |a d1 | are usually small even for large Y ν and T ν . We checked that it is due to the cancelation between the Y ν contribution and the T ν contribution. We also note that there exist samples which correspond to small Y ν and T ν , but are excluded by the XENON-1T experiment. We checked that these samples correspond to a quite large |a d2 | with a d2 /a d1 > 0.
In Fig.6, we project the samples in Fig.4 on a d 2 − λ N plane with the colors indicating the values of T λ N . The left panel and the right panel correspond to Z 12 Z 13 > 0 case and Z 12 Z 13 < 0 case respectively. This figure indicates that for Z 12 Z 13 > 0 case, T λ N prefers to be negative, while for Z 12 Z 13 < 0 case, it tends to positive. In any case, the effect of T λ N is to cancel the λ N contribution to a d2 . This can be understood by the formula where we used the approximation Cν 1ν1 h 2 Cν 1ν1 s and Eq.(2.29). We remind that it is due to the cancelation, λ N as large as 0.3 is still allowed by the XENON-1T experiment. We also remind that the allowed values of λ N and T λ N by the XENON-1T experiment depend on our choice of m h 2 .

LHC constraints on the model
In this section, we examine the constraints from the direct searches for electroweakinos at the LHC on the samples considered in last section. Since Br(χ 0 1,2 →ν 1 ν τ ) = Br(χ 0 1,2 → ν R 1 ν τ ) 50%, Br(χ ± 1 →ν 1 τ ± ) = Br(χ ± 1 →ν R 1 τ ± ) 50% for the samples andν R 1 is longlived at colliders due to its nearly degeneracy withν 1 in mass, we consider the Mono-jet signal from the processes pp →χ ±  [118][119][120] and 13-TeV LHC by ATLAS collaboration [121], all of which have been encoded in the package CheckMATE [122][123][124]. The common requirements of the analyses are: (1) an energetic jet with p T > 100GeV and possible existence of one additional softer jet with ∆φ(j 1 , j 2 ) < 2.5 to suppress large QCD dijet background; (2) large missing energy, typically E miss T > 150GeV; (3) vetoing any event with isolated leptons. With regard to the signal of two hadronic τ s plus E miss T , the strongest limit comes from the analyses of the direct Chargino/Neutralino production at the 8-TeV LHC by   Table 3. Detailed information about the analysis of 2τ + E miss T signal in [125] for six parameter points. SR * stands for the SR with the largest expected sensitivity, and is the net cut efficiency of the signal events.
ATLAS and CMS collaborations [125,126] and the 13-TeV LHC by ATLAS collaboration [127]. As far as our case (i.e. fixed µ at 120GeV) is concerned, the analysis in [125] imposes stronger constraint than that in [127], which can be learned from figure 7 in [127]. The underlying reason is that the analysis in [127] focuses on heavy Chargino case, which requires more energetic jets and larger missing energy than the former. Moreover, we note that the constraint of the analysis in [126] is similar to that in [125] for mχ± 1 < 200GeV, which can be learned by comparing figure 5 in [126] with figure 7 in [125]. So in this work, we only consider the analysis in [125] on the 2τ + E miss T signal. We implement this analysis in the package CheckMATE with the corresponding validation presented in appendix.
To study these signals, we first use the package SARAH [108][109][110] to generate the model files of the ISS-NMSSM in UFO format [128]. Then we use the simulation tools MadGraph/MadEvent [129,130] to generate the parton level events of the processes with Pythia6 [131] for parton fragmentation and hadronization, and Delphes [132] for the fast simulation of the ATLAS or CMS detector. Finally we use the improved CheckMATE to implement the cut selections of the analyses.
For each of the analyses, we consider the signal region (SR) with the largest expected sensitivity for a given mν 1 6 , and calculate its R value defined by R ≡ S/S OBS 95 , where S stands for the number of signal events in the SR with the statistical uncertainty considered and S OBS 95 denotes the observed limit at 95% confidence level for the SR. For the signal which corresponds to several experimental analyses, we select the largest R among the analyses, denoted by R max hereafter, to parameterize the capability of the LHC in exploring the parameter point. If R max is larger than unity, the point is excluded and otherwise it is allowed. In Fig.7, we present our results of R max for the Mono-jet signal (green line) and the 2τ + E T miss signal (red line) respectively. This figure indicates that the constraints from the Mono-jet signal is very weak, and R max reaches only 0.01 in best case. The underlying reason is that the experimental analyses require a relatively large E miss T , which can not be satisfied for most of the events. By contrast, R max for the 2τ + E T miss signal increases monotonously with the enlarged mass splitting betweenχ ± 1 andν 1 , and for mν 1 55GeV it exceeds unity. In order to understand the features of the analysis on the 2τ + E T miss signal, we choose six parameter points and provide in Table 3 more information about the analysis in [125]. As can be seen from this table, the cut efficiency is quite large for mν 1 30GeV, reaching about 19%, and it drops quickly with the increase of mν 1 to 0.6% for mν 1 84GeV.

Phenomenology of the ISS-NMSSM
In the ISS-NMSSM, the DM candidate may be the lightest sneutrino or the lightest neutralino. In this section, we only briefly sum up the phenomenology of the former case. Some of our viewpoints may be applied to the latter case, which will enrich the well studied phenomenology of the NMSSM.
In the ISS-NMSSM, the impact of the sneutrino DM on the phenomenology is reflected in following aspects: • Relaxing greatly the parameter space of the NMSSM, and meanwhile maintaining the naturalness of the model. As we pointed out in [17], so far the DD experiments have put very strong constraints on the natural NMSSM, and consequently it is not easy to get parameter points coinciding with the constraints from DD experiments. In the ISS-NMSSM, however,ν 1 can serve as a viable DM candidate if Min(m h 1 , m A 1 ) < 2mχ0 1 or if the lightest higgsino corresponds toχ 0 1 , which have been illustrated before. These conditions can be easily satisfied in the ISS-NMSSM, and consequently new features in comparison with the original NMSSM may appear in Higgs physics as well as in sparticle physics. Taking the parameter point in Table 2 as an example, we found that it can not get the proper relic density and meanwhile predict an unacceptable large SD cross section for DM-nucleon scattering in the framework of the NMSSM [134]. In the ISS-NMSSM, however, it becomes a phenomenologically viable point.
• Existence of relatively light particles, such as non-standard neutrinos and light Higgs bosons, which, beside exhibiting themselves at colliders by exotic signals, may serve as the decay products of the Higgs bosons and sparticles. This feature makes the search for new particles at colliders quite complicated. For examples, we find from the samples that the non-standard neutrinos may be as light as 30GeV. In this case, a left-handed slepton may decay dominantly into one of the neutrinos plus a higgsinodominated neutralino by the neutrino Yukawa interaction. As a result, the signature of the slepton is distinct from that in the NMSSM. Moreover, since the neutrino has a small left-handed neutrino component, it may be produced in association with one active neutrino at the LHC [90], or with one lepton [135], or in pairs [136][137][138].
Obviously, how to detect these signals is an open question.
• Existence of new interactions which alter the properties of the particles in the NMSSM, and also induce new contribution to some observables. For example, the neutrino Yukawa interaction in the ISS-NMSSM can not only change the decay modes of the higgsinos and the left-handed sleptons in the NMSSM, but also contribute to Higgs boson masses [54,57,67], lepton flavor violating processes [45,47,48,50,58,60,62] as well as muon anomalous magnetic momentum [72,74].
Due to these aspects, the phenomenology of the ISS-NMSSM is quite rich, and may be different from that of the NMSSM. As far as sparticles are concerned, the speciality of the ISS-NMSSM comes from the fact that the couplings ofν 1 with the other particles are usually suppressed, and meanwhile it carries a certain lepton flavor number if Y ν is flavor diagonal. As a result, heavy sparticles will not decay directly intoν 1 , but instead they first decay into a relatively light sparticle with stronger couplings [60,139,140]. This lengthened decay chain makes the decay products of the parent sparticle quite model dependent. For example, if the sleptons are lighter than the higgsinos, the signature of the higgsinos usually corresponds to multi-lepton final state [141,142], which is different from the final states discussed in last section. We remind that in principle Y ν may be flavor non-diagonal, and consequentlyν 1 will not have a definite lepton flavor number any more. This further complicates the sparticle decays.

Conclusions
Given the increasing tension between naturalness and the DD experiments for customary neutralino DM candidate in supersymmetric theories, we discuss the feasibility that the lightest sneutrino acts as a DM candidate to alleviate the tension. For this end, we assume certain symmetries, and extend the field content of the NMSSM in an economical way to incorporate the inverse seesaw mechanism into the framework for neutrino mass. We point out that the resulting theory called ISS-NMSSM not only inherits all the merits of the NMSSM and the seesaw mechanism, but also exhibits new features in both DM physics and sparticle phenomenology. Especially by choosing the sneutrino as DM candidate, we find by analytic formulae that the DM-nucleon scattering rate is usually suppressed in comparison with the neutralino DM in the MSSM, and consequently the constraints from the DD experiments are no longer strong. We also find that the interactions of the sneutrino with the singlet Higgs field alone can account for the measured relic density, and meanwhile predict acceptable cross sections for both direct and indirect DM search experiments. We show these features numerically in physical parameter space, which is obtained by fixing the parameters in the NMSSM sector and scanning the parameters in the sneutrino sector with various experimental constraints (including the LHC search for 2τ + E miss T and Mono-jet+E miss T signals) considered. Finally, we also briefly discuss the phenomenology of the ISS-NMSSM, and point out that it is quite rich and distinct from that of the NMSSM. Given that the LHC experiments have not probed any signals of sparticles, the ISS-NMSSM may deserve a comprehensive study in near future.
Before we end this work, we'd like to compare briefly the ISS-NMSSM with the Type-I seesaw extension of the NMSSM proposed in [44]. In the Type-I seesaw extension, only right-handed neutrino fields are introduced to generate neutrino mass, and the corresponding neutrino Yukawa couplings are of O(10 −6 ), which is at same order as the electron Yukawa coupling in the SM, given that the masses for the right-handed neutrinos are about 1TeV. In both models, the singlet Higgs field plays an important role in various aspects, including generating the higgsino mass and the heavy neutrino masses dynamically, mediating the transition betweenν 1 pair and higgsino pair to keep them in thermal bath in early Universe, acting as DM annihilation final state or mediating DM annihilations, as well as affecting DM-nucleon scattering rate. Consequently both models can yield in certain parameter space thermal DM and a sneutrino-nucleon scattering cross section compatible with DD limits of the recent XENON-1T experiment. On the other hand, the essential difference of the two models comes from following two aspects. One is that in order to accommodate the experimental data for the neutrino oscillations, the electroweak precision measurements and the lepton-flavor violations, one can choose in the ISS-NMSSM a flavor-blind neutrino Yukawa couplings by encoding all the flavor structures into the small lepton-number violating parameter µ X as indicated in Eq.(2.12). This will make the non-unitary limitation mentioned in Eq.(2.8) easily satisfied. By contrast, there is no such freedom in the type-I seesaw extension, and one has to rely on the neutrino Yukawa couplings to account for all the experimental data. So we conclude that the ISS-NMSSM provides more theoretical flexibility in accommodating the data and at same time much richer phenomenology at colliders [138]. The other different comes from the signature of the heavy neutrinos [143]. In the Type-I seesaw extension, due to the Majorana nature of the heavy neutrinos, its associated production with one lepton at the LHC usually results in same-sign di-lepton signal, while in the ISS-NMSSM due to the pseudo-Dirac nature of the neutrinos, the process usually leads to tri-lepton signals. A more dedicated comparison of the two models will be carried out in our forthcoming work.

Acknowledgement
We thank Prof. Xiaojun Bi and Yufeng Zhou for helpful discussion about dark matter indirect detection experiments. This work is supported by the National Natural Science Foundation of China (NNSFC) under grant No. 11575053.

A Appendix
In this section, we validate our code for all SRs in [125]. We work in the MSSM, and consider four cases which correspond toχ ± 1χ ∓ 1 andχ ± 1χ 0 2 productions with bothχ ± 1 andχ 0 2 being wino-dominated,χ ± 1χ ∓ 1 production withχ ± 1 being wino-dominated,τ RτR production andτ LτL production respectively. For each validation, we generate 10000 events in the way introduced in Section 4. Our results are presented in Table 4, 5, 6 and 7 respectively. These tables indicate that we can reproduce the ATLAS results for case 1-3 at 20% level, and case 4 at 30%.  Table 4. Validation of theχ ± 1χ ∓ 1 andχ ± 1χ 0 2 production processes at the 8-TeV LHC by assuming mχ± 1 = mχ0 2 and mτ = mν = (mχ± 1 + mχ0 2 )/2. R ATLAS in the table is the result obtained by ATLAS collaboration, which is taken from the exclusion line of Fig.7a in [125]. SR * and R have same meanings as those in Table 3, and Diff ≡ (R − R ATLAS )/R ATLAS , which parameterizes the deviation of our calculation from its corresponding ATLAS result.  Table 6. Similar to Table 4, but for theτ RτR production process with the corresponding ATLAS results plotted in Fig.8a of [125].
(m τ L , mχ0  Table 7. Similar to Table 4, but for theτ LτL production process with the corresponding ATLAS results plotted in Fig.8b of [125].