A Bayesian analysis of sneutrino DM in the NMSSM with Type-I seesaw mechanism

In the Next-to-Minimal Supersymmetric Standard Model (NMSSM) with extra heavy neutrino superfields, neutrino can acquire its mass via a seesaw mechanism and sneutrino may act as a viable dark matter (DM) candidate. Given the strong tension between the naturalness for $Z$ boson mass and the DM direct detection experiments for customary neutralino DM candidate, we augment the NMSSM with Type-I seesaw mechanism, which is the simplest extension of the theory to predict neutrino mass, and study the scenarios of sneutrino DM. We construct likelihood function with B-physics measurements, LHC Higgs data, DM relic density and its direct and indirect search limits, and perform a comprehensive scan over the parameter space of the theory by Nested Sampling method. We adopt both Bayesian and frequentist statistical quantities to illustrate the favored parameter space of the scenarios, the DM annihilation mechanism as well as the features of DM-nucleon scattering. We find that the scenarios are viable over broad parameter regions,and the DM usually co-annihilated with the Higgsinos to get the measured relic density. Interestingly, our results indicate that a rather low Higgsino mass, $\mu \lesssim 250 {\rm GeV}$, is preferred in the scan, and the DM-nucleon scattering rate is naturally suppressed to coincide with the recent XENON-1T results. Other issues, such as the LHC search for the Higgsinos, are also addressed.


I. INTRODUCTION
A large number of cosmological and astrophysical observations have firmly established the existence of nonbaryonic dark matter (DM) [1][2][3][4]. Among the possible candidates, the weakly interactive massive particles (WIMPs) are most attractive since they naturally lead to right DM abundance [5,6]. So far this type of DM candidates is still compatible with the more and more stringent constraints from DM direct and indirect search experiments, which was recently emphasized in [7]. As the most popular ultraviolet-complete beyond Standard Model (BSM), the minimal supersymmetric Standard Model (MSSM) predicts two WIMP-like DM candidates in the form of sneutrino [8] or neutralino [9] when R-parity is imposed. For the left-handed sneutrino as the lightest supersymmetric particle (LSP), its interaction with the Z boson predicts a very small relic abundance relative to its measured value as well as an unacceptably large DMnucleon scattering rate, which was firstly noticed by T. Falk et al. [10] after considering the Heidelberg-Moscow DM direct detection (DD) experiment [11] and later emphasized in [12]. Consequently, the neutralino DM as the solely viable WIMP candidate has been intensively studied over the past decades. However, with the rapid progress in DM DD experiments (such as PandaX-II [13,14], LUX [15], and XENON-1T [16,17]) in recent years, it was found that the candidate became disfavored by the experiments [18][19][20] assuming that it is fully responsible for the measured DM relic density and that the Higgsino mass μ is of Oð10 2 GeVÞ, which is favored by Z boson mass. This situation motivates us to consider DM physics in extended MSSM.
Besides the DM puzzle, the non-vanishing neutrino mass is another firm evidence on the existence of new physics [21]. Seesaw mechanism is the most popular way to generate the mass, and depending on the introduction of heavy neutrino fields, several variants of this mechanism, such as type-I [22], -II [23], -III [24] and inverse [25,26] seesaw, have been proposed. Among these variants, the type-I mechanism is the most economical one where only a right-handed neutrino field is introduced. In the simplest supersymmetric realization of this mechanism, namely the MSSM with type-I mechanism, a pure right-handed sneutrino LSP [27][28][29] or a mixed left-and right-handed sneutrino LSP [12,30,31] may act as a viable DM candidate. For the former case, the coupling of the candidate to ordinary matter is extremely suppressed either by neutrino Yukawa couplings or by the mass scale of the right-handed neutrino. As a result, its self-and coannihilation cross sections are so tiny that it has to be nonthermal to avoid an overclosed universe [27,28]. For the latter case, a significant chiral mixture of the sneutrinos requires an unconventional supersymmetry breaking mechanism [30]. Furthermore, since the couplings of the DM candidate with the SM particles are determined by the mixing, it is difficult to predict simultaneously the right DM abundance and a suppressed DM-nucleon scattering rate required by the DM DD experiments [32][33][34][35]. 1 These facts reveal that the DM physics in type-I MSSM is also unsatisfactory. This situation is also applied to the inverse seesaw extension of the MSSM [37].
The situation may be changed greatly if one embeds the seesaw mechanism into the NMSSM, which, as one of most economical extensions of the MSSM, is characterized by predicting one gauge singlet Higgs superfieldŜ [38]. It has long been known that the fieldŜ plays an extraordinary role in the model: solving the μ problem of the MSSM [38], enhancing the theoretical prediction about the mass of the SM-like Higgs boson [39][40][41] as well as enriching the phenomenology of the NMSSM (see for example [42][43][44][45][46][47]). In this context we stress that, in the seesaw extension of the NMSSM, it is also responsible for heavy neutrino mass and the annihilation of sneutrino DM [48][49][50], and consequently makes the sneutrino DM compatible with various measurements. The underlying reason for the capability is that the newly introduced heavy neutrino fields in the extension are singlets under the gauge group of the SM model, so they can couple directly withŜ. We also stress that the seesaw extension is essential not only to generate neutrino mass, but also to enrich greatly the phenomenology of the NMSSM given the very strong constraint of the recent DM DD experiments on neutralino DM in the NMSSM [51].
In our previous work [48], we augmented the NMSSM with inverse seesaw mechanism by introducing two types of gauge singlet chiral superfieldsν R andX, which have lepton number −1 and 1 respectively, and sketched the features of sneutrino DM. We found that, due to the assignment of the superfield's charge under the SM gauge group, the scalar component fields ofν R ,X andŜ compose a secluded DM sector, which can account for the measured DM relic abundance, and also be testable by future DM indirect detect experiments and LHC experiments. Since this sector communicates with the SM sector mainly through the small singlet-doublet Higgs mixing, the DM-nucleon scattering rate is naturally suppressed, which is consistent with current DM direct search results. We note that these features should be applied to the type-I seesaw extension of the NMSSM due to the similarities of the two theoretical frameworks. 2 We also note that, in comparison with the inverse seesaw extension, the sneutrino sector in the type-I seesaw extension involves less parameters, and thus is easier to be fully explored in practice. So as a preliminary work in our series of studies on sneutrino DM in different seesaw extensions of the NMSSM, we focus on the type-I extension, and perform a rather sophisticated scan over the vast parameter space of the model by the nested sampling method [52]. The relevant likelihood function is constructed from LHC Higgs data, B-physics measurements, DM relic density and its direct and indirect search limits, and statistic quantities are used to analyze the scan results. As far as we know, such an analysis of the type-I seesaw extended NMSSM has not been done before, and new insights about the model are obtained. For example, we find that for most samples in the 1σ regions of the posterior probability distribution function (PDF) on the plane of the sneutrino DM mass m˜ν 1 versus the Higgsino mass μ, the two masses are nearly degenerate. This is because the singlet Higgs field can mediate the transition between the DM pair and the Higgsino pair, which implies that the DM 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 the DM during freeze-out, and consequently the Higgsinos played an important role in determining DM relic density [53] (in literature such a phenomenon was called coannihilation [54]). As a result, even for very weak couplings of the DM with SM particles, the DM may still reach the correct relic density by coannihilating with the Higgsino-dominated particles. Such a possibility is not discussed in the 1 Numerically speaking, the mixing angle should satisfy sin θ˜ν ∼ 0.02 for m˜ν ¼ 100 GeV to predict the right DM relic density by the Z boson mediated annihilation, which corresponds to the scattering rate at the order of 10 −45 cm −2 [34]. Such a rate has been excluded by the latest XENON-1T experiment, which, on the other side, limits sin θ˜ν < 0.01 by the recent calculation in [35]. Moreover, we find that the correlation between the relic density and the scattering rate is underestimated in Fig. 1 of [36]. 2 Generally speaking, the property of the sneutrino DM in the type-I extension differs from that in the inverse seesaw extension only in two aspects. One is that in the type-I extension, the sneutrino DM is a roughly pure right-handed sneutrino state, while in the inverse seesaw extension, it is a mixture ofν L ,ν R and X with the last two components being dominated. Consequently, the DM physics in the inverse seesaw extension is more complex [48][49][50]. The other is that in both the extensions, CP-even and CP-odd sneutrino states split in mass due to the presence of lepton number violating interactions. In the type-I extension the splitting may be quite large, while in the other model, it tends to be smaller than about 1 GeV. literatures. We also find that in many cases, the model does not need fine tuning to be consistent with the DM DD limits on the DM-nucleon scattering rate. This is a great advantage of the theory in light of the tightness of the limits, but it is not emphasized in the literatures.
The type-I extension of the NMSSM was firstly proposed in [55], and its DM physics was sketched in [56,57]. Since then a lot of works appeared to study the phenomenology of the model [58][59][60][61][62][63][64][65][66][67]. For example, the spectral features of the γ-ray from DM annihilations were investigated in [58,60,61,64,65], and the Higgs physics was discussed in [59,62]. We note that some of these studies focused on the parameter regions which predict a relatively large DM-nucleon scattering rate [56,57]. These regions have now been excluded by the DM DD experiments, and thus the corresponding results are out of date. We also note that some other works were based on multicomponent DM assumption, so they used the upper limit of the DM relic abundance as the criterion for parameter selection [58,60,61,64]. Obviously, the conclusions obtained in this way are less definite than those with single DM candidate assumption. Given the incompleteness of the research in this field and also the great improvements of experimental limits on DM property in recent years, we are encouraged to carry out a comprehensive study on the key features of the model in this work. This paper is organized as follows. In Sec. II, we introduce briefly the basics of the NMSSM with type-I seesaw extension and the property of the sneutrino DM. In Sec. III, we perform a comprehensive scan over the vast parameter space of the model by considering various experimental measurements, and adopt statistic quantities to show the favored parameter space, DM annihilation mechanisms as well as the features of DM-nucleon scattering for the cases that the sneutrino DM is CP-even and CPodd respectively. Collider constraints on sneutrino DM scenarios are discussed in Sec. IV, and conclusions are presented in Sec. V.

