# Singlino-dominated dark matter in general NMSSM

Junjie Caoa, Demin Lia, Jingwei Liana, Yuanfang Yuea, and Haijing Zhoua

###### Abstract

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.

##  1 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}{\rm 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}{\rm 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,
$\tilde{\chi}_{1}^{0}$, 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_{\tilde{\chi}_{1}^{0}}\lesssim 1~{}{\rm 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~{}{\rm 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
gauge-singlet superfield and is another popular realization of supersymmetry,
and to consider a Singlino-dominated neutralino as a DM candidate. The theory
introduces the inputs $\lambda$ and $\kappa$ 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
$\lambda$ and the Higgsino mass $\mu_{tot}$ (see discussions in this paper).
When $\lambda\gtrsim 0.3$ and $\mu_{tot}\lesssim 500~{}{\rm 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
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to t\bar{t},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|\kappa|\simeq\lambda$ with $\lambda\lesssim 0.1$ [9,
14]. It is worth adding that, with the improvement of the experimental
sensitivity in detecting DM, the small upper bound of $\lambda$ 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, $|\kappa|$ must be less than
$\lambda/2$ [12]. Consequently, it should be small after considering the DM
experiments’ substantial limitation on $\lambda$. 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|\kappa|/\lambda$, so $|\kappa|$ may be much larger than
$\lambda$ to predict the Singlino-dominated neutralino as the LSP. The
Singlino-Higgsino complex system has four independent parameters, which may
take $\lambda$, $\mu_{tot}$, $m_{\tilde{\chi}_{1}^{0}}$, and $\kappa$. This
characteristic contrasts with that of the $Z_{3}$-NMSSM, which only contains
three input parameters, i.e., $\lambda$, $\mu_{tot}$, and any one of
$m_{\tilde{\chi}_{1}^{0}}$ and $\kappa$. 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 $\kappa$ determines the interactions among the singlet-
dominated particles, the Singlino-dominated DM can achieve the correct
abundance by the process $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}A_{s}$ through adjusting the value of $\kappa$.

  * •

When the parameter $\lambda$ 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 $\lambda$ and $\kappa$ may play different roles in the
DM physics of the GNMSSM. Similar to the $Z_{3}$-NMSSM, the parameters
$\lambda$, $\mu_{tot}$, and $m_{\tilde{\chi}_{1}^{0}}$ 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 $\kappa$ 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 [17, 18, 19, 20, 21, 13, 9, 10, 14, 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.

##  2 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
massive-charge Higgs and gauginos.

###  2.1 General NMSSM

The GNMSSM augments the MSSM by a gauge singlet superfield $\hat{S}$ that does
not carry any leptonic or baryonic number. Thus, its Higgs sector contains
$\hat{S}$ and two $SU(2)_{L}$ doublet superfields,
$\hat{H}_{u}=(\hat{H}_{u}^{+},\hat{H}_{u}^{0})$ and
$\hat{H}_{d}=(\hat{H}_{d}^{0},\hat{H}_{d}^{-})$. The general form of its
superpotential is [12]

$\displaystyle W_{\rm GNMSSM}$ $\displaystyle=$ $\displaystyle W_{\rm
Yukawa}+\lambda\hat{S}\hat{H_{u}}\cdot\hat{H_{d}}+\frac{1}{3}\kappa\hat{S}^{3}+\mu\hat{H_{u}}\cdot\hat{H_{d}}+\frac{1}{2}\nu\hat{S}^{2}+\xi\hat{S},$
(2.1)

where $W_{\rm Yukawa}$ is the Yukawa term for quarks and leptons, and
$\lambda$ and $\kappa$ are dimensionless coefficients of the trilinear terms
invariant under the $Z_{3}$ symmetry. The terms characterized by the bilinear
mass parameters $\mu$ and $\nu$ and the singlet tadpole parameter $\xi$ break
the $Z_{3}$ symmetry. Without the loss of generality, one can eliminate the
linear term in $\hat{S}$ by displacing the field by a specific constant $c$,
i.e., $\hat{S}=\hat{S}^{\prime}+c$, so that only the bilinear mass terms in
$\hat{S}^{\prime}$ remain to violate the symmetry. As suggested by the studies
in [24, 25, 26, 27], the $\mu$ and $\nu$’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 $\chi$ 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 $\mu$-term is
connected with $\chi$ via gravitino mass $m_{3/2}$, i.e.,
$\mu=\frac{3}{2}m_{3/2}\chi$. It is the only $\rm 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 model111It
must be emphasized here that, contrary to studies in the inflation-inspired
model [31, 33], we concern ourselves with the $\mu$ extension’s generality
irrespective of its origin, and therefore the inflation explanation and its
corresponding gravitino DM scenario are not discussed. and mainly for
simplifying our analysis, a specific GNMSSM is investigated in which all
dimensional parameters in the superpotential are zero except $\mu$, labeled as
$\mu$-extended NMSSM ($\rm\mu NMSSM$) hereafter. The superpotential and
corresponding soft supersymmetry-breaking Lagrangian are reduced to the
following forms:

$\displaystyle W_{\rm\mu NMSSM}$ $\displaystyle=$ $\displaystyle W_{\rm
Yukawa}+(\lambda\hat{S}+\mu)\hat{H_{u}}\cdot\hat{H_{d}}+\frac{1}{3}\kappa\hat{S}^{3},$
$\displaystyle-\mathcal{L}_{soft}$ $\displaystyle=$
$\displaystyle\Bigg{[}A_{\lambda}\lambda SH_{u}\cdot
H_{d}+\frac{1}{3}A_{\kappa}\kappa S^{3}+B_{\mu}\mu H_{u}\cdot
H_{d}+h.c.\Bigg{]}$ (2.2)
$\displaystyle+m^{2}_{H_{u}}|H_{u}|^{2}+m^{2}_{H_{d}}|H_{d}|^{2}+m^{2}_{s}|S|^{2},$

where $H_{u}$, $H_{d}$, and $S$ denote the Higgs superfields’ scalar
components. Throughout this work, the $B_{\mu}$ term is set to be zero since
it plays a minor role in this study and the supersymmetry-breaking mass
parameters, $m^{2}_{H_{u}}$, $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 $\lambda$, $\kappa$, $A_{\lambda}$,
$A_{\kappa}$, $\mu$, and the Higgs fields’ vacuum expectation values (vevs)
$\left\langle H_{u}^{0}\right\rangle=v_{u}/\sqrt{2}$, $\left\langle
H_{d}^{0}\right\rangle=v_{d}/\sqrt{2}$, and $\left\langle
S\right\rangle=v_{s}/\sqrt{2}$ with $v=\sqrt{v_{u}^{2}+v_{d}^{2}}\simeq
246~{}\mathrm{GeV}$. The ratio of the vevs $\tan{\beta}\equiv v_{u}/v_{d}$ and
$\mu_{eff}\equiv\lambda v_{s}/\sqrt{2}$ are chosen as input parameters.

Concerning the Higgs sector, it is convenient to work with the field
combinations $H_{\rm SM}\equiv\sin\beta Re[H_{u}^{0}]+\cos\beta
Re[H_{d}^{0}]$, $H_{\rm NSM}\equiv\cos\beta Re[H_{u}^{0}]-\sin\beta
Re[H_{d}^{0}]$, and $A_{\rm NSM}\equiv\cos\beta Im[H_{u}^{0}]-\sin\beta
Im[H_{d}^{0}]$ [32, 34], where $H_{\rm SM}$ represents the Higgs field in the
SM with the vev $v/\sqrt{2}$, $H_{\rm NSM}$ denotes the other CP-even doublet
Higgs field with a vanishing vev, and $A_{\rm NSM}$ corresponds to the CP-odd
Higgs boson in the MSSM. In the bases ($H_{\rm NSM}$, $H_{\rm SM}$, $Re[S]$),
the elements of the CP-even Higgs fields’ squared mass matrix are formulated
as [33]

$\displaystyle{\cal M}^{2}_{S,11}$ $\displaystyle=$
$\displaystyle\frac{2\mu_{eff}(\lambda
A_{\lambda}+\kappa\mu_{eff})}{\lambda\sin
2\beta}+\frac{1}{2}(2m_{Z}^{2}-\lambda^{2}v^{2})\sin^{2}2\beta,$
$\displaystyle{\cal M}^{2}_{S,12}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}(2m_{Z}^{2}-\lambda^{2}v^{2})\sin 4\beta,$
$\displaystyle{\cal M}^{2}_{S,13}$ $\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2}}(\lambda A_{\lambda}+2\kappa\mu_{eff})v\cos
2\beta,$ $\displaystyle{\cal M}^{2}_{S,22}$ $\displaystyle=$ $\displaystyle
m_{Z}^{2}\cos^{2}2\beta+\frac{1}{2}\lambda^{2}v^{2}\sin^{2}2\beta,$
$\displaystyle{\cal M}^{2}_{S,23}$ $\displaystyle=$
$\displaystyle\frac{v}{\sqrt{2}}\left[2\lambda\mu_{eff}+2\lambda\mu-(\lambda
A_{\lambda}+2\kappa\mu_{eff})\sin 2\beta\right],$ $\displaystyle{\cal
M}^{2}_{S,33}$ $\displaystyle=$ $\displaystyle\frac{\lambda A_{\lambda}\sin
2\beta}{4\mu_{eff}}\lambda v^{2}+\frac{\mu_{eff}}{\lambda}(\kappa
A_{\kappa}+\frac{4\kappa^{2}\mu_{eff}}{\lambda})-\frac{\lambda\mu}{2\mu_{eff}}\lambda
v^{2},$ (2.3)

and in the bases ($A_{\rm NSM}$, $Im[S]$), those for CP-odd Higgs fields are
given by

$\displaystyle{\cal M}^{2}_{P,11}$ $\displaystyle=$
$\displaystyle\frac{2\mu_{eff}(\lambda
A_{\lambda}+\kappa\mu_{eff})}{\lambda\sin 2\beta},$ $\displaystyle{\cal
M}^{2}_{P,22}$ $\displaystyle=$ $\displaystyle\frac{(\lambda
A_{\lambda}+4\kappa\mu_{eff})\sin 2\beta}{4\mu_{eff}}\lambda
v^{2}-\frac{3\mu_{eff}}{\lambda}\kappa
A_{\kappa}-\frac{\lambda\mu}{2\mu_{eff}}\lambda v^{2},$ $\displaystyle{\cal
M}^{2}_{P,12}$ $\displaystyle=$ $\displaystyle\frac{v}{\sqrt{2}}(\lambda
A_{\lambda}-2\kappa\mu_{eff}).$ (2.4)

Three CP-even Higgs mass eigenstates $h_{i}=\\{h,H,h_{\rm s}\\}$ and two CP-
odd Higgs mass eigenstates $a_{i}=\\{A_{H},A_{\rm s}\\}$ are then obtained by

$\displaystyle h_{i}$ $\displaystyle=$ $\displaystyle V_{h_{i}}^{\rm
NSM}H_{\rm NSM}+V_{h_{i}}^{\rm SM}H_{\rm SM}+V_{h_{i}}^{\rm S}Re[S],$
$\displaystyle a_{i}$ $\displaystyle=$ $\displaystyle U_{a_{i}}^{\rm
NSM}A_{\rm NSM}+U_{a_{i}}^{\rm S}Im[S],$ (2.5)

