Singlino-dominated dark matter in general NMSSM

The general Next-to-Minimal Supersymmetric Standard Model (NMSSM) describes the singlino-dominated dark-matter (DM) property by four independent parameters: singlet-doublet Higgs coupling coefficient $\lambda$, Higgsino mass $\mu_{tot}$, DM mass $m_{\tilde{\chi}_1^0}$, and singlet Higgs self-coupling coefficient $\kappa$. The first three parameters strongly influence the DM-nucleon scattering rate, while $\kappa$ usually affects the scattering only slightly. This characteristic implies that singlet-dominated particles may form a secluded DM sector. Under such a theoretical structure, the DM achieves the correct abundance by annihilating into a pair of singlet-dominated Higgs bosons by adjusting $\kappa$'s value. Its scattering with nucleons is suppressed when $\lambda v/\mu_{tot}$ is small. This speculation is verified by sophisticated scanning of the theory's parameter space with various experiment constraints considered. In addition, the Bayesian evidence of the general NMSSM and that of $Z_3$-NMSSM is computed. It is found that, at the cost of introducing one additional parameter, the former is approximately $3.3 \times 10^3$ times the latter. This result corresponds to Jeffrey's scale of 8.05 and implies that the considered experiments strongly prefer the general NMSSM to the $Z_3$-NMSSM.