II. NMSSM WITH TYPE-I MECHANISM
In this section we first recapitulate the basics of the NMSSM with type-I mechanism, including its Lagrangian and Higgs sector, then we concentrate on sneutrino sector by analyzing sneutrino mass matrix, the annihilation channels of sneutrino DM, and its scattering with nucleon.
A. The Lagrangian of the model As the simplest seesaw extension of the NMSSM, the NMSSM with type-I seesaw mechanism introduces three generations of right-handed neutrino fields to generate neutrino mass. Consequently, the extension differs from the NMSSM only in neutrino/sneutrino sector, and the sneutrino DM as the lightest supersymmetric state in this sector couples mainly to Higgs bosons. With the field content presented in Table I, the relevant superpotential and soft breaking terms are [56,57] where W F denotes the superpotential of the MSSM without μ term, and a Z 3 symmetry is introduced to forbid the appearance of any dimensional parameters in W. In the above formulas, the coefficients λ and κ parametrize the interactions among the Higgs fields, Y ν andλ ν are neutrino Yukawa couplings with flavor index omitted, m i (i ¼ H u ; H d ; Á Á Á) denote soft breaking masses, and A i (i ¼ λ; κ; Á Á Á) are soft breaking coefficients for trilinear terms.

B. Higgs sector
Since the soft breaking squared masses m 2 H u , m 2 H d , and m 2 S are related with the vacuum expectation values of the fields H u , H d , and S, , by the minimization conditions of the Higgs potential after the electroweak symmetry breaking [38], it is conventional to take λ, κ, tan β ≡ v u =v d , A λ , A κ , and μ ≡ λv s = ffiffi ffi 2 p as theoretical input parameters in Higgs sector. The elements of the squared mass matrix for CP-even Higgs fields in the basis (S 1 ≡ cos βRe½H 0 u − sin βRe½H 0 d , S 2 ≡ sin βRe½H 0 u þ cos βRe½H 0 d , S 3 ≡ Re½S) are then given by [41] As a result, the model predicts three CP-even Higgs mass eigenstates h i with i ¼ 1, 2, 3 and two CP-odd mass eigenstates A 1 and A 2 , which are the mixtures of the real and imaginary parts of the fields H 0 u , H 0 d , and S, respectively. Throughout this paper, we label these eigenstates in an ascending mass order, i.e., m h 1 < m h 2 < m h 3 and So far the measurement on the property of the 125 GeV Higgs boson at the LHC indicates that it is quite SM-like, which restricts the mixing of the S 2 field with the other fields to be small. This implies from the definition of the S 1 and S 2 fields that the SM-like Higgs boson is Re½H 0 u dominated if tan β ≫ 1, and the heavy neutral doublet dominated mass eigenstates are mainly composed by H 0 d field. In Sec. III and Appendices, we will discuss in detail the implication of the Higgs physics on the input parameters. We remind the reader that it is also popular to adopt the basis (Re½H 0 d , Re½H 0 u , Re½S) in studying the property of the CP-even Higgs bosons since the couplings of the basis with SM fermions are simple.
The model also predicts a pair of charged Higgs H AE ,  [68]. These doubletlike states couple with the sneutrino DM via the interaction λλ ν H u · H dν Ã Rν Ã R þ h:c:, which is induced by the F-term of the Lagrangian. In the limit m H AE → ∞, only the SM-like Higgs boson plays a role in DM physics, but since its interaction with sneutrino pair is proportional to λλ ν v d , its effect is usually unimportant.
As for the singlet-dominated Higgs bosons, collider constraints on them are rather weak and consequently they may be light. These states couple with sneutrino pair with three or four scalar interaction induced by thē λ νŝνν term in the superpotential and its soft breaking term. Consequently, they can act as the annihilation product of the sneutrino DM or mediate the annihilation, and thus play an important role in DM annihilation.

C. Sneutrino sector
In the NMSSM with type-I seesaw extension, the active neutrino mass matrix is given by ν v S denoting the heavy neutrino mass matrix [55]. Since the mass scale of the active neutrinos is ∼0.1 eV, the magnitude of Y ν should be about 10 −6 for the scale of M around 100 GeV. In order to reproduce neutrino oscillation data, m ν must be flavor nondiagonal, which can be realized by assuming that the Yukawa coupling Y ν is nondiagonal, whileλ ν is diagonal. If one further assumes that the soft breaking parameters in sneutrino sector, such as m˜l,mν andĀ λ ν , are flavor diagonal, the flavor mixings of sneutrinos are then extremely suppressed by the off-diagonal elements of Y ν . In this case, it is enough to only consider one generation case in studying the properties of the sneutrino DM, which is what we will do. In the following, we use the symbols λ ν , A λ ν , and mν to denote the 33 element ofλ ν ,Ā λ ν , andmν respectively.
After decomposing sneutrino fields into CP-even and CP-odd parts: one can write down the sneutrino mass matrix in the basis (ν L1 ,ν R1 ,ν L2 ,ν R2 ) as follows (2.5) If all the parameters in the matrix are real, namely there is no CP violation, the real and imaginary parts of the sneutrino fields will not mix and the mass term can be split into two parts for CP-even (CP-odd) states, and the minus signs in the matrix are for the CP-odd states. From this formula, one can learn that the chiral mixings of the sneutrinos are proportional to Y ν , and hence can be ignored safely. So sneutrino mass eigenstate coincides with chiral state. In our study, the sneutrino DM corresponds to the lightest right-handed sneutrino. One can also learn that the mass splitting between the CP-even and CP-odd righthanded states is given by Δm 2 ≡ m 2 even − m 2 odd ¼ 4m 2 RR , which implies that the CP-even stateν R1 is lighter than the CP-odd stateν R2 if m 2 RR < 0 and vice versa. This implies that the sneutrino DM may be either CP-even or CP-odd. In this work, we consider both the possibilities.
Once the form of the sneutrino DMν 1 is given, one can determine its coupling strength with Higgs bosons. For example, for a CP-evenν 1 we have where Z ij with i, j ¼ 1, 2, 3 (Z 0 mn with m, n ¼ 1, 2) are the elements of the matrix to diagonalize the CP-even Higgs mass matrix in the basis (Re½H 0 d , Re½H 0 u , Re½S) [the CPodd Higgs mass matrix in Eq. (2.3)], and for a CP-oddν 1 , its coupling strengths can be obtained from the expressions by the substitution λ ν → −λ ν . As is expected, the first coupling is suppressed by a factor λλ ν cos β if h i as the SMlike Higgs boson is ReðH 0 u Þ dominant, and all the couplings may be moderately large if the Higgs bosons are singlet dominant.

D. DM relic density
It is well known that the WIMP's relic abundance is related to its thermal averaged annihilation cross section at the time of freeze-out [9]. In order to obtain the WIMP's abundance one should solve the following Boltzmann equation where g Ã is the effective number of degrees of freedom at thermal equilibrium, M p is the Plank mass, Y and Y eq are the relic abundance and the thermal equilibrium abundance respectively, and hσvi is WIMP's relativistic thermal averaged annihilation cross section with v denoting the relative velocity between the annihilating particles. hσvi is related to the particle physics model by [69] hσvi ¼ where g i is the number of degrees of freedom, σ ij;kl is the cross section for annihilation of a pair of particles with masses m i and m j into SM particles k and l, p ij is the momentum of incoming particles in their center of mass frame with squared total energy s, and K i (i ¼ 1, 2) are modified Bessel functions. The present day abundance is obtained by integrating Eq. (2.9) from T ¼ ∞ to T ¼ T 0 , where T 0 is the temperature of the Universe today. 3 And WIMP relic density can be written as [69] Ωh 2 ¼ 2.742 × 10 8 M WIMP GeV YðT 0 Þ: ð2:11Þ In the NMSSM with type-I seesaw mechanism, possible annihilation channels of the sneutrino DM include [56,57] (1)ν 1ν1 → VV Ã , VS, ff with V, S and f denoting a vector boson (W or Z), a Higgs boson and a SM fermion, respectively. This kind of annihilations proceed via s-channel exchange of a CP-even Higgs boson. (2)ν 1ν1 → SS Ã via s-channel Higgs exchange, t=uchannel sneutrino exchange, and relevant scalar quartic couplings. (3)ν 1ν1 → ν RνR via s-channel Higgs exchange and t=uchannel neutralino exchange.
denoting a right-handed sneutrino with an opposite CP number to that ofν 1 , and X ð0Þ and Y ð0Þ denoting any possible light state. These annihilation channels are important in determining the relic density only when the CP-even and CP-odd states are nearly degenerate in mass. (5)ν 1H → XY andHH 0 → X 0 Y 0 withH andH 0 being any Higgsino dominated neutralino or Higgsino dominated chargino. These annihilation channels are called coannihilation [53,54], and they become important if the mass splitting between Higgsino and ν 1 is less than about 10%. The expressions of σv for some of the channels are presented in [57]. One can learn from them that the parameters in sneutrino sector, such as λ ν , A ν , and m 2 ν , as well as the parameters in Higgs sector are involved in the annihilations. 4