where $V$ and $U$ are unitary matrices used to diagonalize
$\mathcal{M}^{2}_{S}$ and $\mathcal{M}^{2}_{P}$, respectively. The model also
predicts a pair of charged Higgs bosons, $H^{\pm}=\cos\beta
H_{u}^{\pm}+\sin\beta H_{d}^{\pm}$, and their masses take the following simple
form:

$\displaystyle m^{2}_{H_{\pm}}=\frac{2\mu_{\rm eff}}{\sin
2\beta}\left(\frac{\kappa}{\lambda}\mu_{\rm
eff}+A_{\lambda}\right)+m^{2}_{W}-\lambda^{2}v^{2}.$ (2.6)

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
restricted $\sqrt{\left(V_{h}^{\rm NSM}\right)^{2}+\left(V_{h}^{\rm
S}\right)^{2}}\lesssim 0.1$ and $|V_{h}^{\rm SM}|\sim 1$. Furthermore, from
theoretical points of view, its mass may be significantly affected by the
interaction $\lambda\,\hat{s}\,\hat{H}_{u}\cdot\,\hat{H}_{d}$, the singlet-
doublet Higgs mixing as well as the radiative correction from top/stop loops
[35, 36, 34].

  * •

$H$ and $A_{H}$ are $H_{\rm NSM}$\- and $A_{\rm 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}}\gtrsim 200~{}{\rm 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 DM-nucleon scattering discussed
below in specific situations [40], its contribution to the scattering is
usually far below current DD experimental sensitivities222The $H$-mediated
contribution to the scattering is proportional to $\tan^{2}\beta/m_{H}^{4}$,
and its typical size is $10^{-49}~{}{\rm cm^{2}}$ for $\tan\beta\leq 5$ (see,
e.g., Figure 6 in [5]). With the increase of $\tan\beta$, 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}~{}{\rm cm^{2}}$..
Thus, only the case where $H$ and $A_{H}$ are above $1~{}{\rm 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_{\lambda}$ a significant value, e.g., $A_{\lambda}=2~{}{\rm 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]:

$\displaystyle m_{h_{s}}^{2}$ $\displaystyle\simeq$ $\displaystyle{\cal
M}^{2}_{S,33}-\frac{{\cal M}^{4}_{S,13}}{{\cal M}^{2}_{S,11}-{\cal
M}^{2}_{S,33}},\quad m_{A_{s}}^{2}\simeq{\cal M}^{2}_{P,22}-\frac{{\cal
M}^{4}_{P,12}}{{\cal M}^{2}_{P,11}-{\cal M}^{2}_{P,22}},$ (2.7)
$\displaystyle\frac{V_{h}^{\rm S}}{V_{h}^{\rm SM}}$ $\displaystyle\simeq$
$\displaystyle\frac{{\cal M}^{2}_{S,23}}{m_{h}^{2}-{\cal M}^{2}_{S,33}},\quad
V_{h}^{\rm NSM}\sim 0,\quad V_{h}^{\rm SM}\simeq\sqrt{1+\left(\frac{V_{h}^{\rm
S}}{V_{h}^{\rm SM}}\right)^{2}}\sim 1,$ $\displaystyle\frac{V_{h_{s}}^{\rm
SM}}{V_{h_{s}}^{\rm S}}$ $\displaystyle\simeq$ $\displaystyle\frac{{\cal
M}^{2}_{S,23}}{m_{h_{s}}^{2}-m_{h}^{2}},\quad V_{h_{s}}^{\rm NSM}\sim 0,\quad
V_{h_{s}}^{\rm S}\simeq\sqrt{1+\left(\frac{V_{h_{s}}^{\rm SM}}{V_{h_{s}}^{\rm
S}}\right)^{2}}\sim 1,\quad U_{A_{s}}^{\rm S}\simeq 1.$

It is emphasized that the mass matrices $\mathcal{M}^{2}_{S}$ and
$\mathcal{M}^{2}_{P}$ differ from their corresponding ones in the
$Z_{3}$-NMSSM by additional $\mu$-induced contributions. These contributions
are sometimes vital in this study, e.g., a positively considerable $\mu$ 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 field $\tilde{B}^{0}$, a Wino field
$\tilde{W}^{0}$, Higgsino fields $\tilde{H}_{d}^{0}$ and $\tilde{H}_{u}^{0}$,
and a Singlino field $\tilde{S}^{0}$. In the bases
$\psi^{0}=(-i\tilde{B}^{0},-i\tilde{W}^{0},\tilde{H}_{d}^{0},\tilde{H}_{u}^{0},\tilde{S}^{0})$,
the symmetric neutralino mass matrix is given by [12]

${\cal
M}=\left(\begin{array}[]{ccccc}M_{1}&0&-m_{Z}\sin\theta_{W}\cos\beta&m_{Z}\sin\theta_{W}\sin\beta&0\\\
&M_{2}&m_{Z}\cos\theta_{W}\cos\beta&-m_{Z}\cos\theta_{W}\sin\beta&0\\\
&&0&-\mu_{tot}&-\frac{1}{\sqrt{2}}\lambda v\sin\beta\\\
&&&0&-\frac{1}{\sqrt{2}}\lambda v\cos\beta\\\
&&&&\frac{2\kappa}{\lambda}\mu_{eff}\end{array}\right),$ (2.8)

where $M_{1}$ and $M_{2}$ are gaugino soft breaking masses and
$\mu_{tot}\equiv\mu_{eff}+\mu$. It can be diagonalized by a rotation matrix
$N$, and, subsequently, the mass eigenstates labeled in mass-ascending order
are

$\displaystyle\tilde{\chi}_{i}^{0}=N_{i1}\psi^{0}_{1}+N_{i2}\psi^{0}_{2}+N_{i3}\psi^{0}_{3}+N_{i4}\psi^{0}_{4}+N_{i5}\psi^{0}_{5}.$
(2.9)

Evidently, $N_{i3}$ and $N_{i4}$ characterize the $\tilde{H}_{d}^{0}$ and
$\tilde{H}_{u}^{0}$ components in $\tilde{\chi}_{i}^{0}$, and $N_{i5}$ denotes
the Singlino component. We add that the Higgsino mass and Singlino mass in the
GNMSSM are determined by $\mu_{tot}$ and $\mu_{eff}$, respectively. As a
result, $|2\kappa/\lambda|$ 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 the case of very massive gauginos and
$|m_{\tilde{\chi}_{1}^{0}}^{2}-\mu_{tot}^{2}|\gg\lambda^{2}v^{2}$, the
following approximations for the Singlino-dominated $\tilde{\chi}_{1}^{0}$
[42, 43, 44] are obtained:

$\displaystyle\frac{2\kappa}{\lambda}\mu_{eff}$ $\displaystyle\simeq$
$\displaystyle
m_{\tilde{\chi}_{1}^{0}}-\frac{1}{2}\frac{\lambda^{2}v^{2}(m_{\tilde{\chi}_{1}^{0}}-\mu_{tot}\sin
2\beta)}{m_{\tilde{\chi}_{1}^{0}}^{2}-\mu_{tot}^{2}},\quad N_{11}\sim 0,\quad
N_{12}\sim 0,$ (2.10) $\displaystyle\frac{N_{13}}{N_{15}}$ $\displaystyle=$
$\displaystyle\frac{\lambda v}{\sqrt{2}\mu_{\rm
tot}}\frac{(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot})\sin\beta-\cos\beta}{1-\left(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot}\right)^{2}},\quad\quad\frac{N_{14}}{N_{15}}=\frac{\lambda
v}{\sqrt{2}\mu_{\rm tot}}\frac{(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot})\cos\beta-\sin\beta}{1-\left(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot}\right)^{2}},$ $\displaystyle N_{15}^{2}$ $\displaystyle\simeq$
$\displaystyle\left(1+\frac{N^{2}_{13}}{N^{2}_{15}}+\frac{N^{2}_{14}}{N^{2}_{15}}\right)^{-1}$
$\displaystyle=$
$\displaystyle\frac{\left[1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot})^{2}\right]^{2}}{\left[(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot})^{2}-2(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm tot})\sin
2\beta+1\right]\left(\frac{\lambda v}{\sqrt{2}\mu_{\rm
tot}}\right)^{2}+\left[1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot})^{2}\right]^{2}}.$

These approximations imply the relation

$\displaystyle\frac{Z_{h}}{Z_{s}}$ $\displaystyle=$
$\displaystyle\left(\frac{\lambda v}{\sqrt{2}\mu_{\rm
tot}}\right)^{\\!\\!2}\frac{\left(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot}\right)^{2}-2{(m_{\tilde{\chi}_{1}^{0}}}/{\mu_{\rm tot}})\sin
2\beta+1}{\left[1-\left({m_{\tilde{\chi}_{1}^{0}}}/{\mu_{\rm
tot}}\right)^{2}\right]^{2}},$ (2.11)

where $Z_{h}\equiv N_{13}^{2}+N_{14}^{2}$ and $Z_{s}\equiv N_{15}^{2}$
represent the Higgsino and Singlino fractions in $\tilde{\chi}_{1}^{0}$,
respectively, and $Z_{s}>0.5$ for the Singlino-dominated
$\tilde{\chi}_{1}^{0}$. They also show that, for fixed $\tan\beta$, the
Higgsino compositions in $\tilde{\chi}_{1}^{0}$ depend only on $\lambda$,
$\mu_{tot}$, and $m_{\tilde{\chi}_{1}^{0}}$. In studying the
$\tilde{\chi}_{1}^{0}$’s properties, it is usually convenient to take the
three parameters and $\kappa$ as theoretical inputs.

In this work, the following interactions are crucial:

$\displaystyle C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h_{i}}$
$\displaystyle=$ $\displaystyle V_{h_{i}}^{\rm
SM}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}H_{\rm SM}}+V_{h_{i}}^{\rm
NSM}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}H_{\rm NSM}}+V_{h_{i}}^{\rm
S}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}Re[S]}$ (2.12)
$\displaystyle\simeq$ $\displaystyle\sqrt{2}\lambda\left[V_{h_{i}}^{\rm
SM}N_{15}\left(N_{13}\sin\beta+N_{14}\cos\beta\right)+V_{h_{i}}^{\rm
NSM}N_{15}\left(N_{13}\cos\beta-N_{14}\sin\beta\right)\right]$
$\displaystyle+\sqrt{2}V_{h_{i}}^{\rm S}\left(\lambda{N_{13}}{N_{14}}-\kappa
N_{15}^{2}\right)\,,$ $\displaystyle
C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}a_{i}}$ $\displaystyle=$
$\displaystyle U_{a_{i}}^{\rm
NSM}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}A_{\rm NSM}}+U_{a_{i}}^{\rm
S}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}Im[S]}$ (2.13)
$\displaystyle\simeq$ $\displaystyle-\sqrt{2}U_{a_{i}}^{\rm NSM}\lambda
N_{15}(N_{13}\cos\beta+N_{14}\sin\beta)+\sqrt{2}U_{a_{i}}^{\rm S}(\lambda
N_{13}N_{14}-\kappa N_{15}^{2}),$ $\displaystyle
C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}Z}$ $\displaystyle=$
$\displaystyle\frac{m_{Z}}{v}(N_{13}^{2}-N_{14}^{2}),\quad\quad
C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}G^{0}}\simeq-\sqrt{2}\lambda
N_{15}(N_{13}\sin\beta-N_{14}\cos\beta),$ (2.14) $\displaystyle
C_{h_{i}A_{s}A_{s}}$ $\displaystyle\simeq$ $\displaystyle\lambda
v\,V_{h_{i}}^{\rm SM}(\lambda+\kappa\sin 2\beta)+\lambda\kappa vV_{h_{i}}^{\rm
NSM}\cos 2\beta+\sqrt{2}V_{h_{i}}^{S}(\sqrt{2}\kappa^{2}v_{s}-\kappa
A_{\kappa}).$ (2.15)