Introduction
With the smooth development of various dark-matter (DM) detection experiments in the past decades, DM research has gradually become a hot topic. Massive Weakly Interacting Particles (WIMPs) are considered the most promising DM candidates because they can naturally predict the abundance measured by the Planck experiment in [1,2]. This phenomenon is called the WIMP Miracle in the literature [3,4]. One distinctive feature of the candidates is that the spin-independent (SI) cross-section of their scattering with nucleons is approximately 10 −45 cm 2 since they are usually supposed to couple to the Standard Model (SM) particles through weak interactions [5]. Nevertheless, the DM direct detection experiments give a different answer; that is, so far, the XENON-1T and PandaX-II experiments have restricted the cross-section below the order of 10 −46 cm 2 [6][7][8], which implies that the interaction between DM and nucleons is at most feeble [9]. This fact reveals a general conclusion that simple WIMP theories face significant experimental challenges.
In popular supersymmetric theories, the lightest neutralino,χ 0 1 , as the lightest supersymmetric particle (LSP), is a stable WIMP in theories' natural parameter space and usually acts as a DM candidate. It affects both the signal of supersymmetric particles at colliders and the evolution of the universe, so it has been the focus of DM physics over the past several decades [3]. With the increasing sensitivity of DM direct detection experiments, this candidate is becoming increasingly more difficult to be consistent with both the detection experiments and the abundance naturally [9]. Taking the Minimal Supersymmetric Standard Model (MSSM) as an example, the neutralino with mχ0 1 1 TeV must be Bino-dominated in its component to produce correct relic abundance. In this case, the SI and spin-dependent (SD) cross-sections rely on Higgsino mass in different ways. If the Higgsinos are moderately light, at least one of them would be considerable in size [10]. As a result, the Higgsinos must be heavier than approximately 300 GeV after considering current XENON-1T experimental constraints on the cross-sections and approximately 600 GeV supposing no DM signal detected in the future LZ experiments [10]. Additionally, a global fit of the MSSM to various experimental data known in 2017 preferred the Higgsinos to be heavier than approximately 350 GeV at a 95% confidence level [11]. This conclusion should be significantly improved since both DM detection experiments and the LHC search for supersymmetry have progressed rapidly in recent years. These facts reveal that the theory already has undergone a relatively significant fine-tuning to predict Z boson mass. With further development of DM detection experiments, the tuning will become even more significant.
The dilemma of the MSSM inspired our study of the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [12], which extends the MSSM by one gaugesinglet superfield and is another popular realization of supersymmetry, and to consider a Singlino-dominated neutralino as a DM candidate. The theory introduces the inputs λ and κ to parameterize singlet-doublet Higgs interactions and singlet Higgs self-interactions. Thus, the couplings between the DM and nucleus depend on the Higgsino composition of the DM, which, for fixed DM mass, is determined by the dimensionless parameter λ and the Higgsino mass µ tot (see discussions in this paper). When λ 0.3 and µ tot 500 GeV, the scattering cross-sections of the DM and nucleus tend to be too large to be consistent with the XENON-1T constraints under the premise of the correct DM abundance [13]. This situation remains valid for any DM mass when the constraints from the LHC search for supersymmetry are included [9,14], and it has a great impact on DM abundance [10]. For example, the annihilationsχ 0 1χ 0 1 → tt, hA s , h s A s were believed to play crucial roles in determining the abundance [5,15], where t, h, h s , and A s denote the top quark, SM-like Higgs boson, and singlet-dominated CP-even and CP-odd Higgs bosons, respectively. After considering the XENON-1T constraints, none of them, however, can be fully responsible for the correct abundance. Instead, DM prefered to co-annihilate with the Higgsinos to obtain the density, which corresponds to a correlated parameter space of 2|κ| λ with λ 0.1 [9,14]. It is worth adding that, with the improvement of the experimental sensitivity in detecting DM, the small upper bound of λ will be further reduced.
It is noted that the above conclusions about the Singlino-dominated DM are established in the NMSSM with a Z 3 symmetry. In such a framework, to ensure that the neutralino is the LSP, |κ| must be less than λ/2 [12]. Consequently, it should be small after considering the DM experiments' substantial limitation on λ. This conclusion, however, does not hold in the general NMSSM (GNMSSM) due to the emergence of Z 3 -violating terms. Specifically, the ratio of Singlino and Higgsino masses is no longer 2|κ|/λ, so |κ| may be much larger than λ to predict the Singlinodominated neutralino as the LSP. The Singlino-Higgsino complex system has four independent parameters, which may take λ, µ tot , mχ0 1 , and κ. This characteristic contrasts with that of the Z 3 -NMSSM, which only contains three input parameters, i.e., λ, µ tot , and any one of mχ0 1 and κ. An important application of these differences is that the singlet-dominated particles may form a secluded DM sector [16], which has the following salient features.
• Since the parameter κ determines the interactions among the singlet-dominated particles, the Singlino-dominated DM can achieve the correct abundance by the processχ 0 1χ 0 1 → h s A s through adjusting the value of κ.
• When the parameter λ is small, the secluded sector communicates with the SM sector only through the weak singlet-doublet Higgs mixing. In this case, the interaction between the DM and nucleus is naturally feeble.
It is emphasized that λ and κ may play different roles in the DM physics of the GNMSSM. Similar to the Z 3 -NMSSM, the parameters λ, µ tot , and mχ0 1 mainly affect the coupling between the DM and nucleons, so they have been strongly restricted by DM direct detection (DD) experiments. In contrast, the singlet fields' self-interactions can be entirely responsible for the DM density, where the parameter κ plays a crucial role. Owing to this characteristic, the GNMSSM has a broad parameter space consistent with the current DM experimental results. This fact has been neglected in the literature. Considering the dilemma of the Z 3 -NMSSM in DM physics, the DM physics of the GNMSSM are scrutinized herein. Our results are quite different from those of previous studies on Singlino-dominated DM, which worked in the framework of the Z 3 -NMSSM [9,10,13,14,[17][18][19][20][21][22][23].
The rest of this paper is organized as follows. First, the key features of the GNMSSM are introduced in Section II. In particular, an economic GNMSSM realization is considered. The research strategy and numerical results are described in Section III to show the characteristics of DM physics. Finally, conclusions are drawn in Section IV.

Theoretical preliminaries
In this section, the GNMSSM's basics are reviewed. The Singlino-dominated DM's properties will be elucidated in detail by analytic formulae for the case of massivecharge Higgs and gauginos.

