# Dark Matter Physics in General NMSSM

Lei Menga, Junjie Caoa,∗††*Corresponding author., and Shenshen Yanga

###### Abstract

In the General Next-to-Minimal Supersymmetric Standard Model (GNMSSM), singlet
particles may form a secluded sector of dark matter (DM), in which Singlino-
like DM could achieve the observed relic abundance through various channels
such as $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}h_{s},A_{s}A_{s},h_{s}A_{s}$, where $h_{s}$ and $A_{s}$ represent
singlet-dominated CP-even and CP-odd Higgs bosons. We provide analytical
formulas for both the spin-independent and spin-dependent cross sections of
Singlino DM scattering with nucleons, illustrating their dependence on the
model’s parameters in a clear manner. We also present analytic expressions for
the annihilation cross sections of these three important channels. Based on
these preparations, we conducted Bayesian analyses of the GNMSSM and concluded
that the theory significantly favored Singlino-dominated DM over Bino-like DM
across a much broader range of parameters. The combined results from our
numerical analyses and the formulas distinctly highlight crucial aspects of DM
physics within the GNMSSM.

##  1 Introduction

Although the presence of non-baryonic dark matter (DM) in the cosmos is
compelling from the perspective of astrophysics and cosmological observations
and constitutes a pivotal constituent of the prevailing standard model of
cosmology, its precise nature remains undetermined. Weakly Interacting Massive
Particles (WIMPs) are taken as the favored Cold Dark Matter (CDM) candidates
in the standard $\Lambda$CDM model. These large-mass, stable particles are
favored because they can be produced in thermal equilibrium conditions, and as
a consequence, the DM relic density measured by the Planck experiment can be
naturally predicted, referred to as the ’WIMP Miracle’ [1, 2, 3, 4]. The WIMPs
are typically assumed to couple to Standard Model (SM) particles through weak
interactions, resulting in spin-independent (SI) and spin-dependent (SD)
scattering cross sections with nucleons around $10^{-45}$ and $10^{-39}\ {\rm
cm}^{2}$, respectively [5]. This coincidence has spurred extensive direct [6,
7, 8, 9], indirect [10, 11, 12], and collider searches for WIMPs [13, 14]. In
particular, recent direct detection experiments like LUX-ZEPLIN (LZ) have
achieved unprecedented precision in probing both SI and SD DM-nucleon
scatterings but found no evidence of DM, restricting the scattering cross
sections to be less than the order of $10^{-47}$ and $10^{-42}\ {\rm cm}^{2}$,
respectively [9]. This implies that the interaction between DM and nucleons is
at most feeble, posing significant experimental challenges to simple WIMP
theories [15]. Consequently, there is growing interest in exploring scenarios
that go beyond the traditional thermal WIMP paradigm.

WIMPs can naturally arise in the supersymmetric extension of the SM in
particle physics. Furthermore, supersymmetric theories possess numerous
intriguing features, such as the absence of quadratic divergences, gauge
coupling unification, and distinctive signatures of supersymmetric particles
produced in present and future colliders. In the framework of the Minimal
Supersymmetric Standard Model (MSSM), recognized as one of the most economical
supersymmetric theories, the neutral electroweakino sector encompasses the
mixing of Bino (denoted as $\tilde{B}$), Wino (referred to as $\tilde{W}$),
and neutral Higgsinos (designated as $\tilde{H}_{u}$ and $\tilde{H}_{d}$),
giving rise to the formation of four distinct neutralinos [16, 17, 18]. Among
these, the lightest mass eigenstate, denoted as $\tilde{\chi}_{1}^{0}$,
usually corresponds to the lightest supersymmetric particle (LSP) and acts as
a feasible candidate of DM. It prefers to be $\tilde{B}$-dominated and co-
annihilate with the $\tilde{W}$-like electroweakinos to achieve the measured
relic abundance. As indicated by the research in Ref. [19], the LZ experiment
alone imposes a lower bound on the Higgsino mass, $\mu\gtrsim 380$ GeV, in the
absence of strong cancellations between different contributions to the SI
scattering. Although such a large $\mu$ may be naturally generated by the
well-known Giudice–Masiero mechanism in the gravity-mediated supersymmetry
breaking scenario [20], it leads to severe fine-tuning problems in precisely
predicting the mass of the $Z$ boson, after considering the Large Hadron
Collider (LHC) Higgs discovery and non-observations of any DM and
supersymmetric particle signals in experiments when the theory runs down from
an infrared high-energy scale to the electroweak scale [21, 22, 23].

The dilemma faced by the MSSM inspired us to extend this model. The Next-to-
Minimal Supersymmetric Standard Model with a $\mathbb{Z}_{3}$ symmetry
($\mathbb{Z}_{3}$-NMSSM) introduces a gauge singlet superfield $\hat{S}$ on
top of the MSSM [24, 25]. In this framework, the $\mu$-parameter of the MSSM
is dynamically generated once the scalar component of the superfield $\hat{S}$
acquires a vacuum expectation value (vev) below approximately 1 TeV. This
places it naturally at the electroweak scale. The neutralino sector includes
an additional fermionic partner of $\hat{S}$, known as the Singlino. Both the
Bino-dominated (which is the case in most physical scenarios) and Singlino-
dominated neutralinos are viable DM candidates [26, 27, 28, 5, 29, 30, 31, 32,
33, 34, 15, 35, 36, 37, 38, 39, 40]. Primarily, the Bino-dominated DM
candidate differs from the MSSM given the former’s ability to co-annihilate
with a Singlino-dominated neutralino in achieving the measured abundance.
Nonetheless, this scenario is confined to an exceedingly limited parameter
space defined by specific conditions: $|2\kappa\mu/\lambda|\simeq|M_{1}|$,
moderately large $\lambda$ and $\kappa$, and $\mu$ values exceeding 300 GeV
[5, 15]. This space predicts a significant mixing of the SM Higgs field with
the singlet Higgs field and is thus disfavored by the Higgs property
measurements at the LHC. The characteristics of Singlino-dominated DM are
contingent upon the values of $\lambda$, $\tan\beta$,
$m_{\tilde{\chi}_{1}^{0}}$, and $\mu_{\rm eff}\equiv\lambda v_{s}/\sqrt{2}$
[41], where $\lambda$ denotes the singlet–doublet Higgs Yukawa coupling within
the superpotential and $m_{\tilde{\chi}_{1}^{0}}\simeq 2\kappa\mu/\lambda$
requires $2|\kappa|/\lambda<1$. Considering that both SI and SD DM-nucleon
scattering cross sections scale with $\lambda^{4}$ in the regime of heavy
singlet Higgs bosons and consequently, the LZ experiment generally sets an
upper limit of $\lambda\lesssim 0.1$, the DM is likely to achieve the observed
abundance through one of two main mechanisms [29]: co-annihilating with
Higgsino-dominated electroweakinos or undergoing resonant annihilations
facilitated by singlet-dominated CP-even or CP-odd Higgs bosons. The former
scenario is only viable within a narrow parameter space characterized by
$2|\kappa|/\lambda\simeq 1$, $\lambda<0.1$, and $\mu<400$ GeV, where $\kappa$
represents the self-coupling coefficient of the singlet Higgs field [41]. To
achieve the measured density, the latter scenario requires
$2|m_{\tilde{\chi}_{1}^{0}}|$ to be close to the scalar mass. Thus, all these
scenarios correspond to an exceptionally constrained parameter space,
resulting in a significant suppression of Bayesian evidence111Bayesian
evidence is a fundamental concept in Bayesian statistics [42]. It is typically
denoted as $Z(D|M)\equiv\int{P(D|O(M,\Theta))P(\Theta|M)\prod d\Theta_{i}}$,
where $P(\Theta|M)$ represents the prior probability density function 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
observed values $O$. The Bayesian evidence considers both the theoretical
predictions $O(M,\Theta)$ and the experimental data $D$, and its computation
involves marginalizing the model’s parameters, which requires integrating the
likelihood function and prior probability over all possible parameter values.
In Bayesian statistics, evidence is used to compare different models and
determine which one is more plausible for explaining the observed data. While
a higher value of $Z$ suggests that the respective model is more likely to
agree with the data, an extremely small $Z$ means that the theory requires its
parameters to be fine-tuned for the same..

The situation becomes different in the General NMSSM (GNMSSM), where both
$\lambda$ and $\kappa$ can play distinct roles in DM physics [43, 44].
Specifically, the relationship that holds for Singlino DM in the
$\mathbb{Z}_{3}$-NMSSM, i.e., $2|\kappa|/\lambda<1$, no longer applies,
implying that $|\kappa|$ can be significantly larger than $\lambda$. In this
case, particles predominantly comprising the singlet sector, including the
Singlino-dominated DM and singlet-dominated Higgs bosons, can form an isolated
DM sector [45], where the DM attains the correct relic abundance by tuning the
value of $\kappa$ and through processes such as
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}A_{s},h_{s}h_{s},A_{s}A_{s}$. Here, $h_{s}$ and $A_{s}$ denote singlet-
dominated CP-even and CP-odd Higgs bosons, respectively. Similar to the
situation of $\mathbb{Z}_{3}$-NMSSM, the SI DM-nucleon scattering rate is
proportional to $\lambda^{2}\kappa^{2}$ for moderately light $h_{s}$ and
$\lambda^{4}$ for very massive $h_{s}$. It is thus suppressed to satisfy the
LZ restrictions for a small $\lambda$. As a result, the vacuum of the scalar
potential in the GNMSSM becomes more stable than that of the MSSM [46, 47].
Additionally, given the singlet nature of the DM and the complex mass
hierarchy, the decay chains of heavy supersymmetric particles in the GNMSSM
are lengthened compared to the prediction of the MSSM, which causes their
detection at the LHC to be rather tricky. These characteristics allow a more
extensive parameter space of the GNMSSM to be consistent with current
experimental results.

Based on the challenges faced by the MSSM and $\mathbb{Z}_{3}$-NMSSM, an
extensive investigation of the DM physics in the GNMSSM was performed in this
study. The rest of this paper is organized as follows. Section 2 outlines the
properties of the GNMSSM, Section 3 describes the research strategy and
numerical results to elucidate the characteristics of DM physics, and Section
4 provides a summary of the research findings.

##  2 Theoretical preliminaries

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 given by [24]

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

where $W_{\rm Yukawa}$ is the MSSM superpotential containing quark and lepton
Yukawa couplings, and the dimensionless coupling coefficients $\lambda$ and
$\kappa$ parameterize the interactions between the Higgs fields, the same as
those of the $\mathbb{Z}_{3}$-NMSSM. The bilinear mass parameters $\mu$ and
$\mu^{\prime}$ and the singlet tadpole parameter $\xi$ depict
$\mathbb{Z}_{3}$-symmetry-violating effects, which are advantageous for
solving the tadpole problem [48, 24] and the cosmological domain-wall problem
of the $\mathbb{Z}_{3}$-NMSSM [49, 50, 51]. Given that one of these parameters
can be eliminated by shifting the $\hat{S}$ field by a particular constant and
redefining the other parameters [52], we set $\xi$ to be zero without losing
the generality of this study. As suggested by previous studies [49, 53, 54,
52, 55], the natural values of the electroweak order for $\mu$ and
$\mu^{\prime}$ may arise from breaking the fundamental discrete $R$-symmetry,
$\mathbb{Z}^{R}_{4}$ or $\mathbb{Z}^{R}_{8}$, at high energy scales. They can
significantly change the properties of neutral Higgs bosons and neutralinos
and predict richer phenomenology than the $\mathbb{Z}_{3}$-NMSSM, which is the
focus of this study.

###  2.1 Higgs sector of GNMSSM

The soft supersymmetry-breaking Lagrangian for the Higgs fields in the GNMSSM
is given by

$\displaystyle-\mathcal{L}_{soft}=$ $\displaystyle\Bigg{[}\lambda
A_{\lambda}SH_{u}\cdot H_{d}+\frac{1}{3}\kappa
A_{\kappa}S^{3}+m_{3}^{2}H_{u}\cdot
H_{d}+\frac{1}{2}{m_{S}^{\prime}}^{2}S^{2}+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 scalar components of the Higgs
superfields, and $m^{2}_{H_{u}}$, $m^{2}_{H_{d}}$, and $m^{2}_{S}$ are their
supersymmetry-breaking masses. After the electroweak symmetry breaking, the
neutral Higgs fields acquire non-zero vevs:

$\displaystyle\left\langle
H_{u}^{0}\right\rangle=v_{u}/\sqrt{2},~{}~{}\left\langle
H_{d}^{0}\right\rangle=v_{d}/\sqrt{2},~{}~{}\left\langle
S\right\rangle=v_{s}/\sqrt{2},$ (2.3)

with $v=\sqrt{v_{u}^{2}+v_{d}^{2}}\simeq 246~{}\mathrm{GeV}$, and the three
masses are determined in solving the conditional equations to minimize the
scalar potential. The Higgs sector is then described by ten free parameters:

$\displaystyle\tan\beta,~{}\lambda,~{}\kappa,~{}v_{s},~{}A_{\lambda},~{}A_{\kappa},~{}\mu,~{}\mu^{\prime},~{}m_{3}^{2},~{}{m_{\text{S}}^{\prime}}^{2},$
(2.4)

where $\tan\beta$ is defined by $\tan\beta\equiv v_{u}/v_{d}$.

In revealing the properties of Higgs physics, it is customary to work with the
field combinations of $H_{\rm SM}\equiv\sin\beta{\rm
Re}(H_{u}^{0})+\cos\beta{\rm Re}(H_{d}^{0})$, $H_{\rm NSM}\equiv\cos\beta{\rm
Re}(H_{u}^{0})-\sin\beta{\rm Re}(H_{d}^{0})$, and $A_{\rm
NSM}\equiv\cos\beta{\rm Im}(H_{u}^{0})-\sin\beta{\rm Im}(H_{d}^{0})$, where
$H_{\rm SM}$ denotes the SM Higgs field, and $H_{\rm NSM}$ and $A_{\rm NSM}$
represent the extra doublet Higgs fields [56, 57]. In the bases $\left(H_{\rm
NSM},H_{\rm SM},{\rm Re}[S]\right)$, the mass matrix of the CP-even Higgs
fields takes the following form [24, 57]

$\displaystyle{\cal M}^{2}_{S,11}$ $\displaystyle=$
$\displaystyle\frac{\lambda v_{s}(\sqrt{2}A_{\lambda}+\kappa
v_{s}+\sqrt{2}\mu^{\prime})+2m_{3}^{2}}{\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,\quad{\cal
M}^{2}_{S,13}=-\frac{\lambda v}{\sqrt{2}}(A_{\lambda}+\sqrt{2}\kappa
v_{s}+\mu^{\prime})\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{\lambda v}{\sqrt{2}}\left[(\sqrt{2}\lambda
v_{s}+2\mu)-(A_{\lambda}+\sqrt{2}\kappa v_{s}+\mu^{\prime})\sin
2\beta\right],$ $\displaystyle{\cal M}^{2}_{S,33}$ $\displaystyle=$
$\displaystyle\frac{(A_{\lambda}+\mu^{\prime})\sin
2\beta}{2\sqrt{2}v_{s}}\lambda v^{2}+\frac{\kappa
v_{s}}{\sqrt{2}}(A_{\kappa}+2\sqrt{2}\kappa
v_{s}+3\mu^{\prime})-\frac{\mu}{\sqrt{2}v_{s}}\lambda v^{2}.\quad$ (2.5)

In the bases $\left(A_{\rm NSM},{\rm Im}(S)\right)$, the elements of the CP-
odd Higgs matrix are given by

$\displaystyle{\cal M}^{2}_{P,11}$ $\displaystyle=$
$\displaystyle\frac{\lambda v_{s}(\sqrt{2}A_{\lambda}+\kappa
v_{s}+\sqrt{2}\mu^{\prime})+2m_{3}^{2}}{\sin 2\beta},\quad{\cal
M}^{2}_{P,12}=\frac{\lambda v}{\sqrt{2}}(A_{\lambda}-\sqrt{2}\kappa
v_{s}-\mu^{\prime}),$ $\displaystyle{\cal M}^{2}_{P,22}$ $\displaystyle=$
$\displaystyle\frac{(A_{\lambda}+2\sqrt{2}\kappa v_{s}+\mu^{\prime})\sin
2\beta}{2\sqrt{2}v_{s}}\lambda v^{2}-\frac{\kappa
v_{s}}{\sqrt{2}}(3A_{\kappa}+\mu^{\prime})-\frac{\mu}{\sqrt{2}v_{s}}\lambda
v^{2}-2m_{S}^{\prime\ 2}.\quad\quad$ (2.6)

Three CP-even mass eigenstates $h_{i}=\\{h,H,h_{s}\\}$ and two CP-odd mass
eigenstates $a_{i}=\\{A_{H},A_{s}\\}$ are acquired through unitary rotations
of $V$ and $V_{P}$ to diagonalize ${\cal{M}}_{S}^{2}$ and ${\cal{M}}_{P}^{2}$,
respectively, leading to

$\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 V_{P,a_{i}}^{\rm
NSM}A_{\rm NSM}+V_{P,a_{i}}^{\rm S}Im[S].$ (2.7)

Among these states, $h$ corresponds to the scalar discovered at the LHC, $H$
and $A_{H}$ represent heavy doublet-dominated Higgs bosons, and $h_{s}$ and
$A_{s}$ denote singlet-dominated states. These states are also labeled as
$h_{i}$ (i=1,2,3) and $A_{j}$ (j=1,2) in ascending mass orders for
convenience, i.e., $m_{h_{1}}<m_{h_{2}}<m_{h_{3}}$ and $m_{A_{1}}<m_{A_{2}}$
in this study. Therefore, $h\equiv h_{1}$ and $m_{h_{s}}>m_{h}$ for the
$h_{1}$ scenario, and $h\equiv h_{2}$ and $m_{h}>m_{h_{s}}$ for the $h_{2}$
scenario. The model also predicts a pair of charged Higgs bosons,
$H^{\pm}=\cos\beta H_{u}^{\pm}+\sin\beta H_{d}^{\pm}$, with their masses given
by

$\displaystyle m^{2}_{H^{\pm}}$ $\displaystyle=$ $\displaystyle\frac{\lambda
v_{s}(\sqrt{2}A_{\lambda}+\kappa v_{s}+\sqrt{2}\mu^{\prime})+2m_{3}^{2}}{\sin
2\beta}+m^{2}_{W}-\frac{1}{2}\lambda^{2}v^{2}.$ (2.8)

Note that charged Higgs bosons $H^{\pm}$ degenerate with the CP-even doublet
scalar $H$ and the CP-odd scalar $A_{H}$ in mass.

So far, the LHC experiments have performed intensive searches for the
additional Higgs bosons $H$, $A_{H}$, $H^{\pm}$, $h_{s}$, and $A_{s}$ and set
limits on their properties, such as their masses and couplings [58, 59]. We
note that the parameters $A_{\kappa}$, $\mu$, $\mu^{\prime}$, $m_{3}^{2}$, and
$m_{S}^{\prime\ 2}$ are not directly related to experimental measurements.
This motivated us to use the masses of CP-odd heavy doublet Higgs field,
$CP$-even and -odd singlet Higgs fields, and Higgsino and Singlino fields,
denoted as $m_{A}\equiv\sqrt{{\cal M}^{2}_{P,11}}$, $m_{B}\equiv\sqrt{{\cal
M}^{2}_{S,33}}$, $m_{C}\equiv\sqrt{{\cal M}^{2}_{P,22}}$,
$\mu_{\text{tot}}\equiv\mu_{\text{eff}}+\mu$, and
$m_{N}\equiv\frac{2\kappa}{\lambda}\mu_{\rm eff}+\mu^{\prime}$, respectively,
as theoretical inputs. The former set of parameters are then given by

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\mu_{tot}-\frac{\lambda}{\sqrt{2}}v_{s},\quad\mu^{\prime}=m_{N}-\sqrt{2}\kappa
v_{s},\quad m^{2}_{3}=\frac{m^{2}_{A}\sin{2\beta}}{2}-\lambda
v_{s}(\frac{\kappa
v_{s}}{2}+\frac{\mu^{\prime}}{\sqrt{2}}+\frac{A_{\lambda}}{\sqrt{2}}),$
$\displaystyle\kappa A_{\kappa}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}m^{2}_{B}}{v_{s}}+\frac{\lambda\mu
v^{2}}{v_{s}^{2}}-\frac{\lambda(A_{\lambda}+\mu^{\prime})v^{2}\sin{2\beta}}{2v_{s}^{2}}-2\sqrt{2}\kappa^{2}v_{s}-3\kappa\mu^{\prime},$
$\displaystyle m_{S}^{\prime 2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left[m^{2}_{C}+\frac{\lambda\mu
v^{2}}{\sqrt{2}v_{s}}+\frac{\kappa
v_{s}(3A_{\kappa}+\mu^{\prime})}{\sqrt{2}}-\frac{\lambda(A_{\lambda}+2\sqrt{2}\kappa
v_{s}+\mu^{\prime})v^{2}\sin{2\beta}}{2\sqrt{2}v_{s}}\right],$ (2.9)

and Eqs. (2.5), (2.6), and (2.8) take the following simplified forms:

$\displaystyle{\cal M}^{2}_{S,11}$ $\displaystyle=$ $\displaystyle
m_{A}^{2}+\frac{1}{2}(2m_{Z}^{2}-\lambda^{2}v^{2})\sin^{2}2\beta,\quad{\cal
M}^{2}_{S,12}=-\frac{1}{4}(2m_{Z}^{2}-\lambda^{2}v^{2})\sin 4\beta,$
$\displaystyle{\cal M}^{2}_{S,13}$ $\displaystyle=$
$\displaystyle-\frac{\lambda v}{\sqrt{2}}(A_{\lambda}+m_{N})\cos
2\beta,\quad{\cal
M}^{2}_{S,22}=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{\lambda
v}{\sqrt{2}}\left[2\mu_{tot}-(A_{\lambda}+m_{N})\sin 2\beta\right],\quad{\cal
M}^{2}_{S,33}=m_{B}^{2},\quad{\cal M}^{2}_{P,11}=m_{A}^{2},$
$\displaystyle\quad{\cal M}^{2}_{P,22}$ $\displaystyle=$ $\displaystyle
m_{C}^{2},\quad{\cal M}^{2}_{P,12}=\frac{\lambda
v}{\sqrt{2}}(A_{\lambda}-m_{N}),\quad
m^{2}_{H^{\pm}}=m_{A}^{2}+m^{2}_{W}-\frac{1}{2}\lambda^{2}v^{2}.$ (2.10)

In the case of very massive charged Higgs bosons, the following approximations
hold [5]:

$\displaystyle m_{h_{s}}^{2}$ $\displaystyle\simeq$ $\displaystyle
m_{B}^{2}-\frac{{\cal M}^{4}_{S,13}}{m_{A}^{2}-m_{B}^{2}},\quad
m_{A_{s}}^{2}\simeq m_{C}^{2}-\frac{{\cal
M}^{4}_{P,12}}{m_{A}^{2}-m_{C}^{2}},\quad\frac{V_{P,A_{s}}^{\rm
NSM}}{V_{P,A_{s}}^{\rm S}}=\frac{{\cal
M}^{2}_{P,12}}{m_{A_{s}}^{2}-m_{A}^{2}}\simeq 0,$
$\displaystyle\frac{V_{h}^{\rm S}}{V_{h}^{\rm SM}}$ $\displaystyle\simeq$
$\displaystyle\frac{{\cal M}^{2}_{S,23}}{m_{h}^{2}-m_{B}^{2}},\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.$ (2.11)

Evidently, the singlet fields will decouple from the doublet Higgs fields in
the limit of $\lambda\to 0$, and the field masses $m_{B}$, $m_{C}$, and
$m_{N}$ can be regarded as physical particle masses to a good approximation.
Additionally, the singlet masses are independent and may all take small values
since they are weakly constrained by current experiments.

###  2.2 Neutralino sector of GNMSSM

The mixing between fermionic partners for neutral Higgs bosons and gauginos
results in five neutralinos and two charginos, denoted as
$\tilde{\chi}_{i}^{0}$ ($i=1,\dots,5$) and $\tilde{\chi}_{i}^{\pm}$ ($i=1,2$),
respectively. In the gauge eigenstate bases
$\psi^{0}=\left(-i\tilde{B},-i\tilde{W},\tilde{H}_{d}^{0},\tilde{H}_{u}^{0},\tilde{S}\right)$,
the symmetric neutralino mass matrix is given by [24]

$M_{\tilde{\chi}^{0}}=\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_{\text{tot}}&-\frac{1}{\sqrt{2}}\lambda v\sin\beta\\\
&&&0&-\frac{1}{\sqrt{2}}\lambda v\cos\beta\\\
&&&&m_{\text{N}}\end{array}\right),$ (2.12)