These expressions indicate that

$\displaystyle
C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h_{s}}\sim-\sqrt{2}\kappa,\quad
C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}A_{s}}\sim-\sqrt{2}\kappa,\quad
C_{h_{s}A_{s}A_{s}}\sim\sqrt{2}(\sqrt{2}\kappa^{2}v_{s}-\kappa A_{\kappa}),$
(2.16) $\displaystyle C_{hA_{s}A_{s}}\sim\lambda v(\lambda+\kappa\sin
2\beta),\quad C_{HA_{s}A_{s}}\sim\lambda\kappa
v,\quad|C_{h_{s}A_{s}A_{s}}|\gg|C_{hA_{s}A_{s}}|,|C_{HA_{s}A_{s}}|.$

and the other interactions are suppressed.

###  2.2 DM Relic Density

In the GNMSSM, Singlino-dominated DM could achieve the measured abundance by
the processes $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}A_{s},t\bar{t}$, or its co-annihilation with one or more slightly heavier
states, such as the next-to-lightest neutralino $\tilde{\chi}_{2}^{0}$, the
lightest chargino $\tilde{\chi}_{1}^{\pm}$, or sleptons. Under a non-
relativistic approximation, the thermally averaged cross-section of DM pair
annihilation can be expanded as [3, 45]

$\left\langle\sigma_{A}v\right\rangle=a+b\left\langle
v^{2}\right\rangle+\mathcal{O}(\left\langle v^{4}\right\rangle)\approx
a+6\frac{b}{x}$ (2.17)

where $a$ and $b$ correspond to the s- and p-wave contributions of the
process, respectively, and $x\equiv m_{\tilde{\chi}_{1}^{0}}/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]

$\displaystyle\Omega
h^{2}=0.12\left(\frac{80}{g_{*}}\right)^{1/2}\left(\frac{x_{F}}{25}\right)\left(\frac{2.3\times
10^{-26}\mathrm{cm^{3}/s}}{\langle\sigma v\rangle_{x_{F}}}\right)\;,$ (2.18)

where $g_{*}\sim 80$ is the number of effectively relativistic degrees of
freedom, and $x_{F}\sim 25$ is obtained by solving the freeze-out equation in
[45].

The key features of these annihilations are shown in the following.

1\. $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to 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,
$\langle\sigma v\rangle_{x_{F}}^{h_{s}A_{s}}$ is given by [5, 45]

$\displaystyle\left\langle\sigma v\right\rangle_{x_{F}}^{h_{s}A_{s}}\simeq$
$\displaystyle\frac{1}{64\pi
m_{\tilde{\chi}_{1}^{0}}^{2}}\left\\{\left[1-\frac{\left(m_{h_{s}}+m_{A_{s}}\right)^{2}}{4m_{\tilde{\chi}_{1}^{0}}^{2}}\right]\left[1-\frac{\left(m_{h_{s}}-m_{A_{s}}\right)^{2}}{4m_{\tilde{\chi}_{1}^{0}}^{2}}\right]\right\\}^{1/2}|{\cal{A}}_{s}+{\cal{A}}_{t}|^{2},$

where the $s$\- and $t$-channel contributions are approximated by

$\displaystyle{\cal{A}}_{s}$ $\displaystyle\simeq$
$\displaystyle\frac{-2m_{\tilde{\chi}_{1}^{0}}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}A_{s}}C_{h_{s}A_{s}A_{s}}}{m_{A_{s}}^{2}-4m_{\tilde{\chi}_{1}^{0}}^{2}},$
$\displaystyle{\cal{A}}_{t}$ $\displaystyle\simeq$
$\displaystyle-2C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h_{s}}\,C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}A_{s}}\left[1+\frac{2m_{A_{s}}^{2}}{4m_{\tilde{\chi}_{1}^{0}}^{2}-\left(m_{h_{s}}^{2}+m_{A_{s}}^{2}\right)}\right],$
(2.19)

respectively, if there are no resonant contributions. With the formulae from
Eqs.(2.12)–(2.16), one can learn that, once the involved particles’ masses are
fixed, the density is mainly determined by the parameter $\kappa$. The authors
of Ref. [5] estimated the value of $\kappa$ to predict the density, and
concluded that

$\displaystyle|\kappa|\sim
0.15\times\left(\frac{m_{\tilde{\chi}_{1}^{0}}}{300~{}{\rm GeV}}\right)^{1/2}$
(2.20)

in the case of $|{\cal{A}}_{t}|\gg|{\cal{A}}_{s}|$.

2\. $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to t\bar{t}$

This annihilation is mediated by the s-channel exchange of $Z$ and Higgs
bosons. 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, $\langle\sigma
v\rangle_{x_{F}}^{t\bar{t}}$ is approximated by [5]

$\displaystyle\langle\sigma v\rangle_{x_{F}}^{t\bar{t}}\sim 2\times
10^{-26}\,\frac{{\rm cm}^{3}}{{\rm
s}}\left(\frac{\left|C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}G^{0}}\right|}{0.1}\right)^{2}\left(\frac{m_{\tilde{\chi}_{1}^{0}}}{300\,{\rm
GeV}}\right)^{-2},$ (2.21)

where

$\displaystyle C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}G^{0}}$
$\displaystyle\simeq$ $\displaystyle-\sqrt{2}\lambda
N_{15}\left(N_{13}s_{\beta}-N_{14}c_{\beta}\right)\simeq\frac{\lambda^{2}v}{\mu_{tot}}\frac{Z_{s}(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})\cos
2\beta}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}}.$ (2.22)

These formulae reflect that a moderately large $\lambda$ 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
$\tilde{\chi}_{1}^{0}$ and the Higgsino-dominated electroweakinos is less than
approximately $10\%$ [45, 46]. In this case, $\sigma$ in Eq.(2.18) should be
replaced by $\sigma_{eff}$ given by

$\displaystyle\sigma_{eff}$
$\displaystyle=\sum^{N}_{a,b}\sigma_{ab}\frac{g_{a}g_{b}}{g^{2}_{eff}}(1+\Delta_{a})^{3/2}(1+\Delta_{b})^{3/2}\times
exp[-x(\Delta_{a}+\Delta_{b})].$ (2.23)

In this formula, $\sigma_{ab}$ is the cross-section of the annihilation $ab\to
XY$ with
$a,b=\tilde{\chi}_{1}^{0},\tilde{\chi}_{2}^{0},\tilde{\chi}_{1}^{\pm},\cdots$
and $X,Y$ denoting any possible final states,
$\Delta_{i}\equiv(m_{i}-m_{\tilde{\chi}_{1}^{0}})/m_{\tilde{\chi}_{1}^{0}}$
($i=a,b$) represents the mass splitting between $\tilde{\chi}_{1}^{0}$ and the
initial state $i$, and $x\equiv m_{\tilde{\chi}_{1}^{0}}/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_{eff}$ is defined by

$\displaystyle
g_{eff}\equiv\sum^{N}_{a}g_{a}(1+\Delta_{a})^{3/2}exp(-x\Delta_{a}).$ (2.24)

It is evident that $\sigma_{eff}$ decreases monotonously as the co-
annihilating particles depart from $\tilde{\chi}_{1}^{0}$ in mass.

Finally, the annihilations $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
hA_{s},h_{s}h_{s}$ are briefly discussed. The former process proceeds in a
similar way as $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$. From
the expressions of $C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h}$ and
$C_{hA_{s}A_{s}}$ and the fact that
$C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h_{s}}\gg
C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h}$ and $C_{h_{s}A_{s}A_{s}}\gg
C_{hA_{s}A_{s}}$ in most cases, one can conclude that the abundance is
significantly affected only when $\lambda$ and $V_{h}^{S}$ are tremendously
large. This condition, however, has been tightly limited by the XENON-1T and
LHC experiments. Process $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}h_{s}$ proceeds by the s-channel exchange of the CP-even Higgs bosons and
the t-channel exchange of the neutralinos. It is a p-wave process [47], so its
effect on the abundance becomes crucial only when $\kappa$ is large and
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ is phase-space-
suppressed. This situation happens only in some corners of the theory’s
parameter space.

###  2.3 DM-nucleon scattering cross-sections

In the limit of very massive squarks, the SD scattering of
$\tilde{\chi}_{1}^{0}$ with nucleons is mediated by a $Z$ boson. The
scattering cross-section is then given by [43, 44]

$\displaystyle\sigma_{\tilde{\chi}_{1}^{0}-N}^{\rm SD}$ $\displaystyle\simeq$
$\displaystyle C_{N}\times\left(\frac{N_{13}^{2}-N_{14}^{2}}{0.1}\right)^{2},$
(2.25)

where $C_{p}\simeq 4\times 10^{-4}\rm~{}pb$ for protons and $C_{n}\simeq
3.1\times 10^{-4}\rm~{}pb$ for neutrons, and

$\displaystyle N_{13}^{2}-N_{14}^{2}=\big{(}\frac{\lambda
v}{\sqrt{2}\mu_{tot}}\big{)}^{2}\,\frac{Z_{s}\cos
2\beta}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}},$ (2.26)

from the formulae in Eq.(2.10). These expressions indicate that the SD cross-
section is proportional to $(\lambda v/\mu_{tot})^{4}$ and it increases as
$m_{\tilde{\chi}_{1}^{0}}$ approaches from below to $\mu_{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]

$\displaystyle\sigma^{SI}_{\tilde{\chi}^{0}_{1}-{N}}=\frac{4\mu_{r}^{2}}{\pi}|f^{(N)}|^{2},$
(2.27)

where $N=p,n$ denotes proton and neutron, respectively, and $\mu_{r}$ is the
reduced mass of $\tilde{\chi}_{1}^{0}$ and nucleon given by $\mu_{r}\equiv
m_{N}m_{\tilde{\chi}^{0}_{1}}/(m_{N}+m_{\tilde{\chi}^{0}_{1}})$. The effective
couplings $f^{(N)}$ are

$\displaystyle f^{(N)}=\sum_{h_{i}=\\{h,H,h_{\rm
s}\\}}\\!f^{(N)}_{h_{i}}=\sum_{h_{i}=\\{h,H,h_{\rm
s}\\}}\\!\frac{C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h_{i}}C_{NNh_{i}}}{2m^{2}_{h_{i}}},$
(2.28)

where $C_{NNh_{i}}$ represents the $h_{i}$-nucleon coupling, which is given by