E. DM direct detection
Sinceν 1 in this work is a right-handed scalar with definite CP and lepton numbers, its scattering with nucleon N (N ¼ p, n) proceeds only by exchanging CP-even Higgs bosons. In the nonrelativistic limit, this process is described by an effective operator Lν 1 N ¼ f Nν1ν1ψ N ψ N with the coefficient f N given by [72] Consequently, the spindependent cross section for the scattering vanishes, and the spin independent cross section is given by [72] Þ is the reduced mass of the nucleon with mν 1 , and ξ i with i ¼ 1, 2, 3 are defined by to facilitate our analysis. Obviously, ξ i represents the h i contribution to the cross section. In our numerical calculation of the DM-proton scattering rate σ SĨ ν 1 −p , we use the default setting of the package MICROMEGAS [69,73,74] for the nucleon form factors, σ πN ¼ 34 MeV and σ 0 ¼ 42 MeV, and obtain F ðpÞ u ≃ 0.15 and F ðpÞ d ≃ 0.14. 5 In this case, Eq. (2.14) can be approximated by ð2:16Þ 3 Generally speaking, it is accurate enough to integrate the Boltzmann equation from T ¼ m DM to T ¼ T 0 in getting the DM abundance [70]. In the code MICROMEGAS, it actually starts the integration from the temperature T 1 defined by . This temperature is moderately higher than the freeze-out temperature T f ≃ m DM =25, which is defined by the equation YðT f Þ ¼ 2.5Y eq ðT f Þ [70]. 4 We note that the works [56,57] failed to consider the annihilation channels (4) and (5). Especially, the channel (5) is testified to be most important by our following study. We also note that in calculating σv of the channelν 1ν1 → A i A j , the authors neglected the contribution from the t=u-channel sneutrino exchange. 5 Note that different choices of the pion-nucleon sigma term σ πN and σ 0 can induce an uncertainty of Oð10%Þ on F From this approximation and also the expression of Cν 1ν1 h i in Eq. (2.7), one can get following important features about the scatting of the sneutrino DM with nucleon: (i) For each of the h i contributions, it depends not only on the parameters in Higgs sector, but also on the parameters in sneutrino sector such as λ ν and A λ ν . This feature lets the theory have a great degree of freedom to adjust the contribution size so that the severe cancellation among the contributions can be easily achieved. By contrast, in the NMSSM with neutralino as a DM candidate, the contribution depends on, beside the DM mass, only the parameters in Higgs sector, and consequently it is not easy to reach the blind spot for DM-nucleon scattering due to the tight constraints on the Higgs parameters from LHC experiments [79][80][81][82]. (ii) Each of the h i contributions can be suppressed. To be more specific, the ReðH 0 d Þ dominated Higgs is usually at TeV scale, so its contribution is suppressed by the squared mass; the ReðH 0 u Þ dominated scalar corresponds to the SM-like Higgs boson, and its coupling withν 1 can be suppressed by cos β and/or λ ν , or by the accidental cancellation among different terms in Cν 1ν1 h i . In most cases, the contribution from the singlet-dominated scalar is most important, but such a contribution is obviously suppressed by the doublet-singlet mixing in the scalar. We emphasize that these features make the theory compatible with the strong constraints from the DM DD experiments in broad parameter spaces. In order to parametrize the degree of the cancellation among the contributions, we define the fine tuning quantity Δ FT as Obviously, Δ FT ∼ 1 is the most ideal case for any DM theory in predicting an experimentally allowed σ SĨ ν 1 −p , and the larger Δ FT becomes, the more unnatural the theory is.

III. NUMERICAL RESULTS
In our study of the sneutrino DM scenarios, we utilize the package SARAH-4.11.0 [83][84][85] to build the model, the codes SPHENO-4.0.3 [86] and FLAVORKIT [87] to generate the particle spectrum and compute low energy flavor observables, respectively, and the package MICROMEGAS 4.3.4 [69,73,74] to calculate DM observables by assuming that the lightest sneutrino is the sole DM candidate in the universe. We also consider the bounds from the direct searches for Higgs bosons at LEP, Tevatron and LHC by the packages HIGGSBOUNDS-5.0.0 [88] and HIGGSSIGNAL-2.0.0 [89]. Note that in calculating the radiative correction to the mass spectrum of the Higgs bosons, the code SPHENO-4.0.3 only includes full one-and two-loop effects using a diagrammatic approach with vanishing external momenta [86]. This leaves an uncertainty less than about 3 GeV for the SM-like Higgs boson mass.