where $M_{1}$ and $M_{2}$ are gaugino soft-breaking masses,
$s_{W}=\sin\theta_{W}$, and $c_{w}=\cos\theta_{W}$. We use the rotation matrix
$N$ to diagonalize this mass matrix and label the resulting mass eigenstates
by an ascending mass order

$\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.13)

where the matrix elements $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 say
$\tilde{\chi}_{1}^{0}$ is Singlino-dominated if $N_{15}^{2}>0.5$.

In the case that the gauginos are very massive and
$\mu_{\text{tot}}^{2}-m_{N}^{2}\gg\lambda^{2}v^{2}$, the following
approximations for the Singlino-dominated $m_{\tilde{\chi}_{1}^{0}}$ hold [60,
61, 62]:

$\displaystyle m_{\tilde{\chi}_{1}^{0}}$ $\displaystyle\simeq$ $\displaystyle
m_{N}+\frac{1}{2}\frac{\lambda^{2}v^{2}(m_{\tilde{\chi}_{1}^{0}}-\mu_{\text{tot}}\sin
2\beta)}{m_{\tilde{\chi}_{1}^{0}}^{2}-\mu_{\text{tot}}^{2}}\simeq m_{N},\quad
N_{11}\sim 0,\quad N_{12}\sim 0,$ (2.14) $\displaystyle\frac{N_{13}}{N_{15}}$
$\displaystyle=$ $\displaystyle\frac{\lambda v}{\sqrt{2}\mu_{\rm
tot}}\frac{(m_{\tilde{\chi}_{1}^{0}}/\mu_{\text{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_{\text{tot}}}\frac{(m_{\tilde{\chi}_{1}^{0}}/\mu_{\text{tot}})\cos\beta-\sin\beta}{1-\left(m_{\tilde{\chi}_{1}^{0}}/\mu_{\text{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_{\text{tot}})^{2}\right]^{2}}{\left[(m_{\tilde{\chi}_{1}^{0}}/{\text{tot}})^{2}-2(m_{\tilde{\chi}_{1}^{0}}/\mu_{\text{tot}})\sin
2\beta+1\right]\left(\frac{\lambda
v}{\sqrt{2}\mu_{\text{tot}}}\right)^{2}+\left[1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{\text{tot}})^{2}\right]^{2}}.$

The ratio of the Higgsino fraction in $\tilde{\chi}_{1}^{0}$, denoted as
$Z_{h}\equiv N_{13}^{2}+N_{14}^{2}$, to the Singlino fractions $Z_{s}\equiv
N_{15}^{2}$ is acquired by

$\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.15)

These approximations reveal that the characteristics of $\tilde{\chi}_{1}^{0}$
are represented by five independent parameters: $\tan\beta$, $\lambda$,
$\kappa$, $\mu_{\text{tot}}$, and $m_{\tilde{\chi}_{1}^{0}}$, where $\kappa$
quantifies the self-interactions of the singlet fields. They reflect that both
$\lambda$ and $\mu_{\text{tot}}$ play a significant role in determining
$Z_{h}$; specifically, a smaller $\lambda$ or a larger $\mu_{\rm tot}$ can
suppress $Z_{h}$. Notably, if one chooses $m_{N}$, or equivalently
$m_{\tilde{\chi}_{1}^{0}}$, as a theoretical input, $\kappa$ does not
explicitly affect $m_{\tilde{\chi}_{1}^{0}}$, and it is independent of
$\lambda$ in influencing the other properties of $\tilde{\chi}_{1}^{0}$. This
feature distinguishes it from the Singlino-dominated $\tilde{\chi}_{1}^{0}$ in
the $Z_{3}$-NMSSM, whose properties are determined by four parameters, i.e.,
$\tan\beta$, $\lambda$, $\mu_{\text{eff}}$, and
$m_{\tilde{\chi}_{1}^{0}}\simeq 2\kappa\mu_{\text{eff}}/\lambda$, with
$\kappa$ satisfying $2|\kappa|<\lambda$ [41].

We point out the difference between the GNMSSM and $\mu$-extended NMSSM
discussed in Ref. [43]: in the former, the singlet masses $m_{B}$, $m_{C}$,
and $m_{N}$ and the Higgsino mass $\mu_{\rm tot}$ are independent, while in
the latter, they are correlated.

###  2.3 DM physics of GNMSSM

In the GNMSSM, the following interactions are pertinent to the DM physics:

$\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.16)
$\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 V_{P,a_{i}}^{\rm
NSM}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}A_{\rm NSM}}+V_{P,a_{i}}^{\rm
S}C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}Im[S]}$ (2.17)
$\displaystyle\simeq$ $\displaystyle-\sqrt{2}V_{P,a_{i}}^{\rm NSM}\lambda
N_{15}(N_{13}\cos\beta+N_{14}\sin\beta)+\sqrt{2}V_{P,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.18) $\displaystyle
C_{h_{i}h_{s}h_{s}}$ $\displaystyle\simeq$ $\displaystyle\lambda
v\,V_{h_{i}}^{\rm SM}V_{h_{s}}^{\rm S}V_{h_{s}}^{\rm S}(\lambda-\kappa\sin
2\beta)-\lambda\kappa vV_{h_{i}}^{\rm NSM}V_{h_{s}}^{\rm S}V_{h_{s}}^{\rm
S}\cos 2\beta$ (2.19) $\displaystyle+\sqrt{2}V_{h_{i}}^{\rm S}V_{h_{s}}^{\rm
S}V_{h_{s}}^{\rm S}(3\sqrt{2}\kappa^{2}v_{s}+\kappa
A_{\kappa}+3\kappa\mu^{\prime}),$ $\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$
(2.20) $\displaystyle+\sqrt{2}V_{h_{i}}^{S}(\sqrt{2}\kappa^{2}v_{s}-\kappa
A_{\kappa}+\kappa\mu^{\prime}).$

In the small $\lambda$ case, we conclude that

$\displaystyle
C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h_{s}}\simeq-\sqrt{2}\kappa,\quad
C_{h_{s}h_{s}h_{s}}\simeq\sqrt{2}(3\sqrt{2}\kappa^{2}v_{s}+\kappa
A_{\kappa}+3\kappa\mu^{\prime}),$ $\displaystyle
C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}A_{s}}\simeq-\sqrt{2}\kappa,\quad
C_{h_{s}A_{s}A_{s}}\simeq\sqrt{2}(\sqrt{2}\kappa^{2}v_{s}-\kappa
A_{\kappa}+\kappa\mu^{\prime}),$ (2.21)

where $A_{\kappa}$ and $\mu^{\prime}$ are determined by $\lambda$, $\kappa$,
$\tan\beta$, $v_{s}$, $A_{\lambda}$, $\mu_{\rm tot}$, $m_{N}$, and $m_{B}$
from Eq. (2.9), and the other interactions, such as
$C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h}$,
$C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}Z}$, $C_{hh_{s}h_{s}}$, and
$C_{hA_{s}A_{s}}$, are suppressed by $\lambda$.

####  2.3.1 DM relic density

The thermally averaged cross section for WIMP DM pair annihilation under a
non-relativistic approximation and in the absence of co-annihilation can be
expanded as [3, 63]

$\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.22)

where $a$ corresponds to the $s$-wave contribution at a zero relative
velocity, $b$ encompasses contributions from both $s$-wave and $p$-wave
processes, and $x\equiv m_{\tilde{\chi}_{1}^{0}}/T$, with $T$ denoting the
temperature of the early Universe. After integrating the density function from
the freeze-out temperature $x_{F}=m/T_{F}$ to infinity, the present thermal
relic abundance is given by [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.23)

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

As indicated by the following study, the GNMSSM predicts a Singlino-dominated
DM in most cases and achieves the observed relic abundance primarily through
the following channels:

  1. 1.

$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$

This channel, which can be the most important annihilation of
$\tilde{\chi}_{1}^{0}$ when $m_{\tilde{\chi}_{1}^{0}}>m_{h_{s}}$ and
$m_{A_{s}}>2m_{\tilde{\chi}_{1}^{0}}-m_{h_{s}}$, occurs by the $s$-channel
exchange of CP-even Higgs bosons and $t$-channel exchange of neutralinos. The
full analytic expression of $\langle\sigma v\rangle_{x_{F}}^{h_{s}h_{s}}$ is
very complex [64]. However, this can be significantly simplified in the case
of minor $\lambda$ paired with significant $|\kappa|$, where only the
contributions from the exchanges of $h_{s}$ and $\tilde{\chi}_{1}^{0}$ are
crucial. Correspondingly, $\langle\sigma v\rangle_{x_{F}}^{h_{s}h_{s}}$ is
approximated by [65]

$\displaystyle\left\langle\sigma v\right\rangle_{x_{F}}^{h_{s}h_{s}}$
$\displaystyle\simeq$ $\displaystyle\frac{v^{2}_{F}}{192\pi
m_{\tilde{\chi}_{1}^{0}}^{2}}\sqrt{1-\frac{m_{h_{s}}^{2}}{m_{\tilde{\chi}_{1}^{0}}^{2}}}\times\left\\{\frac{8C_{h_{s}h_{s}h_{s}}C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h_{s}}^{3}m_{\tilde{\chi}_{1}^{0}}^{3}(2m_{h_{s}}^{2}-5m_{\tilde{\chi}_{1}^{0}}^{2})}{(m_{h_{s}}^{2}-4m_{\tilde{\chi}_{1}^{0}}^{2})(m_{h_{s}}^{2}-2m_{\tilde{\chi}_{1}^{0}}^{2})}\right.$
(2.24)
$\displaystyle\left.+\frac{3C_{h_{s}h_{s}h_{s}}^{2}C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h_{s}}^{2}m_{\tilde{\chi}_{1}^{0}}^{2}}{(m_{h_{s}}^{2}-4m_{\tilde{\chi}_{1}^{0}}^{2})^{2}}+\frac{16C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h_{s}}^{4}(9m_{\tilde{\chi}_{1}^{0}}^{8}-8m_{\tilde{\chi}_{1}^{0}}^{6}m_{h_{s}}^{2}+2m_{h_{s}}^{8})}{(m_{h_{s}}^{2}-2m_{\tilde{\chi}_{1}^{0}}^{2})^{4}}\right\\}.\quad\quad\quad$

This formula together with Eq. (2.21) indicate that $\left\langle\sigma
v\right\rangle_{x_{F}}^{h_{s}h_{s}}$ depends not only on $\kappa$,
$m_{h_{s}}$, and $m_{\tilde{\chi}_{1}^{0}}$, but also on $\lambda$,
$\tan\beta$, $v_{s}$, $\mu_{\rm tot}$, and $A_{\lambda}$ via
$C_{h_{s}h_{s}h_{s}}$. In particular, $\kappa$ plays a critical role in
determining the magnitude of $\left\langle\sigma
v\right\rangle_{x_{F}}^{h_{s}h_{s}}$, since most of its dominant contributions
are proportional to $\kappa^{4}$. If $C_{h_{s}h_{s}h_{s}}$ is sufficiently
small so that the s-channel contribution can be neglected, we have

$\displaystyle\left\langle\sigma v\right\rangle_{x_{F}}^{h_{s}h_{s}}$
$\displaystyle\simeq$ $\displaystyle\frac{3v_{F}^{2}\kappa^{4}}{16\pi
m_{\tilde{\chi}_{1}^{0}}^{2}},$ (2.25)

if $m_{\tilde{\chi}_{1}^{0}}^{2}\gg m_{h_{s}}^{2}$, which implies that

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

allowing the observed density to be acquired.

  2. 2.

$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$

This channel can be the dominant annihilation of $\tilde{\chi}_{1}^{0}$ if
$m_{\tilde{\chi}_{1}^{0}}>m_{A_{s}}$ and
$m_{h_{s}}>2m_{\tilde{\chi}_{1}^{0}}-m_{A_{s}}$. It proceeds in a way similar
to $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$, and
$\left\langle\sigma v\right\rangle_{x_{F}}^{A_{s}A_{s}}$ is given by [65]