General NMSSM
The GNMSSM augments the MSSM by a gauge singlet superfieldŜ that does not carry any leptonic or baryonic number. Thus, its Higgs sector containsŜ and two SU (2) L doublet superfields, . The general form of its superpotential is [12] where W Yukawa is the Yukawa term for quarks and leptons, and λ and κ are dimensionless coefficients of the trilinear terms invariant under the Z 3 symmetry. The terms characterized by the bilinear mass parameters µ and ν and the singlet tadpole parameter ξ break the Z 3 symmetry. Without the loss of generality, one can eliminate the linear term inŜ by displacing the field by a specific constant c, i.e.,Ŝ =Ŝ +c, so that only the bilinear mass terms inŜ remain to violate the symmetry. As suggested by the studies in [24][25][26][27], the µ and ν's natural smallness could stem from the breaking of an underlying discrete R symmetry, Z R 4 or Z R 8 , at a low-energy scale. Additionally, it is noticeable that a non-minimal coupling χ of a Higgs bilinear to gravity was introduced to implement a superconformal symmetry in the GNMSSM [28,29], and this coupling could drive inflation in the early universe [30]. In this case, the extra µ-term is connected with χ via gravitino mass m 3/2 , i.e., µ = 3 2 m 3/2 χ. It is the only Z 3 -breaking term in superpotential due to the superconformal symmetry breaking at the Planck scale [31]. In this study, partly motivated by the inflation-inspired model 1 and mainly for simplifying our analysis, a specific GNMSSM is investigated in which all dimensional parameters in the superpotential are zero except µ, labeled as µ-extended NMSSM (µNMSSM) hereafter. The superpotential and corresponding soft supersymmetry-breaking Lagrangian are reduced to the following forms: where H u , H d , and S denote the Higgs superfields' scalar components. Throughout this work, the B µ term is set to be zero since it plays a minor role in this study and the supersymmetry-breaking mass parameters, m 2 Hu , m 2 H d , and m 2 s , are fixed by solving the conditional equations to minimize the scalar potential. The Higgs sector is then described by parameters λ, κ, A λ , A κ , µ, and the Higgs fields' vacuum expectation values (vevs) H 0 and in the bases (A NSM , Im[S]), those for CP-odd Higgs fields are given by where V and U are unitary matrices used to diagonalize M 2 S and M 2 P , respectively. The model also predicts a pair of charged Higgs bosons, H ± = cos βH ± u + sin βH ± d , and their masses take the following simple form: In this study, the following properties of the Higgs bosons are utilized: • h corresponds to the scalar discovered at the LHC. The LHC Higgs data have 0.1 and |V SM h | ∼ 1. Furthermore, from theoretical points of view, its mass may be significantly affected by the interaction λŝĤ u ·Ĥ d , the singlet-doublet Higgs mixing as well as the radiative correction from top/stop loops [34][35][36].
• H and A H are H NSM -and A NSM -dominated states. They are approximately degenerate in mass with the charged Higgs bosons if they are significantly heavier than v. The LHC search for extra Higgs bosons, the measured property of h, and the indirect constraints from B-physics prefer m H , m A H 200 GeV [37] when the alignment condition is satisfied [38,39]. It is worth noting that, although H of several hundreds of GeV may play a significant role in the DMnucleon scattering discussed below in specific situations [40], its contribution to the scattering is usually far below current DD experimental sensitivities 2 . Thus, only the case where H and A H are above 1 TeV is considered in this paper to simplify the analysis of DM phenomenology. This case is a theoretical hypothesis, and it can be realized by simply setting A λ a significant value, e.g., A λ = 2 TeV in the following numerical study. It does not change essentially the conclusions of this paper.
• h s and A s are mainly composed of the singlet field. They may be moderately light (e.g., several tens of GeV) without conflicting with any collider constraints. As will be shown below, they play important roles in DM physics.
• In the case of massive charged Higgs bosons, the following approximations hold [5]: It is emphasized that the mass matrices M 2 S and M 2 P differ from their corresponding ones in the Z 3 -NMSSM by additional µ-induced contributions. These contributions 2 The H-mediated contribution to the scattering is proportional to tan 2 β/m 4 H , and its typical size is 10 −49 cm 2 for tan β ≤ 5 (see, e.g., Figure 6 in [5]). With the increase of tan β, the lower bound of m H increases rapidly to satisfy the constraints from the LHC search for extra Higgs bosons (see, e.g., Figure 2 in [41] for the MSSM results.). Therefore, the contribution is difficult to reach 10 −47 cm 2 .
are sometimes vital in this study, e.g., a positively considerable µ can significantly reduce the mass of h s and A s so that they may act as the annihilation product of moderately light DM.
The neutralino sector comprises a Bino fieldB 0 , a Wino fieldW 0 , Higgsino fields H 0 d andH 0 u , and a Singlino fieldS 0 . In the bases , the symmetric neutralino mass matrix is given by [12] where M 1 and M 2 are gaugino soft breaking masses and µ tot ≡ µ ef f + µ. It can be diagonalized by a rotation matrix N , and, subsequently, the mass eigenstates labeled in mass-ascending order arẽ Evidently, N i3 and N i4 characterize theH 0 d andH 0 u components inχ 0 i , and N i5 denotes the Singlino component. We add that the Higgsino mass and Singlino mass in the GNMSSM are determined by µ tot and µ ef f , respectively. As a result, |2κ/λ| may be much larger than 1 in obtaining Singlino-dominated DM, which is the starting point of this work. In addition, there is one more input parameter to describe the Singlino-Higgsino mixing system than in the case of the Z 3 -NMSSM. These characteristics distinguish the GNMSSM from the Z 3 -NMSSM, which, as introduced previously, has important implications in DM physics.
In this work, the following interactions are crucial: These expressions indicate that and the other interactions are suppressed.