A. Scan strategy
The previous discussion indicates that only the parameters in Higgs and sneutrino sectors are involved in the DM physics. We perform a sophisticated scan over these parameters in following ranges: by setting A b ¼ A t and the other less important parameters in Table II. All the parameters are defined at the scale Q ¼ 1 TeV. Our interest in the parameter space is due to following considerations: (i) Since the masses of the heavy doublet-dominated Higgs bosons are usually at TeV scale, which is favored by direct searches for extra Higgs boson at the LHC, their effects on the DM physics are decoupled and thus become less important. We fix the parameter A λ , which is closely related with m H AE [38], at 2 TeV as an example to simplify the calculation given that the scan is very time-consuming for our computer clusters. Consequently, we find that the Higgs bosons are heavier than about 1 TeV for almost all samples. We will discuss the impact of the other choices of A λ on the results in Appendices of this work. (ii) It is well known that the radiative correction from top/stop and bottom/sbottom loops to the Higgs mass spectrum plays an important role for SUSY to coincide with relevant experimental measurements at the LHC. In our calculation, we include such an varying A t over a broad region, jA t j ≤ 5 TeV. We remind that within the region, the color and charge symmetries of the theory remain unbroken [90,91], and the LHC search for third generation squarks does not impose any constraints. (iii) The upper bounds of the parameters λ, κ, λ ν , and tan β coincide with the perturbativity of the theory up to Planck scale [38]. (iv) Because the Higgsino mass μ is directly related to Z boson mass, naturalness prefers μ ∼ Oð10 2 GeVÞ [38]. So we require 100 GeV ≤ μ ≤ 300 GeV in this work, where the lower bound comes from the LEP search for chargino and neutralinos, and the upper bound is imposed by hand. We will discuss the effect of a wider range of μ on the results in Appendix B. (v) Since the sneutrino DM must be lighter than the Higgsino, its soft breaking mass m˜ν is therefore upper bounded by about 300 GeV, and jA λ ν j ≲ 1 TeV is favored by naturalness. (vi) In order to get correct EWSB and meanwhile predict the masses of the singlet dominated scalars around 100 GeV, jA κ j cannot be excessively large from naturalness argument (see Appendix A). (vii) We require the other dimensional parameters, such as M i with i ¼ 1, 2, 3 and m˜l, sufficiently large so that their prediction on sparticle spectrum is consistent with the results of the direct search for sparticles at the LHC. In order to make the conclusions obtained in this work as complete as possible, we adopt the MultiNest algorithm introduced in [52] in the scan, which is implemented in our code EASYSCAN_HEP [92], with totally more than 2 × 10 8 physical samples computed. 6 The output of the scan includes the Bayesian evidence defined by where PðΘjMÞ is called prior probability density function (PDF) for the parameters Θ ¼ ðΘ 1 ; Θ 2 ; Á Á ÁÞ in the model M, and PðDjOðM; ΘÞÞ ≡ LðΘÞ is the likelihood function for the theoretical predictions on the observables O confronted with their experimentally measured values D.
Computationally, the evidence is an average likelihood, and it depends on the priors of the model's parameters. For different scenarios in one theory, the larger Z is, the more readily the corresponding scenario agrees with the data. The output also includes weighted and unweighted parameter samples which are subject to the posterior PDF PðΘjM; DÞ. This PDF is given by and it reflects the state of our knowledge about the parameters Θ given the experimental data D, or alternatively speaking, the updated prior PDF after considering the impact of the experimental data. This quantity may be sensitive to the shape of the prior, but the sensitivity can be counterbalanced by sufficient data and thus lost in certain conditions [94]. Obviously, one can infer from the distribution of the samples the underlying physics of the model.
In the scan, we take flat distribution for all the parameters in Eq. (3.1). 7 We will discuss the influence of different prior PDFs on our results in Appendix B and verify that the flat distribution can predict the right posterior distributions. The likelihood function we construct contains where each contribution on the right side of the equation is given as follows: (i) For the likelihood function of the Higgs searches at colliders, we set where χ 2 comes from the fit of the theoretical prediction on the property of the SM-like Higgs boson to relevant LHC data with its value calculated by the code HIGGSSIGNAL [89], and A 2 reflects whether the parameter point is allowed or excluded by the direct searches for extra Higgs at colliders. In practice, we take a total (theoretical and experimental) uncertainty of 3 GeV for the SM-like Higgs boson mass in calculating the χ 2 with the package HIGGSSIGNAL, and set by the output of the code HIGGSBOUNDS [88] either A 2 ¼ 0 for the case of experimentally allowed at 95% confidence level (C.L.) or A 2 ¼ 100 for the other cases.
Since the Bayesian evidence for the scenario where the lightest Higgs boson h 1 acts as the 6 To be more specific, we performed five independent scans for each case mentioned in the text. In each scan, we set the nlive parameter in MultiNest (which denotes the number of active or live points used to determine the isolikelihood contour in each iteration [52,93]) at 2000. 7 Dimensional parameters usually span wide ranges, and if one sets them flat distributed, the efficiency of the scan is not high in general. Due to this reason, special treatments of such kind of parameters were frequently adopted in literatures [95,96]. In Appendix B of this work, we will show that the prior setting of simple log distribution for A κ can overemphasize the low jA κ j region in the scan, and consequently deform the shape of posterior distributions. SM-like Higgs boson is much larger than that with the next lightest Higgs boson h 2 as the SM-like Higgs boson [40,41] (see the discussion at the end of this work), we focus on the scenario where h 1 corresponds to the Higgs boson discovered at the LHC, and present the results according to the CP property of the sneutrino DMν 1 .
(ii) For the second, third, and fourth contributions, i.e., the likelihood functions about the measurements of BrðB s → μ þ μ − Þ, BrðB s → X s γÞ and DM relic density Ων 1 , they are Gaussian distributed, i.e., where O th denotes the theoretical prediction of the observable O, O exp represents its experimental central value and σ is the total (including both theoretical and experimental) uncertainty. We take the experimental information about the B physics measurements from latest particle data book [97] and the relic density Ων 1 from [98]. (iii) For the likelihood function of DM DD experiments L DD , we take a Gaussian form with a mean value of zero [99]: where σ stands for DM-nucleon scattering rate, and δ σ is evaluated by δ 2 σ ¼ UL 2 σ =1.64 2 þ ð0.2σÞ 2 with UL σ denoting the upper limit of the latest XEN-ON1T results on σ at 90% C.L. [17], and 0.2σ parametrizing theoretical uncertainties. (iv) For the likelihood function of DM indirect search results from dwarf galaxies L ID , we use the data of Fermi-LAT collaboration [100], and adopt the likelihood function proposed in [101,102]. About the likelihood function in Eq. (3.3), we have three more explanations. The first is that we do not consider the constraint on line signal of γ-ray from Fermi-LAT data since it is rather weak [103]. The second is that we do not include the likelihood function of the LHC search for Higgsinos since the involved Monte Carlo simulation is time consuming, and meanwhile the most favored parameter region satisfy the constraint automatically. We will address this issue in Sec. IV. And the last point is that during the scan, one usually encounters the nonphysical situation where the squared mass of any scalar particle is negative orν 1 is not the LSP. In this case, we set the likelihood function to be sufficiently small, e.g, e −100 .
Given the posterior PDF and the likelihood function, one can obtain statistic quantities such as marginal posterior PDF and profile likelihood function (PL). The marginal posterior PDF for a given set of parameters (Θ A ; Θ B ; …) is defined by integrating the posterior PDF PðΘjM; DÞ in Eq. (3.2) over the rest model parameters. For example, the one-dimensional (1D) and two-dimensional (2D) marginal posterior PDFs are given by In practice, these PDFs are calculated by the sum of weighted samples in a chain with user-defined bins, and their densities as a function of Θ A and ðΘ A ; Θ B Þ respectively reflect the preference of the samples obtained in the scan. On the other hand, the frequentist PL is defined as the largest likelihood value in a certain parameter space. Take the 1D and 2D PLs as an example, we get them by the procedure Obviously, PL reflects the preference of a theory on the parameter space, and for a given point on Θ A − Θ B plane, the value of LðΘ A ; Θ B Þ represents the capability of the point in the theory to account for experimental data by varying the other parameters. In the following, we display our results by these quantities. We also use other statistic quantities, such as 1σ and 2σ credible regions (CRs) for the marginal PDF, 1σ and 2σ confidence intervals (CIs) for the PL, and posterior mean, median and mode, to illustrate the features of the sneutrino DM scenario. The definition of these quantities can be found in the Appendix of [104], and we use the package SUPERPLOT [104] with kernel density estimation to get them. Note that all the quantities depend on the parameter space in Eq. (3.1), which is inspired by the physics of the theory.

B. Favored parameter regions
In this subsection, we discuss the favored parameter space of the scan for the case that the sneutrino DM is a CP-even scalar. In Fig. 1, we show the 1D marginal PDFs and PLs for the parameters λ, κ, tan β, A κ , μ, A t , λ ν , A λ ν , and mν respectively. We also present 2D results in Fig. 2 on κ − λ, tan β − λ, and μ − λ planes with color bar representing marginal posterior PDF for the top panels and PL for the bottom panels. The 1σ and 2σ CRs, the 1σ and 2σ CIs together with the best fit point, posterior means, medians, and modes are also shown in these panels. From Fig. 1, one can learn following features: (i) The marginal PDF of λ is peaked around 0.2 and its 1σ CR corresponds to the region 0.15 ≲ λ ≲ 0.37. The underlying reason for this conclusion is that the parameter λ affects both the mass and the couplings of the SM-like Higgs boson with the other SM particles [41]. A large λ tends to enhance the singlet component in the SM-like Higgs boson (see the expression of M 2 23 in Eq. (2.2) of Sec. II and also the discussion in Appendix A), and is thus disfavored by the LHC Higgs data. 8 On the other hand, a small λ is also tightly limited since it can enhance theν 1ν1 h i coupling [see Eq. (2.7) in Sec. II], and consequently the rate of the DM-nucleon scattering.
FIG. 1. One-dimensional marginal posterior PDFs and PLs for the parameters λ, κ, tan β, A κ , μ, λ ν , A λ ν and m˜ν respectively. Credible regions, confidence intervals and other statistic quantities are also presented. 8 In fact, as one can learn from the results in [41], the parameter space of the NMSSM is rather limited for the large λ case to predict a SM-like Higgs boson with mass around 125 GeV.
In the second part of Appendix A, we show the impact of different experimental data on various marginal PDFs and PLs with the results presented from Figs. 14 and 15 and summarized in Table IV. The marginal PDF of λ in Fig. 14 reveals that it is mainly determined by the Higgs and DM observables, which verifies our conjectures. (ii) The posterior PDFs of κ and A κ are maximized around −0.6 and 200 GeV respectively, and the PDFs show an approximate reflection symmetry.
Since the underlying reason for the behavior is somewhat complicated, we will discuss it in an elaborated way in Appendixes A and B. (iii) A relatively large tan β, i.e., 10 ≲ tan β ≲ 38 for 1σ CR, is preferred by the marginal PDF. This is due to at least two facts. One is that, since λ is usually less than about 0.37 so that its contribution to the squared mass of the SM-like Higgs boson is sub-dominant, a large tan β is helpful to enhance the tree level MSSM prediction of the mass [38]. The other is that the heavy doublet dominated Higgs bosons in our model contribute to the process B s → μ þ μ − , and the effect can be enhanced by a factor of tan 6 β [105]. So a too large tan β is not favored by B physics. We note that these two facts are actually reflected in Fig. 14, where we investigate the effect of various observables on the shape of the PDF and PL for tan β.
(iv) A moderately large μ ranging from about 170 GeV to 280 GeV for 1σ CR is favored by the posterior PDF. The reason is twofold. One is that due to the specific choice of A λ and the third generation squark masses, Higgs physics and B physics observables are insensitive to this parameter, and consequently μ is evenly distributed in the range from about 100 GeV to 300 GeV (We will discuss this issue in detail in Appendix A, B and C). The other reason is that since the sneutrino DM tends to coannihilate with the Higgsinos to get its right relic density, which requires their mass splitting to be less than about 10% (see Fig. 3 below), a large μ provides a broader parameter space for the annihilation.
In Appendix C of this work, we will show that any modification of the scan range in Eq. (3.1) or the setting of A λ may alter the marginal PDF of μ, and in some cases moderately low values of μ are preferred.
(v) jA t j ≳ 2 TeV is preferred to enhance the SM-like Higgs boson mass by large radiative correction from stop loops [41]. Meanwhile a positive A t is more favored by B physics since stop-chargino loops can mediate the transition of bottom quark to strange quark. These reasons are illustrated in the third row of Fig. 15. (vi) For each input parameter, the range covered by the 1σ CR is significantly less than that of the 1σ CI. This difference is caused by the fact that the marginal PDF defined in Eq. (3.7) of this section is determined not only by the likelihood function, but also by the parameter space. By contrast, the PL is affected only by the likelihood function.
In practice, one usually encounters the phenomenon that with the increase of the setting nlive in the MultiNest algorithm, the predicted PLs increase so that they become more or less constants over the whole scanned parameter range (see for example the κ ∼ 0 region of the PL in the third row of Fig. 16 in Appendix B). The reason is that, with a low setting of nlive, the scan is not elaborate enough so that it misses the fine tuning cases where the theoretical prediction of a ceratin physical observable approaches to its measured value by strong cancellation of different contributions or by other subtle mechanisms. This situation can be improved by increasing the sampling in the scan, 9 but it cannot change the fact that the case is rare, and hence the corresponding posterior PDF is suppressed. In a word, the volume effect usually makes the marginal PDFs and corresponding PLs quite different. Figure 2 shows the correlation of λ with the other parameters. From this figure, one can see that the 1σ and 2σ CIs on tan β − λ plane roughly overlap. One reason for this is that the observables considered in this work, especially the SM-like Higgs boson mass, are sensitive to the two parameters, and so is the likelihood function in Eq. (3.3) of this section. Consequently, any shift of the parameters can alter the likelihood value significantly. One can also see that the 2σ CIs in the figure are usually isolated. This is because the marginal PDFs in these regions are relatively small so that fewer samples in the scan concentrate on the regions. We checked that, with the increase of the nlive, the regions will be connected and both the 1σ and 2σ CIs usually expand since the improved scan can survey more fine tuning cases. In the second part of Appendix C in this work, we show the impact of different nlives on the CIs.
Before we end this subsection, we want to emphasize two points. One is that the posterior PDFs and the PLs are determined by the total likelihood function, not by any individual contribution. For a good theory, the optimal parameter space for one component of the likelihood function should be compatible with those for the other components of the function so that the Bayesian evidence of the theory is not suppressed. In Appendix A, we show the impacts of different experimental measurements on the distribution of the parameters in Higgs sector. By comparing the figures presented in this Appendix, one can learn that the Higgs observables play an important role in FIG. 3. Similar to Fig. 2, but projected on m A 1 − m˜ν 1 , m h 2 − m˜ν 1 , and μ − m˜ν 1 planes with the red dashed lines corresponding to the cases of m A 1 ¼ m˜ν 1 , m h 2 ¼ 2m˜ν 1 , and μ ¼ m˜ν 1 , respectively. determining most of the results, and the other observables are concordant in selecting favorable parameter space of the Type-I seesaw extension. This implies that once the mass spectrum in the Higgs sector is determined, one can adjust the parameters in the sneutrino sector to be consistent with the DM measurements. This is a simplified way to get good parameter points in the model. The other is that the posterior marginal PDFs are affected by the choice of the prior distribution for the input parameters, and they may be quite different if one compares the results of two distinct prior PDFs with insufficient sampling in the scan. In Appendix B, we investigate the difference induced by the flat distribution adopted in this work and the log distribution usually used in the scan by Markov chain method. Our study shows that although the predictions of the log distribution are stable with the increase of the setting nlive, the prior distribution overemphasizes the low jA κ j and low μ regions, and thus deforms the shape of the posterior PDFs from what the underlying physics predicts. By contrast, the flat distribution is able to predict the right results, but at the cost of a long time calculation by clusters when one sets a large nlive.