$\displaystyle
C_{NNh_{i}}=-\frac{m_{N}}{v}\left[F^{(N)}_{d}\left(V_{h_{i}}^{\rm
SM}-\tan\beta V_{h_{i}}^{\rm NSM}\right)+F^{(N)}_{u}\left(V_{h_{i}}^{\rm
SM}+\frac{1}{\tan\beta}V_{h_{i}}^{\rm NSM}\right)\right]\,,$ (2.29)

with $F^{(N)}_{d}\equiv f^{(N)}_{d}+f^{(N)}_{s}+\frac{2}{27}f^{(N)}_{G}$ and
$F^{(N)}_{u}\equiv f^{(N)}_{u}+\frac{4}{27}f^{(N)}_{G}$. The nucleon form
factors $f^{(N)}_{q}=m_{N}^{-1}\left<N|m_{q}q\bar{q}|N\right>$ for $q=u,d,s$
represent the normalized light quark contribution to nucleon mass, and
$f^{(N)}_{G}=1-\sum_{q=u,d,s}f^{(N)}_{q}$ denotes other heavy quarks’ mass
fraction in the nucleon [50, 51]. For the default settings of the package
micrOMEGAs about $f_{q}^{N}$ [52, 53, 54], $F_{u}^{p}\simeq F_{u}^{n}\simeq
0.15$ and $F_{d}^{p}\simeq F_{d}^{n}\simeq 0.13$. Consequently, one can
conclude that
$\sigma^{SI}_{\tilde{\chi}^{0}_{1}-p}\simeq\sigma^{SI}_{\tilde{\chi}^{0}_{1}-n}$.

In the case of very massive charged Higgs bosons, $V_{h}^{\rm NSM}\sim 0$,
$V_{h_{s}}^{\rm NSM}\sim 0$, and the $H$-mediated contribution can be safely
neglected since it is further suppressed by a factor $1/m_{H}^{4}$.
Consequently, the SI cross-section is given by [14]

$\displaystyle\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-N}$ $\displaystyle\simeq$
$\displaystyle 5\times 10^{-45}\,{\rm
cm^{2}}~{}\left(\frac{\cal{A}}{0.1}\right)^{2},$ (2.30)

where

$\displaystyle{\cal{A}}=\left(\frac{125\rm GeV}{m_{h}}\right)^{2}V_{h}^{\rm
SM}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h}+\left(\frac{125\rm
GeV}{m_{h_{s}}}\right)^{2}V_{h_{s}}^{\rm
SM}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h_{s}}~{}.$ (2.31)

Furthermore, using the formulae from Eqs.(2.10)–(2.15), one may conclude that

$\displaystyle C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h}$
$\displaystyle\simeq$ $\displaystyle\frac{\mu_{tot}}{v}\,\big{(}\frac{\lambda
v}{\mu_{tot}}\big{)}^{2}\,\frac{Z_{s}V_{h}^{\rm
SM}(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot}-\sin
2\beta)}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}}+\frac{\lambda}{2\sqrt{2}}\big{(}\frac{\lambda
v}{\mu_{tot}}\big{)}^{2}\frac{Z_{s}V_{h}^{\rm S}\sin
2\beta}{\big{[}1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}\big{]}}$ (2.32)
$\displaystyle-\sqrt{2}\kappa Z_{s}V_{h}^{\rm S}\left[1+\big{(}\frac{\lambda
v}{\sqrt{2}\mu_{tot}}\big{)}^{2}\frac{1}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}}\frac{\mu_{eff}}{\mu_{tot}}\right],$
$\displaystyle C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h_{s}}$
$\displaystyle\simeq$ $\displaystyle\frac{\mu_{tot}}{v}\,\big{(}\frac{\lambda
v}{\mu_{tot}}\big{)}^{2}\,\frac{Z_{s}V_{h_{s}}^{\rm
SM}(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot}-\sin
2\beta)}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}}+\frac{\lambda}{2\sqrt{2}}\big{(}\frac{\lambda
v}{\mu_{tot}}\big{)}^{2}\frac{Z_{s}V_{h_{s}}^{\rm S}\sin
2\beta}{\big{[}1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}\big{]}}$ (2.33)
$\displaystyle-\sqrt{2}\kappa Z_{s}V_{h_{s}}^{\rm
S}\left[1+\big{(}\frac{\lambda
v}{\sqrt{2}\mu_{tot}}\big{)}^{2}\frac{1}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}}\frac{\mu_{eff}}{\mu_{tot}}\right].$

These formulae indicate that, if the LHC-discovered scalar is pure SM-like,
i.e., $V_{h}^{\rm SM}=1$, $V_{h}^{\rm S}=0$, and $V_{h_{s}}^{\rm SM}=0$, the
expression of ${\cal{A}}$ is simplified as

$\displaystyle{\cal{A}}=\left(\frac{125\rm
GeV}{m_{h}}\right)^{2}\frac{\mu_{tot}}{v}\,\big{(}\frac{\lambda
v}{\mu_{tot}}\big{)}^{2}\,\frac{Z_{s}(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot}-\sin
2\beta)}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{tot})^{2}}.$ (2.34)