$\displaystyle\left\langle\sigma v\right\rangle_{x_{F}}^{A_{s}A_{s}}$
$\displaystyle\simeq$ $\displaystyle\frac{v^{2}_{F}}{128\pi
m_{\tilde{\chi}_{1}^{0}}^{2}}\sqrt{1-\frac{m_{A_{s}}^{2}}{m_{\tilde{\chi}_{1}^{0}}^{2}}}\times$
(2.27)
$\displaystyle\left\\{\frac{32}{3}\frac{C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}A_{s}}^{4}m_{\tilde{\chi}_{1}^{0}}^{4}(m_{A_{s}}^{2}-m_{\tilde{\chi}_{1}^{0}}^{2})^{2}}{(m_{A_{s}}^{2}-2m_{\tilde{\chi}_{1}^{0}}^{2})^{4}}+\frac{4C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h_{s}}^{2}C_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}A_{s}}^{2}m_{\tilde{\chi}_{1}^{0}}^{2}m_{h_{s}}^{2}}{(m_{h_{s}}^{2}-4m_{\tilde{\chi}_{1}^{0}}^{2})^{2}+m_{h_{s}}^{2}\Gamma_{h_{s}}^{2}}\right\\}.\quad\quad\quad$

Compared with $\left\langle\sigma v\right\rangle_{x_{F}}^{h_{s}h_{s}}$,
$\left\langle\sigma v\right\rangle_{x_{F}}^{A_{s}A_{s}}$ relies on one more
parameter, $m_{A_{s}}$, and for $m_{h_{s}}\simeq 2m_{\tilde{\chi}_{1}^{0}}$,
it can be resonantly enhanced. As a result, it possesses more complex features
than $\left\langle\sigma v\right\rangle_{x_{F}}^{h_{s}h_{s}}$. Similar to the
previous discussion, if we neglect the s-channel contribution and assume
$m_{\tilde{\chi}_{1}^{0}}^{2}\gg m_{A_{s}}^{2}$, we conclude that

$\displaystyle\left\langle\sigma v\right\rangle_{x_{F}}^{A_{s}A_{s}}$
$\displaystyle\simeq$ $\displaystyle\frac{v_{F}^{2}\kappa^{4}}{48\pi
m_{\tilde{\chi}_{1}^{0}}^{2}}$ (2.28)

and

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

to predict the observed density.

  3. 3.

$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$

This process occurs if $2m_{\tilde{\chi}_{1}^{0}}>m_{h_{s}}+m_{A_{s}}$ and
proceeds through the $s$-channel exchange of $Z$ and CP-odd Higgs bosons and
the $t$-channel exchange of neutralinos. Compared with the previous two
annihilations, the full analytic expression of $\langle\sigma
v\rangle_{x_{F}}^{h_{s}A_{s}}$ is much more complex [64]. However, given that
the contributions from the exchanges of $A_{s}$ and $\tilde{\chi}_{1}^{0}$ are
dominant for the small $\lambda$ and sizable $|\kappa|$ case, it can be
simplified as follows [5, 63]

$\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.30)

respectively, assuming that $h_{s}$ is not exceptionally light to forbid
$A_{s}$-mediated resonant annihilation. Furthermore, from the formula of
$C_{h_{s}A_{s}A_{s}}$ in Eq. (2.21), one can infer that if $|\sqrt{2}\kappa
v_{s}-A_{\kappa}+\mu^{\prime}|\ll(4m_{\tilde{\chi}_{1}^{0}}^{2}-m_{A_{s}}^{2})/m_{\tilde{\chi}_{1}^{0}}$,
$|{\cal{A}}_{t}|$ is much larger than $|{\cal{A}}_{s}|$, reflecting that the
$t$-channel contribution prevails. In this case, $\left\langle\sigma
v\right\rangle_{x_{F}}^{h_{s}A_{s}}$ is simplified as

$\displaystyle\left\langle\sigma v\right\rangle_{x_{F}}^{h_{s}A_{s}}$
$\displaystyle\simeq$ $\displaystyle\frac{\kappa^{4}}{4\pi
m_{\tilde{\chi}_{1}^{0}}^{2}},$ (2.31)

if $m_{\tilde{\chi}_{1}^{0}}$ is much larger than $m_{A_{s}}$ and $m_{h_{s}}$,
and

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

to predict the observed density.

Notably, once the process $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}A_{s}$ is kinematically allowed, at least one of the annihilations
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$ and
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$ is open. Since the
former channel is dominated by s-wave contributions while the latter ones are
p-wave processes, $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$
dominantly contributes to $\left\langle\sigma v\right\rangle_{x_{F}}$.
Additionally, the annihilation $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
hA_{s}$, which proceeds in a manner similar to
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$, may provide
substantial contributions to $\left\langle\sigma v\right\rangle_{x_{F}}$ when
$\lambda$ and $V^{S}_{h}$ take exceptionally large values. Such a possibility,
however, has been rigorously bounded by the results of the LZ experiment and
the measurements of Higgs property at the LHC.

To achieve the observed relic abundance, the Singlino-dominated DM can also
co-annihilate with electroweakinos. Corresponding processes include
$\tilde{\chi}_{i}\tilde{\chi}_{i}$, $\tilde{\chi}_{i}\tilde{\chi}_{j}$,
$\tilde{\chi}_{j}\tilde{\chi}_{j}\to XX^{\prime}$, where $\tilde{\chi}_{i}$
and $\tilde{\chi}_{j}$ denote either the LSP or next-to-lightest
supersymmetric particle (NLSP) and $XX^{\prime}$ represents SM particles.
Basically, these annihilations influence the abundance only when the mass
difference between $\tilde{\chi}_{1}^{0}$ and its co-annihilation partner is
less than approximately $10\%$ [63, 66], which necessitates the tuning of the
two theoretically independent masses and thus leads to the suppression of
Bayesian evidence.

####  2.3.2 DM-nucleon scattering

In the regime of heavy squarks, the SI scattering of DM with nucleons
predominantly arises from the $t$-channel exchange of $CP$-even Higgs bosons.
The cross section takes the following form [61, 67]

$\displaystyle\sigma^{\rm SI}_{N}=\frac{4\mu_{r}^{2}}{\pi}|f^{N}|^{2},\quad
f^{N}=\sum_{i}^{3}f^{N}_{h_{i}}=\sum_{i}^{3}\frac{C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h_{i}}C_{NNh_{i}}}{2m^{2}_{h_{i}}},$
(2.33)

where $N=p,n$ denotes a proton ($p$) or neutron ($n$), and $\mu_{r}\equiv
m_{N}m_{\tilde{\chi}^{0}_{1}}/(m_{N}+m_{\tilde{\chi}^{0}_{1}})$ is the reduced
mass of the DM-nucleon composite. $C_{NNh_{i}}$ represents the strength of
Higgs coupling to the nucleon, 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]\,.\quad$ (2.34)

Here, $F^{N}_{d}$ and $F^{N}_{u}$ are defined as $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}$, where the nucleon form factors
$f^{(N)}_{q}\equiv m_{N}^{-1}\left<N|m_{q}q\bar{q}|N\right>$ for $q=u,d,s$
denote the normalized light quark contribution to the nucleon mass and
$f^{(N)}_{G}\equiv 1-\sum_{q=u,d,s}f^{(N)}_{q}$ represents other heavy quarks’
mass fraction in the nucleon [68, 69]. For the default settings of the package
micrOMEGAs for $f_{q}^{(N)}$ [70, 71, 72], $F_{u}^{p}\simeq F_{u}^{n}\simeq
0.15$, $F_{d}^{p}\simeq F_{d}^{n}\simeq 0.13$, and therefore, $\sigma^{\rm
SI}_{\tilde{\chi}^{0}_{1}-p}\simeq\sigma^{\rm SI}_{\tilde{\chi}^{0}_{1}-n}$.

Because the $H$-mediated contribution to $\sigma^{\rm SI}_{N}$ usually plays a
minor role222This contribution is proportional to $\tan^{2}\beta/m_{H}^{4}$
and is typically less than $10^{-49}~{}{\rm cm^{2}}$ for $\tan\beta\leq 5$, as
indicated by Figure 6 of [5]. With the further increase of $\tan\beta$,
$m_{H}$ is stringently limited by the LHC search for extra Higgs bosons (see,
e.g., Figure 22 of [73] for the MSSM results.). Consequently, the contribution
has difficulty reaching $10^{-47}~{}{\rm cm^{2}}$. Therefore, barring the
blind spots discussed in [74], it is less crucial., we focus on the
contribution from the $t$-channel exchange of the SM-like Higgs boson $h$ and
the singlet Higgs boson $h_{s}$. Using the formulas from Eqs.(2.14)–(2.20), we
obtain the couplings $C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h}$ and
$C_{\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}h_{s}}$ for the Singlino-dominated
$\tilde{\chi}_{1}^{0}$ as follows:

$\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.35)
$\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.36)
$\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],$

and the expression of the SI cross section is [43]

$\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.37)

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.38)

Given the approximations of $V_{h}^{\rm S}$ and $V_{h_{s}}^{\rm SM}$ in Eq.
(2.11), these formulas suggest that the primary contribution to $\cal{A}$ in
the series expansion of $\lambda$ is proportional to $\lambda\kappa$ if the
singlet–doublet mixing is substantial and the contribution mediated by $h$
significantly cancels out that mediated by $h_{s}$. This dynamic is critical
only when $h_{s}$ is moderately light. In contrast, if $h_{s}$ is tremendously
massive and consequently $h$ is purely SM like (i.e., $V_{h}^{\rm SM}\simeq
1$, $V_{h}^{\rm S}\simeq 0$, and $V_{h_{s}}^{\rm SM}\simeq 0$), ${\cal{A}}$’s
expression is significantly simplified:

$\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.39)

indicating that ${\cal{A}}$ is proportional to $\lambda^{4}$. Either a small
$\lambda$ or the correlation $m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm tot}\simeq\sin
2\beta$, known as the blind-spot condition of the NMSSM [75], is required to
suppress the SI cross section.

Regarding the spin-dependent (SD) scattering cross section, only the
$t$-channel exchange of a $Z$ boson contributes in the limit of very massive
squarks. It is approximated by [61, 60]

$\displaystyle\sigma_{\tilde{\chi}_{1}^{0}-N}^{\rm SD}$ $\displaystyle\simeq$
$\displaystyle C_{N}\times 10^{-4}~{}{\rm
pb}\times\left(\frac{N_{13}^{2}-N_{14}^{2}}{0.1}\right)^{2},$ (2.40)
$\displaystyle\simeq$ $\displaystyle C_{N}\times 10^{-2}~{}{\rm
pb}\times\left(\frac{\lambda v}{\sqrt{2}\mu_{\rm
tot}}\right)^{4}\left(\frac{N_{15}^{2}\cos{2\beta}}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot})^{2}}\right)^{2},$ (2.41)

where $C_{N}\simeq 4.0$ for protons and $C_{N}\simeq 3.1$ for neutrons. This
formula indicates that the SD cross section is proportional to $(\lambda
v/\mu_{tot})^{4}$ and increases as $m_{\tilde{\chi}_{1}^{0}}$ approaches
$\mu_{tot}$ from below. It is the same as the expression of the SD scattering
cross section for Bino-dominated DM,

$\displaystyle\sigma_{\tilde{B}_{1}^{0}-N}^{\rm SD}$ $\displaystyle\simeq$
$\displaystyle C_{N}\times 10^{-2}~{}{\rm
pb}\times\left(\frac{v}{\sqrt{2}\mu_{\rm
tot}}\right)^{4}\left(\frac{\cos{2\beta}}{1-(m_{\tilde{\chi}_{1}^{0}}/\mu_{\rm
tot})^{2}}\right)^{2},$ (2.42)

except for an additional factor $\lambda^{4}N_{15}^{2}$.

In summary, the aforementioned formulas suggest the plausible existence of a
secluded DM sector, composed of singlet states [45]. Specifically, the
Singlino-dominated DM achieves the observed relic abundance primarily through
the processes $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}A_{s},h_{s}h_{s},A_{s}A_{s}$, or $h_{s}/A_{s}$-funnel by tuning $\kappa$.
Since the DM sector communicates with the SM sector mainly by the weak
singlet–doublet Higgs mixing, the DM-nucleon scattering is suppressed by
$\lambda$.

##  3 Numerical results

Here the sampling strategy of this study is introduced and the features of DM
physics in the GNMSSM are demonstrated. The model file of the GNMSSM was
generated by the package SARAH 4.14.3 [76, 77, 78, 79], the particle spectrum
and low-energy flavor measurements were calculated using the SPheno 4.0.4 [80,
81] and FlavorKit [82] code, respectively, and DM observables were computed
with the package micrOMEGAs 5.0.4 [83, 84, 85, 86, 87, 88]. We analyzed the
acquired samples using the marginal posterior probability density function
(PDF) in Bayesian inference and the profile likelihood (PL) in frequentist
statistics [89]333In frequentist statistics, the one-dimensional PL refers to
the maximum likelihood value in a specific parameter space [89]. For a given
set of input parameters
$\Theta\equiv\left(\Theta_{1},\Theta_{2},\cdots\right)$, the one-dimensional
PL, denoted as
$\mathcal{L}\left(\Theta_{A}\right)=\max_{\Theta_{1},\cdots,\Theta_{A-1},\Theta_{A+1},\cdots}\mathcal{L}(\Theta),$
(3.1) is obtained by maximizing the likelihood function while varying the
other parameters. At a given point $\Theta_{A}$,
$\mathcal{L}\left(\Theta_{A}\right)$ reflects the capability of that point
within the theory to explain the experimental data. The one-dimensional
marginal posterior PDF is derived through the integration of the posterior PDF
from the Bayesian theorem with respect to the remaining inputs of the model,
i.e.,  $P\left(\Theta_{A}\right)=\int P(\Theta)d\Theta_{1}d\Theta_{2}\cdots
d\Theta_{A-1}d\Theta_{A+1}\cdots\cdots.$ (3.2) The PL reflects the data’s
preference for the parameter space, while the posterior PDF represents the
preference of the samples obtained during the scan..

###  3.1 Research strategy

Table 1: In this study, the parameter space was explored for the $h_{1}$
scenario, assuming all the inputs were flatly distributed in the prior since
they have clear physical meanings in the small $\lambda$ case. We assumed that
the soft trilinear coefficients for the third-generation squarks were equal,
i.e., $A_{t}=A_{b}$, for simplicity and let them vary because they could
significantly influence the Standard Model (SM)-like Higgs boson mass by
radiative corrections. We also fixed $M_{3}=3~{}{\rm TeV}$ and selected a
shared value of $2~{}{\rm TeV}$ for the other unmentioned dimensional
parameters to be consistent with the Large Hadron Collider (LHC) search for
new physics, noting that they are not crucial to this study. We define all
these parameters at the renormalization scale $Q_{inp}=1~{}{\rm TeV}$. The
$h_{2}$ scenario corresponds to the same parameter space as that of the
$h_{1}$ scenario except that $1~{}{\rm GeV}\leq m_{B}\leq 130~{}{\rm GeV}$.

Parameter | Prior | Range | Parameter | Prior | Range
---|---|---|---|---|---
$\lambda$ | Flat |  $0$–$0.75$ | $\kappa$ | Flat |  $-0.75$–$0.75$
$\tan\beta$ | Flat |  $1$–$60$ | $v_{s}/{\rm TeV}$ | Flat |  $10^{-3}$–$1.0$
$\mu_{\rm tot}/{\rm TeV}$ | Flat |  $0.1$–$1.0$ | $m_{N}/{\rm TeV}$ | Flat |  $-1.0$–$1.0$
$m_{B}/{\rm TeV}$ | Flat |  $0.1$–$1.0$ | $m_{C}/{\rm TeV}$ | Flat |  $10^{-3}$–$1.0$
$M_{1}/{\rm TeV}$ | Flat |  $-1.0$–$1.0$ | $M_{2}/{\rm TeV}$ | Flat |  $0.1$–$1.0$
$A_{t}/{\rm TeV}$ | Flat |  $-5.0$–$5.0$ | $A_{\lambda}/{\rm TeV}$ | Flat |  $-2.5$–$2.5$

We started by constructing a likelihood function that incorporated the DM
relic abundance, results from the LZ experiment [9], and other relevant
constraints on the GNMSSM. The likelihood function is given by