DM Relic Density
In the GNMSSM, Singlino-dominated DM could achieve the measured abundance by the processesχ 0 1χ 0 1 → h s A s , tt, or its co-annihilation with one or more slightly heavier states, such as the next-to-lightest neutralinoχ 0 2 , the lightest charginoχ ± 1 , or sleptons. Under a non-relativistic approximation, the thermally averaged crosssection of DM pair annihilation can be expanded as [3,45] where a and b correspond to the s-and p-wave contributions of the process, respectively, and x ≡ mχ0 1 /T with T denoting the temperature in the early Universe.
After integrating the density function from x F = m/T F at freeze-out temperature to infinity, the present relic density can be expressed as [5] Ωh 2 = 0.12 80 where g * ∼ 80 is the number of effectively relativistic degrees of freedom, and x F ∼ 25 is obtained by solving the freeze-out equation in [45].
The key features of these annihilations are shown in the following. 1.χ 0 1χ 0 1 → h s A s This process proceeds through the s-channel exchange of Z and CP-odd Higgs bosons, and the t-channel exchange of neutralinos. Consequently, σv hsAs x F is given by [5,45] σv hsAs where the s-and t-channel contributions are approximated by Owing to the largeness of top quark mass, its amplitude is mainly contributed by the longitudinal polarization of the Z boson when the mediators of the Higgs bosons are far off-shell. As a result, σv tt x F is approximated by [5] σv tt (2.22) These formulae reflect that a moderately large λ is needed to account for the density. As will be shown below, this situation usually leads to sizable DM-nucleon scattering rates.
3. Co-annihilation with Higgsino-dominated electroweakinos This annihilation affects the abundance when the mass splitting betweenχ 0 1 and the Higgsino-dominated electroweakinos is less than approximately 10% [45,46]. In this case, σ in Eq.(2.18) should be replaced by σ ef f given by In this formula, σ ab is the cross-section of the annihilation ab → XY with a, b = χ 0 1 ,χ 0 2 ,χ ± 1 , · · · and X, Y denoting any possible final states, ∆ i ≡ (m i − mχ0 1 )/mχ0 1 (i = a, b) represents the mass splitting betweenχ 0 1 and the initial state i, and x ≡ mχ0 1 /T . g i denotes the internal degrees of freedom for the initial state i, which is 2 for a Majorana fermion and 4 for a Dirac fermion, and g ef f is defined by (

2.24)
It is evident that σ ef f decreases monotonously as the co-annihilating particles depart fromχ 0 1 in mass. Finally, the annihilationsχ 0 It is a p-wave process [47], so its effect on the abundance becomes crucial only when κ is large andχ 0 1χ 0 1 → h s A s is phase-space-suppressed. This situation happens only in some corners of the theory's parameter space.

DM-nucleon scattering cross-sections
In the limit of very massive squarks, the SD scattering ofχ 0 1 with nucleons is mediated by a Z boson. The scattering cross-section is then given by [43,44] σ SD where C p 4 × 10 −4 pb for protons and C n 3.1 × 10 −4 pb for neutrons, and from the formulae in Eq.(2.10). These expressions indicate that the SD cross-section is proportional to (λv/µ tot ) 4 and it increases as mχ0 1 approaches from below to µ tot . The SI scattering is mainly contributed by the t-channel exchange of the CP-even Higgs bosons. Correspondingly, the cross-section is given by [48,49] σ SĨ χ 0 27) where N = p, n denotes proton and neutron, respectively, and µ r is the reduced mass ofχ 0 1 and nucleon given by µ r ≡ m N mχ0 1 /(m N + mχ0 1 ). The effective couplings f (N ) are 28) where C N N h i represents the h i -nucleon coupling, which is given by Consequently, the SI cross-section is given by [14] σ SĨ (2.33) These formulae indicate that, if the LHC-discovered scalar is pure SM-like, i.e., V SM h = 1, V S h = 0, and V SM hs = 0, the expression of A is simplified as