C. DM annihilation mechanisms
In this subsection, we study the annihilation mechanisms of the CP-even sneutrino DM. For this purpose, we project the posterior PDF and PL on m A 1 − mν 1 , m h 2 − mν 1 , and μ − mν 1 planes in Fig. 3, and use the red dashed lines to denote the cases of m A 1 ¼ mν 1 , m h 2 ¼ 2mν 1 , and μ ¼ mν 1 , respectively. We also show in Fig. 4  (i) For the samples in 1σ CRs, the singlet dominated A 1 and h 2 are heavier than about 300 GeV. (ii) A small portion of samples satisfy m A 1 ≲ mν 1 , which implies that the DM may annihilate dominantly into A 1 A 1 final state. Since this annihilation channel is mainly a s-wave process, hσvi 0 ≃ hσvi FO ≃ 3 × 10 −26 cm 3 s −1 with hσvi 0 denoting current DM annihilation rate and hσvi FO representing the rate at freeze-out temperature, and consequently the annihilation is tightly limited by the Fermi-LAT data from dwarf galaxies [106]. As was pointed in our previous work [48], such a constraint can be avoided by the forbidden annihilation proposed in [54] for m˜ν 1 ≲ 100 GeV 10 or simply by requiring mν 1 > 100 GeV where the constraints become loose [106]. Our results verify this point. Note that for the case considered here, the sneutrino DM and the singlet-dominated A 1 compose a self-contained DM sector, which interacts with the SM sector only via the small Higgs-doublet component of A 1 , and the corresponding likelihood function may be quite large, which is shown by the PL distribution on m A 1 − mν 1 plane. This is a typical hidden or secluded DM scenario [107]. Note that this case is rare since it requires specific kinematics, i.e., A 1 should be slightly heavier thanν 1 . which implies hσvi FO > hσvi 0 if the h 2 -mediated contribution is dominant in the annihilation [54].
(iv) For all samples in the 1σ CR on the μ − mν 1 plane, μ and mν 1 are approximately equal. In this case, the DM and Higgsinos can be in thermal equilibrium in early Universe before DM freeze-out due to the transition between the DM pair and the Higgsino pair, which is mediated by the singlet-dominated Higgs bosons. As a result,ν 1 can reach its correct relic density by coannihilating with the Higgsinodominated neutralinos, although its couplings with the SM particles are tiny. This situation is quite similar to the sneutrino DM scenario in the inverse seesaw extension of the NMSSM [48]. Note that a heavyν 1 corresponds to a larger parameter space for the annihilation, so the posterior PDF is maximized at μ ≃ 300 GeV. In Appendix C of this work, we will further discuss the situation of the marginal PDF for μ > 300 GeV case. (v) The lightest CP-odd sneutrinoν 1;I may roughly degenerate in mass withν 1 . In this case,ν 1 can coannihilate withν 1;I to reach the right relic density by exchanging the CP-odd scalar A 1 . In particular, the relation m A 1 ≃ m˜ν 1 þ m˜ν 1;I may hold when 55 GeV ≤ν 1 ≤ 190 GeV, which means that the coannihilation is mediated by a resonant A 1 .

(vi) From the PL distributions in Figs. 3 and 4, one can
infer that the DM may also annihilate by a resonant h 1 , but the probability of the occurrence is small. (vii) We emphasize that these results are based on flat prior PDF of the input parameters, which has no bias on these parameters before the scan. If one focuses on low jA κ j region, e.g., by choosing a log prior PDF of jA κ j, a light A 1 becomes preferred, and consequently more samples below the red line m A 1 ¼ mν 1 can be obtained. We will discuss this issue in Appendix B. In summary, Figs. 3 and 4 reveal the fact that the singlet dominated Higgs boson may either serve as the DM annihilation product, or mediate the DM annihilations in the type-I extension of the NMSSM, and consequently, the lightest sneutrino in the theory can act as a viable DM candidate. This situation is quite different from the type-I extension of the MSSM. Moreover, we note that thẽ ν 1 -Higgsino coannihilation was neglected in previous studies, but from our study it is actually the most important annihilation mechanism from both the marginal PDF and PL distributions.