$\displaystyle\mathcal{L}$ $\displaystyle\equiv$
$\displaystyle\mathcal{L}_{\Omega h^{2}}\times\mathcal{L}_{\rm
LZ}\times\mathcal{L}_{\rm Const},$ $\displaystyle\mathcal{L}_{\Omega h^{2}}$
$\displaystyle=$ $\displaystyle\exp\left[-\frac{1}{2}\left(\frac{\Omega
h^{2}-0.120}{0.012}\right)^{2}\right],\quad\mathcal{L}_{\rm
LZ}=\exp\left[-\frac{\sigma^{\rm
SI}_{\tilde{\chi}_{1}^{0}-p}}{2\delta_{\sigma}^{2}}\right],$
$\displaystyle\mathcal{L}_{\rm Const}$ $\displaystyle=$
$\displaystyle\left\\{\begin{aligned} &1&&{\rm if\ satisfying\ all\
experimental\ constraints}\\\ &\exp\left[-100\right]&&{\rm
otherwise}\end{aligned}\right..$ (3.3)

In this formulation, we assumed the relic abundance $\Omega h^{2}$ followed a
Gaussian distribution. We used its central value from the Planck experiment
[2] with a theoretical uncertainty of $10\%$ in its calculation. Regarding the
likelihood function of the LZ experiment, $\mathcal{L}_{\rm LZ}$, we modeled
it with a Gaussian distribution centered at zero and defined
$\delta_{\sigma}^{2}={\rm UL}_{\sigma}^{2}/1.64^{2}+(0.2\sigma)^{2}$, where
${\rm UL}_{\sigma}$ refers to the upper limit of the LZ results on the DM-
nucleon SI scattering rate at a $90\%$ C. L., and $0.2\sigma$ accounts for
theoretical uncertainties [90]. $\mathcal{L}_{\text{Const}}$ includes the
following experimental constraints:

  * •

Higgs data fit: The properties of the SM-like Higgs boson $h$ must align with
the LHC Higgs data at a $95\%$ confidence level. This condition was examined
using the HiggsSignal 2.6.2code, requiring the fit’s $p$-value to exceed 0.05
[91, 92, 93, 94].

  * •

Extra Higgs searches: This aspect was tested using the HiggsBounds 5.10.2 code
to conduct a comprehensive search for additional Higgs bosons beyond the SM at
LEP, Tevatron, and the LHC [95, 96, 97, 98, 99].

  * •

Indirect DM searches: The Fermi-LAT collaboration has made years of
observations of dwarf galaxies, limiting the annihilation cross section as a
function of the DM mass444see website: www-glast.stanford.edu/pub data/1048.
We employed the likelihood function proposed in Ref. [100, 101] to implement
this constraint.

  * •

$B$-physics observables: The branching ratios of $B_{s}\to\mu^{+}\mu^{-}$ and
$B\to X_{s}\gamma$ should remain consistent with their experimental
measurements at the $2\sigma$ level [102].

  * •

Vacuum stability: The vacuum state of the scalar potential should be either
stable or long lived [47]. This criterion was rigorously examined using
VevaciousPlusPlus [103], the C++ version of Vevacious [104]. Note that its
implications on the GNMSSM were comprehensively addressed in a recent study
[47].

Table 2: Experimental analyses included in the package SModelS-2.1.1.

Name | Scenario | Final State | $\texttt{Luminosity}(\texttt{fb}^{\texttt{-1}})$
---|---|---|---
| CMS-SUS-17-010 [105]
---
| $\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}\rightarrow
W^{\pm}\tilde{\chi}_{1}^{0}W^{\mp}\tilde{\chi}_{1}^{0}$
---
$\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}\rightarrow\nu\tilde{\ell}\ell\tilde{\nu}$
2$\ell$ \+ $E_{\rm T}^{\rm miss}$ | 35.9
| CMS-SUS-17-009 [106]
---
$\tilde{\ell}\tilde{\ell}\rightarrow\ell\tilde{\chi}_{1}^{0}\ell\tilde{\chi}_{1}^{0}$ | 2$\ell$ \+ $E_{\rm T}^{\rm miss}$ | 35.9
| CMS-SUS-17-004 [107]
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow Wh(Z)\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | n$\ell$(n>=0) + nj(n>=0) + $E_{\rm T}^{\rm miss}$ | 35.9
| CMS-SUS-16-045 [108]
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow W^{\pm}\tilde{\chi}_{1}^{0}h\tilde{\chi}_{1}^{0}$ | 1$\ell$ 2b + $E_{\rm T}^{\rm miss}$ | 35.9
| CMS-SUSY-16-039 [109]
---
|
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow\ell\tilde{\nu}\ell\tilde{\ell}$
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow\tilde{\tau}\nu\tilde{\ell}\ell$
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow\tilde{\tau}\nu\tilde{\tau}\tau$
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow
WZ\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow
WH\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$
n$\ell(n\textgreater{}0)$($\tau$) + $E_{\rm T}^{\rm miss}$ | 35.9
| CMS-SUS-16-034 [110]
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow W\tilde{\chi}_{1}^{0}Z(h)\tilde{\chi}_{1}^{0}$ | n$\ell$(n>=2) + nj(n>=1) $E_{\rm T}^{\rm miss}$ | 35.9
| ATLAS-1803-02762 [111]
---
| $\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow
WZ\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow\nu\tilde{\ell}l\tilde{\ell}$
$\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}\rightarrow\nu\tilde{\ell}\nu\tilde{\ell}$
$\tilde{\ell}\tilde{\ell}\rightarrow\ell\tilde{\chi}_{1}^{0}\ell\tilde{\chi}_{1}^{0}$
n$\ell$ (n>=2) + $E_{\rm T}^{\rm miss}$ | 36.1
| ATLAS-1812-09432 [112]
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow Wh\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | n$\ell$ (n>=0) + nj(n>=0) + nb(n>=0) + n$\gamma$(n>=0) + $E_{\rm T}^{\rm miss}$ | 36.1
| ATLAS-1806-02293 [113]
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow WZ\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | n$\ell$(n>=2) + nj(n>=0) + $E_{T}^{miss}$ | 36.1
| ATLAS-1912-08479 [114]
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow W(\rightarrow l\nu)\tilde{\chi}_{1}^{0}Z(\rightarrow\ell\ell)\tilde{\chi}_{1}^{0}$ | 3$\ell$ \+ $E_{\rm T}^{\rm miss}$ | 139
| ATLAS-1908-08215 [115]
---
|
$\tilde{\ell}\tilde{\ell}\rightarrow\ell\tilde{\chi}_{1}^{0}\ell\tilde{\chi}_{1}^{0}$
---
$\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}$
2$\ell$ \+ $E_{\rm T}^{\rm miss}$ | 139
| ATLAS-1909-09226 [116]
---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow Wh\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | 1$\ell$ \+ h($\bm{\rightarrow}$ bb) + $E_{\rm T}^{\rm miss}$ | 139

Table 3: Experimental analyses of the electroweakino production processes
considered in this study, categorized by the topologies of the
electroweakino’s signal.

Scenario | Final State | Name
---|---|---
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow WZ\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | $n\ell(n\geq 2)+nj(n\geq 0)+\text{E}_{\text{T}}^{\text{miss}}$ |  CMS-SUS-20-001($137fb^{-1}$) [117]
ATLAS-2106-01676($139fb^{-1}$) [118]
CMS-SUS-17-004($35.9fb^{-1}$) [107]
CMS-SUS-16-039($35.9fb^{-1}$) [109]
ATLAS-1803-02762($36.1fb^{-1}$) [111]
ATLAS-1806-02293($36.1fb^{-1}$) [113]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow\ell\tilde{\nu}\ell\tilde{\ell}$ | $n\ell(n=3)+\text{E}_{\text{T}}^{\text{miss}}$ |  CMS-SUS-16-039($35.9fb^{-1}$) [109]
ATLAS-1803-02762($36.1fb^{-1}$) [111]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow\tilde{\tau}\nu\ell\tilde{\ell}$ | $2\ell+1\tau+\text{E}_{\text{T}}^{\text{miss}}$ |  CMS-SUS-16-039($35.9fb^{-1}$) [109]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow\tilde{\tau}\nu\tilde{\tau}\tau$ | $3\tau+\text{E}_{\text{T}}^{\text{miss}}$ |  CMS-SUS-16-039($35.9fb^{-1}$) [109]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow Wh\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | $n\ell(n\geq 1)+nb(n\geq 0)+nj(n\geq 0)+\text{E}_{\text{T}}^{\text{miss}}$ |  ATLAS-1909-09226($139fb^{-1}$) [119]
CMS-SUS-17-004($35.9fb^{-1}$) [107]
CMS-SUS-16-039($35.9fb^{-1}$) [109]
ATLAS-1812-09432($36.1fb^{-1}$)[112]
CMS-SUS-16-034($35.9fb^{-1}$)[110]
CMS-SUS-16-045($35.9fb^{-1}$) [108]
$\tilde{\chi}_{1}^{\mp}\tilde{\chi}_{1}^{\pm}\rightarrow WW\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | $2\ell+\text{E}_{\text{T}}^{\text{miss}}$ |  ATLAS-1908-08215($139fb^{-1}$) [115]
CMS-SUS-17-010($35.9fb^{-1}$) [105]
$\tilde{\chi}_{1}^{\mp}\tilde{\chi}_{1}^{\pm}\rightarrow 2\tilde{\ell}\nu(\tilde{\nu}\ell)$ | $2\ell+\text{E}_{\text{T}}^{\text{miss}}$ |  ATLAS-1908-08215($139fb^{-1}$) [115]
CMS-SUS-17-010($35.9fb^{-1}$) [105]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow ZW\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | $2j(\text{large})+\text{E}_{\text{T}}^{\text{miss}}$ | ATLAS-2108-07586($139fb^{-1}$) [120]
$\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}\rightarrow
WW\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow(h/Z)W\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | $j(\text{large})+b(\text{large})+\text{E}_{\text{T}}^{\text{miss}}$ | ATLAS-2108-07586($139fb^{-1}$) [120]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{3}^{0}\rightarrow(h/Z)Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\mp}\rightarrow h/ZW\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0},\tilde{\chi}_{1}^{0}\rightarrow\gamma/Z\tilde{G}$ | $2\gamma+n\ell(n\geq 0)+nb(n\geq 0)+nj(n\geq 0)+\text{E}_{\text{T}}^{\text{miss}}$ | ATLAS-1802-03158($36.1fb^{-1}$) [121]
$\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}\rightarrow
WW\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0},\tilde{\chi}_{1}^{0}\rightarrow\gamma/Z\tilde{G}$
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow ZW\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0},\tilde{\chi}_{1}^{0}\rightarrow h/Z\tilde{G}$ | $n\ell(n\geq 4)+\text{E}_{\text{T}}^{\text{miss}}$ | ATLAS-2103-11684($139fb^{-1}$) [122]
$\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}\rightarrow
WW\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0},\tilde{\chi}_{1}^{0}\rightarrow
h/Z\tilde{G}$
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0}\rightarrow
Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0},\tilde{\chi}_{1}^{0}\rightarrow
h/Z\tilde{G}$
$\tilde{\chi}_{1}^{\mp}\tilde{\chi}_{1}^{0}\rightarrow
W\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0},\tilde{\chi}_{1}^{0}\rightarrow
h/Z\tilde{G}$
$\tilde{\chi}_{i}^{0,\pm}\tilde{\chi}_{j}^{0,\mp}\rightarrow\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}+\chi_{soft}\rightarrow ZZ/H\tilde{G}\tilde{G}$ | $n\ell(n\geq 2)+nb(n\geq 0)+nj(n\geq 0)+\text{E}_{\text{T}}^{\text{miss}}$ |  CMS-SUS-16-039($35.9fb^{-1}$) [109]
CMS-SUS-17-004($35.9fb^{-1}$) [107]
CMS-SUS-20-001($137fb^{-1}$) [117]
$\tilde{\chi}_{i}^{0,\pm}\tilde{\chi}_{j}^{0,\mp}\rightarrow\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}+\chi_{soft}\rightarrow HH\tilde{G}\tilde{G}$ | $n\ell(n\geq 2)+nb(n\geq 0)+nj(n\geq 0)+\text{E}_{\text{T}}^{\text{miss}}$ |  CMS-SUS-16-039($35.9fb^{-1}$) [109]
CMS-SUS-17-004($35.9fb^{-1}$) [107]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow W^{*}Z^{*}\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | $3\ell+\text{E}_{\text{T}}^{\text{miss}}$ |  ATLAS-2106-01676($139fb^{-1}$) [118]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}\rightarrow Z^{*}W^{*}\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ | $2\ell+nj(n\geq 0)+\text{E}_{\text{T}}^{\text{miss}}$ |  ATLAS-1911-12606($139fb^{-1}$) [123]
ATLAS-1712-08119($36.1fb^{-1}$) [124]
|  CMS-SUS-16-048($35.9fb^{-1}$) [125]
$\tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{\pm}+\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{\mp}+\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{1}^{0}$ | $2\ell+nj(n\geq 0)+\text{E}_{\text{T}}^{\text{miss}}$ |  ATLAS-1911-12606($139fb^{-1}$) [123]
ATLAS-1712-08119($36.1fb^{-1}$) [124]
CMS-SUS-16-048($35.9fb^{-1}$) [125]

Next, we performed sophisticated scans over the parameter space in Table 1 for
the $h_{1}$ and $h_{2}$ scenario using the MultiNest algorithm with
${\it{nlive}}=8000$ [126]555Nlive in the MultiNest method represents the
number of active or live points to determine the iso-likelihood contour in
each iteration [126, 127]. The larger it is, the more detailed the scan
results will be in surveying the parameter space. and the likelihood function
in Eq. (3.3). Regarding the acquired samples, we were particularly interested
in those that can explain the relic abundance at the $2\sigma$ level, satisfy
the LZ bounds, and be consistent with all the experimental constraints. We
will analyze their features using statistical quantities and illustrate the
differences of the two scenarios.

Finally, we assessed the alignment of the samples with the results of the LHC
search for electroweakinos. To reduce the calculation time, we initially
utilized the program SModelS-2.1.1 [128], which encoded various event-
selection efficiencies by the topologies of the electroweakino signals listed
in Table 2 to implement the task. Given that the exclusion capability of
SModelS-2.1.1 on the samples is limited by its database and strict working
prerequisites, we further surveyed the surviving samples by simulating the
analyses listed in Table 3. We finished this task by the following procedures.
We concentrated on the electroweakino production processes

$\begin{split}pp\to\tilde{\chi}_{i}^{0}\tilde{\chi}_{j}^{\pm}&,\quad
i=2,3,4,5;\quad j=1,2\\\
pp\to\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{\mp}&,\quad i,j=1,2\\\
pp\to\tilde{\chi}_{i}^{0}\tilde{\chi}_{j}^{0}&,\quad i,j=2,3,4,5\end{split}$
(3.4)

and calculated their cross sections to next-to-leading order using the program
Prospino2 [129]. Subsequently, we generated $10^{5}$ events for these
processes by MadGraph_aMC@NLO [130, 131] and furnished their parton shower and
hadronization by the program PYTHIA8 [132]. Detector simulations were
implemented with the program Delphes [133]. Finally, we fed the event files
into the package CheckMATE-2.0.29 [134, 135, 136] to calculate the $R$ value
defined by $R\equiv max\\{S_{i}/S_{i,obs}^{95}\\}$ for all the involved
analyses, where $S_{i}$ represents the simulated event number of the $i$-th
signal region (SR), and $S_{i,obs}^{95}$ is the corresponding $95\%$
confidence level upper limit. $R>1$ means that the sample is experimentally
excluded if the involved uncertainties are neglected [137], while $R<1$
indicates that it is consistent with the experimental analyses.

###  3.2 Global features

Table 4: Primary annihilation channels for Singlino-dominated dark matter (DM)
in the $h_{1}$ scenario and their marginal posterior probability density
functions (PDFs) normalized to the total Bayesian evidence of the scenario.
The second line takes into account all 21097 samples acquired from the scan,
while the third line only includes those further satisfying the LHC
constraints, which resulted in smaller PDFs than their corresponding ones in
the second line. Regarding the co-annihilation in this table, the annihilation
partners were Higgsino-like electroweakinos in most cases and Wino-like
electroweakinos in a few cases, as indicated by the samples.

$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ | $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$ | Co-annihilation | $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$
---|---|---|---
$58.1\%$ | $34.8\%$ | $2.5\%$ | $2.3\%$
$58.0\%$ | $32.7\%$ | $2.4\%$ | $2.3\%$

Table 5: Similar to Table 4, but presenting the $h_{2}$ scenario results. The
channel $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to q\bar{q}$ in this table
proceeded primarily through exchanging a resonant $h_{s}$, $A_{s}$, or $h$.

$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$ | $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ | Co-annihilation | $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to q\bar{q}$
---|---|---|---
$42.3\%$ | $12.6\%$ | $5.9\%$ | $2.1\%$
$18.3\%$ | $9.7\%$ | $4.1\%$ | $1.4\%$

Figure 1:  One-dimensional PLs (red line) and marginal posterior PDFs (black
line) of $\lambda$, $\kappa$, and $\tan\beta$ before implementing the LHC
constraints. The left panels show the results of the $h_{1}$ scenario, while
the right panels depict those of the $h_{2}$ scenario. Statistical measures
such as the confidence interval and credible region were introduced in Ref.
[89].

Figure 2: Same as Fig. 1, but for the distributions of
$m_{\tilde{\chi}_{1}^{0}}$ and $\mu_{\rm tot}$.

By analyzing the samples acquired by the scans, we made the following
observations:

  * •

We obtained 21229 samples for the $h_{1}$ scenario and 20511 samples for the
$h_{2}$ scenario. These sample sets accounted for the relic abundance at the
$2\sigma$ level, remained consistent with the LZ bounds, and satisfied all the
experimental constraints encoded in $\mathcal{L}_{\rm Const}$. The Bayesian
evidences were $\ln Z_{h_{1}}=-11.28\pm 0.037$ and $\ln Z_{h_{2}}=-12.12\pm
0.038$. The Jeffrey’s scale, defined as the ratio of $\ln Z_{h_{1}}$ to $\ln
Z_{h_{2}}$, was equal to 0.84; thus, the $h_{1}$ scenario was better suited to
account for experimental results than the $h_{2}$ scenario [42]. The primary
reason was that the $h_{2}$ scenario introduced additional non-trivial
theoretical constraints, specifically $m_{B}\lesssim 130~{}{\rm GeV}$ and
$\lambda\mu_{\rm tot}\lesssim 30~{}{\rm GeV}$, based on the Higgs data fit of
this analysis, leading to a more stringent limitation on its parameter space.

  * •

The DM candidate may have been dominated by the Singlino or Bino field in its
components. In the $h_{1}$ scenario, 21097 out of 21229 samples predicted a
Singlino-dominated $\tilde{\chi}_{1}^{0}$, contributing to $99.3\%$ of the
Bayesian evidence. In contrast, 13113 out of 20511 samples in the $h_{2}$
scenario predicted a Singlino-dominated $\tilde{\chi}_{1}^{0}$, accounting for
approximately $64.8\%$ of the Bayesian evidence.

  * •

After including constraints from the LHC search for supersymmetry, the sample
numbers were reduced to 20709 and 12604 for the scenarios $h_{1}$ and $h_{2}$,
respectively. Correspondingly, the Bayesian evidences decreased by $2.4\%$ and
$38.8\%$, respectively. These results indicated that the LHC constraint had
significant impacts on the DM physics of the $h_{2}$ scenario, but hardly
affected that of the $h_{1}$ scenario.

  * •