To suppress the SI cross-section, either $\lambda$ is small or
$m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm tot}\simeq\sin 2\beta$, 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_{h_{s}}^{\rm SM}>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 $\lambda$ and $\mu_{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 $\kappa$ may also
affect the cross-section significantly, but it occurs only when
$\kappa\gg\lambda$ and $V_{h}^{\rm S}$ is sizable. This situation is different
from that of the SD cross-section, which is insensitive to $\kappa$ in all
cases.

##  3 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 $\mu$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].

###  3.1 Research Strategy

Parameter | Value | Parameter | Value | Parameter | Value
---|---|---|---|---|---
$A_{\lambda}$ | $2\mathrm{TeV}$ | $m_{\tilde{q}}$ | $2\mathrm{TeV}$ | $m_{\tilde{l}}$ | $2\mathrm{TeV}$
$M_{1}$ | $3\mathrm{TeV}$ | $M_{2}$ | $3\mathrm{TeV}$ | $M_{3}$ | $3\mathrm{TeV}$
Table 1: Setting of different fields’ soft-breaking parameters. $A_{\lambda}$,
$m_{\tilde{q}}$, $m_{\tilde{l}}$, and $M_{i}$ (i=1,2,3) are for Higgs,
squarks, sleptons, and different gaugino fields, respectively.

The following likelihood function is taken:

$\displaystyle\mathcal{L}=\mathcal{L}_{\Omega
h^{2}}\times\mathcal{L}_{DD}\times\mathcal{L}_{IDD}\times\mathcal{L}_{Higgs}\times\mathcal{L}_{B}$
(3.1)

where $\mathcal{L}_{\Omega h^{2}}$ is a Gaussian likelihood function for the
measured DM abundance, $\mathcal{L}_{DD}$ encodes the upper bounds on the SI
and SD cross-sections from the XENON-1T experiments [6, 61], and
$\mathcal{L}_{IDD}$ incorporates the DM indirect search information of the
Fermi-LAT experiment [63, 62]. $\mathcal{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, $\mathcal{L}_{B}$ takes a Gaussian form to
describe the measured branching ratios of $B\to X_{s}\gamma$ and
$B^{0}_{s}\to\mu^{+}\mu^{-}$ [66]. The concrete expression of $\mathcal{L}$
was given in our previous works [59, 60].

The following parameter space is scanned:

$\displaystyle 0<\lambda\leq 0.70,\quad|\kappa|\leq 0.70,\quad 1\leq{\rm
tan}\beta\leq 60,\quad|A_{t}|\leq 5~{}{\rm TeV},$ (3.2) $\displaystyle
0<\mu_{eff}\leq 500~{}{\rm GeV},\quad 100~{}{\rm GeV}\leq|\mu_{tot}|\leq
500~{}{\rm GeV},\quad|A_{\kappa}|\leq 1000~{}{\rm GeV},$

by assuming the input parameters are flat distributed. Upper bounds are
imposed on $\lambda$ and $\kappa$ to keep the theory perturbative up to the
Grand Unification scale. $\mu_{eff}$, $\mu_{tot}$, and $A_{\kappa}$ 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 $|\mu_{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
$\mu$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 [74, 55, 75, 76, 77, 78]. We
consider $h\equiv h_{1}$ and $h\equiv 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.

|  $h\equiv h_{1}$ scenario |  $h\equiv h_{2}$ scenario
---|---|---
Parameters | 1 $\sigma$ | 2 $\sigma$ | 1 $\sigma$ | 2 $\sigma$
$\lambda$ | 0.04 $\sim$ 0.24 | 0.02 $\sim$ 0.41 | 0.05 $\sim$ 0.18 | 0.03 $\sim$ 0.31
$\kappa$ | -0.25 $\sim$ 0.31 | -0.31 $\sim$ 0.41 | -0.17 $\sim$ 0.22 | -0.30 $\sim$ 0.30
$\tan{\beta}$ | 5.9 $\sim$ 23.3 | 4.7 $\sim$ 49.3 | 5.7 $\sim$ 12.1 | 5.1 $\sim$ 20.2
$\mu$ | 245 $\sim$ 420 | 163 $\sim$ 467 | 150 $\sim$ 319 | 51 $\sim$ 408
$\mu_{eff}$ | 29 $\sim$ 118 | 14 $\sim$ 253 | 33 $\sim$ 148 | 18 $\sim$ 298
$\mu_{tot}$ | 325 $\sim$ 475 | 255 $\sim$ 495 | 225 $\sim$ 425 | 136 $\sim$ 475
$A_{\kappa}$ | -375 $\sim$ 128 | -576 $\sim$ 397 | -652 $\sim$ 339 | -776 $\sim$ 534
$A_{t}$ | 3148 $\sim$ 4248 | 2849 $\sim$ 4448 | 2651 $\sim$ 4150 | -3845 $\sim$ 4550
$m_{\tilde{\chi}_{1}^{0}}$ | 240 $\sim$ 407 | 182 $\sim$ 466 | 181 $\sim$ 346 | 118 $\sim$ 404
Table 2: One-dimensional 1$\sigma$ and 2$\sigma$ credible regions for input
parameters in $h\equiv h_{1}$ and $h\equiv h_{2}$ scenarios, where $h$,
$h_{1}$, and $h_{2}$ denote the SM-like Higgs boson and the lightest and next-
to-lightest CP-even Higgs boson, respectively. Dimensional parameters are in
units of GeV. $h\equiv h_{1}$ scenario: lnZ = -65.79 $\pm$ 0.046
---
$\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;A_{s}$ | $\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;t\;\bar{t}$ | $\rm\quad\;\,\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;h_{s}$ | Co-annihilation
$88\%$ | $8\%$ | $3\%$ | $0.7\%$
$h\equiv h_{2}$ scenario: lnZ = -68.23 $\pm$ 0.051
$\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;A_{s}$ | $\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;t\;\bar{t}$ | Co-annihilation |  $h$-funnel
$76\%$ | $12\%$ | $11.6\%$ | $0.3\%$
Table 3: Dominant annihilation channels and their normalized posterior
probabilities for $h\equiv h_{1}$ and $h\equiv 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.

###  3.2 Posterior probability density

In Bayesian statistics, the one-dimensional marginal posterior PDF $P(\theta)$
with $\theta$ being an input parameter is defined by integrating the posterior
PDF over the other input parameters. This decides the one-dimensional
$1\sigma$ ($2\sigma$) 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\sigma$ and $2\sigma$ credible regions for
the input parameters are determined. The results are presented in Table 2.
Compared with the scanned region, one can learn that some parameters’ favored
range is reduced significantly. For example, $\lambda$ and $|\kappa|$ prefer
the range $\lambda\lesssim 0.4$ and $|\kappa|\lesssim 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
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ is the dominant
channel for most cases, and $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
t\bar{t}$ is another popular channel. In contrast, the co-annihilation is
vital only in rare cases in the $h\equiv 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 $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$, the co-
annihilation, and the $h$-funnel can be fully responsible for the DM
abundance. Compared with the first mechanism for the abundance, the latter two
need a stricter condition, i.e., $|m_{\tilde{\chi}_{1}^{0}}|\simeq|\mu_{tot}|$
and $m_{\tilde{\chi}_{1}^{0}}\simeq m_{h}/2$, respectively. Therefore, their
Bayesian evidence is relatively small333Bayesian evidence is usually defined
by $Z(D|M)\equiv\int{P(D|O(M,\Theta))P(\Theta|M)\prod d\Theta_{i}}$, where
$P(\Theta|M)$ is the prior PDF of input parameters
$\Theta=(\Theta_{1},\Theta_{2},\cdots)$ in model $M$ and
$P(D|O(M,\Theta))\equiv\mathcal{L}(\Theta)$ is the likelihood function for
observables $O$, which considers both theoretical predictions, $O(M,\Theta)$,
and experimental data, $D$. Computationally, the evidence represents an
average likelihood in the surveyed parameter space. It depends on the priors
of the model’s parameters. For different scenarios in a theory, the larger $Z$
is, the more readily it is for the corresponding scenario to agree with the
data. The details of Bayesian evidence and its applications were presented in,
e.g., [79]..

  * •

Owing to the tight constraints of the XENON-1T experiments on $\lambda$,
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to t\bar{t}$ contributes to the
total annihilation cross-section at the freeze-out temperature by less than
$50\%$. This fact implies that, although
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to t\bar{t}$ 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, $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}h_{s}$ contributes to the abundance by, at most, $40\%$. This situation
is quite similar to that of $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
t\bar{t}$. Furthermore, since $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}h_{s}$ dominates over $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}A_{s}$ only when the strict conditions of
$m_{\tilde{\chi}_{1}^{0}}>m_{h_{s}}$,
$2|m_{\tilde{\chi}_{1}^{0}}|\lesssim(m_{h_{s}}+m_{A_{s}})$ and
$|\kappa|\gtrsim 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\pm 0.046$ for the $h\equiv h_{1}$
scenario and $\ln Z=-68.21\pm 0.051$ for the $h\equiv h_{2}$ scenario. Since
Jeffreys’ scale defined by $\delta\ln Z\equiv\ln(Z_{h_{1}}/Z_{h_{2}})$ [80] is
$2.42$, it is concluded that the considered experiments prefer the $h\equiv
h_{1}$ scenario slightly to the $h\equiv h_{2}$ scenario 444Given two
scenarios to be compared together, Jeffreys’ scale provides a calibrated
spectrum of significance for the relative strength of the scenarios [80]. For
the application of the Jeffreys’ scale in particle physics, see, for e.g.,
[81]. . One reason for this conclusion is that the latter scenario needs
tuning of its parameter space to explain the $125{\rm GeV}$ Higgs data [34].
Also compared in the present paper is the Bayesian evidence of the $\mu\neq 0$
and $\mu=0$ cases (note that the $\mu=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 $\mu\neq 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 $\lambda$ and $\kappa$ 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 $\mu$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 $\mu$ contributes to the Higgs and
neutralino mass matrices, and, consequently, the DM physics differs
significantly.

Figure 1: Two-dimensional profile likelihood of function $\mathcal{L}_{\rm
DM}$ in Eq. (3.3) for $h\equiv h_{1}$ scenario, which is projected onto
$|\kappa|-\lambda$ and $|\mu_{tot}|-m_{\tilde{\chi}_{1}^{0}}$ planes of
$\mu$NMSSM. Since $\chi^{2}_{\rm Best}\simeq 0$ for the best point (marked
with a star symbol in figure), the boundaries for $1\sigma$ (red solid line)
and $2\sigma$ (white solid line) confidence intervals correspond to
$\chi^{2}\simeq 2.3$ and $\chi^{2}\simeq 6.18$, respectively. Dotted-dashed
line represents relation $2{\kappa}=\lambda$. Points around
$m_{\tilde{\chi}_{1}^{0}}\simeq|\mu_{tot}|$ on right-hand panel obtained
correct abundance by co-annihilation. This figure shows $\mu$NMSSM’s parameter
region that can well explain DM experiments.

Figure 2: Same as Figure 1, but for profile likelihood projected onto
$\sigma_{\tilde{\chi}_{1}^{0}-p}^{\rm SI}-m_{\tilde{\chi}_{1}^{0}}$ and
$\sigma_{\tilde{\chi}_{1}^{0}-n}^{\rm SD}-m_{\tilde{\chi}_{1}^{0}}$ plane. Red
curve denotes latest XENON-1T’s bound on SI and SD cross-section.

###  3.3 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

$\displaystyle\mathcal{L}_{\rm
DM}(\Theta_{A},\Theta_{B})=\mathop{\max}_{\Theta_{1},\cdots,\Theta_{A-1},\Theta_{A+1},\cdots,\Theta_{B-1},\Theta_{B+1},\cdots}\left\\{\mathcal{L}_{\Omega
h^{2}}\times\mathcal{L}_{DD}\times\mathcal{L}_{IDD}\right\\}(\Theta),$ (3.3)

where $\Theta=(\Theta_{1}\equiv\lambda,\Theta_{2}\equiv\kappa,\cdots)$ are the
input parameters of $\mathcal{L}_{i}$ ($i=\Omega h^{2},\cdots$), and the
maximization is realized by varying the parameters other than $\Theta_{A}$ and
$\Theta_{B}$. The $1\sigma$ ($2\sigma$) confidence interval on the
$\Theta_{A}-\Theta_{B}$ plane is determined by the criteria of
$(\delta\chi^{2}\equiv\chi^{2}-\chi^{2}_{Best})\leq 2.3$ ($\delta\chi^{2}\leq
6.18$), where $\chi^{2}\equiv-2\ln{\mathcal{L}_{\rm
DM}}(\Theta_{A},\Theta_{B})$ and $\chi^{2}_{Best}$ denotes the best point’s
$\chi^{2}$. Two subtleties about $\mathcal{L}_{\rm DM}(\Theta_{A},\Theta_{B})$
must be clarified. One is that it is completely different from the PL of the
total likelihood function $\mathcal{L}$ in Eq. (3.1). Specifically, the
$\chi^{2}$ of $\mathcal{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 $\mathcal{L}(\Theta_{A},\Theta_{B})$ is roughly
equivalent to performing a Higgs data fit on the $\Theta_{A}-\Theta_{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 $\mathcal{L}$ strongly disfavors these samples, they should be
neglected. Therefore, only the following samples are considered in studying
$\mathcal{L}_{\rm DM}(\Theta_{A},\Theta_{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.

  * •

It accounts for the measurement of $Br(B\to X_{s}\gamma)$ and
$Br(B_{s}^{0}\to\mu^{+}\mu^{-})$ at $2\sigma$ level [66].

Figure 3: Same as Figure 1, but for $h\equiv h_{2}$ scenario.

Figure 4: Same as Figure 2, but for $h\equiv h_{2}$ scenario.

Using the refined samples, different 2D PL maps are shown in Figures 1 to 4.
Figure 1 plots the credible interval (CI) on $|\kappa|-\lambda$ and
$|\mu_{tot}|-m_{\tilde{\chi}_{1}^{0}}$ plane for the $h\equiv h_{1}$ scenario.
This plot indicates that the theory can explain the DM experiments quite well
when $\lambda\lesssim 0.5$, $0.05\lesssim|\kappa|\lesssim 0.5$,
$m_{\tilde{\chi}_{1}^{0}}\gtrsim 140~{}{\rm GeV}$, and $|\mu_{tot}|\gtrsim
160~{}{\rm GeV}$. The best point locates at approximately $\lambda=0.05$,
$\kappa=0.28$, $m_{\tilde{\chi}_{1}^{0}}=350~{}{\rm GeV}$, and
$\mu_{tot}=450~{}{\rm GeV}$. Additionally, $2|\kappa|/\lambda$ can be much
larger than 1 to keep $\tilde{\chi}_{1}^{0}$ Singlino-dominated, which is one
distinct feature of the $\mu$NMSSM. Figure 2 shows the $1\sigma$ and $2\sigma$
CIs on the $\sigma_{\tilde{\chi}_{1}^{0}-p}^{\rm SI}-m_{\tilde{\chi}_{1}^{0}}$
and $\sigma_{\tilde{\chi}_{1}^{0}-n}^{\rm SD}-m_{\tilde{\chi}_{1}^{0}}$
planes. This plot indicates that the SI and SD cross-sections can be as low as
$10^{-51}~{}{\rm cm^{2}}$ and $10^{-46}~{}{\rm cm^{2}}$, respectively, over a
broad mass range of $\tilde{\chi}_{1}^{0}$. This conclusion comes from the
fact that $\lambda$ can be small and the ratio
$m_{\tilde{\chi}_{1}^{0}}/\mu_{tot}$ in Eqs.(2.32)-(2.32) can be much less
than 1. Figures 3 and 4 are similar to Figures 1 and 2, except for the
$h\equiv h_{2}$ scenario. Compared with the $h\equiv h_{1}$ scenario, this
scenario has the following different characteristics.

  * •

The CI on $|\kappa|-\lambda$ plane shrinks significantly.

  * •

$m_{\tilde{\chi}_{1}^{0}}$ as low as $100~{}{\rm GeV}$ can explain the
experiments.

  * •

The best point locates at approximately $\kappa=0.2$ and
$m_{\tilde{\chi}_{1}^{0}}=250~{}{\rm GeV}$.

In addition to the fact that the parameter space of the $h\equiv 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 $\tilde{\chi}_{1}^{0}$ can make the annihilation
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to 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.

###  3.4 Constraints from LHC search for electroweakinos

Figure 5: Constraint from latest ATLAS analysis of di-lepton signal at 13-TeV
LHC, which is shown on $|\mu_{tot}|-m_{\tilde{\chi}_{1}^{0}}$ plane. Analysis
is concentrated on compressed mass spectra [86] and excludes the yellow-shaded
region at $95\%$ confidence level when $\lambda=0.01$ and $\tan\beta=10$ are
fixed.

|  $\rm h\equiv h_{1}$ Scenario |  $\rm h\equiv h_{2}$ Scenario
---|---|---
| Point $\rm I$ | Point $\rm II$ | Point $\rm III$ | Point $\rm IV$
$\lambda$ | 0.026 | 0.227 | 0.038 | 0.077
$\kappa$ | 0.236 | 0.328 | -0.212 | 0.199
$\tan{\beta}$ | 9.750 | 6.698 | 9.454 | 8.092
$\mu~{}(\rm GeV)$ | 322.6 | 243.8 | 226.7 | 291.3
$\mu_{\rm tot}~{}(\rm GeV)$ | 335.0 | 331.9 | 248.0 | 338.2
$A_{\kappa}~{}(\rm GeV)$ | -310.5 | -238.1 | 21.58 | -415.8
$A_{t}~{}(\rm GeV)$ | 4244 | 4112 | 2999 | 2959
$m_{h_{1}}~{}(\rm GeV)$ | 125.2 | 125.2 | 125.1 | 93.67
$m_{h_{2}}~{}(\rm GeV)$ | 134.8 | 203.9 | 233.3 | 125.2
$m_{h_{3}}~{}(\rm GeV)$ | 835.3 | 1273 | 794.6 | 1053
$m_{A_{s}}~{}(\rm GeV)$ | 326.9 | 310.6 | 86.62 | 388.4
$m_{\tilde{\chi}_{1}^{0}}~{}(\rm GeV)$ | 228.6 | 249.2 | 236.7 | 255.4
$m_{\tilde{\chi}_{2}^{0}}~{}(\rm GeV)$ | 343.5 | 343.2 | 254.5 | 386.7
$m_{\tilde{\chi}_{3}^{0}}~{}(\rm GeV)$ | 345.4 | 345.5 | 258.0 | 348.7
$m_{\tilde{\chi}_{1}^{\pm}}~{}(\rm GeV)$ | 344.6 | 340.7 | 255.7 | 387.8
$V_{h}^{\rm NSM},V_{h}^{\rm SM},V_{h}^{\rm S}$ | -0.20, -0.94, 0.30 | -0.30, -0.95, -0.06 | -0.21, -0.98, 0.02 | 0.01, 0.99, 0.17
$V_{h_{s}}^{\rm NSM},V_{h_{s}}^{\rm SM},V_{h_{s}}^{\rm S}$ | -0.05, -0.29, -0.96 | -0.05, 0.08, -0.99 | -0.02, -0.01, -0.99 | 0.01, -0.17, 0.99
$N_{13},N_{14},N_{15}$ | 0.01, -0.02, 0.99 | 0.14, -0.21, 0.97 | -0.17, -0.18, -0.97 | 0.04, -0.07, 0.99
$N_{23},N_{24},N_{25}$ | -0.71, 0.71, 0.02 | 0.71, 0.71, 0.05 | -0.71, 0.71, 0.05 | -0.71, 0.70, 0.08
$N_{33},N_{34},N_{35}$ | 0.71, 0.71, 0.01 | -0.69, 0.68, 0.25 | -0.69, -0.69, -0.25 | 0.71, 0.71, 0.02
$\Omega h^{2}$ | 0.125 | 0.125 | 0.131 | 0.124
$\sigma^{SI}_{\tilde{\chi}_{1}^{0}-p}(\rm cm^{2})$ | $2.83\times 10^{-48}$ | $2.59\times 10^{-47}$ | $1.01\times 10^{-47}$ | $6.36\times 10^{-48}$
$\sigma^{SD}_{\tilde{\chi}_{1}^{0}-n}(\rm cm^{2})$ | $2.47\times 10^{-46}$ | $1.95\times 10^{-41}$ | $5.49\times 10^{-43}$ | $2.50\times 10^{-44}$
Annihilation processes |  |
---
$\;\;\;71\%\rm\;\;\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;A_{s}$
$\;\;\;20\%\rm\;\;\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h\;\;\,A_{s}$
$\;\;7.4\%\rm\;\;\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;h_{s}$
$1.3\%\rm\;\;\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h\;h_{s}$
……
|
---
$34\%\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;t\;\bar{t}$
$\quad
33\%\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;A_{s}$
$\quad
27\%\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;h_{s}$
$\;\;1.4\%\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h\;A_{s}$
$\qquad\
1.1\%\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;W^{+}\;W^{-}$
……
|
---
$5.4\%\rm\quad\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\;\to\;d\;\bar{u}$
$5.4\%\rm\quad\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\;\to\;s\;\bar{c}$
$4.0\%\rm\quad\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\;\to\;b\;\bar{t}$
$4.0\%\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{\pm}\;\to\;s\;\bar{c}$
$4.0\%\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{\pm}\,\to\;d\;\bar{u}$
……
| $96\%\rm\quad\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;A_{s}$
---
$3.4\%\rm\;\;\;\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\;\to\;h_{s}\;h_{s}$
$\tilde{\chi}_{2}^{0}$ Decay channels | $100\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;Z$ | $100\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;Z$ |  |
---
$16\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;d\;\bar{d}\;(\,d=d,\;s\,)$
$12\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;u\;\bar{u}(\;u=u,\;c\;)$
$10\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;b\;\bar{b}$
$7.4\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;\nu\;\bar{\nu}\;(\,\nu=\nu_{e},\;\nu_{\mu},\;\nu_{\tau}\,)$
$3.7\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;l\;\bar{l}(\;l=e,\;\mu\;)$
$3.5\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;\tau\;\bar{\tau}\;$
| $57\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;Z$
---
$43\%\rm\quad\tilde{\chi}_{2}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;h_{s}$
$\tilde{\chi}_{3}^{0}$ Decay channels | $100\%\rm\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;Z$ | $100\%\rm\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;Z$ |  |
---
$11\%\rm\quad\;\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;d\;\bar{d}\;(\,d=d,\;s\,)$
$8.3\%\rm\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;u\;\bar{u}$
$7.9\%\rm\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;c\;\bar{c}$
$6.8\%\rm\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;b\;\bar{b}$
$5.6\%\rm\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{3}^{0}\;b\;\bar{b}$
$5.0\%\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;\nu\;\bar{\nu}\;(\,\nu=\nu_{e},\;\nu_{\mu},\;\nu_{\tau}\,)$
……
| $99.7\%\rm\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;Z$
---
$0.3\%\rm\quad\tilde{\chi}_{3}^{0}\;\to\;\tilde{\chi}_{1}^{0}\;h_{s}$
$\tilde{\chi}_{1}^{\pm}$ Decay channels | $~{}100\%\rm\quad\tilde{\chi}_{1}^{\pm}\;\to\;\tilde{\chi}_{1}^{0}\;W^{\pm}$ | $~{}100\%\rm\quad\tilde{\chi}_{1}^{\pm}\;\to\;\tilde{\chi}_{1}^{0}\;W^{\pm}$ |  |
---
$~{}33\%\rm\quad\tilde{\chi}_{2}^{\pm}\;\to\;\tilde{\chi}_{1}^{0}\;u\;\bar{d}(\,u\;\bar{d}=u\;\bar{d},\;c\;\bar{s}\,)$
$~{}11\%\rm\quad\tilde{\chi}_{1}^{\pm}\;\to\;\tilde{\chi}_{1}^{0}\;\nu\;\bar{l}\;(\,\nu\;\bar{l}=\nu_{e}\;\bar{e},\;\nu_{\nu}\;\bar{\mu},\;\nu_{\tau}\;\bar{\tau}\,)$
$~{}100\%\rm\quad\tilde{\chi}_{1}^{\pm}\;\to\;\tilde{\chi}_{1}^{0}\;W^{\pm}$

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 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 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_{\tilde{\chi}_{1}^{0}}$ in the simplified model.
Recently, the ATLAS collaboration analyzed $139{\rm fb^{-1}}$ proton-proton
collision data collected at the LHC with $\sqrt{s}=13{\rm TeV}$ and concluded
that the LHC had already explored the region with the Wino mass up to
approximately $700{\rm GeV}$ and $m_{\tilde{\chi}_{1}^{0}}$ up to $300{\rm
GeV}$ [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_{\tilde{\chi}_{1}^{0}}\lesssim 100~{}{\rm GeV}$555A similar conclusion was
obtained as shown in Figure 3 of [10], where the results on the upper panel of
Figure 8 in [85] are re-interpreted in terms of the Higgsino pair production
process. The latter result was obtained by a combined analysis of the multi-
lepton signal by the CMS collaboration. Its exclusion capability on the Wino
mass versus DM mass plane is roughly identical to that of the analyses in [82,
83, 84]. . As a specific application, the $h$-funnel region in the $h\equiv
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.9{\rm fb^{-1}}$
data [85].

Above discussion reveals that the analyses in [82, 83, 84, 85] have little
restriction on the theory since $\tilde{\chi}_{1}^{0}$ 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 $139{\rm fb^{-1}}$
data of the LHC. Events with missing transverse momentum and two same-flavor,
oppositely charged, low-transverse-momentum leptons are selected, and are
further categorized by the presence of hadronic activity from initial-state
radiation. The analysis was reproduced in the present work by using the
simulation tools<EMAIL_ADDRESS>[87, 88] to generate the parton-level
events of the processes
$pp\to\tilde{\chi}_{2}^{0}\tilde{\chi}_{3}^{0}(j),\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}(j),\tilde{\chi}_{3}^{0}\tilde{\chi}_{1}^{\pm}(j)$,
and $\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}(j)$; Pythia-8.2 [89] for
parton fragmentation and hadronization; Delphes-3.4.2 [90] for fast simulation
of the performance of the ATLAS detector; and CheckMATE-2.0.26 [91, 92, 93] to
implement the analysis’ cut selections. The validation of the method is
provided in our recent work [14].

In Figure 5, the excluded region is shown on the
$|\mu_{tot}|-m_{\tilde{\chi}_{1}^{0}}$ plane. $\lambda=0.01$, $\tan\beta=10$,
and $M_{1}=M_{2}=3~{}{\rm TeV}$ are fixed, and the procedure depicted in [59]
is followed to plot the figure. The result indicates that the analysis can
explore $m_{\tilde{\chi}_{1}^{0}}$ up to $230~{}{\rm 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 $\tilde{\chi}_{1}^{0}$ 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\equiv h_{1}$ scenario and the last one is from the
$h\equiv h_{2}$ scenario. They all predict $m_{\tilde{\chi}_{1}^{0}}\gtrsim
230{\rm GeV}$ 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 $\tilde{\chi}_{i}^{0}$ ($i=2,3$) and
$\tilde{\chi}_{1}^{0}$ is larger than any of the neutral Higgs boson masses,
$\tilde{\chi}_{i}^{0}$ will decay with sizable branching ratios into the Higgs
bosons. Otherwise, it will decay into a real or virtual $Z$ boson.

  * •

Although the process $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$
is mainly responsible for the abundance in most cases, its role diminishes
when $m_{h_{s}}+m_{A_{s}}$ approaches $2m_{\tilde{\chi}_{1}^{0}}$ from below.
This situation is demonstrated by Points 1 and 2.

  * •

Since $V_{h}^{S}=0.3$ is sizable, $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
hA_{s}$ for Point 1 also plays a crucial role in determining the abundance.

  * •

Since $\lambda=0.227$ is sizable, $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$
and $\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ for Point 2 are relatively
large.

  * •

Since $m_{\tilde{\chi}_{1}^{0}}/\mu_{tot}=0.95$ approaches 1, $\sigma^{\rm
SI}_{\tilde{\chi}_{1}^{0}-p}$ and $\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$
for Point 3 are relatively large even though $\lambda$ is small.

  * •

Point 4 predicts a significant cancellation between the $h$\- and
$h_{s}$-mediated contributions to $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$.

##  4 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
$\mu\hat{H}_{u}\cdot\hat{H}_{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: $\lambda$, $\mu_{tot}$,
$m_{\tilde{\chi}_{1}^{0}}$, and $\kappa$. The first three parameters strongly
influence the DM-nucleon scattering rate, while $\kappa$ 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
$\kappa$’s value. Its scattering with nucleons is suppressed when $\lambda
v/\mu_{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 $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ in most cases, and
the SI and SD cross-sections may be as low as $10^{-51}{\rm cm^{2}}$ and
$10^{-46}{\rm 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_{\tilde{\chi}_{1}^{0}}$
and $\kappa$ are no longer independent once one takes $\lambda$ and
$\mu_{tot}$ as inputs to study the DM’s property. Second, $|\kappa|$ must be
less than $\lambda/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 $\mu$-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 $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to 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
$\lambda\simeq 2|\kappa|$ with $\lambda\lesssim 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.

## Acknowledgements

We thank Di Zhang for useful discussions about the LHC search for
supersymmetry. This work is supported by the National Natural Science
Foundation of China (NNSFC) under Grant Nos. 11575053 and 12075076.

## References

  * [1] Planck collaboration, P. A. R. Ade et al., _Planck 2015 results. XIII. Cosmological parameters_ , _Astron. Astrophys._ 594 (2016) A13, arXiv:1502.01589 [hep-ph].
  * [2] N. Aghanim et al. [Planck Collaboration], _Planck 2018 results. VI. Cosmological parameters_ , arXiv:1807.06209 [hep-ph].
  * [3] G. Jungman, M. Kamionkowski and K. Griest, _Supersymmetric dark matter_ , _Physics Reports_ 333 (1995) 167.
  * [4] G. Bertone, D. Hooper and J. Silk, _Particle dark matter: Evidence, candidates and constraints_ , Phys. Rept. 405, 279-390 (2005) doi:10.1016/j.physrep.2004.08.031, arXiv:hep-ph/0404175 [hep-ph].
  * [5] S. Baum, M. Carena, N. R. Shah and C. E. Wagner, _Higgs portals for thermal Dark Matter. EFT perspectives and the NMSSM_ , _JHEP_ 04 (2018) 069, arXiv:1712.09873 [hep-ph].
  * [6] XENON collaboration, E. Aprile et al., _Dark Matter Search Results from a One Ton-Year Exposure of XENON1T_ , _Phys. Rev. Lett._ 121 (2018) 111302, arXiv:1805.12562 [hep-ph].
  * [7] Q. Wang et al. [PandaX-II Collaboration], _Results of dark matter search using the full PandaX-II exposure_ , Chin. Phys. C 44, no. 12, 125001 (2020) arXiv:2007.15469 [astro-ph.CO].
  * [8] X. Cui et al. [PandaX-II Collaboration], _Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment_ , Phys. Rev. Lett. 119, no. 18, 181302 (2017) arXiv:1708.06917 [astro-ph.CO].
  * [9] J. Cao, Y. He, L. Shang, Y. Zhang and P. Zhu, _Current status of a natural NMSSM in light of LHC 13 TeV data and XENON-1T results_ , _Phys. Rev._ D99 (2019) 075020, arXiv:1810.09143[hep-ph].
  * [10] J. Cao, L. Meng, Y. Yue, H. Zhou and P. Zhu, _Suppressing the scattering of wimp dark matter and nucleons in supersymmetric theories_ , _Phys. Rev. D_ 101 (2020) 075003, arXiv:1910.14317 [hep-ph]
  * [11] E. Bagnaschi, K. Sakurai, M. Borsato, O. Buchmueller, M. Citron, J. C. Costa, A. De Roeck, M. J. Dolan, J. R. Ellis and H. Flächer, et al., _Likelihood Analysis of the pMSSM11 in Light of LHC 13-TeV Data_ ,Eur. Phys. J. C 78, no.3, 256 (2018) doi:10.1140/epjc/s10052-018-5697-0, arXiv:1710.11091 [hep-ph].
  * [12] U. Ellwanger, C. Hugonie and A. M. Teixeira, _The Next-to-Minimal Supersymmetric Standard Model_ , _Phys. Rept._ 496 (2010) 1 arXiv:0910.1785 [hep-ph].
  * [13] J. Cao, Y. He, L. Shang, W. Su, P. Wu and Y. Zhang, _Strong constraints of LUX-2016 results on the natural NMSSM_ , JHEP 10, 136 (2016) doi:10.1007/JHEP10(2016)136, arXiv:1609.00204 [hep-ph].
  * [14] H. Zhou, J. Cao, J. Lian, Y. Yue and D. Zhang, _Singlino-dominated dark matter in $Z_{3}$-NMSSM_, arXiv:2102.05309 [hep-ph].
  * [15] J. Cao, L. Shang, P. Wu, J. M. Yang and Y. Zhang, _Interpreting the galactic center gamma-ray excess in the NMSSM_ , _JHEP_ 10 (2015) 030, arXiv:1506.06471 [hep-ph].
  * [16] M. Pospelov, A. Ritz and M. B. Voloshin, _Secluded WIMP Dark Matter_ , _Phys. Lett. B_ 662 (2008) 53 arXiv:0711.4866 [hep-ph].
  * [17] D. Das, U. Ellwanger and A. M. Teixeira, _Modified Signals for Supersymmetry in the NMSSM with a Singlino-like LSP_ , JHEP 1204, 067 (2012) arXiv:1202.5244 [hep-ph].
  * [18] U. Ellwanger and A. M. Teixeira, _NMSSM with a singlino LSP: possible challenges for searches for supersymmetry at the LHC_ , JHEP 1410, 113 (2014) arXiv:1406.7221 [hep-ph].
  * [19] U. Ellwanger, _Present Status and Future Tests of the Higgsino-Singlino Sector in the NMSSM_ , JHEP 1702, 051 (2017) arXiv:1612.06574 [hep-ph].
  * [20] U. Ellwanger and C. Hugonie, _The higgsino-singlino sector of the NMSSM: combined constraints from dark matter and the LHC_ , Eur. Phys. J. C 78, no. 9, 735 (2018) arXiv:1806.09478 [hep-ph].
  * [21] Q. F. Xiang, X. J. Bi, P. F. Yin and Z. H. Yu, _Searching for Singlino-Higgsino Dark Matter in the NMSSM_ , Phys. Rev. D 94 (2016) no.5, 055031 arXiv:1606.02149 [hep-ph].
  * [22] W. Abdallah, A. Chatterjee and A. Datta, _Revisiting singlino dark matter of the natural $Z_{3}$-symmetric NMSSM in the light of LHC_, JHEP 1909, 095 (2019) arXiv:1907.06270 [hep-ph].
  * [23] M. Guchait and A. Roy, Phys. Rev. D 102, no.7, 075023 (2020) doi:10.1103/PhysRevD.102.075023, arXiv:2005.05190 [hep-ph].
  * [24] M. L. Hyun, R. Stuart, R. Michael, G. R. Graham, S. Roland, S.-H. Kai et al., _A unique $Z^{R}_{4}$ symmetry for the MSSM_, arXiv:1009.0905 [hep-ph].
  * [25] M. L. Hyun, R. Stuart, R. Michael, G. R. Graham, S. Roland, S.-H. Kai et al., _Discrete R symmetries for the MSSM and its singlet extensions_ , arXiv:1102.3595 [hep-ph].
  * [26] G. R. Graham, S.-H. Kai and S. Florian, _The fine-tuning of the generalised NMSSM_ , arXiv:1108.1284 [hep-ph].
  * [27] G. R. Graham, S.-H. Kai and S. Florian, _The generalised NMSSM at one loop: fine tuning and phenomenology_ , arXiv:1205.1509 [hep-ph].
  * [28] F. Sergio, K. Renata, L. Andrei, M. Alessio and V. P. Antoine, _Jordan Frame Supergravity and Inflation in NMSSM_ , arXiv:1004.0712 [hep-ph].
  * [29] F. Sergio, K. Renata, L. Andrei, M. Alessio and V. P. Antoine, _Superconformal Symmetry, NMSSM, and Inflation_ , arXiv:1008.2942 [hep-ph].
  * [30] B. E. Martin and R. T. J. D, _Inflation with Non-minimal Gravitational Couplings in Supergravity_ , arXiv:0912.2718 [hep-ph].
  * [31] W. G. Hollik, S. Liebler, G. Moortgat-Pick, S. Pa?ehr and G. Weiglein, _Phenomenology of the inflation-inspired NMSSM at the electroweak scale_ , _Eur. Phys. J. C_ 79 (2019) 75, arXiv:1809.07371 [hep-ph].
  * [32] D. J. Miller, R. Nevzorov and P. M. Zerwas, _The higgs sector of the next-to-minimal supersymmetric standard model_ , _Nuclear Physics B_ 681 (2004) 3, arXiv:hep-ph/0304049
  * [33] W. G. Hollik, C. Li, G. Moortgat-Pick and S. Paasch, _Phenomenology of a Supersymmetric Model Inspired by Inflation_ , arXiv:2004.14852 [hep-ph].
  * [34] J. J. Cao, Z. X. Heng, J. M. Yang, Y. M. Zhang and J. Y. Zhu, _A SM-like Higgs near 125 GeV in low energy SUSY: a comparative study for MSSM and NMSSM_ , JHEP 1203, 086 (2012) doi:10.1007/JHEP03(2012)086, arXiv:1202.5821 [hep-ph].
  * [35] L. J. Hall, D. Pinner and J. T. Ruderman, _A Natural SUSY Higgs Near 126 GeV_ , JHEP 1204, 131 (2012) doi:10.1007/JHEP04(2012)131, arXiv:1112.2703 [hep-ph].
  * [36] U. Ellwanger, _A Higgs boson near 125 GeV with enhanced di-photon signal in the NMSSM_ , JHEP 1203, 044 (2012) doi:10.1007/JHEP03(2012)044, arXiv:1112.3548 [hep-ph].
  * [37] S. Baum, N. R. Shah and K. Freese, JHEP 04, 011 (2019) doi:10.1007/JHEP04(2019)011, arXiv:1901.02332 [hep-ph].
  * [38] M. Carena, I. Low, N. R. Shah and C. E. M. Wagner, JHEP 04, 015 (2014) doi:10.1007/JHEP04(2014)015, arXiv:1310.02248 [hep-ph]. [arXiv:1310.2248 [hep-ph]].
  * [39] M. Carena, H. E. Haber, I. Low, N. R. Shah and C. E. M. Wagner, Phys. Rev. D 93, no.3, 035013 (2016) doi:10.1103/PhysRevD.93.035013, arXiv:1510.09137 [hep-ph].
  * [40] P. Huang and C. E. M. Wagner, Phys. Rev. D 90, no.1, 015018 (2014) doi:10.1103/PhysRevD.90.015018, arXiv:1404.0392 [hep-ph].
  * [41] G. Aad et al. [ATLAS], Phys. Rev. Lett. 125, no.5, 051801 (2020) doi:10.1103/PhysRevLett.125.051801, arXiv:2002.12223 [hep-ex].
  * [42] C. Cheung, M. Papucci, D. Sanford, N. R. Shah and K. M. Zurek, _NMSSM Interpretation of the Galactic Center Excess_ , Phys. Rev. D 90, no.7, 075011 (2014) doi:10.1103/PhysRevD.90.075011, arXiv:1406.6372 [hep-ph].
  * [43] M. Badziak, M. Olechowski and P. Szczerbiak, _Blind spots for neutralino dark matter in the NMSSM_ , _JHEP_ 03 (2016) 179, arXiv:1512.02472 [hep-ph].
  * [44] M. Badziak, M. Olechowski and P. Szczerbiak, _Spin-dependent constraints on blind spots for thermal singlino-higgsino dark matter with(out) light singlets_ , _Journal of High Energy Physics_ 2017 (2017) ,arXiv:1705.00227 [hep-ph]
  * [45] K. Griest and D. Seckel, _Three exceptions in the calculation of relic abundances_ , _Physical Review D_ 43 (1991) 3191.
  * [46] M. J. Baker et al., _The Coannihilation Codex_ , JHEP 1512, 120 (2015) , arXiv:1510.03434 [hep-ph].
  * [47] K. Griest, M. Kamionkowski and M. S. Turner, _Supersymmetric Dark Matter Above the W Mass_ Phys. Rev. D 41, 3565 (1990). .
  * [48] M. Badziak, M. Olechowski and P. Szczerbiak, _Blind spots for neutralinos in nmssm with light singlet scalar_ , arXiv:1601.00768 [hep-ph].
  * [49] A. Pierce, N. R. Shah and K. Freese, _Neutralino dark matter with light staus_ , 2013, arXiv:1309.7351 [hep-ph].
  * [50] M. Drees and M. Nojiri, _Neutralino - nucleon scattering revisited_ , Phys. Rev. D 48, 3483 (1993) , arXiv:hep-ph/9307208.
  * [51] M. Drees and M. M. Nojiri, _New contributions to coherent neutralino - nucleus scattering_ , Phys. Rev. D 47, 4226 (1993) , arXiv:hep-ph/9210272.
  * [52] G. Bélanger, F. Boudjema, A. Pukhov and A. Semenov, _Dark matter direct detection rate in a generic model with micromegas2.2_ , _Computer Physics Communications_ 180 (2009) 747–767.
  * [53] J. M. Alarcon, J. Martin Camalich and J. A. Oller, _Phys. Rev. D_ 85, 051503 (2014), arXiv:1110.3797[hep-ph].
  * [54] J. M. Alarcon, L. S. Geng, J. Martin Camalich and J. A. Oller, _Phys. Lett. B_ 730, 342-346 (2014), arXiv:1209.2870[hep-ph].
  * [55] G. Belanger, F. Boudjema, C. Hugonie, A. Pukhov and A. Semenov, _Relic density of dark matter in the NMSSM_ , _JCAP_ 09 (2005) 001, arXiv:hep-ph/0505142.
  * [56] F. Feroz, M. P. Hobson and M. Bridges, _MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics_ , Mon. Not. Roy. Astron. Soc. 398, 1601 (2009) , arXiv:0809.3437 [astro-ph].
  * [57] F. Feroz, M. P. Hobson, E. Cameron and A. N. Pettitt, _Importance Nested Sampling and the MultiNest Algorithm_ , arXiv:1306.2144 [astro-ph.IM].
  * [58] A. Fowlie and M. H. Bardsley, _Superplot: a graphical interface for plotting and analysing MultiNest output_ , Eur. Phys. J. Plus 131, no. 11, 391 (2016) , arXiv:1603.00555 [physics.data-an].
  * [59] J. Cao, J. Li, Y. Pan, L. Shang, Y. Yue and D. Zhang, _Bayesian analysis of sneutrino dark matter in the NMSSM with a type-I seesaw mechanism_ , _Phys. Rev. D_ 99 (2019) 115033, arXiv:1807.03762[hep-ph].
  * [60] J. Cao, Y. He, Y. Pan, Y. Yue, H. Zhou and P. Zhu, _Impact of leptonic unitarity and dark matter direct detection experiments on the NMSSM with inverse seesaw mechanism_ , JHEP 2012, 023 (2020) , arXiv:1903.01124 [hep-ph].
  * [61] XENON collaboration, E. Aprile et al., _Constraining the spin-dependent WIMP-nucleon cross sections with XENON1T_ , _Phys. Rev. Lett._ 122 (2019) 141301, arXiv:1902.03234 [hep-ph].
  * [62] For Fermi-LAT data analysis overview, see https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/overview.html.
  * [63] M. Ackermann et al. [Fermi-LAT Collaboration], _Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data_ , Phys. Rev. Lett. 115, no. 23, 231301 (2015) , arXiv:1503.02641 [astro-ph.HE].
  * [64] P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak and G. Weiglein, _Probing the Standard Model with Higgs signal rates from the Tevatron, the LHC and a future ILC_ , JHEP 1411, 039 (2014) , arXiv:1403.1582 [hep-ph].
  * [65] P. Bechtle, S. Heinemeyer, O. Stal, T. Stefaniak and G. Weiglein, _Applying Exclusion Likelihoods from LHC Searches to Extended Higgs Sectors_ , Eur. Phys. J. C 75, no. 9, 421 (2015) , arXiv:1507.06706 [hep-ph].
  * [66] M. Tanabashi et al. [Particle Data Group], _Review of Particle Physics_ , Phys. Rev. D 98 (2018) no.3, 030001
  * [67] F. Staub, _SARAH_ , arXiv:0806.0538 [hep-ph].
  * [68] F. Staub, _SARAH 3.2: Dirac Gauginos, UFO output, and more_ ,_Computer Physics Communications 184 (2013) pp. 1792-1809_ (2013) , arXiv:1207.0906 [hep-ph].
  * [69] F. Staub, _SARAH4: A tool for (not only SUSY) model builders_ ,_Computer Physics Communications 185 (2014) pp. 1773-1790_ (2014) , arXiv:1309.7223 [hep-ph].
  * [70] F. Staub, _Exploring new models in all detail with sarah_ ,_Advances in High Energy Physics_ (2015) , arXiv:1503.04200 [hep-ph].
  * [71] W. Porod, _SPheno, a program for calculating supersymmetric spectra, SUSY particle decays and SUSY particle production at e+ e- colliders_ ,_Comput. Phys. Commun._ 153 (2003) 275, arXiv:hep-ph/0301101.
  * [72] W. Porod and F. Staub, _Spheno 3.1: extensions including flavour, cp-phases and models beyond the mssm_ ,_Computer Physics Communications_ 183 (2012) 2458–2469, arXiv:1104.1573 [hep-ph].
  * [73] W. Porod, F. Staub and A. Vicente, _A Flavor Kit for BSM models_ , Eur. Phys. J. C 74, no. 8, 2992 (2014) doi:10.1140/epjc/s10052-014-2992-2, arXiv:1405.1434 [hep-ph].
  * [74] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _MicrOMEGAs: A Program for calculating the relic density in the MSSM_ ,_Comput. Phys. Commun._ 149 (2002) 103, arXiv:hep-ph/0112278.
  * [75] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _MicrOMEGAs 2.0: A Program to calculate the relic density of dark matter in a generic model_ ,_Comput. Phys. Commun._ 176 (2007) 367, arXiv:hep-ph/0607059.
  * [76] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _micromegas : a tool for dark matter studies_ , 2010, arXiv:1005.4133 [hep-ph]
  * [77] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _micrOMEGAs $\\_$3: A program for calculating dark matter observables_,_Comput. Phys. Commun._ 185 (2014) 960, arXiv:1305.0237[hep-ph].
  * [78] D. Barducci, G. Belanger, J. Bernon, F. Boudjema, J. D. Silva, S. Kraml et al., _Collider limits on new physics within micromegas4.3_ , 2017, arXiv:1606.03834 [hep-ph].
  * [79] P. Gregory, Bayesian Logical Data Analysis for the Physical Sciences. Cambridge University Press, 2005.
  * [80] H. Jeffreys (1961). The Theory of Probability (3rd ed.). Oxford. p. 432.
  * [81] F. Feroz, B. C. Allanach, M. Hobson, S. S. AbdusSalam, R. Trotta and A. M. Weber, _Bayesian Selection of sign(mu) within mSUGRA in Global Fits Including WMAP5 Results_ , JHEP 0810, 064 (2008),arXiv:0807.4512 [hep-ph].
  * [82] G. Aad et al. [ATLAS Collaboration], _Search for electroweak production of charginos and sleptons decaying into final states with two leptons and missing transverse momentum in $\sqrt{s}=13$ TeV $pp$ collisions using the ATLAS detector_, Eur. Phys. J. C 80, no. 2, 123 (2020) [arXiv:1908.08215 [hep-ex]].
  * [83] G. Aad et al. [ATLAS Collaboration], _Search for chargino-neutralino production with mass splittings near the electroweak scale in three-lepton final states in $\sqrt{s}$=13 TeV $pp$ collisions with the ATLAS detector_, Phys. Rev. D 101, no. 7, 072001 (2020) [arXiv:1912.08479 [hep-ex]].
  * [84] The ATLAS collaboration [ATLAS Collaboration], _Search for chargino-neutralino pair production in final states with three leptons and missing transverse momentum in $\sqrt{s}$ = 13 TeV p-p collisions with the ATLAS detector_, ATLAS-CONF-2020-015.
  * [85] A. M. Sirunyan et al. [CMS Collaboration], _Combined search for electroweak production of charginos and neutralinos in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, JHEP 1803 (2018) 160 [arXiv:1801.03957 [hep-ex]].
  * [86] G. Aad et al. [ATLAS Collaboration], _Searches for electroweak production of supersymmetric particles with compressed mass spectra in $\sqrt{s}=$ 13 TeV $pp$ collisions with the ATLAS detector_, Phys. Rev. D 101, no. 5, 052005 (2020) [arXiv:1911.12606 [hep-ex]].
  * [87] J. Alwall et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , JHEP 1407, 079 (2014) [arXiv:1405.0301 [hep-ph]].
  * [88] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, _MadGraph 5 : Going Beyond_ , JHEP 1106, 128 (2011) [arXiv:1106.0522 [hep-ph]].
  * [89] T. Sjostrand, S. Mrenna and P. Z. Skands, _PYTHIA 6.4 Physics and Manual_ , JHEP 0605, 026 (2006) [hep-ph/0603175].
  * [90] J. de Favereau et al. [DELPHES 3 Collaboration], _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , JHEP 1402, 057 (2014) [arXiv:1307.6346 [hep-ex]].
  * [91] D. Dercks, N. Desai, J. S. Kim, K. Rolbiecki, J. Tattersall and T. Weber, _CheckMATE 2: From the model to the limit_ , arXiv:1611.09856 [hep-ph].
  * [92] J. S. Kim, D. Schmeier, J. Tattersall and K. Rolbiecki, _A framework to create customised LHC analyses within CheckMATE_ , Comput. Phys. Commun. 196, 535 (2015) [arXiv:1503.01123 [hep-ph]].
  * [93] M. Drees, H. Dreiner, D. Schmeier, J. Tattersall and J. S. Kim, _CheckMATE: Confronting your Favourite New Physics Model with LHC Data_ , Comput. Phys. Commun. 187, 227 (2015) [arXiv:1312.2591 [hep-ph]].