Parameter Value Parameter Value Parameter Value
(2.34) To suppress the SI cross-section, either λ is small or mχ0 1 /µ tot sin 2β, which is called the blind-spot condition of the NMSSM [43,48]. However, if h s contains a sizable SM Higgs field component, e.g., V SM hs > 0.1, the h s -mediated contribution may be crucial when h s is much lighter than h. This may drastically cancel the h-mediated contribution. In either case, the SI cross-section is usually dominated by the h-mediated contribution, and it is sensitive to the parameters λ and µ tot . Note that this feature still holds when the H-mediated contribution is not negligibly small. Therefore, the DM direct detection experiments can tightly limit these two parameters. Finally, it is noticeable that κ may also affect the cross-section significantly, but it occurs only when κ λ and V S h is sizable. This situation is different from that of the SD cross-section, which is insensitive to κ in all cases.

Numerical Results
In this section, the characteristics of the Singlino-dominated DM are demonstrated by numerical results. A likelihood function is constructed for the nest sampling algorithm adopted in this work to guide sophisticated scans over the µNMSSM's parameter space [56,57]. The scanned samples are studied by statistical quantities such as the posterior probability density function (PDF) in Bayesian statistics and the profile likelihood (PL) in frequency statistics. The former reflects the preference of the obtained samples to the input parameters, and the latter, however, represents the theory's capability to explain experimental data. These statistical quantities have been introduced concisely in [58] and are applied to scrutinize the phenomenology of the NMSSM's seesaw extensions [10,59,60].

Research Strategy
The following likelihood function is taken: where L Ωh 2 is a Gaussian likelihood function for the measured DM abundance, L DD encodes the upper bounds on the SI and SD cross-sections from the XENON-1T experiments [6,61], and L IDD incorporates the DM indirect search information of the Fermi-LAT experiment [62,63]. L Higgs includes a fit of the h's property to corresponding data of the LHC by the code HiggsSignal-2.2.3 [64]. It also includes the constraints from the direct search for extra Higgs bosons at the LEP, Tevatron, and LHC, which are implemented by the code HiggsBounds-5.3.2 [65]. Finally, L B takes a Gaussian form to describe the measured branching ratios of B → X s γ and B 0 s → µ + µ − [66]. The concrete expression of L was given in our previous works [59,60].
The following parameter space is scanned: 0 < µ ef f ≤ 500 GeV, 100 GeV ≤ |µ tot | ≤ 500 GeV, |A κ | ≤ 1000 GeV, by assuming the input parameters are flat distributed. Upper bounds are imposed on λ and κ to keep the theory perturbative up to the Grand Unification scale. µ ef f , µ tot , and A κ are restricted within relatively narrow ranges to naturally break the electroweak symmetry, and the parameter A t is varied to include the critical radiative correction of the top-stop loops to Higgs mass. In addition, a lower bound is set on |µ tot | by the LEP search for electroweakinos [66]. For the other less crucial parameters, their values are fixed as in Table 1.
The package SARAH-4.14.3 [67][68][69][70] was utilized to build the µNMSSM's model file, and the codes SPheno-4.0.3 [71,72] and FlavorKit [73] were used to generate the particle spectrum and calculate low-energy flavor observables, respectively. The Singlino-dominated neutralino is assumed to be the sole DM candidate in the Universe and the observables in DM physics are computed with the package MicrOMEGAs 5.0.4 [55,[74][75][76][77][78]. We consider h ≡ h 1 and h ≡ h 2 scenarios, where h 1 and h 2 denote the lightest and the next-to-lightest Higgs bosons. Details of each obtained sample are reported for further analysis. For example, all processes that affect the abundance are listed and their fractions to the total DM annihilation cross-section at the freeze-out temperature are presented. This information helps identify the dominant contribution to the abundance.