As discussed in Section 2.3.1, the relative importance between the
annihilation channels $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}h_{s},A_{s}A_{s},h_{s}A_{s}$ for the Singlino-dominated DM depended on
the mass spectrum and couplings of the samples. We have categorized the
samples based on their dominant annihilation processes and evaluated their
impact on the Bayesian evidence. Tables 4 and 5 detail their contributions to
the total Bayesian evidence for the $h_{1}$ and $h_{2}$ scenario,
respectively, both before and after incorporating the LHC constraints. These
tables reveal that most samples in the $h_{1}$ scenario predominantly
annihilated through $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$,
aligning with the observed relic abundance. In contrast, the $h_{2}$ scenario
showed that a substantial fraction of the samples primarily annihilated via
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$, facing stringent
restrictions from the LHC search for the supersymmetry.

  * •

In the $h_{2}$ scenario, a significant fraction of samples predicted a Bino-
dominated $\tilde{\chi}_{1}^{0}$, accounting for $35.2\%$ and $28.0\%$ of the
total Bayesian evidence before and after taking the LHC constraints into
consideration, respectively. Such a DM achieved the observed relic abundance
primarily through co-annihilating with Wino-like electroweakinos. It preferred
a negative value of $M_{1}$ to suppress the SI DM-nucleon scattering by
canceling different contributions [15]. The LZ experiment alone set a lower
bound of approximately $380~{}{\rm GeV}$ on $\mu_{\rm tot}$ [19], and the LHC
search for the compressed spectrum of supersymmetry by leptonic signals
required $|m_{\tilde{\chi}_{1}^{0}}|\gtrsim 220~{}{\rm GeV}$ [118]. Our
results confirmed these characteristics.

These aforementioned conclusions indicated that the GNMSSM preferred a
Singlino-dominated DM, instead of a Bino-dominated DM, across a wide range of
parameters. This feature distinguished it from the MSSM and
$\mathbb{Z}_{3}$-NMSSM.

Figure 3: Same as Fig. 2, but for the results of the $h_{2}$ scenario after
considering the LHC constraints.

In the analysis discussed below, we studied the statistical distributions of
the input parameters that determined the property of $\tilde{\chi}_{1}^{0}$.
We plotted one-dimensional PLs and marginal posterior PDFs of $\lambda$,
$\kappa$, and $\tan\beta$ in Fig. 1 and those of $m_{\tilde{\chi}_{1}^{0}}$
and $\mu_{\rm tot}$ in Fig. 2. The left panels in both figures show the
results of the $h_{1}$ scenario without considering the LHC constraints, while
the right panels depict those of the $h_{2}$ scenario. From these plots,
several insights can be inferred:

  * •

$\lambda$’s PL and posterior PDF distributions indicated that $\lambda\lesssim
0.09$ for almost all samples in the $h_{1}$ scenario and $\lambda\lesssim
0.05$ for almost all samples in the $h_{2}$ scenario. These ranges are notably
narrower than those listed in Table 1, and the small values of $\lambda$ are a
distinctive feature of the GNMSSM. This trend primarily stemmed from the LZ
constraint discussed in Section 2: in scenarios where $\tilde{\chi}_{1}^{0}$
is Singlino dominated, the main contribution to $\sigma^{\rm
SI}_{\tilde{\chi}_{1}^{0}-N}$ is proportional to $\lambda^{2}\kappa^{2}$ for
substantial singlet–doublet Higgs mixing and $\lambda^{4}$ in the heavy
$h_{s}$ limit. In contrast, when $\tilde{\chi}_{1}^{0}$ is Bino dominated, as
significant in the $h_{2}$ scenario, the LZ constraint required $\mu_{\rm
tot}$ to exceed $380~{}{\rm GeV}$, and the Higgs data fit limited
$\lambda\mu_{\rm tot}$ in ${\cal{M}}_{\rm S,23}^{2}$ to approximately
$30~{}{\rm GeV}$ or less. In any scenario, a small $\lambda$ is preferred.

  * •

The value of $\kappa$ varied widely from $-0.75$ to $0.75$. Notably, in the
$h_{1}$ scenario, there were gaps around $\kappa=0$ in the PL and posterior
PDF distributions, while they were continuous in the $h_{2}$ scenario. This
distinction arose because the Singlino component predominantly contributed to
the DM candidate for nearly all samples in the $h_{1}$ scenario, necessitating
a non-zero $\kappa$ to account for the observed relic abundance. In contrast,
in the $h_{2}$ scenario, a significant portion of samples featured Bino-
dominated DM, whose properties were less influenced by $\kappa$.

  * •

For both the $h_{1}$ and $h_{2}$ scenarios, $\tan\beta$ spanned a broad range
from 5 to 60. The distributions did not show a significant preference for
specific regions, suggesting that the DM physics are not sensitive to
$\tan\beta$.

  * •

The PL distributions suggest that $\mu_{\rm tot}$ can range from $100~{}{\rm
GeV}$ to $1~{}{\rm TeV}$ in both scenarios. However, $\mu_{\rm tot}$’s
posterior PDF distributions indicated a preference for smaller $\mu_{\rm tot}$
in the $h_{2}$ scenario. This preference was driven by the constraints from
the Higgs data fit, which limited $\lambda\mu_{\rm tot}$ in ${\cal{M}}_{\rm
S,23}^{2}$ to be no more than approximately $30~{}{\rm GeV}$, a restriction
not applicable to the $h_{1}$ scenario.

  * •

Similarly, a relatively light DM was favored in the $h_{2}$ scenario, stemming
from the condition $|m_{\tilde{\chi}_{1}^{0}}|<\mu_{\rm tot}$ and the observed
preference for a smaller $\mu_{\rm tot}$ in this scenario. Additionally, a
light DM is self-consistent with the annihilation
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$, where $h_{s}$ could
be as light as several GeV in the $h_{2}$ scenario.

We also investigated the impact of the LHC search for supersymmetry on the
distributions. In the $h_{1}$ scenario, we observed that it only shifted the
lower bounds of $m_{\tilde{\chi}_{1}^{0}}$ and $\mu_{\rm tot}$ in the PL
distribution from $34$ and $105~{}{\rm GeV}$ to approximately $48$ and
$165~{}{\rm GeV}$, respectively. In contrast, in the $h_{2}$ scenario, apart
from altering the lower bounds on $m_{\tilde{\chi}_{1}^{0}}$ and $\mu_{\rm
tot}$ from $14$ and $100~{}{\rm GeV}$ to approximately $40$ and $155~{}{\rm
GeV}$, respectively, the LHC constraint significantly modified their posterior
PDFs by shifting the peaks of their distribution toward higher mass values.
The updated results are shown in Fig. 3 compared to their corresponding ones
in Fig. 2. This phenomenon arose due to the crucial nature of the LHC
constraints when the DM and Higgsinos had moderately low masses. This effect
was less pronounced on other distributions of the $h_{2}$ scenario.

Figure 4:  Two-dimensional distributions of the marginal posterior PDF for the
likelihood $\mathcal{L}$ in Eq. 3.3, projected onto the
$|m_{\tilde{\chi}_{1}^{0}}|-\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ and
$|m_{\tilde{\chi}_{1}^{0}}|-\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-n}$ planes.
The left panels show the results of the $h_{1}$ scenario, while the right
panels show those of the $h_{2}$ scenario. The solid and dashed lines surround
the $1\sigma$ and $2\sigma$ credible regions, respectively, and the samples
falling within these regions contributed $65.3\%$ and $95.4\%$ of the total
Bayesian evidence [89]. All the considered samples were consistent with the
LHC constraint.

Figure 5: Scattering plots of the samples that primarily annihilated through
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s},A_{s}A_{s},h_{s}A_{s}$
channels to achieve the observed relic abundance, projected onto the
$|m_{\tilde{\chi}_{1}^{0}}|-m_{h_{s}}$,
$|m_{\tilde{\chi}_{1}^{0}}|-m_{A_{s}}$, and
$|m_{\tilde{\chi}_{1}^{0}}|-m_{h_{s}}$ planes, respectively. The color bars on
the left panels denote the values of $|\kappa|$, while those on the right
panels represent the values of $\ln\lambda$.

Figure 6: Similar to Fig. 5, except that the results are shown on the
$|m_{\tilde{\chi}_{1}^{0}}|-\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ (left
panels) and $|m_{\tilde{\chi}_{1}^{0}}|-\mu_{\rm tot}$ (right panels) planes.
Samples marked by gray were excluded by the LHC constraint, while those in
blue were consistent with all the experimental constraints.

Finally, we projected the samples passing the LHC constraint onto the
$|m_{\tilde{\chi}_{1}^{0}}|-\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ and
$|m_{\tilde{\chi}_{1}^{0}}|-\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ planes
to obtain two-dimensional posterior PDF maps, which are shown in Fig. 4. This
figure reveals that the samples in the $2\sigma$ credible region predicted
$1.0\times 10^{-49}\lesssim\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}/{\rm
cm^{2}}\lesssim 2.5\times 10^{-46}$ and $1.0\times 10^{-50}\lesssim\sigma^{\rm
SD}_{\tilde{\chi}_{1}^{0}-n}/{\rm cm^{2}}\lesssim 3.2\times 10^{-43}$ in the
$h_{1}$ scenario. In contrast, analyzing the $h_{2}$ scenario yielded
$5.0\times 10^{-49}\lesssim\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}/{\rm
cm^{2}}\lesssim 2.0\times 10^{-46}$ and $1.0\times 10^{-52}\lesssim\sigma^{\rm
SD}_{\tilde{\chi}_{1}^{0}-n}/{\rm cm^{2}}\lesssim 1.0\times 10^{-43}$ for
Singlino-like DM and $3.0\times 10^{-49}\lesssim\sigma^{\rm
SI}_{\tilde{\chi}_{1}^{0}-p}/{\rm cm^{2}}\lesssim 2.0\times 10^{-46}$,
$1.0\times 10^{-43}\lesssim\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}/{\rm
cm^{2}}\lesssim 5.0\times 10^{-42}$ for Bino-like DM. These results showed
significant difference of $\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ for
Singlino- and Bino-like DM, which was explained in Eqs. (2.41) and (2.42).
Additionally, they revealed that with the future increase in experimental
sensitivities to $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$, a significant
portion of the samples would become disfavored. Numerically speaking, we
observed that the Bayesian evidences decreased by $46\%$ and $62\%$ for the
$h_{1}$ and $h_{2}$ scenarios, respectively, if the sensitivities were
improved by a factor of 5. Consequently, the upper bounds of $\lambda$ were
reduced to approximately $0.05$ and $0.02$, respectively.

###  3.3 Three annihilation channels

As introduced in Section 2.3, the relic abundance and nucleon scattering cross
sections of the Singlino-dominated DM depend not only on parameters within the
neutralino sector, including $\lambda$, $\kappa$, $\tan\beta$, $\mu_{\rm
tot}$, and $m_{\tilde{\chi}_{1}^{0}}$, but also on those within the Higgs
sector, such as $v_{s}$, $A_{\lambda}$, $m_{B}$ (or equivalently $m_{h_{s}}$),
and $m_{C}$ ($m_{A_{s}}$). Consequently, the characteristics of the DM are
complex. Nevertheless, valuable insights can still be obtained from scattering
plots derived from our sample sets. We categorized the samples based on their
predominant annihilation processes for this study’s purposes. We particularly
focus on samples where the primary annihilation mechanisms were
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$,
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$, and
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$, and label them as
Type-I, Type-II, and Type-III samples, respectively. Their respective counts
were approximately 8000, 8000, and 16000 samples obtained by the scans, which
created a challenge in plot representations and complicated our analyses due
to their sheer volume. To address this issue, we applied a method introduced
in Ref. [19] to reduce the sample density for the first two types of samples,
while for the latter type, we controlled its sample size by implementing a
mass constraint of $0.8m_{h_{s}}\lesssim m_{A_{s}}\lesssim 1.2m_{h_{s}}$ to
simplify our analysis. By comparing plots derived from the full sample sets
with those from the reduced sample sets, we observed that the latter
efficiently captured essential features of the annihilations without losing
the generality of our findings.

Point P1 | Point P2
---|---
$\lambda$ | 0.012 | $m_{h_{s}}$ | 269.8 GeV | $\lambda$ | 0.003 | $m_{h_{s}}$ | 119.3 GeV
$\kappa$ | -0.40 | $m_{h}$ | 125.6 GeV | $\kappa$ | 0.56 | $m_{h}$ | 124.9 GeV
$\tan{\beta}$ | 16.58 | $m_{H}$ | 1826 GeV | $\tan{\beta}$ | 11.75 | $m_{H}$ | 1871 GeV
$v_{s}$ | 692.3 GeV | $m_{A_{s}}$ | 737.5 GeV | $v_{s}$ | 609.2 GeV | $m_{A_{s}}$ | 735.6 GeV
$\mu_{\rm tot}$ | 390.4 GeV | $m_{A_{H}}$ | 1826 GeV | $\mu_{\rm tot}$ | 405.2 GeV | $m_{A_{H}}$ | 1871 GeV
$M_{1}$ | -747.5 GeV | $m_{H^{\pm}}$ | 1827 GeV | $M_{1}$ | -718.8 GeV | $m_{H^{\pm}}$ | 1872 GeV
$M_{2}$ | 822.5 GeV | $m_{\tilde{\chi}_{1}^{0}}$ | 324.2 GeV | $M_{2}$ | 614.6 GeV | $m_{\tilde{\chi}_{1}^{0}}$ | 339 GeV
$A_{t}$ | -3983 GeV | $m_{\tilde{\chi}_{2}^{0}}$ | 393.2 GeV | $A_{t}$ | -4233 GeV | $m_{\tilde{\chi}_{2}^{0}}$ | 400.5 GeV
$A_{\lambda}$ | 1654 GeV | $m_{\tilde{\chi}_{3}^{0}}$ | -400 GeV | $A_{\lambda}$ | 1562 GeV | $m_{\tilde{\chi}_{3}^{0}}$ | -414.6 GeV
$m_{A}$ | 2000 GeV | $m_{\tilde{\chi}_{4}^{0}}$ | -755.3 GeV | $m_{A}$ | 2000 GeV | $m_{\tilde{\chi}_{4}^{0}}$ | 660.3 GeV
$m_{B}$ | 280.8 GeV | $m_{\tilde{\chi}_{5}^{0}}$ | 865.1 GeV | $m_{B}$ | 109.5 GeV | $m_{\tilde{\chi}_{5}^{0}}$ | -726.7 GeV
$m_{N}$ | 326.3 GeV | $m_{\tilde{\chi}_{1}^{\pm}}$ | 395 GeV | $m_{N}$ | 345.1 GeV | $m_{\tilde{\chi}_{1}^{\pm}}$ | 402.5 GeV
$\Omega h^{2}$ | 0.1199 | $m_{\tilde{\chi}_{2}^{\pm}}$ | 865.3 GeV | $\Omega h^{2}$ | 0.12 | $m_{\tilde{\chi}_{2}^{\pm}}$ | 660.6 GeV
$\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ | 1.73$\times 10^{-46}{\rm~{}cm^{2}}$ | $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ | 2.24$\times 10^{-47}{\rm~{}cm^{2}}$ | $\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ | 7.42$\times 10^{-49}{\rm~{}cm^{2}}$ | $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ | 7.19$\times 10^{-47}{\rm~{}cm^{2}}$
$V_{h_{s}}^{\rm S},~{}V_{h_{s}}^{\rm SM},~{}V_{h}^{\rm S},~{}V_{h}^{\rm SM}$ | -0.999, -0.019, 0.019, -0.999 | $V_{h_{s}}^{\rm S},~{}V_{h_{s}}^{\rm SM},~{}V_{h}^{\rm S},~{}V_{h}^{\rm SM}$ | 0.982, -0.190, 0.190, 0.982
$N_{11},~{}N_{12},~{}N_{13},~{}N_{14},~{}N_{15}$ | 0.001, 0.002, -0.012, 0.015, -0.999 | $N_{11},~{}N_{12},~{}N_{13},~{}N_{14},~{}N_{15}$ | 0.001, 0.001, -0.003, 0.004, -0.999
$N_{21},~{}N_{22},~{}N_{23},~{}N_{24},~{}N_{25}$ | -0.028, -0.126, 0.706, -0.696, -0.020 | $N_{21},~{}N_{22},~{}N_{23},~{}N_{24},~{}N_{25}$ | -0.028, -0.239, 0.696, -0.676, -0.005
$N_{31},~{}N_{32},~{}N_{33},~{}N_{34},~{}N_{35}$ | -0.081, -0.042, -0.704, -0.704, -0.002 | $N_{31},~{}N_{32},~{}N_{33},~{}N_{34},~{}N_{35}$ | -0.090, -0.048, -0.703, -0.704, -0.001
$N_{41},~{}N_{42},~{}N_{43},~{}N_{44},~{}N_{45}$ | 0.996, -0.004, -0.037, -0.078, -0.001 | $N_{41},~{}N_{42},~{}N_{43},~{}N_{44},~{}N_{45}$ | 0.007, -0.970, -0.136, 0.202, 0.001
$N_{51},~{}N_{52},~{}N_{53},~{}N_{54},~{}N_{55}$ | 0.003, -0.991, -0.060, 0.118, 0.001 | $N_{51},~{}N_{52},~{}N_{53},~{}N_{54},~{}N_{55}$ | 0.996, -0.005, -0.043, -0.084, -0.001
annihilations | Fractions [%] | annihilations | Fractions [%]
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$ | 99.9 | $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}/h_{s}h$ | 92.7/7.02
Decays | Branching ratios [%] | Decays | Branching ratios [%]
$\tilde{\chi}^{0}_{2}\to\tilde{\chi}^{0}_{1}Z^{\ast}$ | 100 | $\tilde{\chi}^{0}_{2}\to\tilde{\chi}^{0}_{1}Z^{\ast}$ | 100
$\tilde{\chi}^{0}_{3}\to\tilde{\chi}^{0}_{1}Z^{\ast}$ | 96.5 | $\tilde{\chi}^{0}_{3}\to\tilde{\chi}^{\mp}_{1}(W^{\pm})^{\ast}/\tilde{\chi}^{0}_{2}Z^{\ast}/\tilde{\chi}^{0}_{1}Z^{\ast}$ | 51.8/40.7/8.5
$\tilde{\chi}^{0}_{4}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{3}h/\tilde{\chi}^{0}_{2}Z/\tilde{\chi}^{0}_{3}Z$ | 50/22.5/21.5/3 | $\tilde{\chi}^{0}_{4}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{3}Z/\tilde{\chi}^{0}_{2}h$ | 55.8/22.5/19.2
$\tilde{\chi}^{0}_{5}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{3}Z/\tilde{\chi}^{0}_{2}h$ | 51.7/22.2/21 | $\tilde{\chi}^{0}_{5}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{5}h/\tilde{\chi}^{0}_{2}Z$ | 50.4/21.8/21.6
$\tilde{\chi}^{+}_{1}\to\tilde{\chi}^{0}_{1}(W^{+})^{\ast}$ | 100 | $\tilde{\chi}^{+}_{1}\to\tilde{\chi}^{0}_{1}(W^{+})^{\ast}$ | 100
$\tilde{\chi}^{+}_{2}\to\tilde{\chi}^{0}_{2}W^{+}/\tilde{\chi}^{0}_{3}W^{+}/\tilde{\chi}^{+}_{1}Z/\tilde{\chi}^{+}_{1}h/$ | 25.4/25.4/25.1/23.7 | $\tilde{\chi}^{+}_{2}\to\tilde{\chi}^{0}_{2}W^{+}/\tilde{\chi}^{+}_{1}Z/\tilde{\chi}^{0}_{3}W^{+}/\tilde{\chi}^{+}_{1}h$ | 27.3/25.8/25.6/20.5
$h_{s}\to W^{+}W^{-}/ZZ/hh{b}$ | 56.4/23/20.4 | $h_{s}\to b\bar{b}/WW^{\ast}/gg/\tau^{+}\tau^{-}$ | 62.9/16.2/9.1/7.2
$h\to b\bar{b}/WW^{\ast}/gg/\tau^{+}\tau^{-}$ | 53.7/26.6/8.2/6.2 | $h\to b\bar{b}/WW^{\ast}/gg/\tau^{+}\tau^{-}/c\bar{c}$ | 54.9/25.2/8.3/6.4/2.5
$H\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{0}_{5}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{5}/\tau^{+}\tau^{-}$ | 49.3/26.2/13.7/6/4.9 | $H\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{0}_{4}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{4}/\tilde{\chi}^{+}_{1}\tilde{\chi}^{-}_{1}$ | 46.6/15.9/12.4/7/3.7
$A_{H}\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{5}$ | 38.5/26.3/13.5 | $A_{H}\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{4}/\tilde{\chi}^{+}_{1}\tilde{\chi}^{-}_{1}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{5}$ | 42.9/14.1/12.4/5.8/2.9
$H^{+}\to t\bar{b}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{+}_{2}/\tilde{\chi}^{0}_{5}\tilde{\chi}^{+}_{5}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{+}_{2}/\tilde{\chi}^{0}_{4}\tilde{\chi}^{+}_{1}$ | 27.2/20.3/19.7/19.3/6.8 | $H^{+}\to\tilde{\chi}^{0}_{2}\tilde{\chi}^{+}_{2}/\tilde{\chi}^{0}_{4}\tilde{\chi}^{+}_{1}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{+}_{2}/t\bar{b}/\tilde{\chi}^{0}_{5}\tilde{\chi}^{+}_{1}$ | 25.9/24.4/22.9/14.8/6.3
$R$ value: 0.16 | Signal Region: SR-incWh-SFOS-09 in Ref. [118] |  $R$ value: 0.33 | Signal Region: SR-G05 in Ref. [120]