D. DM-nucleon scattering
In Sec. II, we have shown by analytic formulas the natural suppression of the cross section forν 1 -nucleon scattering. In this subsection, we present relevant numerical results. In Fig. 5, we plot the marginal posterior PDFs and PLs on ξ 2 − ξ 1 planes (top panels) and σν 1 −p − mν 1 planes (bottom panels). The top panels indicate that the typical magnitudes of ξ i (i ¼ 1; 2) are 10 −11 GeV −3 , and both the CRs and CIs are not symmetric under the sign exchange of ξ i . Correspondingly, the cross section forν 1 -nucleon scattering is usually less than 10 −46 cm 2 with its marginal PDF maximized around 4 × 10 −47 cm 2 , and may be as low as 10 −49 cm 2 for a significant part of samples (see bottom panels).
From Eqs. (2.7) and (2.16) in Sec. II and the favored parameter space in Fig. 2, one can learn that the scattering rate is sensitive to the parameters λ ν and A λ ν , so we concentrate on these two parameters in following discussion. In the upper panels of Fig. 6, we show the posterior PDF and PL on A λ ν − λ ν plane, and in the lower panels we plot the distribution of the fine tuning parameter Δ FT defined in Eq. (2.17) of Sec. II on the same plane. The upper panels indicate that both the 1σ posterior PDF and 1σ PL spread a broad region on the plane, reflecting that the theory can readily accommodate the tight constraints of the XENON1T-2018 experiment. The lower panels, on the other hand, show that the theory does not need any fine tuning to survive the tight constraints.
In a similar way, one may discuss the physics about a CP-odd sneutrino DM. The results are presented in Fig. 7 for 1D posterior PDFs and PLs of the input parameters, and in Figs. 8-12 for 2D distributions which are analogy with Figs. 2-6. By comparing each panel for the CP-odd case with its corresponding one for the CP-even case, one can learn that the two cases are somewhat similar except that the former favors positive κ and A λ ν and a negative A κ , while the latter prefers opposite signs. Such a difference can be understood from following facts: (i) As indicated by the expression of sneutrino mass in Eq. (2.5) of Sec. II, positive κ and A λ ν can predict a CP-odd sneutrino state which is lighter than its corresponding CP-even state. couplings of the CP-odd state are related with those of the CP-even state by the substitution κ → −κ and A λ ν → −A λ ν . Moreover, the analysis in Appendix A reveals that the Higgs spectrum is invariant under the transformation κ → −κ and A κ → −A κ if A λ ≫ κv s . These intrinsic relations imply that, as a rough approximation, the results of the CP-odd case can be obtained from those of the CP-even case by the replacement κ → −κ, A κ → −A κ , and A λ ν → −A λ ν , which is justified by our calculation. (iii) In any case, a negative κA κ is favored by the positiveness of the squared mass for the CP-odd singlet field. Related discussions about the parameters κ and A κ can be found in Appendixes A and B.
Besides, we compute the Bayesian evidence for the CP-odd case, and find ln Z ¼ −50.12 AE 0.08, which is slightly larger than ln Z ¼ −50.39 AE 0.06 for the CP-even case.
Since the Jeffreys' scale defined in [108] is 0.27 and small, we infer that the data has no significant preference of the CP-odd case over the CP-even case. 11 We also show the detailed information of the best fit point for the two cases in Table III. We find that the χ 2 s of the best points are roughly FIG. 7. Same as Fig. 1, but for the case that the DM candidateν 1 is CP-odd. 11 Given two scenarios to be compared beside one another, the Jeffreys' scale presents a calibrated spectrum of significance for the relative strength between the Bayesian evidences of the scenarios [108]. For the application of Jeffreys' scale in particle physics, see for example [109]. equal, χ 2 min ≃ 74.6, with most of its contribution coming from the Higgs sector, which contains 74 observables in the peak-centered method of the package HIGGSSIGNAL [89].