Posterior probability density
In Bayesian statistics, the one-dimensional marginal posterior PDF P (θ) with θ being an input parameter is defined by integrating the posterior PDF over the other input parameters. This decides the one-dimensional 1σ (2σ) credible region by requiring that the region's posterior probability accounts for 68% (95%) of the total probability [58]. In this study, each sample's posterior probability density is calculated by the nest sampling algorithm and the 1σ and 2σ credible regions for the input parameters are determined. The results are presented in Table 2 h ≡ h 1 scenario: lnZ = -65.79 ± 0.046 11.6% 0.3% Table 3. Dominant annihilation channels and their normalized posterior probabilities for h ≡ h 1 and h ≡ h 2 scenarios. In obtaining the values in this table, each sample's most critical channel for the abundance was identified and sequentially used to classify the samples. The posterior probability densities of the same type of samples were then summed.
region, one can learn that some parameters' favored range is reduced significantly. For example, λ and |κ| prefer the range λ 0.4 and |κ| 0.3 because, as mentioned above, it is more easily consistent with various DM measurements. In Table 3, dominant annihilation channels for the obtained samples and their posterior probability are listed. This table shows thatχ 0 1χ 0 1 → h s A s is the dominant channel for most cases, andχ 0 1χ 0 1 → tt is another popular channel. In contrast, the co-annihilation is vital only in rare cases in the h ≡ h 1 scenario. The characteristics of these channels were introduced in the preceding section. In the following, only several new features after elaborated analysis of the samples are listed.
• Any ofχ 0 1χ 0 1 → h s A s , the co-annihilation, and the h-funnel can be fully re-sponsible for the DM abundance. Compared with the first mechanism for the abundance, the latter two need a stricter condition, i.e., |mχ0 1 | |µ tot | and mχ0 1 m h /2, respectively. Therefore, their Bayesian evidence is relatively small 3 .
• Owing to the tight constraints of the XENON-1T experiments on λ,χ 0 1χ 0 1 → tt contributes to the total annihilation cross-section at the freeze-out temperature by less than 50%. This fact implies that, althoughχ 0 1χ 0 1 → tt is the most considerable contribution in some cases, it must be aided by other processes to achieve the correct abundance.
• Owing to its p-wave nature,χ 0 1 > m hs , 2|mχ0 1 | (m hs + m As ) and |κ| 0.3 are satisfied, its evidence tends to be small. The next object of focus is the Bayesian evidence for different scenarios. Our calculation shows that ln Z = −65.79 ± 0.046 for the h ≡ h 1 scenario and ln Z = −68.21 ± 0.051 for the h ≡ h 2 scenario. Since Jeffreys' scale defined by δ ln Z ≡ ln(Z h 1 /Z h 2 ) [80] is 2.42, it is concluded that the considered experiments prefer the h ≡ h 1 scenario slightly to the h ≡ h 2 scenario 4 . One reason for this conclusion is that the latter scenario needs tuning of its parameter space to explain the 125GeV Higgs data [34]. Also compared in the present paper is the Bayesian evidence of the µ = 0 and µ = 0 cases (note that the µ = 0 case corresponds to the Z 3 -NMSSM), and it is found that the Jeffrey's-scale result is 8.05. This result implies that the experiments strongly favor the µ = 0 case because it is readily consistent with the DM experiments. Realizing that the evidence relies on the prior probability assumptions, it is recalculated by assuming that λ and κ are log-distributed and the other input parameters are flat-distributed. It is found that Jeffrey's scale changes only slightly after recalculation.
Finally, the experimentally favored parameter space of the µNMSSM is compared with that of the Z 3 -NMSSM, which was obtained in our recent work [14]. It is  found that they are quite different. The reason is, as mentioned before, that the additional parameter µ contributes to the Higgs and neutralino mass matrices, and, consequently, the DM physics differs significantly.