Table 6: Benchmark points where the Singlino-dominated DM achieved the
observed relic abundance primarily through the process
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$. Both points
satisfied all experimental constraints. In the left and right parts of this
table, $h_{1}$ and $h_{2}$ were predicted to be the SM-like Higgs boson,
respectively.

Point P3 | Point P4
---|---
$\lambda$ | 0.006 | $m_{h_{s}}$ | 818.5 GeV | $\lambda$ | 0.002 | $m_{h_{s}}$ | 98.97 GeV
$\kappa$ | 0.45 | $m_{h}$ | 124.7 GeV | $\kappa$ | 0.13 | $m_{h}$ | 125.6 GeV
$\tan{\beta}$ | 26.64 | $m_{H}$ | 1670 GeV | $\tan{\beta}$ | 25.25 | $m_{H}$ | 1730 GeV
$v_{s}$ | 886.7 GeV | $m_{A_{s}}$ | 119.5 GeV | $v_{s}$ | 285.8 GeV | $m_{A_{s}}$ | 50.57 GeV
$\mu_{\rm tot}$ | 569 GeV | $m_{A_{H}}$ | 1670 GeV | $\mu_{\rm tot}$ | 428.2 GeV | $m_{A_{H}}$ | 1730 GeV
$M_{1}$ | 692.6 GeV | $m_{H^{\pm}}$ | 1665 GeV | $M_{1}$ | 732.2 GeV | $m_{H^{\pm}}$ | 1728 GeV
$M_{2}$ | 303.9 GeV | $m_{\tilde{\chi}_{1}^{0}}$ | 217.5 GeV | $M_{2}$ | 656.5 GeV | $m_{\tilde{\chi}_{1}^{0}}$ | -65.31 GeV
$A_{t}$ | -2748 GeV | $m_{\tilde{\chi}_{2}^{0}}$ | 314.8 GeV | $A_{t}$ | -3168 GeV | $m_{\tilde{\chi}_{2}^{0}}$ | 422.4 GeV
$A_{\lambda}$ | -1307 GeV | $m_{\tilde{\chi}_{3}^{0}}$ | 583.5 GeV | $A_{\lambda}$ | -731.2 GeV | $m_{\tilde{\chi}_{3}^{0}}$ | -441.9 GeV
$m_{A}$ | 2000 GeV | $m_{\tilde{\chi}_{4}^{0}}$ | -584.1 GeV | $m_{A}$ | 2000 GeV | $m_{\tilde{\chi}_{4}^{0}}$ | 699.8 GeV
$m_{B}$ | 834.7 GeV | $m_{\tilde{\chi}_{5}^{0}}$ | 706.1 GeV | $m_{B}$ | 101 GeV | $m_{\tilde{\chi}_{5}^{0}}$ | 742 GeV
$m_{N}$ | 217.5 GeV | $m_{\tilde{\chi}_{1}^{\pm}}$ | 315.2 GeV | $m_{N}$ | -65.37 GeV | $m_{\tilde{\chi}_{1}^{\pm}}$ | 428.5 GeV
$\Omega h^{2}$ | 0.12 | $m_{\tilde{\chi}_{2}^{\pm}}$ | 596.2 GeV | $\Omega h^{2}$ | 0.1119 | $m_{\tilde{\chi}_{2}^{\pm}}$ | 701.9 GeV
$\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ | 5.37$\times 10^{-49}{\rm~{}cm^{2}}$ | $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-n}$ | 7.86$\times 10^{-49}{\rm~{}cm^{2}}$ | $\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ | 6.11$\times 10^{-51}{\rm~{}cm^{2}}$ | $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-n}$ | 1.03$\times 10^{-47}{\rm~{}cm^{2}}$
$V_{h_{s}}^{\rm S},~{}V_{h_{s}}^{\rm SM},~{}V_{h}^{\rm S},~{}V_{h}^{\rm SM}$ | 0.999, 0.002, -0.002, 0.999 | $V_{h_{s}}^{\rm S},~{}V_{h_{s}}^{\rm SM},~{}V_{h}^{\rm S},~{}V_{h}^{\rm SM}$ | -0.999, 0.043, -0.043, -0.999
$N_{11},~{}N_{12},~{}N_{13},~{}N_{14},~{}N_{15}$ | 0.001, -0.002, 0.001, -0.002, 0.999 | $N_{11},~{}N_{12},~{}N_{13},~{}N_{14},~{}N_{15}$ | -0.001, 0.001, 0.001, 0.001, -0.999
$N_{21},~{}N_{22},~{}N_{23},~{}N_{24},~{}N_{25}$ | 0.013, -0.975, 0.192, -0.108, -0.002 | $N_{21},~{}N_{22},~{}N_{23},~{}N_{24},~{}N_{25}$ | 0.098, -0.212, 0.699, -0.676, -0.001
$N_{31},~{}N_{32},~{}N_{33},~{}N_{34},~{}N_{35}$ | 0.264, 0.208, 0.664, -0.668, -0.002 | $N_{31},~{}N_{32},~{}N_{33},~{}N_{34},~{}N_{35}$ | -0.025, 0.048, 0.703, 0.709, 0.001
$N_{41},~{}N_{42},~{}N_{43},~{}N_{44},~{}N_{45}$ | -0.023, 0.059, 0.703, 0.708, 0.001 | $N_{41},~{}N_{42},~{}N_{43},~{}N_{44},~{}N_{45}$ | 0.173, 0.965, 0.108, -0.166, -0.001
$N_{51},~{}N_{52},~{}N_{53},~{}N_{54},~{}N_{55}$ | 0.964, -0.042, -0.167, 0.201, 0.001 | $N_{51},~{}N_{52},~{}N_{53},~{}N_{54},~{}N_{55}$ | 0.980, -0.148, -0.071, 0.115, 0.001
annihilations | Fractions [%] | annihilations | Fractions [%]
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$ | 99.8 | $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$ | 99.9
Decays | Branching ratios [%] | Decays | Branching ratios [%]
$\tilde{\chi}^{0}_{2}\to\tilde{\chi}^{0}_{1}Z$ | 100 | $\tilde{\chi}^{0}_{2}\to\tilde{\chi}^{0}_{1}Z/\tilde{\chi}^{0}_{1}h$ | 62.2/37.2
$\tilde{\chi}^{0}_{3}\to\tilde{\chi}^{\mp}_{1}W^{\pm}/\tilde{\chi}^{0}_{2}h/\tilde{\chi}^{0}_{2}Z$ | 69.5/27.7/2.8 | $\tilde{\chi}^{0}_{3}\to\tilde{\chi}^{0}_{1}h/\tilde{\chi}^{0}_{1}Z/\tilde{\chi}^{0}_{2}Z^{\ast}/\tilde{\chi}^{\mp}_{1}(W^{\pm})^{\ast}$ | 45.3/29.5/14.6/5.8
$\tilde{\chi}^{0}_{4}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{2}Z/\tilde{\chi}^{0}_{2}h$ | 67.8/30.1/1.7 | $\tilde{\chi}^{0}_{4}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{3}Z/\tilde{\chi}^{0}_{2}h$ | 65.1/17.7/15.5
$\tilde{\chi}^{0}_{5}\to\tilde{\chi}^{\pm}_{2}W^{\mp}/\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{4}Z/\tilde{\chi}^{0}_{2}Z$ | 37.1/29.7/19.4/1.8 | $\tilde{\chi}^{0}_{5}\to\tilde{\chi}^{0}_{3}Z/\tilde{\chi}^{0}_{2}h/\tilde{\chi}^{\pm}_{2}W^{\mp}$ | 38.4/34.8/22.2
$\tilde{\chi}^{+}_{1}\to\tilde{\chi}^{0}_{1}W^{+}$ | 100 | $\tilde{\chi}^{+}_{1}\to\tilde{\chi}^{0}_{1}W^{+}$ | 100
$\tilde{\chi}^{+}_{2}\to\tilde{\chi}^{0}_{2}W^{+}/\tilde{\chi}^{+}_{1}Z/\tilde{\chi}^{+}_{1}h$ | 35.6/34.5/29.9 | $\tilde{\chi}^{+}_{2}\to\tilde{\chi}^{0}_{2}W^{+}/\tilde{\chi}^{+}_{1}Z/\tilde{\chi}^{0}_{3}W^{+}/\tilde{\chi}^{+}_{1}h$ | 26.4/26/25.4/22.1
$h_{s}\to A_{s}A_{s}/\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}$ | 84.4/15.6 | $h_{s}\to b\bar{b}/\tau^{+}\tau^{-}/gg/c\bar{c}$ | 79.6/8.9/7.2/3.4
$h\to b\bar{b}/WW^{\ast}/gg/\tau^{+}\tau^{-}$ | 55.1/24.9/8.4/6.4 | $h\to b\bar{b}/A_{s}A_{s}/WW^{\ast}/gg/\tau^{+}\tau^{-}$ | 39.3/26.6/19.6/6.1/4.6
$H\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{4}/\tau^{+}\tau^{-}/\tilde{\chi}^{0}_{5}\tilde{\chi}^{0}_{5}$ | 39.8/30/9.3/7.7/5.7 | $H\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{0}_{4}\tau^{+}\tau^{-}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{0}_{5}$ | 41.2/28.6/7.8/7.8/4.5
$A_{H}\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{3}/\tau^{+}\tau^{-}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{4}$ | 39.8/28.9/8.3/7.7/5.8 | $A_{H}\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tau^{+}\tau^{-}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{4}$ | 41.3/26.8/7.8/6.9
$H^{+}\to t\bar{b}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{+}_{1}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{+}_{2}/\tilde{\chi}^{0}_{4}\tilde{\chi}^{+}_{1}$ | 39.5/16.4/16.3/15.4 | $H^{+}\to t\bar{b}/\tilde{\chi}^{0}_{4}\tilde{\chi}^{+}_{1}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{+}_{2}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{+}_{2}/\tau^{+}\nu_{\tau}$ | 41/17.6/14.2/15.3/8.6
$R$ value: 0.79 | Signal Region: SR-incWZ-03 in Ref. [118] |  $R$ value: 0.85 | Signal Region: SRG07-0j-mll in Ref. [120]

Table 7: Same as Table 6 except that the Singlino-dominated DM achieved the
observed relic abundance primarily through the process
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$.

Point P5 | Point P6
---|---
$\lambda$ | 0.026 | $m_{h_{s}}$ | 462 GeV | $\lambda$ | 0.001 | $m_{h_{s}}$ | 76.85 GeV
$\kappa$ | -0.22 | $m_{h}$ | 124.9 GeV | $\kappa$ | 0.18 | $m_{h}$ | 125.2 GeV
$\tan{\beta}$ | 42.42 | $m_{H}$ | 2389 GeV | $\tan{\beta}$ | 36.19 | $m_{H}$ | 2148 GeV
$v_{s}$ | 630.1 GeV | $m_{A_{s}}$ | 263.4 GeV | $v_{s}$ | 200.7 GeV | $m_{A_{s}}$ | 282.8 GeV
$\mu_{\rm tot}$ | 514.5 GeV | $m_{A_{H}}$ | 2389 GeV | $\mu_{\rm tot}$ | 259 GeV | $m_{A_{H}}$ | 2148 GeV
$M_{1}$ | 855.2 GeV | $m_{H^{\pm}}$ | 2410 GeV | $M_{1}$ | 704.2 GeV | $m_{H^{\pm}}$ | 2159 GeV
$M_{2}$ | 822.6 GeV | $m_{\tilde{\chi}_{1}^{0}}$ | -389.8 GeV | $M_{2}$ | 621.4 GeV | $m_{\tilde{\chi}_{1}^{0}}$ | 209.1 GeV
$A_{t}$ | 2492 GeV | $m_{\tilde{\chi}_{2}^{0}}$ | 515 GeV | $A_{t}$ | 2424 GeV | $m_{\tilde{\chi}_{2}^{0}}$ | 256.7 GeV
$A_{\lambda}$ | 1123 GeV | $m_{\tilde{\chi}_{3}^{0}}$ | -529.8 GeV | $A_{\lambda}$ | 1501 GeV | $m_{\tilde{\chi}_{3}^{0}}$ | -270.8 GeV
$m_{A}$ | 2000 GeV | $m_{\tilde{\chi}_{4}^{0}}$ | 858.4 GeV | $m_{A}$ | 2000 GeV | $m_{\tilde{\chi}_{4}^{0}}$ | 661.8 GeV
$m_{B}$ | 465.2 GeV | $m_{\tilde{\chi}_{5}^{0}}$ | 871.9 GeV | $m_{B}$ | 76.24 GeV | $m_{\tilde{\chi}_{5}^{0}}$ | 712.3 GeV
$m_{N}$ | -390.2 GeV | $m_{\tilde{\chi}_{1}^{\pm}}$ | 520.2 GeV | $m_{N}$ | 209.5 GeV | $m_{\tilde{\chi}_{1}^{\pm}}$ | 262.1 GeV
$\Omega h^{2}$ | 0.12 | $m_{\tilde{\chi}_{2}^{\pm}}$ | 867 GeV | $\Omega h^{2}$ | 0.12 | $m_{\tilde{\chi}_{2}^{\pm}}$ | 662.9 GeV
$\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ | 8.24$\times 10^{-46}{\rm~{}cm^{2}}$ | $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ | 1.9$\times 10^{-47}{\rm~{}cm^{2}}$ | $\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ | 9.89$\times 10^{-50}{\rm~{}cm^{2}}$ | $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ | 3.76$\times 10^{-48}{\rm~{}cm^{2}}$
$V_{h_{s}}^{\rm S},~{}V_{h_{s}}^{\rm SM},~{}V_{h}^{\rm S},~{}V_{h}^{\rm SM}$ | 0.999, 0.022, -0.022, 0.999 | $V_{h_{s}}^{\rm S},~{}V_{h_{s}}^{\rm SM},~{}V_{h}^{\rm S},~{}V_{h}^{\rm SM}$ | -0.999, 0.009, -0.009, -0.999
$N_{11},~{}N_{12},~{}N_{13},~{}N_{14},~{}N_{15}$ | 0.001, -0.001, -0.014, -0.019, 0.999 | $N_{11},~{}N_{12},~{}N_{13},~{}N_{14},~{}N_{15}$ | -0.001, 0.001, -0.002, 0.002, -0.999
$N_{21},~{}N_{22},~{}N_{23},~{}N_{24},~{}N_{25}$ | -0.089, 0.163, -0.702, 0.687, 0.003 | $N_{21},~{}N_{22},~{}N_{23},~{}N_{24},~{}N_{25}$ | 0.068, -0.141, 0.711, -0.685, -0.003
$N_{31},~{}N_{32},~{}N_{33},~{}N_{34},~{}N_{35}$ | 0.022, -0.040, -0.704, -0.708, -0.023 | $N_{31},~{}N_{32},~{}N_{33},~{}N_{34},~{}N_{35}$ | 0.031, -0.060, -0.700, -0.711, -0.001
$N_{41},~{}N_{42},~{}N_{43},~{}N_{44},~{}N_{45}$ | -0.564, -0.821, -0.046, 0.074, 0.001 | $N_{41},~{}N_{42},~{}N_{43},~{}N_{44},~{}N_{45}$ | 0.112, 0.983, 0.056, -0.133, -0.001
$N_{51},~{}N_{52},~{}N_{53},~{}N_{54},~{}N_{55}$ | 0.821, -0.545, -0.089, 0.145, 0.001 | $N_{51},~{}N_{52},~{}N_{53},~{}N_{54},~{}N_{55}$ | 0.991, -0.099, -0.034, 0.085, 0.001
annihilations | Fractions [%] | annihilations | Fractions [%]
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ | 99.6 | $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}/h_{s}h_{s}$ | 97.1/2.87
Decays | Branching ratios [%] | Decays | Branching ratios [%]
$\tilde{\chi}^{0}_{2}\to\tilde{\chi}^{0}_{1}Z$ | 100 | $\tilde{\chi}^{0}_{2}\to\tilde{\chi}^{0}_{1}Z^{\ast}$ | 100
$\tilde{\chi}^{0}_{3}\to\tilde{\chi}^{0}_{1}h/\tilde{\chi}^{0}_{1}Z$ | 92.5/7.38 | $\tilde{\chi}^{0}_{3}\to\tilde{\chi}^{0}_{2}Z^{\ast}/\tilde{\chi}^{\mp}_{1}(W^{\pm})^{\ast}$ | 79.5/20.4
$\tilde{\chi}^{0}_{4}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{3}Z/\tilde{\chi}^{0}_{2}h$ | 84.7/7.27/6.9 | $\tilde{\chi}^{0}_{4}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{3}Z/\tilde{\chi}^{0}_{2}h$ | 59.6/18.3/16
$\tilde{\chi}^{0}_{5}\to\tilde{\chi}^{0}_{3}Z/\tilde{\chi}^{0}_{2}h/\tilde{\chi}^{0}_{2}Z/\tilde{\chi}^{0}_{3}h$ | 47.9/45.1/3.12/2.39 | $\tilde{\chi}^{0}_{5}\to\tilde{\chi}^{\pm}_{1}W^{\mp}/\tilde{\chi}^{0}_{3}Z/\tilde{\chi}^{0}_{2}h/\tilde{\chi}^{0}_{2}Z/\tilde{\chi}^{0}_{3}h$ | 31.2/30.1/27.5/5.8/5.2
$\tilde{\chi}^{+}_{1}\to\tilde{\chi}^{0}_{1}(W^{+})^{\ast}$ | 100 | $\tilde{\chi}^{+}_{1}\to\tilde{\chi}^{0}_{1}(W^{+})^{\ast}\ /\tilde{\chi}^{0}_{2}(W^{+})^{\ast}$ | 53.2/46.8
$\tilde{\chi}^{+}_{2}\to\tilde{\chi}^{+}_{1}Z/\tilde{\chi}^{0}_{2}W^{+}/\tilde{\chi}^{0}_{3}W^{+}/\tilde{\chi}^{+}_{1}h$ | 25.5/25.4/25.3/23.6 | $\tilde{\chi}^{+}_{2}\to\tilde{\chi}^{0}_{3}W^{+}/\tilde{\chi}^{+}_{1}Z/\tilde{\chi}^{0}_{2}W^{+}/\tilde{\chi}^{+}_{1}h$ | 26/25.4/25/23.2
$h_{s}\to W^{+}W^{-}/hh/ZZ/b\bar{b}$ | 44.5/23/20/12.3 | $h_{s}\to b\bar{b}/gg/\tau^{+}\tau^{-}$ | 73/11/8
$h\to b\bar{b}/WW^{\ast}/gg/\tau^{+}\tau^{-}/c\bar{c}$ | 54.7/25.4/8.4/6.3/2.5 | $h\to b\bar{b}/WW^{\ast}/gg/\tau^{+}\tau^{-}$ | 54.2/25.9/8.4/6.3
$H\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tau^{+}\tau^{-}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{0}_{5}$ | 55.1/18.6/12.6/5.9 | $H\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tau^{+}\tau^{-}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{0}_{4}$ | 49.6/23.7/9.5/6
$A_{H}\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tau^{+}\tau^{-}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{5}$ | 55.1/18/12.6/5.5 | $A_{H}\to b\bar{b}/\tilde{\chi}^{\pm}_{1}\tilde{\chi}^{\mp}_{2}/\tau^{+}\tau^{-}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{4}$ | 49.6/23.5/9.5/5.7
$H^{+}\to t\bar{b}/\tau^{+}\nu_{\tau}/\tilde{\chi}^{0}_{4}\tilde{\chi}^{+}_{1}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{+}_{2}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{+}_{2}$ | 54.4/13.7/12.5/9.6/9.2 | $H^{+}\to t\bar{b}/\tilde{\chi}^{0}_{4}\tilde{\chi}^{+}_{1}/\tilde{\chi}^{0}_{2}\tilde{\chi}^{+}_{2}/\tilde{\chi}^{0}_{3}\tilde{\chi}^{+}_{2}$ | 49/13.5/12.5/11.7
$R$ value: 0.10 | Signal Region: 3LI in Ref. [118] |  $R$ value: 0.74 | Signal Region: SR-WZoff-high-njd in Ref. [120]