IV. CONSTRAINTS FROM THE LHC EXPERIMENTS
In supersymmetric theories, moderately light Higgsinos are favored by naturalness. So they are expected to be copiously produced at the LHC, and the detection of their signals can significantly limit the theory [36,67,110,111]. In the NMSSM with Type-I seesaw mechanism, the neutral Higgsinos mix with the Singlino, the fermionic component field of the singlet superfieldŜ, to form mass eigenstates called neutralinos. Due to the tininess of the Yukawa coupling Y ν , the Higgsino-dominated neutralinos may have a much stronger coupling with theν 1 ν R state, which is FIG. 10. Similar to Fig. 4, but for the case that the DM candidateν 1 is CP-odd. induced by the λ νŝνν term in the superpotential, than with theν 1 ν L state, which comes from the Y νl ·Ĥ uν interaction. The decay product of the neutralinos is then complicated by the decay chainχ 0 →ν 1 ν R with ν R → W ðÃÞ l, Z ðÃÞ ν L , h ðÃÞ ν L [112]. 12 On the other hand, this situation is simplified greatly for the Higgsino-dominated chargino, which may decay mainly by the channelχ AE →ν 1 τ AE if the LSPν 1 carries τ flavor, and meanwhile the magnitude of the 33 element in the matrix Y ν is much larger than that of the other elements. 13 So in this subsection we investigate the constraints of the ATLAS analysis at 8 TeV LHC on the signal of two hadronic τs plus E miss T [115], which arises from the process pp →χ AE 1χ ∓ 1 → 2τ þ E miss T in the type-I extension. We note that recently both the ATLAS and CMS collaboration updated their analysis on the 2τ þ E miss T signal with the data collected at 13 TeV LHC [113,[116][117][118]. As a comparison, we also study the constraints from the renewed ATLAS analysis [116].
In our calculation, we use the simulation tools MADGRAPH/MADEVENT [119,120] to generate the parton level events of the processes, PYTHIA6 [121] for parton fragmentation and hadronization, DELPHES [122] for fast simulation of the performance of the ATLAS detector, and CHECKMATE [123][124][125] to implement the cut selections of the analysis. The validation on the implementation of the analysis in CHECKMATE was provided in our work [48,126].
The procedure to get the constraints is as follows: we first determine the signal region (SR) with the largest expected sensitivity for a given parameter point ðmν 1 ; m˜χAE 1 Þ (see footnote 6 in [48] for more details), then we 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. Obviously, R represents the capability of the LHC in exploring the point. R > 1 implies that the point is excluded, or else it is allowed. In  Fig. 6, but for the case that the DM candidateν 1 is CP-odd. 12 If the neutralino decayχ 0 →ν 1 ν R is kinematically forbidden, χ 0 has to decay intoν 1 ν L assuming no mass splitting among the Higgsino-dominated particles. In this case, the experimental constraints on the monojet signal from the neutralino pair production is rather weak, which has been shown in our previous work [48]. 13 We remind that although the effect of the neutrino Yukawa coupling Y ν on sneutrino mass can be safely neglected, the nondiagonality of Y ν can induce flavor changing decays χ AE 1 →ν 1 e,ν 1 μ, which has been tightly limited by the latest CMS search for sleptons if the mass splitting betweenχ AE 1 andν 1 is sizable [113]. From theoretical point of view, these decays can be suppressed greatly if the 33 element of Y ν is much larger than the other elements. In practice, one can easily obtain the hierarchy structure by adopting the Casas-Ibarra parametrization of Y ν and scanning randomly the angles of the orthogonal matrix R and the right-handed neutrino masses involved in the parametrization [114]. learn that R increases monotonously with the enlarged mass splitting Δm ≡ m˜χAE 1 − mν 1 for fixed m˜χAE 1 , and R ¼ 1 corresponds to Δm ¼ 65 GeV, 70 GeV, 85 GeV, 150 GeV, respectively, for the m˜χAE 1 s. The underlying physics is that a large mass splitting tends to enhance the cut efficiency, and the number of the signal event S depends on both the cut efficiency and the production rate σðpp →χ AE 1χ ∓ 1 Þ, which decreases as the chargino becomes heavy [127]. The right panel shows similar features, except that the R values are significantly smaller than corresponding ones at 8 TeV LHC. This reflects that the analysis at the 8 TeV LHC is more powerful in limiting the model than that at the 13 TeV LHC when m˜χAE 1 ≲ 260 GeV. The difference comes from the fact that the updated analysis focuses on heavy chargino case, which requires more energetic jets and larger missing energy than the original analysis.
With respect to the results in Fig. 13, two points should be noted. One is that we have assumed Brðχ AE 1 →ν 1 τÞ ¼ 100% in getting the figure. But in practice,ν 1 may be eitherν 1;R orν 1;I , and the two channelsχ AE 1 →ν 1;R τ,ν 1;I τ can be kinematically accessible simultaneously. Moreover, the chargino may decay first intoχ 0 1 jj 0 =π AE after considering the mass splitting among the Higgsino-dominated particles [67]. So in case of Brðχ AE 1 →ν 1 τÞ < 1, the R value in Fig. 13 should be rescaled by the factor of Br 2 ðχ AE 1 →ν 1 τÞ, which can weaken the constraint. The other is that, from the plots on the μ −ν 1 plane in Figs. 3 and 9, the mass splitting Δm ¼ μ − mν 1 is of Oð10 GeVÞ for most samples. In this case, the LHC searches for the 2τ signal actually have no exclusion ability.

V. CONCLUSIONS
Given the strong tension between the naturalness for Z boson mass and the DD experiments for customary neutralino DM candidate in minimal supersymmetric theories, it is essential to explore the DM physics in any extension of the MSSM or the NMSSM. In this work, we augment the NMSSM with type-I seesaw mechanism,  In getting this figure, we choose m˜χAE ¼ 120 GeV, 160 GeV, 200 GeV, 240 GeVas benchmark points, and calculate the cross sections of the process at next-to-leading order by the code PROSPINO [127]. We also assume Brðχ AE 1 →ν 1 τÞ ¼ 100%. If the assumption is not satisfied, the R value should be rescaled by the factor Br 2 ðχ AE 1 →ν 1 τÞ.
which is the simplest extension to reconcile the neutrino nonzero masses and neutrinos oscillation experiment, and carry out a comprehensive study on sneutrino DM physics. The highlight of the theory is that 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 the DM 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 contributing to DM-nucleon scattering rate. Moreover, since sneutrino DM does not interact with SM particles directly, the scattering of the sneutrino DM with nucleon can be suppressed in a natural way and by several mechanisms so that the tension is alleviated greatly even after considering the latest XENON1T results. In order to illustrate these features, we carry out a sophisticated scan over the vast parameter space by Nested Sampling method, and adopt both Bayesian and frequentist statistical quantities to analyze the favored parameter space of different scenarios, the DM annihilation mechanism as well as the behavior of the DM-nucleon scattering confronted with the tight DD constraints. To get the statistical inference, we construct the likelihood function by considering Higgs data, B-physics measurements, DM relic density as well as DM direct detection and indirect detection limits. We obtain the following key conclusions: (1) The model provides a viable sneutrino DM candidate over a broad parameter space. In particular, moderately light Higgsinos, μ ≲ 250 GeV, are allowed for a large portion of the parameter space, which, as a theoretical advantage, can predict Z boson mass in a natural way. To the best of our knowledge, the model may be the minimal framework to accommodate sneutrino as WIMP DM in supersymmetric theory. (2) The DM and the singlet-dominated Higgs bosons can compose a realistic DM sector, which communicates with the SM sector only by the small doublet component of the bosons. This is a typical feature of hidden or secluded DM scenario. (3) In most cases, the DM co-annihilated with the Higgsinos to get its right relic density, which was omitted in previous studies. (4) The sneutrino DM scenarios can satisfy the tight constraints from the recent XENON1T experiment on DM-nucleon scattering without any fine tuning. (5) The sneutrino DM scenarios can naturally survive the constraints from the direct searches for electroweakinos with the final state of 2τ þ E T miss signal at the LHC. Finally, we have more explanations about this work: (i) Throughout our discussion, we do not consider the scenario that the second lightest CP-even Higgs h 2 as the discovered Higgs boson. Compared with this scenario, the situation for h 1 as the SM-like Higgs boson with a CP-even sneutrino DM is much more strongly supported by the experimental data. Numerically speaking, we find that the Jeffreys' scale is about 6.0 (5.9) for the h 2 scenario with a CP-even sneutrino DM (a CP-odd sneutrino DM). There are at least two reasons for the suppression of the evidence for the h 2 scenario. One is that the parameter space to predict m h 2 ≃ 125 GeV is relatively narrow [41]. The other is that a light h 1 can enhance the scattering of the sneutrino DM with nucleon, and in order to satisfy the constraint from the XENON-1T experiment, the DM usually annihilated by a resonant h 1 or h 2 to get its measured relic density. The tuning to get the right density is usually more than 100 [128]. Since the h 2 scenario becomes relatively unimportant after considering the latest XENON-1T experiment, we will present its features elsewhere. (ii) Although we have optimized our computer code for the calculation, it still took us about 0.4 million corehours for Intel I9 7900X CPU to finish the scans. This is challenging for our computer system. (iii) The conclusions listed above may be applied to the inverse seesaw extension of the NMSSM proposed in [48] due to the similarities of the two frameworks, and a careful study of the model is necessary. In this Appendix, we study the impacts of different experimental measurements on the marginal PDFs of the theoretical input parameters in the h 1 scenario with DM being a CP-even sneutrino. We note that in the type-I seesaw extended NMSSM, the neutrino and sneutrino sectors have little effect on the Higgs and B physics observables considered in the text. So in order to study the influence of the observables on the parameters, we first fixm˜ν ¼ 0.4 TeV × 1,λ ν ¼ 0.3 × 1, andĀ λ ν ¼ 0 with 1 denoting unit matrix in flavor space, and take the other unimportant parameters in the NMSSM except A λ from Table II, then we explore following parameter space jκj ≤ 0.7; 1 ≤ tan β ≤ 60; jA κ j ≤ 1 TeV; by the MultiNest algorithm with flat prior PDF for the input parameters.

Preference of Higgs physics on input parameters
In order to get the key features of the parameter space favored by Higgs data, we take L ¼ L Higgs in the scan, and adopt a large nlive, nlive ¼ 24000, since the involved calculation is relatively fast and we want to get the marginal PDFs as precise as possible. We realize that varying A λ may be helpful to improve our understanding on the features, so we choose A λ ¼ 2 TeV, 4 TeV, 6 TeV, and 10 TeV as four benchmark cases in the study. Part of our results, i.e., the posterior PDFs of the parameters λ, κ, tan β, A κ , μ, and A t for A λ ¼ 2 TeV, are presented in some of the right panels in Figs. [16][17][18]. From these results, we have following observations: (1) With A λ becoming larger, the most favored value of λ decreases monotonously from about 0.09 for A λ ¼ 2 TeV to about 0.035 for A λ ¼ 10 TeV. Correspondingly, the most favored value of the product λA λ increases from about 180 GeV to about 350 GeV. The product κA κ is always negative. Its distribution is peaked at κA κ ≃ −40 GeV, and extends to about −400 GeV for 1σ CR. The shape of the distribution is insensitive to the choice of A λ .
(3) With the increase of A λ , the most favored value of tan β moves upward from about 13 for A λ ¼ 2 TeV to about 22 for A λ ¼ 10 TeV. (4) For A λ ¼ 2 TeV, the marginal PDF of μ is evenly distributed in the range from 120 GeV to 300 GeV (see the right panel in the first row of Fig. 18), but with the increase of A λ , the preference of a moderately low μ gradually emerges (see the right panel of Fig. 23). (5) The marginal PDFs of the ratio μ=λ ≡ v s = ffiffi ffi 2 p and jκμj=λ are peaked at 950 GeV and 400 GeV, respectively for A λ ¼ 2 TeV, and these most favored values move to 1200 GeV and 600 GeV, respectively for A λ ¼ 10 TeV. This reflects the general conclusion that the ratios tend to increase with A λ .
(6) As for the ratio which parametrizes the degree of cancellation between different terms in 23 element of the CP-even Higgs mass matrix (see discussion below), our results indicate that its marginal PDF is peaked around 0.4 and its 1σ CR corresponds to 0.2 ≲ R 23 ≲ 1 for A λ ¼ 2 TeV. With the increase of A λ to 10 TeV, these quantities shift to about 1.8 and 1 ≲ R 23 ≲ 4 respectively. Our results also indicate that a small λ usually allows R 23 to deviate from 1 significantly, e.g., R 23 >3 is possible only when λ ≲ 0.15 for A λ ¼ 10 TeV, which is shown by the 2D marginal PDF on R 23 − λ plane. We will explain these behaviors later.
With the definition of R 23 , one can calculate its derivatives to the parameters λ, κ, tan β, μ, and A λ . The expression of these derivatives indicates that R 23 is very sensitive to tan β for the favored values of the parameters given in items 1-4.
(7) The marginal PDF of A t is insensitive to A λ , and its shape is quite similar to that in Fig. 18 for any benchmark cases of A λ .
Without a strong mixing of the singlet field with the other fields, the right sides of the approximations denote the squared masses of the singlet states. Then the positivity of M 2 P;22 requires a negative κA κ . Moreover, from item 2 and 5 one can estimate that the masses of the singlet dominated scalars are usually below 1 TeV for A λ ¼ 2 TeV. (ii) The square of the mass scale for the heavy doublet field, FIG. 18. Same as Fig. 16, but for the distributions of the parameters μ and A t .
is enhanced by a factor 1=ðλ sin 2βÞ. With the increase of A λ , this factor tends to increase and consequently the doublet scalars become heavier. We checked that M A ≳ 2 TeV (M A ≳ 4.5 TeV) in its 1σ CR for A λ ¼ 2 TeV (A λ ¼ 10 TeV). (iii) In the decoupling limit of the heavy doublet field, the effective Higgs potential at electroweak scale contains only the SM Higgs field and the singlet field. In the basis (S 2 , S 3 ), the effective squared mass matrix is if M 2 33 ≫ ð100 GeVÞ 2 . This approximation indicates that the mixing angle can be suppressed either by a small λ=κ, or by the cancellation between 2λμ and ðλA λ þ 2κμÞ sin 2β if R 23 ∼ 1. Such a suppression is favored by the 125 GeV Higgs data of the LHC. From the results in item 1, 2, and 3 or 5 and 6, one can infer that the mixing is small and it is very sensitive to tan β instead of to μ.
Moreover, the natural size of μ can be inferred by the minimization condition of the Higgs potential with respect to the field H 0 u as follows [38] 14 : In fact, we studied the marginal PDF of the expression A λ tan β−κ=λ . We found that it is peaked around 300 GeV for A λ ¼ 10 TeV, and with A λ diminishing, the value usually becomes smaller. in the effective mass matrix, one can learn that terms proportional to κ appear either in the form λκv 2 sin 2β or in the form κ 2 μ 2 v 2 =M 2 A , and obviously all of them are suppressed. As a result, the Higgs sector is approximately invariant under the transformation κ → −κ and A κ → −A κ [41], and that is why the marginal posterior PDFs of κ and A κ have a rough reflection symmetry.

Cumulative effect of various measurements on final results
In this subsection, we perform two independent scans over the parameter space with the likelihood function given by L ¼ L Higgs and L ¼ L Higgs × L BrðB s →μ þ μ − Þ × L BrðB s →X s γÞ respectively, and compare their predictions on the marginal PDFs of the input parameters with those presented in Fig. 1.
Our results are presented in Figs. 14 and 15 for the 1D marginal posterior PDFs and PLs of the parameter λ, κ, tan β, A κ , μ, and A t . Since the features of the PDFs and PLs for the parameters λ, tan β, μ, and A t have been depicted in the text, we in the following concentrate on those of κ and A κ . The main conclusions are (i) As we mentioned before, the marginal distribution of κ is roughly symmetric under the substitution κ → −κ if only the Higgs measurements are considered. Especially when κ approaches zero from either direction, the symmetry becomes more and more accurate (see the left panel of Fig. 14) since the terms which violate the symmetry are proportional to κ. Similarly, the marginal posterior PDF of A κ also exhibits a rough reflection symmetry if only the Note that this formula is based on the assumption of μ 2 ≫ jm 2 H u j; m 2 Z , and it can be used only in estimating the magnitude of μ.
Higgs observables are considered (see the left panel in the first row of Fig. 15).
We point out that a moderately large κ can suppress the mixing between the SM Higgs field and the singlet field, which may explain why the κ distribution is maximized at κ ≃ 0.6. (ii) B physics measurements have little impact on the PDFs and PLs of κ and A κ . The reason is that the parameters κ and A κ mainly affect the property of the singlet field, whose effects on B → X s γ and B s → μ þ μ − are usually negligible. (iii) The marginal PDFs of κ and A κ are modified significantly by the DM physics, and a negative κ is now favored (see the right panels in the second row of Fig. 14 and the first row of Fig. 15). This can be understood from the expression of the sneutrino mass in Sec. II, which shows that it depends on the parameter κ, and the CP-even state as the lightest sneutrino prefers a negative κ. In Table IV, we summarize the dependence of these 1D marginal PDFs and PLs on different experimental measurements. This table indicates that all the PDFs except that of μ are affected greatly by the Higgs data for A λ ¼ 2 TeV, while the PLs (except that of A t ) usually have a weak dependence on the Higgs data.

APPENDIX B: IMPACT OF DIFFERENT PRIOR PDFs ON POSTERIOR PDFs
In this section, we study the difference of the posterior PDFs induced by the choices of the flat prior PDF adopted in the text and the log prior PDF where λ, κ, and λ ν are flat distributed and the other input parameters are log distributed. From the discussion in previous section, we learn that the Higgs data play an important role in determining most of the posterior PDFs. Noting that the calculation of the L Higgs is much faster than that of the total likelihood function in Eq. (3.3) of Sec. III, we first consider the scan in Eq. (A1) of Appendix A with L ¼ L Higgs as an example to illustrate the features of the difference. We remind that even after the simplification, the sampling process is still time consuming for our computer clusters when we take a large value of the nlive, which is needed to test the stability of the results when more and more samplings are considered. We then compare the results calculated by the total likelihood function with different choices of the prior PDFs and also with different settings of nlive. experimental data, we rewrite Eq. (3.7) in Sec. III as follows where πðΘÞ is the prior PDF of the parameter Θ if one assumes that all input parameters are independent in their prior PDFs, ZðΘÞ is the Bayesian evidence for a fixed Θ and obtained by integrating over the rest parameter space, and by contrast Z is the evidence integrated over the whole parameter space. This formula provides an alternative way to calculate the marginal PDF PðΘjDÞ by the following procedures (i) First calculate ZðΘÞ for a series of Θ i ; (ii) Find the largest value of ZðΘÞπðΘÞ among the Θ i ; (iii) Normalize all ZðΘ i ÞπðΘ i Þ to this largest value, and plot the normalized ZðΘÞπðΘÞ as a function of Θ; (iv) As a byproduct, one may also get the normalized ZðΘÞ distribution in a similar way.
Since the evidence ZðΘÞ can be calculated rather precisely and meanwhile its uncertainty can be estimated, the result of this method can serve as a criteria to judge the quality of the PðΘjDÞ obtained in other ways. Moreover, by comparing the shape of the PðΘjDÞ obtained directly by the code MULTINEST with that of the ZðΘÞ distribution, one can judge which factor plays the dominant role in deciding the PðΘjDÞ, the prior PDF of Θ, ZðΘÞ, or the both. In Fig. 19, we show the normalized distribution of ZðA κ Þ (left panel) and ZðμÞ (right panel). These evidences are obtained by the scan in Eq. (A1) of Appendix A with A κ and μ fixed, respectively. They are calculated by the code MULTINEST with a relatively large nlive, nlive ¼ 6000, and consequently the uncertainty of the evidences is below 0.1%. Comparing the distributions in Fig. 19 with corresponding results in Figs. 17 and 18, one can conclude that the flat prior PDF predicts rather precisely the distributions of A κ and μ, even though they are moderately unstable with the increase of nlive, while the log prior PDF overemphasizes low jA κ j region so that the two distributions are affected greatly by the prior PDF.

Differences induced by the prior PDFs in type-I extended NMSSM
As a supplement to the results in the text, we present in this subsection the 1D marginal PDFs and PLs of A κ and μ obtained from the scan in Eq. (3.1) of Sec. III with the log prior PDF for the input parameters and the total likelihood function in Eq. (3.3). In order to show the effect of the setting nlive, we carry out two independent scans with nlive ¼ 2000 and 6000 respectively. We also compare the results with those obtained in a similar way but by the flat prior PDF. The results are shown in Fig. 20 for A κ and Fig. 21 for μ. From these figures, one can learn that (i) For the flat prior PDF, the marginal posterior PDF of A κ extends over a broad range without any strong preference of certain regions, and its relative size in low jA κ j region gradually increases with the increase of nlive. The peak of the μ posterior distribution locates around 280 GeV, and the location moves slowly towards to a lower value of μ with the increase of nlive. (ii) The marginal posterior PDFs of A κ and μ from the log prior PDF are rather stable with the increase of the setting nlive, and they differ greatly from corresponding prediction of the flat prior PDF: the A κ distribution is sharply peaked around A κ ¼ 0 region, and the μ distribution is maximized at about 160 GeV. As we pointed out before, these features are induced mainly by the log prior distributions of A κ and μ. We also compare the difference of the 2D marginal posterior PDFs on μ − mν 1 plane, which is the main conclusion of this work, induced by the two prior PDFs with the setting nlive ¼ 6000. The results are presented in Fig. 22. From this figure, one can learn that the DM prefers to coannihilate with the Higgsinos to get its measured relic density regardless the choice of the prior PDFs, and the difference mainly comes from wether a relatively small μ is favored or not. The fundamental reason of the difference, as we discussed above, is that the log prior PDF overemphasizes low jA κ j region. One can also infer that the log prior PDF is more likely to predict a light A 1 than the flat PDF so that there are more samples which allow the channel ν 1ν1 → A 1 A 1 to happen in DM annihilation.
APPENDIX C: OTHER RELATED ISSUES 1. Subtleness about the marginal PDF of μ From the formula of marginal PDF in Eq. (3.7) of Sec. III and also from the discussion in Appendix A, one can learn that the posterior PDF depends not only on the parameter space considered in the scan, but also on the choice of A λ . Considering that the posterior PDF of μ is relatively sensitive to these settings among the input parameters of the theory, we further study its features in the following by changing the boundary of the space and the value of A λ .
We first repeat the scan in Eq. (A1) of Appendix A with L ¼ L Higgs and A λ ¼ 2 TeV, but this time we restrict A κ within a narrow range −200 GeV ≤ A κ ≤ 200 GeV. This case is more likely to predict a light A 1 , and thus enhance the probability of occurrence for the annihilatioñ ν 1ν1 → A 1 A 1 , which is quite similar to the prediction of the log prior PDF (see the right panel of Fig. 22). As a result, one may expect that a light μ is preferred. This is verified by the left panel of Fig. 23, where we show the marginal PDF and PL of μ for this case.
Next we reset A λ to be 10 TeV and perform the same scan as that in Eq. (A1) of Appendix A with L ¼ L Higgs . We find the Bayesian evidence for the new A λ is reduced by a factor of e 1.37 in comparison with that for A λ ¼ 2 TeV case. This reflects the fact that A λ can change the distribution of the likelihood function L over the parameter space, especially the theory becomes more tuned as one increases A λ to coincide with the Higgs data. Correspondingly, the FIG. 22. Comparison of the marginal posterior PDF on μ − m˜ν 1 plane obtained by the flat prior PDF (left panel) and the log prior PDF (right) with the total likelihood function. We have set nlive ¼ 6000 to get these panels. posterior distribution of μ may be altered significantly. The marginal PDF of μ obtained by this new A λ is presented in the right panel of Fig. 23, where one learns that moderately low values of μ are also favored.
Finally we investigate how large μ is preferred to remain consistent with all the experimental data in A λ ¼ 2 TeV case. For this end, we extend the upper bound of μ in Eq. (3.1) of Sec. III to 1 TeV, and then perform same scan as that in Sec. III A. In Fig. 24, we show the 1D marginal PDF of μ (left panel) and 2D marginal PDF on μ −ν 1 plane (right panel). The left panel indicates that the marginal PDF of μ reaches its peak around μ ≃ 300 GeV, and then falls rapidly with the increase of μ. Two reasons may account for this behavior. One is that although the Higgs data have weak constraints on the parameter μ, the electroweak symmetry breaking disfavors a very large μ, which was discussed in Appendix A. In fact, we once plotted again the right panel of Fig. 19 by allowing 100 GeV ≤ μ ≤ 1000 GeV and found that ZðμÞ decreases monotonously when μ exceeds about 350 GeV. The other reason is that one of the most important roles of Higgsinos in the theory is to coannihilate with sneutrino DM to get its right relic density. Since we have restricted the upper bound of the DM mass around 300 GeV in the scan, there are no strong motivation for the theory to prefer a large μ. This fact is illustrated in the right panel of Fig. 24.