Profile likelihood of DM physics
The object of focus in this subsection is the theory's capability to explain the DM experiments. DM's two-dimensional (2D) PL is defined as  where Θ = (Θ 1 ≡ λ, Θ 2 ≡ κ, · · · ) are the input parameters of L i (i = Ωh 2 , · · · ), and the maximization is realized by varying the parameters other than Θ A and Θ B . The 1σ (2σ) confidence interval on the Θ A − Θ B plane is determined by the criteria of ( Best denotes the best point's χ 2 . Two subtleties about L DM (Θ A , Θ B ) must be clarified. One is that it is completely different from the PL of the total likelihood function L in Eq. (3.1). Specifically, the χ 2 of L is extremely dominated by the Higgs data fit, which contains 102 measurements in the code HiggsSignal-2.2.3 [64]. Thus, studying the distribution of L(Θ A , Θ B ) is roughly equivalent to performing a Higgs data fit on the Θ A − Θ B plane, which has nothing to do with DM physics. The other subtlety is that some of the scanned samples may be in serious conflict with the Higgs data or B-physics measurements, although they coincide well with the DM experiments.
Since L strongly disfavors these samples, they should be neglected. Therefore, only the following samples are considered in studying L DM (Θ A , Θ B ).
• Its p value in the Higgs data fit is larger than 0.05, which implies that it coincides with the data at 95% confidence level.
• It is consistent with the extra Higgs search results implemented in the code HiggsBounds-5.3.2.
Using the refined samples, different 2D PL maps are shown in Figures 1 to 4.  µ tot = 450 GeV. Additionally, 2|κ|/λ can be much larger than 1 to keepχ 0 1 Singlinodominated, which is one distinct feature of the µNMSSM. Figure 2 shows the 1σ and 2σ CIs on the σ SĨ  Figures 1 and 2, except for the h ≡ h 2 scenario. Compared with the h ≡ h 1 scenario, this scenario has the following different characteristics.
• The CI on |κ| − λ plane shrinks significantly. In addition to the fact that the parameter space of the h ≡ h 2 scenario is relatively narrow to fit well with the Higgs data [34], another reason for the difference is that, since h s = h 1 is light, a moderately lightχ 0 1 can make the annihilationχ 0 1χ 0 1 → h s A s occur.
A detailed explanation of the PL figure, including how to plot it and understand it correctly, was presented in our previous works [10,60]. For the sake of brevity, it is not repeated here. Besides, the fact that the marginal posterior PDF and the PL differ significantly in their definitions is worth noting. Therefore, they are complementary to each other in describing the results of the present work.

Constraints from LHC search for electroweakinos
In the GNMSSM, the natural electroweak symmetry breaking tends to a Higgsino mass of several hundreds of GeV. In this case, the LHC will produce Higgsino pair events copiously, and the Higgsino's property is strongly restricted by searching for  Figure 5. Constraint from latest ATLAS analysis of di-lepton signal at 13-TeV LHC, which is shown on |µ tot |−mχ0 1 plane. Analysis is concentrated on compressed mass spectra [86] and excludes the yellow-shaded region at 95% confidence level when λ = 0.01 and tan β = 10 are fixed.
multi-lepton signals. Up to now, experimental analyses in this aspect usually considered Wino pair production and provided the Wino mass bound as a function of mχ0 1 in the simplified model. Recently, the ATLAS collaboration analyzed 139fb −1 proton-proton collision data collected at the LHC with √ s = 13TeV and concluded that the LHC had already explored the region with the Wino mass up to approximately 700GeV and mχ0 1 up to 300GeV [82][83][84]. These analyses were applied to the Higgsino pair production process in the present study by elaborate Monte Carlo simulations and it was found that they can effectively limit the Higgsino's property when mχ0 1 100 GeV 5 . As a specific application, the h-funnel region in the h ≡ h 2 scenario has been excluded by the analyses. This result is consistent with our previous observation in [9] that the area has been ruled out by CMS's search for the multi-lepton signal at the 13-TeV LHC with 35.9fb −1 data [85].
Above discussion reveals that the analyses in [82][83][84][85] have little restriction on the theory sinceχ 0 1 in the GNMSSM is preferred heavier than 100 GeV. This fact motivates us to consider an experimental analysis that can detect high-DM-mass region [86]. The analysis aims to explore compressed mass spectra and utilizes the 139fb −1 data of the LHC. Events with missing transverse momentum and two same-   Table 4. Detailed information of four benchmark points that agree well with all considered experiments. Number before each annihilation process represents its fraction in contributing to total DM annihilation cross-section at freeze-out temperature. Number before each decay denotes its branching ratio.
In Figure 5, the excluded region is shown on the |µ tot | − mχ0 1 plane. λ = 0.01, tan β = 10, and M 1 = M 2 = 3 TeV are fixed, and the procedure depicted in [59] is followed to plot the figure. The result indicates that the analysis can explore mχ0 1 up to 230 GeV and significantly impact the co-annihilation mechanism. In addition, it is noticeable that the excluded region in the present work is broader than that on the upper panel of Figure 14 in [86]. This is becauseχ 0 1 is Singlino-dominated in this work, instead of Higgsino-dominated in [86], so more processes contribute to the signal.
To fully demonstrate the Singlino-dominated DM scenario's characteristics, the details of four benchmark points are provided in Table 4. The first three points belong to the h ≡ h 1 scenario and the last one is from the h ≡ h 2 scenario. They all predict mχ0 1 230GeV and thus are difficult to detect at the LHC. In addition to the characteristics listed in the table, the following features merit emphasizing.
• If the mass splitting betweenχ 0 i (i = 2, 3) andχ 0 1 is larger than any of the neutral Higgs boson masses,χ 0 i will decay with sizable branching ratios into the Higgs bosons. Otherwise, it will decay into a real or virtual Z boson.
• Although the processχ 0 1χ 0 1 → h s A s is mainly responsible for the abundance in most cases, its role diminishes when m hs + m As approaches 2mχ0 1 from below. This situation is demonstrated by Points 1 and 2.
• Since V S h = 0.3 is sizable,χ 0 1χ 0 1 → hA s for Point 1 also plays a crucial role in determining the abundance.

Conclusions
Motivated by the increasingly tight limitation of the direct detection experiments on traditional neutralino DM in the natural MSSM and Z 3 -NMSSM, the Z 3 -NMSSM is extended in the present work by adding a µĤ u ·Ĥ d term in its superpotential, and whether the Singlino-dominated neutralino can act as a feasible DM candidate is studied. Different from the Z 3 -NMSSM, the extended theory describes the neutralino's property by four independent parameters: λ, µ tot , mχ0 1 , and κ. The first three parameters strongly influence the DM-nucleon scattering rate, while κ usually only slightly affects the scattering. This characteristic implies that singlet-dominated particles may form a secluded DM sector. Under such a theoretical structure, the Singlino-dominated neutralino achieves the correct abundance by annihilating into a pair of singlet-dominated Higgs bosons by adjusting κ's value. Its scattering with nucleons is suppressed when λv/µ tot is small.
Our speculations are verified by numerical results. Specifically, a likelihood function containing the current experimental and theoretical situations of DM physics, Higgs physics, and B physics was constructed, which was then utilized to guide sophisticated scans of the theory's parameter space by the nest sampling algorithm. The scanned samples were analyzed by several statistical quantities, such as the marginal posterior PDF and profile likelihood, to reveal the underlying physics. The results show that, in the theory's natural space for electroweak symmetry breaking, the DM obtains the correct abundance byχ 0 1χ 0 1 → h s A s in most cases, and the SI and SD cross-sections may be as low as 10 −51 cm 2 and 10 −46 cm 2 , respectively. Furthermore, the LHC search for electroweakinos restricts the theory very weakly.
It is emphasized herein that the Z 3 -NMSSM differs significantly from the extended theory in at least three aspects. First, mχ0 1 and κ are no longer independent once one takes λ and µ tot as inputs to study the DM's property. Second, |κ| must be less than λ/2 to keep the lightest neutralino Singlino-dominated. Its magnitude should be small after considering the tight constraints from the direct detection experiments. Third, due to the absence of the µ-induced contributions to the Higgs squared mass, the singlet-dominated Higgs particles tend to be more massive. Combined with these facts, the DM cannot obtain the correct abundance byχ 0 1χ 0 1 → h s A s . Instead, it must co-annihilate with the Higgsinos to keep consistent with the DM experimental results, which corresponds to the correlated parameter space λ 2|κ| with λ 0.1. The Bayesian evidence of the two theories was compared and it was found that their Jeffrey's-scale value is 8.05. This result implies that the considered experiments strongly prefer the extended theory. Since it resurrects the Z 3 -NMSSM's broad parameter space that has been experimentally excluded, the extended theory is attractive and worthy of a careful study.