Table 8: Same as Table 6 except that the Singlino-dominated DM achieved the
observed relic abundance primarily through the process
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$.

In Fig. 5, we project three types of samples onto the
$|m_{\tilde{\chi}_{1}^{0}}|-m_{h_{s}}$,
$|m_{\tilde{\chi}_{1}^{0}}|-m_{A_{s}}$, and
$|m_{\tilde{\chi}_{1}^{0}}|-m_{h_{s}}$ planes, respectively. The color bar on
the left panels indicates the value of $|\kappa|$, which is crucial in
determining the relic abundance. On the right panel, the color bar shows the
value of $\ln\lambda$, critical for DM-nucleon scattering cross sections.
Analysis of the left panels yielded the following conclusions:

  * •

Increasing $m_{h_{s}}$ in annihilation processes
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$ and
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ or increasing
$m_{A_{s}}$ in annihilation process
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$ favors a larger
$|\kappa|$ to reach the observed relic abundance at a fixed
$|m_{\tilde{\chi}_{1}^{0}}|$. In particular, $|\kappa|$ reaches its maximum
when the scalar mass closely matches the DM mass due to the phase space
suppression in these annihilations.

  * •

An increase in $|m_{\tilde{\chi}_{1}^{0}}|$ correlates with a preference for a
higher $|\kappa|$ to achieve the observed relic abundance, which is
demonstrated by the approximations in Eqs. (2.26), (2.29), and (2.32) for
simplified cases.

  * •

Among the three annihilations, $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
h_{s}A_{s}$ requires the smallest $|\kappa|$ to achieve the observed relic
abundance due to its s-wave nature, in contrast with the p-wave nature of the
other channels. The comparison of $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to
A_{s}A_{s}$ with $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$
shows that the former prefers a larger $|\kappa|$, as explained by Eqs. (2.26)
and (2.29).

  * •

The s-wave nature of the annihilation
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ with $h_{s}\to
b\bar{b}$ and $A_{s}\to b\bar{b}$ also allows for indirect DM search
experiments to establish a lower bound on $|m_{\tilde{\chi}_{1}^{0}}|$
(approximately $58~{}{\rm GeV}$). This bound is substantially less stringent
than that of the direct annihilation
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to b\bar{b}$, as the power produced
through the former channel is suppressed by its relatively heavy DM mass [138,
139]. In contrast, no experimental limits apply to p-wave annihilations
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$ and
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$, permitting
potentially much lighter DM candidates down to ten GeV.

  * •

The complex dependence of the annihilation cross sections on the model’s
parameters suggests that all values of $\kappa$ ranging from 0.05 to 0.75 are
capable of providing an explanation for the observed relic abundance.

The information provided in the right panels is as follows:

  * •

For Type-I and -III samples, $\lambda$ is typically very small under the
$h_{2}$ scenario due to significant singlet–doublet Higgs mixing, resulting in
$\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ being proportional to
$\lambda^{2}\kappa^{2}$. Therefore, the results of the LZ experiment
stringently limit the magnitude of $\lambda$. However, a significant $\lambda$
may still be possible in rare cases where contributions mediated by $h_{s}$
cancel out those mediated by $h$ [75, 41].

  * •

As $m_{h_{s}}$ and $|m_{\tilde{\chi}_{1}^{0}}|$ increase from $300~{}{\rm
GeV}$ for Type-I and -III samples, larger $\lambda$ values may arise. In such
instances, the singlet–doublet Higgs mixing is very small, resulting in the
dominant contribution to $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ being
proportional to $\lambda^{4}$. Moreover, the contribution mediated by $h_{s}$
is gradually reduced, and simultaneously, the LZ constraint is relaxed, both
allowing for relatively large $\lambda$ values, particularly when
$\tilde{\chi}_{1}^{0}$ and $h_{s}$ are all massive.

  * •

In Type-II samples, the value of $\lambda$ can significantly increase when the
values of $|m_{\tilde{\chi}_{1}^{0}}|$ and $m_{A_{s}}$ are enhanced from their
starting point at $m_{\tilde{\chi}_{1}^{0}}\simeq 10~{}{\rm GeV}$ and
$m_{A_{s}}\simeq 10~{}{\rm GeV}$. Specifically, for this type of sample,
$h_{s}$ is always heavier than $\tilde{\chi}_{1}^{0}$, exceeding $300~{}{\rm
GeV}$ for almost all cases, which results in a minor impact on $\sigma^{\rm
SI}_{\tilde{\chi}_{1}^{0}-p}$. The behavior of $\lambda$ can be explained by
two typical mass spectrum configurations: increasing
$|m_{\tilde{\chi}_{1}^{0}}|$ while keeping $m_{A_{s}}$ constant or fixing the
ratio $m_{A_{s}}/|m_{\tilde{\chi}_{1}^{0}}|$ and increasing
$|m_{\tilde{\chi}_{1}^{0}}|$. In both cases, the lower limit of $m_{h_{s}}$
increases since $m_{h_{s}}>2m_{\tilde{\chi}_{1}^{0}}-m_{A_{s}}$, leading to
the tendency to decrease its contribution to $\sigma^{\rm
SI}_{\tilde{\chi}_{1}^{0}-p}$. Along with the relaxation of the LZ constraint,
this situation allows for larger values of $\lambda$.

  * •

Compared to Type-I and -III samples, a larger portion of Type-II samples
predict a sizable $\lambda$. The prevalence of large $\lambda$ values in Type-
II samples can be attributed to the tendency for a heavy $h_{s}$, resulting in
the dominant contribution to $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ being
proportional to $\lambda^{4}$. Additionally, Type-I samples favor smaller
values of $\lambda$ compared to Type-III samples largely due to a significant
number of samples predicting relatively small $\mu_{\rm tot}$ values.

We have also projected three types of samples onto the
$|m_{\tilde{\chi}_{1}^{0}}|-\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ and
$|m_{\tilde{\chi}_{1}^{0}}|-\mu_{\rm tot}$ planes to create Fig. 6. In this
visualization, gray indicates samples excluded by the LHC constraint, while
blue represents samples in agreement with all experimental constraints. Fig. 6
reveals the following key observations:

  * •

The LZ experiment and the LHC’s search for supersymmetry complement each other
in restricting the parameter space of the GNMSSM.

  * •

Type-II samples tend to predict a smaller $\sigma^{\rm
SI}_{\tilde{\chi}_{1}^{0}-p}$ than the other sample types.

  * •

The LHC constraint is particularly effective in ruling out samples that
predict moderately light DM and Higgsinos, notably enhancing lower limits on
$|m_{\tilde{\chi}_{1}^{0}}|$ from approximately $23$, $14$, and $58~{}{\rm
GeV}$ to approximately $70$, $47$, and $89~{}{\rm GeV}$ for Type-I, -II, and
-III samples, respectively. Additionally, it set a lower bound of
approximately $180~{}{\rm GeV}$ on the Higgsino mass, which is much smaller
than that of the MSSM [19].

  * •

Among three types of samples, the Type-I samples are significantly more
impacted by the LHC constraint due to its inclination toward moderately light
Higgsinos.

In Tables 6, 7, and 8, we present benchmark points for three types of samples
in both $h_{1}$ and $h_{2}$ scenarios. The detailed information of these
points elucidate the underlying physics of Singlino-dominated DM.

##  4 Summary

Traditional supersymmetric models like the MSSM and $\mathbb{Z}_{3}$-NMSSM are
facing increasingly stringent constraints as direct detection experiments for
DM and searches for supersymmetry at the LHC continue to advance rapidly. This
has led to significant fine-tuning of the parameters within these models to
predict viable DM candidates, as evidenced by Refs. [15, 29, 41]. In response
to these challenges, we undertook Bayesian analyses of the GNMSSM to
thoroughly explore its DM physics.

First, we examined the theoretical structure of the GNMSSM and selected a set
of physical parameters as our input criteria. We constructed a comprehensive
likelihood function that incorporated current experimental and theoretical
knowledge on DM physics—including relic density and direct and indirect DM
searches—as well as Higgs and B physics. This likelihood function was then
used to guide detailed scans of the theory’s parameter space using the nested
sampling algorithm. We analyzed the resulting samples using statistical
measures such as the PL and marginal posterior PDF to uncover underlying
physical phenomena. Our findings indicated that the DM candidates in the
GNMSSM may predominantly be either Singlino- or Bino-like DM. Specifically,
the Singlino-like scenario accounted for $99.3\%$ of the Bayesian evidence in
the $h_{1}$ scenario and $64.8\%$ in the $h_{2}$ scenario. Given that the
Bayesian evidence significantly favored the $h_{1}$ scenario over the $h_{2}$
scenario, we concluded that the GNMSSM generally preferred Singlino-dominated
DM across a much broader parameter range. This preference distinctly set the
GNMSSM apart from the MSSM and $\mathbb{Z}_{3}$-NMSSM, highlighting its unique
position within the landscape of supersymmetric models.

Second, we examined the characteristics of Singlino-dominated DM and
identified significant distinctions from the $\mathbb{Z}_{3}$-NMSSM. In the
latter, DM properties were determined by four parameters: $\lambda$,
$\tan\beta$, $\mu_{\rm eff}$, and $m_{\tilde{\chi}_{1}^{0}}\simeq
2\kappa\mu_{\rm eff}/\lambda$. For a Singlino-dominated neutralino to remain
the LSP, $|\kappa|$ had to be smaller than $\lambda/2$. Furthermore, the LZ
constraint on $\lambda$ necessitated $|\kappa|$ to be very small. As a result,
the Singlino DM in the $\mathbb{Z}$-NMSSM needed to co-annihilate with
electroweakinos to match the observed relic abundance. However, such dynamics
differed in the GNMSSM due to the introduction of $\mathbb{Z}_{3}$-violating
terms where the DM’s properties were described through five independent
parameters: $\lambda$, $\kappa$, $\tan\beta$, $\mu_{tot}$, and
$m_{\tilde{\chi}_{1}^{0}}$. Specifically, allowing for decorrelation of
$m_{\tilde{\chi}_{1}^{0}}$ from $\mu_{tot}$ enables $|\kappa|$ to exceed
$\lambda$ when predicting Singlino-dominant DM. This structural difference
significantly impacted the potential creation of a secluded DM sector within
singlet particles [45]. The key characteristics of the secluded sector in our
analysis included the following:

  * •

Annihilation primarily occurred through the channels
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$,
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$, and
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ to achieve the
observed relic abundance. Among the parameters affecting the annihilation
cross sections, $\kappa$ played a crucial role, with
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ preferring the
smallest $|\kappa|$ value compared to other channels. However, due to the
complex dependence of the annihilation cross sections on the model parameters,
the values of $|\kappa|$ ranging from $0.05$ to $0.75$ were all capable of
explaining the observed relic densities.

  * •

The s-wave nature of $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$
resulted in a lower limit on $m_{\tilde{\chi}_{1}^{0}}$ at approximately
$60~{}{\rm GeV}$ from indirect DM search experiments, whereas in other
channels, $m_{\tilde{\chi}_{1}^{0}}$ could decrease to as low as ten GeV.

  * •

The equations governing the interactions between DM and nucleons were
characterized by parameters $\lambda$, $\kappa$, $\tan\beta$, $\mu_{tot}$, and
$m_{\tilde{\chi}_{1}^{0}}$, with $\lambda$ exerting the most significant
influence. It was observed that $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ was
proportional to $\lambda^{2}\kappa^{2}$ in the presence of substantial
singlet–doublet Higgs mixing, and to $\lambda^{4}$ if $h_{s}$ was very
massive. In contrast, $\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ consistently
varied as a function of approximately $\lambda^{4}$. Numerical analyses
indicated that for most samples in the $h_{1}$ scenario, $\lambda\lesssim
0.09$; for the $h_{2}$ scenario, it was suggested that $\lambda\lesssim 0.05$.
Consequently, $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ and $\sigma^{\rm
SD}_{\tilde{\chi}_{1}^{0}-n}$ could be as low as $10^{-49}$ and $10^{-52}\
{\rm cm^{2}}$, respectively. Additionally, among the considered channels,
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$ in the $h_{2}$
scenario faced most stringent constraints from the LZ experiment due to the
significant Higgs mixing and its preference for relatively light Higgsinos.

  * •

For Singlino-like DM in the GNMSSM, although $\sigma^{\rm
SI}_{\tilde{\chi}_{1}^{0}-p}$ could approach the exclusion limits set by the
LZ experiment, $\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ remained well below
these bounds. Conversely, both $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$ and
$\sigma^{\rm SD}_{\tilde{\chi}_{1}^{0}-n}$ could reach exclusion limits for
Bino-like DM. This observation suggested that if SD scattering were observed
in near future, the viability of Singlino-like DM would be significantly
undermined.

  * •

In exceptional cases, Singlino-dominated DM could achieve observed relic
abundance by co-annihilating with either Higgsino- (for most cases) or Wino-
like electroweakinos, with permissible values for $\mu_{\rm tot}$ down to
approximate $200~{}{\rm GeV}$. In contrast, Bino-like DM with a mass at the
electroweak scale could only co-annihilate with Wino-like electroweakinos to
achieve correct abundance, and the LZ experiment set a threshold for $\mu_{\rm
tot}$ at approximately $380~{}{\rm GeV}$.

Finally, we investigated the impact of the LHC’s search for supersymmetry on
DM physics, which was effective in excluding scenarios that predicted
moderately light DM and Higgsinos. Our findings suggested that this constraint
significantly raised the lower limits for $|m_{\tilde{\chi}_{1}^{0}}|$ from
approximately $23$, $14$, and $58~{}{\rm GeV}$ in
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$,
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$, and
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$ to about $70$, $47$,
and $89~{}{\rm GeV}$, respectively. Moreover, a lower limit of around
$180~{}{\rm GeV}$ was imposed on the Higgsino mass, representing a significant
departure from assumptions in the MSSM [19]. Notably, among these channels,
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$ was most influenced
by the LHC constraint due to its preference for moderately light Higgsinos.

This study differed from our previous work cited in Ref. [43] in the following
aspects.

  * •

In contrast to Ref. [43], which focused on the $\mu$-term-extended NMSSM and
investigated the viability of Singlino-dominated DM, Bayesian analyses were
conducted in this study within a general theoretical framework of the NMSSM.
Specifically, we proposed a set of parameters with clear physical meanings as
theoretical inputs and concluded that the theory favors Singlino-dominated DM
across a much wider parameter space compared to Bino-like DM.

  * •

We incorporated the most recent experimental results to constrain the DM
physics of the GNMSSM, including the 2023 LZ search for both SI and SD DM-
nucleon scattering, various LHC searches for supersymmetry, and measurements
of $h$’s couplings at the LHC with $136~{}{\rm fb^{-1}}$ of data. We derived
mass bounds on $\tilde{\chi}_{1}^{0}$ from indirect DM search experiments and
the LHC supersymmetry search.

  * •

We provide analytical formulas for the SI and SD cross sections of DM-nucleon
scattering. These demonstrate that $\sigma^{\rm SI}_{\tilde{\chi}_{1}^{0}-p}$
is proportional to $\lambda^{2}\kappa^{2}$ amid substantial singlet–doublet
Higgs mixing and $\lambda^{4}$ in heavy $h_{s}$ limit, whereas $\sigma^{\rm
SD}_{\tilde{\chi}_{1}^{0}-n}$ consistently scales with $\lambda^{4}$.
Additionally, we present analytic expressions for the annihilation cross
sections of $\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}h_{s}$,
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to A_{s}A_{s}$, and
$\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}\to h_{s}A_{s}$.

Our study delivered valuable insights into the DM physics of the GNMSSM,
reopening the $\mathbb{Z}_{3}$-NMSSM’s extensive parameter space previously
excluded by experiments. We eagerly anticipate further results from DM and LHC
experiments that may provide more clues about the properties of DM and help
pinpoint the correct theoretical framework.

## Acknowledgments

J. Cao thanks Dr. Yuanfang Yue and Jinwei Lian for helpful discussions. This
work was supported by the National Natural Science Foundation of China (NNSFC)
under grant No. 12075076.

## References

  * [1] Planck collaboration, P. A. R. Ade et al., _Planck 2015 results. XIII. Cosmological parameters_ , _Astron. Astrophys._ 594 (2016) A13 [1502.01589].
  * [2] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [3] K. Griest and M. Kamionkowski, _Supersymmetric dark matter_ , _Phys. Rept._ 333 (2000) 167.
  * [4] G. Bertone, D. Hooper and J. Silk, _Particle dark matter: Evidence, candidates and constraints_ , _Phys. Rept._ 405 (2005) 279 [hep-ph/0404175].
  * [5] S. Baum, M. Carena, N. R. Shah and C. E. M. Wagner, _Higgs portals for thermal Dark Matter. EFT perspectives and the NMSSM_ , _JHEP_ 04 (2018) 069 [1712.09873].
  * [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 [1805.12562].
  * [7] PandaX-II collaboration, Q. Wang et al., _Results of dark matter search using the full PandaX-II exposure_ , _Chin. Phys. C_ 44 (2020) 125001 [2007.15469].
  * [8] PandaX-II collaboration, X. Cui et al., _Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment_ , _Phys. Rev. Lett._ 119 (2017) 181302 [1708.06917].
  * [9] LZ collaboration, J. Aalbers et al., _First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment_ , _Phys. Rev. Lett._ 131 (2023) 041002 [2207.03764].
  * [10] PAMELA collaboration, O. Adriani et al., _PAMELA results on the cosmic-ray antiproton flux from 60 MeV to 180 GeV in kinetic energy_ , _Phys. Rev. Lett._ 105 (2010) 121101 [1007.0821].
  * [11] Fermi-LAT collaboration, M. Ackermann et al., _Measurement of separate cosmic-ray electron and positron spectra with the Fermi Large Area Telescope_ , _Phys. Rev. Lett._ 108 (2012) 011103 [1109.0521].
  * [12] AMS 01 collaboration, M. Aguilar et al., _Cosmic-ray positron fraction measurement from 1 to 30-GeV with AMS-01_ , _Phys. Lett. B_ 646 (2007) 145 [astro-ph/0703154].
  * [13] J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait and H.-B. Yu, _Constraints on Light Majorana dark Matter from Colliders_ , _Phys. Lett. B_ 695 (2011) 185 [1005.1286].
  * [14] P. J. Fox, R. Harnik, J. Kopp and Y. Tsai, _Missing Energy Signatures of Dark Matter at the LHC_ , _Phys. Rev. D_ 85 (2012) 056011 [1109.4398].
  * [15] 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 [1910.14317].
  * [16] J. F. Gunion and H. E. Haber, _Higgs Bosons in Supersymmetric Models. 1._ , _Nucl. Phys. B_ 272 (1986) 1.
  * [17] H. E. Haber and G. L. Kane, _The Search for Supersymmetry: Probing Physics Beyond the Standard Model_ , _Phys. Rept._ 117 (1985) 75.
  * [18] A. Djouadi, _The Anatomy of electro-weak symmetry breaking. II. The Higgs bosons in the minimal supersymmetric model_ , _Phys. Rept._ 459 (2008) 1 [hep-ph/0503173].
  * [19] Y. He, L. Meng, Y. Yue and D. Zhang, _Impact of the recent measurement of (g-2) $\mu$, the LHC search for supersymmetry, and the LZ experiment on the minimal supersymmetric standard model_, _Phys. Rev. D_ 108 (2023) 115010 [2303.02360].
  * [20] G. F. Giudice and A. Masiero, _A Natural Solution to the mu Problem in Supergravity Theories_ , _Phys. Lett. B_ 206 (1988) 480.
  * [21] A. Arvanitaki, M. Baryakhtar, X. Huang, K. van Tilburg and G. Villadoro, _The Last Vestiges of Naturalness_ , _JHEP_ 03 (2014) 022 [1309.3568].
  * [22] J. A. Evans, Y. Kats, D. Shih and M. J. Strassler, _Toward Full LHC Coverage of Natural Supersymmetry_ , _JHEP_ 07 (2014) 101 [1310.5758].
  * [23] H. Baer, V. Barger, D. Mickelson and M. Padeffke-Kirkland, _SUSY models under siege: LHC constraints and electroweak fine-tuning_ , _Phys. Rev. D_ 89 (2014) 115019 [1404.2277].
  * [24] U. Ellwanger, C. Hugonie and A. M. Teixeira, _The Next-to-Minimal Supersymmetric Standard Model_ , _Phys. Rept._ 496 (2010) 1 [0910.1785].
  * [25] M. Maniatis, _The Next-to-Minimal Supersymmetric extension of the Standard Model reviewed_ , _Int. J. Mod. Phys. A_ 25 (2010) 3505 [0906.0777].
  * [26] J. Cao, Y. He, L. Shang, W. Su and Y. Zhang, _Natural NMSSM after LHC Run I and the Higgsino dominated dark matter scenario_ , _JHEP_ 08 (2016) 037 [1606.04416].
  * [27] U. Ellwanger, _Present Status and Future Tests of the Higgsino-Singlino Sector in the NMSSM_ , _JHEP_ 02 (2017) 051 [1612.06574].
  * [28] 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) 055031 [1606.02149].
  * [29] 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. D_ 99 (2019) 075020 [1810.09143].
  * [30] 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 (2018) 735 [1806.09478].
  * [31] F. Domingo, J. S. Kim, V. M. Lozano, P. Martin-Ramiro and R. Ruiz de Austri, _Confronting the neutralino and chargino sector of the NMSSM with the multilepton searches at the LHC_ , _Phys. Rev. D_ 101 (2020) 075010 [1812.05186].
  * [32] S. Baum, N. R. Shah and K. Freese, _The NMSSM is within Reach of the LHC: Mass Correlations \ & Decay Signatures_, _JHEP_ 04 (2019) 011 [1901.02332].
  * [33] M. van Beekveld, S. Caron and R. Ruiz de Austri, _The current status of fine-tuning in supersymmetry_ , _JHEP_ 01 (2020) 147 [1906.10706].
  * [34] W. Abdallah, A. Chatterjee and A. Datta, _Revisiting singlino dark matter of the natural $Z_{3}$-symmetric NMSSM in the light of LHC_, _JHEP_ 09 (2019) 095 [1907.06270].
  * [35] M. Guchait and A. Roy, _Light Singlino Dark Matter at the LHC_ , _Phys. Rev. D_ 102 (2020) 075023 [2005.05190].
  * [36] W. Abdallah, A. Datta and S. Roy, _A relatively light, highly bino-like dark matter in the Z 3-symmetric NMSSM and recent LHC searches_, _JHEP_ 04 (2021) 122 [2012.04026].
  * [37] D. Das, U. Ellwanger and A. M. Teixeira, _Modified Signals for Supersymmetry in the NMSSM with a Singlino-like LSP_ , _JHEP_ 04 (2012) 067 [1202.5244].
  * [38] U. Ellwanger and A. M. Teixeira, _NMSSM with a singlino LSP: possible challenges for searches for supersymmetry at the LHC_ , _JHEP_ 10 (2014) 113 [1406.7221].
  * [39] J. Cao, F. Li, J. Lian, Y. Pan and D. Zhang, _Impact of LHC probes of SUSY and recent measurement of (g $-$ 2)μ on $\mathbb{Z}$3-NMSSM_, _Sci. China Phys. Mech. Astron._ 65 (2022) 291012 [2204.04710].
  * [40] J. Cao, L. Meng and Y. Yue, _Electron and muon anomalous magnetic moments in the Z3-NMSSM_ , _Phys. Rev. D_ 108 (2023) 035043 [2306.06854].
  * [41] H. Zhou, J. Cao, J. Lian and D. Zhang, _Singlino-dominated dark matter in Z3-symmetric NMSSM_ , _Phys. Rev. D_ 104 (2021) 015017 [2102.05309].
  * [42] H. Jeffreys, _Theory of probability_ , Oxford Classic Texts in the Physical Sciences. The Clarendon Press, Oxford University Press, New York, 1998.
  * [43] J. Cao, D. Li, J. Lian, Y. Yue and H. Zhou, _Singlino-dominated dark matter in general NMSSM_ , _JHEP_ 06 (2021) 176 [2102.05317].
  * [44] J. Cao, X. Jia, L. Meng, Y. Yue and D. Zhang, _Status of the singlino-dominated dark matter in general Next-to-Minimal Supersymmetric Standard Model_ , _JHEP_ 03 (2023) 198 [2210.08769].
  * [45] M. Pospelov, A. Ritz and M. B. Voloshin, _Secluded WIMP Dark Matter_ , _Phys. Lett. B_ 662 (2008) 53 [0711.4866].
  * [46] 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 [1809.07371].
  * [47] W. G. Hollik, G. Weiglein and J. Wittbrodt, _Impact of Vacuum Stability Constraints on the Phenomenology of Supersymmetric Models_ , _JHEP_ 03 (2019) 109 [1812.04644].
  * [48] U. Ellwanger, _NONRENORMALIZABLE INTERACTIONS FROM SUPERGRAVITY, QUANTUM CORRECTIONS AND EFFECTIVE LOW-ENERGY THEORIES_ , _Phys. Lett. B_ 133 (1983) 187.
  * [49] S. A. Abel, _Destabilizing divergences in the NMSSM_ , _Nucl. Phys. B_ 480 (1996) 55 [hep-ph/9609323].
  * [50] C. F. Kolda, S. Pokorski and N. Polonsky, _Stabilized singlets in supergravity as a source of the mu - parameter_ , _Phys. Rev. Lett._ 80 (1998) 5263 [hep-ph/9803310].
  * [51] C. Panagiotakopoulos and K. Tamvakis, _Stabilized NMSSM without domain walls_ , _Phys. Lett. B_ 446 (1999) 224 [hep-ph/9809475].
  * [52] G. G. Ross and K. Schmidt-Hoberg, _The Fine-Tuning of the Generalised NMSSM_ , _Nucl. Phys. B_ 862 (2012) 710 [1108.1284].
  * [53] H. M. Lee, S. Raby, M. Ratz, G. G. Ross, R. Schieren, K. Schmidt-Hoberg et al., _A unique $\mathbb{Z}_{4}^{R}$ symmetry for the MSSM_, _Phys. Lett. B_ 694 (2011) 491 [1009.0905].
  * [54] H. M. Lee, S. Raby, M. Ratz, G. G. Ross, R. Schieren, K. Schmidt-Hoberg et al., _Discrete R symmetries for the MSSM and its singlet extensions_ , _Nucl. Phys. B_ 850 (2011) 1 [1102.3595].
  * [55] G. G. Ross, K. Schmidt-Hoberg and F. Staub, _The Generalised NMSSM at One Loop: Fine Tuning and Phenomenology_ , _JHEP_ 08 (2012) 074 [1205.1509].
  * [56] 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_ 03 (2012) 086 [1202.5821].
  * [57] D. J. Miller, R. Nevzorov and P. M. Zerwas, _The Higgs sector of the next-to-minimal supersymmetric standard model_ , _Nucl. Phys. B_ 681 (2004) 3 [hep-ph/0304049].
  * [58] ATLAS collaboration, G. Aad et al., _Search for heavy Higgs bosons decaying into two tau leptons with the ATLAS detector using $pp$ collisions at $\sqrt{s}=13$ TeV_, _Phys. Rev. Lett._ 125 (2020) 051801 [2002.12223].
  * [59] ATLAS collaboration, G. Aad et al., _Search for charged Higgs bosons decaying into a top quark and a bottom quark at $\sqrt{\mathrm{s}}$ = 13 TeV with the ATLAS detector_, _JHEP_ 06 (2021) 145 [2102.10076].
  * [60] M. Badziak, M. Olechowski and P. Szczerbiak, _Spin-dependent constraints on blind spots for thermal singlino-higgsino dark matter with(out) light singlets_ , _JHEP_ 07 (2017) 050 [1705.00227].
  * [61] M. Badziak, M. Olechowski and P. Szczerbiak, _Blind spots for neutralinos in NMSSM with light singlet scalar_ , _PoS_ PLANCK2015 (2015) 130 [1601.00768].
  * [62] C. Cheung, M. Papucci, D. Sanford, N. R. Shah and K. M. Zurek, _NMSSM Interpretation of the Galactic Center Excess_ , _Phys. Rev. D_ 90 (2014) 075011 [1406.6372].
  * [63] K. Griest and D. Seckel, _Three exceptions in the calculation of relic abundances_ , _Phys. Rev. D_ 43 (1991) 3191.
  * [64] T. Nihei, L. Roszkowski and R. Ruiz de Austri, _Exact cross-sections for the neutralino WIMP pair annihilation_ , _JHEP_ 03 (2002) 031 [hep-ph/0202009].
  * [65] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre et al., _The waning of the WIMP? A review of models, searches, and constraints_ , _Eur. Phys. J. C_ 78 (2018) 203 [1703.07364].
  * [66] M. J. Baker et al., _The Coannihilation Codex_ , _JHEP_ 12 (2015) 120 [1510.03434].
  * [67] A. Pierce, N. R. Shah and K. Freese, _Neutralino Dark Matter with Light Staus_ , 1309.7351.
  * [68] M. Drees and M. Nojiri, _Neutralino - nucleon scattering revisited_ , _Phys. Rev. D_ 48 (1993) 3483 [hep-ph/9307208].
  * [69] M. Drees and M. M. Nojiri, _New contributions to coherent neutralino - nucleus scattering_ , _Phys. Rev. D_ 47 (1993) 4226 [hep-ph/9210272].
  * [70] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _Dark matter direct detection rate in a generic model with micrOMEGAs 2.2_ , _Comput. Phys. Commun._ 180 (2009) 747 [0803.2360].
  * [71] J. M. Alarcon, J. Martin Camalich and J. A. Oller, _The chiral representation of the $\pi N$ scattering amplitude and the pion-nucleon sigma term_, _Phys. Rev. D_ 85 (2012) 051503 [1110.3797].
  * [72] J. M. Alarcon, L. S. Geng, J. Martin Camalich and J. A. Oller, _The strangeness content of the nucleon from effective field theory and phenomenology_ , _Phys. Lett. B_ 730 (2014) 342 [1209.2870].
  * [73] ATLAS collaboration, G. Aad et al., _Interpretations of the ATLAS measurements of Higgs boson production and decay rates and differential cross-sections in $pp$ collisions at $\sqrt{s}=13$ TeV_, 2402.05742.
  * [74] P. Huang and C. E. M. Wagner, _Blind Spots for neutralino Dark Matter in the MSSM with an intermediate $m_{A}$_, _Phys. Rev. D_ 90 (2014) 015018 [1404.0392].
  * [75] M. Badziak, M. Olechowski and P. Szczerbiak, _Blind spots for neutralino dark matter in the NMSSM_ , _JHEP_ 03 (2016) 179 [1512.02472].
  * [76] F. Staub, _SARAH_ , 0806.0538.
  * [77] F. Staub, _SARAH 3.2: Dirac Gauginos, UFO output, and more_ , _Comput. Phys. Commun._ 184 (2013) 1792 [1207.0906].
  * [78] F. Staub, _SARAH 4 : A tool for (not only SUSY) model builders_ , _Comput. Phys. Commun._ 185 (2014) 1773 [1309.7223].
  * [79] F. Staub, _Exploring new models in all detail with SARAH_ , _Adv. High Energy Phys._ 2015 (2015) 840780 [1503.04200].
  * [80] 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 [hep-ph/0301101].
  * [81] W. Porod and F. Staub, _SPheno 3.1: Extensions including flavour, CP-phases and models beyond the MSSM_ , _Comput. Phys. Commun._ 183 (2012) 2458 [1104.1573].