Impact of recent (g − 2)μ measurement on the light CP-even Higgs scenario in general Next-to-Minimal Supersymmetric Standard Model

The General Next-to-Minimal Supersymmetric Standard Model (GNMSSM) is an attractive theory that is free from the tadpole problem and the domain-wall problem of Z3-NMSSM, and can form an economic secluded dark matter (DM) sector to naturally predict the DM experimental results. It also provides mechanisms to easily and significantly weaken the constraints from the LHC search for supersymmetric particles. These characteristics enable the theory to explain the recently measured muon anomalous magnetic moment, (g − 2)μ, in a broad parameter space that is consistent with all experimental results and at same time keeps the electroweak symmetry breaking natural. This work focuses on a popular scenario of the GNMSSM in which the next-to-lightest CP-even Higgs boson corresponds to the scalar discovered at the Large Hadron Collider (LHC). Both analytic formulae and a sophisticated numerical study show that in order to predict the scenario without significant tunings of relevant parameters, the Higgsino mass μtot ≲ 500 GeV and tan β ≲ 30 are preferred. This character, if combined with the requirement to account for the (g − 2)μ anomaly, will entail some light sparticles and make the LHC constraints very tight. As a result, this scenario can explain the muon anomalous magnetic moment in very narrow corners of its parameter space.


Introduction
The latest measurement of the muon anomalous magnetic moment a µ ≡ (g − 2) µ /2 announced by the Fermilab National Accelerator Laboratory (FNAL) [1] is in full agreement with the Brookhaven National Laboratory (BNL) E821 result [2]. The combined experimental average is given by equation ( In addition, the Run-1 results in Fermilab also imply that a more thorough analysis in future experiments will most probably substantiate the excess of a µ in 5σ discovery level. In recent years, this situation has inspired continuous attention to a µ . In particular, it was widely conjectured that the anomaly may arise from new physics beyond the SM (see, e.g., ref. [24] and the references therein). Among a variety of theories that can account for the anomalous magnetic moment, supersymmetry (SUSY) is especially promising due JHEP03(2022)203 to its elegant structure and natural solutions to many puzzles in the SM, such as the hierarchy problem, the unification of different forces, and the dark matter (DM) mystery [25][26][27][28]. Studies of low-energy supersymmetric models have indicated that the source of the significant deviation can be totally or partially contributed to smuon-neutralino or sneutrino-chargino loop effects [24,. Although SUSY has multiple theoretical advantages, it has been strongly restricted by DM direct detection (DD) experiments, such as XENON-1T [80,81] and PandaX-4T [82,83] experiments, as well as LHC sparticle searches [84][85][86][87][88][89][90][91]. As a consequence, some of its economical realizations become unnatural for electro-weak symmetry breaking in interpreting the anomalous magnetic moment. In the Minimal Supersymmetric Standard Model (MSSM) with R-parity conservation [26,[92][93][94], the lightest neutralino is usually the lightest supersymmetric particle (LSP), and thus a viable DM candidate. In order to fully account for the DM abundance measured by the Planck experiment [95], it must be Bino-dominated when it is lighter than 1 TeV [96]. In this case, the XENON-1T experiment and the LHC experiments prefer that the magnitude of the Higgsino mass parameter, µ, be larger than about 500 GeV in explaining the discrepancy in the 2 σ level [78]. This implies a fine-tuning at the order of 1 % in predicting m Z [97]. This conclusion may be understood from the features of the DM annihilation mechanisms in the MSSM: • In the case that the DM co-annihilates with Wino-dominated particles to obtain the measured abundance, the Higgsino mass prefers to be much larger than the DM mass, which was explained in appendix A of this work. In addition, the SUSY explanation of the deviation together with the results of the LHC search for SUSY can further restrict µ values in a significant way. 1 • In the case that the DM co-annihilates with Higgsino-dominated particles to obtain the measured abundance, the DM DD experiments have required |µ| to be as large as several TeV, because DM-nucleon scattering rates are enhanced by a factor of 1/(1 − m 2 χ 0 1 /µ 2 ) 2 as m 2 χ 0 1 /µ 2 → 1 (see discussions in appendix A).
• In the case that the DM co-annihilates with Sleptons to obtain the measured abundance, the LHC's searches for electroweakinos require Higgsinos to be massive because Wino-and Higgsino-dominated electroweakinos can decay into Sleptons and thus enhance the production rate of lepton signals at the LHC [67].
• In the case that the DM co-annihilates with Squarks or Gluinos to obtain the measured abundance, the LHC's searches for colored sparticles require the DM mass to be heavier than 1 TeV [98]. 1 In carrying out this study, we also investigated the characteristics of MSSM and Z3-NMSSM in a similar way to this work. We found that the co-annihilations with Wino-and Slepton-dominated particles are the main annihilation mechanisms of the bino-dominated DM , i.e., their corresponding Bayesian evidences are the largest in comparison with the other annihilation mechanisms. We also found that the DM preferred to be heavier than about 300 GeV and |µ| 500 GeV when we implemented detailed Monte Carlo simulations for the constraints from the latest LHC searches for electroweakinos. These observations significantly improve the conclusions of [78] and [79], since more LHC analyses were considered carefully.

JHEP03(2022)203
• In the case that the DM obtains the measured abundance by the SM-like Higgs funnel or Z funnel, the LHC's searches for eletroweakinos prefer massive Higgsinos because the DM is relatively light and the LHC's constraints on sparticle mass spectrum are rather strong [99].
• The case in which the DM obtains the measured abundance by the resonance of heavy doublet Higgs bosons is rare. One reason for this is that this case requires significant tuning of SUSY parameters to realize the correlation |mχ0 1 | m A /2, where m A denotes the mass of CP-odd Higgs bosons in MSSM. Another reason is that the LHC's searches for exotic Higgs bosons prefer the bosons to be very massive [100]. Consequently, the DM is also massive.
Next, we consider the Next-to-Minimal Supersymmetric Standard Model with a Z 3 symmetry (Z 3 -NMSSM) [101,102]. This model extends MSSM by a gauge-singlet Higgs superfieldŜ, and has the advantage that either a Bino-dominated (in most physical cases) or a Singlino-dominated neutralino can act as a viable DM candidate [99,[103][104][105][106][107][108][109][110][111][112][113][114]. The Bino-dominated DM candidate differs from the MSSM prediction mainly in that it could coannihilate with a Singlino-dominated neutralino to obtain the measured abundance [106]. This situation, however, occurrs in very narrow parameter space characterized by |2κµ/λ| |M 1 |, moderately large λ and κ, and |µ| 300 GeV [106,112]. In addition, a SUSY µ is insensitive to Yukawa coupling λ because the Singlino field has no mixing with Wino and Bino fields, and it does not couple directly to the muon lepton. As a result, the formulae to calculate a SUSY µ in NMSSM are same as those at the lowest order of the mass-insertion approximation in MSSM [67]. Considering these features, we expected that the Binodominated DM case in the Z 3 -NMSSM and MSSM would not show significant differences in explaining the discrepancy. The properties of the Singlino-dominated DM are determined by λ, the Higgsino mass µ (denoted by µ tot in this work), and tanβ for a given DM mass [115]. A relatively large λ can increase the DM-nucleon scattering cross-sections, and so far, λ 0.3 is disfavored by the XENON-1T experiments [112,115]. This conclusion implies that the traditional DM annihilation channels,χ 0 1χ 0 1 → tt, h s A s , hA s , where t, h, h s , and A s denote the top quark, SM-like Higgs boson, and singlet-dominated CPeven and CP-odd Higgs bosons, respectively, can not be fully responsible for the measured abundance [112]. As a result, the DM is more likely to obtain the abundance by means of the co-annihilation with the Higgsino-dominated particles, which corresponds to a correlated parameter space of 2|κ| λ with λ 0.1. The Bayesian evidence in this case is heavily suppressed owing to the very narrow parameter space, which entails a certain degree of finetuning to meet the DM experiments [115]. In addition, the interpretation of the magnetic moment causes the Z 3 -NMSSM to be further restricted by the updated searches for SUSY at the LHC with 139 f b −1 data. In particular, the region of tan β 30 in figure 7 of [99] has been excluded because both the DM and Higgsinos are relatively light. Such a situation, as we will show below, was frequently encountered in this work.
The dilemma of MSSM and Z 3 -NMSSM inspired us to study the general Next-to-Minimal Supersymmetric Standard Model (GNMSSM) [116]. Unlike Z 3 -NMSSM, GN-MSSM usually predicts the Singlino-dominated neutralino as a viable DM candidate due JHEP03(2022)203 to its following specific theoretical feature: the properties of the Singlino-dominated DM are described by λ, µ tot , mχ0 1 , tanβ, and κ, among which the first four parameters determine the DM couplings to nucleon, and κ mainly dominates the DM couplings to singletdominated Higgs bosons [116]. Consequently, singlet-dominated particlesχ 0 1 , h s , and A s can constitute a secluded DM sector, where the measured DM abundance can be achieved by the h s /A s -mediated resonant annihilation into SM particles or through the annihilation process ofχ 0 1χ 0 1 → h s A s by adjusting the value of κ. Given that this sector interacts with SM matters only through weak singlet-doublet Higgs field mixing, the DM-nucleon scattering rate can be naturally suppressed by λv/µ tot when λ is small [116]. Since the parameters need no significant tuning to be consistent with the constraints from the DM experiments, the corresponding Bayesian evidence is significantly larger than that for the Bino-dominated DM case [116]. Other characteristics of the theory include that, due to the very weak couplings of the Singlino-dominated DM to other sparticles, heavy sparticles initially prefer to decay into next-to-LSP (NLSP) or next-next-to-LSP (NNLSP). As a result, their decay chains are lengthened and their signals become complicated. In addition, the DM as LSP may be moderately heavy, since the annihilationχ 0 1χ 0 1 → h s A s requires m DM > (m hs + m As )/2. These features weaken significantly the limitations from the LHC's searches for SUSY. Specifically, in our recent work, we studied a SUSY µ in a simplified version of GNMSSM, which we called µ-extended NMSSM (µNMSSM) [67]. We found that, by presuming the DM and LHC experiments are satisfied, µNMSSM can explain the discrepancy in a broad parameter space where Higgsinos are lighter than about 500 GeV.
In our previous work [67], we considered only the h 1 scenario, in which the lightest CP-even Higgs boson corresponds to the SM-like Higgs boson discovered at the LHC. A typical feature of NMSSM is that the next lightest CP-even Higgs boson may also act as the SM-like Higgs boson, which has been dubbed the h 2 scenario in the literature. Thus, a full understanding of GNMSSM necessitates the study of the discrepancy in the h 2 scenario. In particular, given that specific configurations of Higgs parameters are needed to predict m h 1 < (m h 2 125 GeV) and h 2 to be SM-like, it is conceivable that the h 2 scenario suffers tighter experimental constraints than the h 1 scenario. 2 This leaves in doubt the idea that the h 2 scenario can explain the discrepancy. As a result, a careful examination of the experimental constraints on the h 2 scenario is needed, which is the focus of this work.
This work is organized as follows. In section 2, we briefly introduce the basics of GNMSSM and the SUSY contribution to the moment. In section 3, we perform a sophisticated scan over the broad parameter space of µNMSSM, and show the features of the theory in explaining the discrepancy. By using specific Monte Carlo simulations, we also comprehensively study the constraints from the LHC's searches for SUSY. In section 4, we concentrate on the GNMSSM, which has much broader parameter space than µNMSSM, and perform a similar study to those in section 3. Lastly, we draw conclusions in section 5.
2 This conclusion may also be understood intuitively as follows: the lightness of h1 and the premise that h1 and h2 are singlet-dominated and SM-like, respectively, lead to the tendency that some parameters in the Higgs sector are relatively small. Consequently, light sparticles are usually predicted. This phenomenon is similar to the well-known fact that the natural result for electroweak-symmetry breaking prefers |µ| 500 GeV [97].

JHEP03(2022)203 2 Theoretical preliminaries
It is well-known that the superpotential of the popular Z 3 -NMSSM is given by [101,102] where the Yukawa terms W Yukawa are the same as those in MSSM, Higgs superfields, and λ, κ are dimensionless couplings coefficient parameterizing the Z 3 -invariant trilinear terms. GNMSSM differs from Z 3 -MSSM in that its superpotential does not respect the Z 3 symmetry, and thus it contains the following most general renormalizable couplings: Historically, the terms characterized by the bilinear mass parameters µ, µ and the singlet tadpole parameter ξ were introduced to solve the tadpole problem [101,117] and the cosmological domain-wall problem of Z 3 -NMSSM [118][119][120], and the ξ-term can be eliminated by shifting theŜ field and redefining the µ parameter [121]. 3 The bilinear terms could stem from an underlying discrete R symmetry, Z R 4 or Z R 8 , after supersymmetry breaking, and be naturally at the electroweak scale [118,[121][122][123][124]. Note that these extra terms can change the properties of Higgs bosons and neutralinos (in comparison with the Z 3 -NMSSM prediction) and significantly alter the phenomenology of the theory. As emphasized in the introduction, this is one of main motivations of this work.

Higgs sector of GNMSSM
Corresponding to the potential in eq. (2.2), the soft-breaking terms of the GNMSSM are given by [101,102] where H u , H d and S denote the scalar components of the Higgs superfields. The softbreaking mass parameters m 2 Hu , m 2 H d and m 2 S can be fixed by solving the conditional equations for minimizing the scalar potential and then expressing them in terms of the vacuum expectation values (vevs) of the scalar fields: GeV. As usual, the ratio of the two Higgs doublet vevs is defined as tan β ≡ v u /v d , and an effective µ-parameter of MSSM is generated by Consequently, the Higgs sector is described by ten free parameters: tan β, µ eff , the Yukawa couplings λ and κ, the soft-breaking trilinear coefficients A λ and A κ , the bilinear mass parameters µ and µ , and their soft-breaking parameters m 2 3 and m 2 S . The GNMSSM predicts three CP-even Higgs bosons h i = {h, H, h s }, two CP-odd Higgs bosons a i = {A H , A s }, and a pair of charged Higgs bosons H ± = cos βH ± u + sin βH ± d . In We also obtain the following approximations as shown in equation (2.10): These formulae reveal the following facts: • Parameters A λ and m 3 mainly determine the heavy Higgs boson masses, and they have little impact on the other Higgs bosons' mass spectrums.
• m hs and m As depend on parameters λ, κ, µ eff , µ tot , A κ and µ . In addition, m As also depends on m S . This implies that, even when λ, κ, µ eff , µ tot , and µ are fixed, m hs and m As can still vary freely by the adjustment of A κ and m S , respectively. This situation is different from that of Z 3 -NMSSM, where µ tot ≡ µ eff , µ = 0, and m S = 0, and consequently, the masses of singlet fields are correlated [99].
• The most important feature is that the latest LHC Higgs data have imposed an upper limit of about 40 GeV on |λµ tot | in the tremendously large tan β limit, since |λµ tot | may induce a sizable V S h . Furthermore, since a small λ is preferred by DM DD experiments in the Singlino-dominated DM case, |µ eff |, |µ tot | and |µ tot /µ eff | in eq. (2.8) are unlikely to be exceedingly large. Otherwise, strong cancellations among the different terms on the right side of eq. (2.8) are needed to predict m hs < 125 GeV, which makes the theory fine-tuned.
Given that too many parameters are involved in the Higgs sector, the h 2 scenario is studied using the following strategy. First, we assume the charged Higgs bosons to be JHEP03(2022)203 very massive by setting A λ = 2 TeV and m 3 = 1 TeV, following the discussion above. Second, we investigate the characteristics of µNMSSM, where µ and m S are taken to be zero. 4 This model contains most of the key features of the GNMSSM [116], and thus has pedagogical significance. Finally, we concentrate on the GNMSSM by treating µ, µ and m S as variables, and investigate its features in explaining the discrepancy.

Neutralino sector of GNMSSM
The neutralino sector in the GNMSSM consists of the mixtures among the Bino fieldB, the Wino fieldW , the Higgsino fieldsH 0 d ,H 0 u and the Singlino fieldS. Its mass matrix in the basis (−iB, −iW ,H 0 d ,H 0 u ,S) takes the following form [101], as shown in equation (2.11): where M 1 and M 2 are gaugino soft-breaking masses, and µ tot represents the Higgsino mass. This matrix can be diagonalized by a rotation matrix N , and subsequently the mass eigenstates are expressed by equation (2.12): whereχ 0 i (i = 1, 2, 3, 4, 5) are labeled in an ascending mass order. N i3 and N i4 characterize theH 0 d andH 0 u components inχ 0 i , and N i5 denotes the Singlino component. In the case of very massive gauginos and |m 2 χ 0 1 − µ 2 tot | λ 2 v 2 , the following approximations are obtained for the Singlino-dominatedχ 0 1 [131][132][133], given by equations (2.13): 4 Note that µNMSSM as the most economical realization of GNMSSM could arise from the Z3-NMSSM when it was embedded into canonical superconformal supergravity in the Jordan frame, and had applications to the inflation in the early universe [126][127][128][129][130]. This is an interesting realization of supersymmetry in particle physics.

JHEP03(2022)203
These approximations indicate that the mass of the Singlino-dominated DM is determined by the parameters λ, κ, µ eff , µ tot , and µ . In particular, λ and κ are two independent parameters in predicting |mχ0 1 | < |µ tot |. This situation is different from that of the Z 3 -NMSSM, where µ ≡ 0, µ tot ≡ µ eff , and consequently, |κ| must be less than λ/2 to predict the Singlino-dominated neutralino as the LSP [101]. They also indicate that, for fixed tan β, the Higgsino compositions inχ 0 1 depend only on λ, µ tot , and mχ0 1 . Therefore, it is convenient to take the three parameters and κ as theoretical inputs in studying theχ 0 1 's properties, where κ determines the interactions among the singlet-dominated particles. This characteristic contrasts with that of the Z 3 -NMSSM, which only needs the three input parameters of λ, µ tot , and any of mχ0 1 or κ to describeχ 0 1 properties [115]. These differences imply that the singlet-dominated particles may form a secluded DM sector [134], which has the following salient features: • The Singlino-dominated DM can achieve the correct abundance by the process χ 0 1χ 0 1 → h s A s , through adjusting the value of κ, or by the h s /A s -mediated resonant annihilation into SM particles.
• Since the secluded sector communicates with the SM sector only through the weak singlet-doublet Higgs mixing, the interaction between the DM and nucleus is naturally feeble when λ is small.
We added that, even when λ, κ, µ eff and µ tot are fixed, m 0 χ 1 can still vary freely through the tuning of µ . In addition to the processχ 0 1χ 0 1 → h s A s and h s /A s resonant annihilation, the DM has other annihilation channels for obtaining the measured abundance [67], e.g., coannihilation with Higgsino-dominated electroweakinos and/or sleptons, and resonant Z/h annihilations. Owing to these features, the GNMSSM has a broad parameter space consistent with the current DM experimental results. As a result, it is the Singlino-dominated LSP, instead of the Bino-dominated LSP, that is most favored to be a viable DM candidate.

Muon g-2 in the GNMSSM
The SUSY source of the muon anomalous magnetic moment a SUSY µ mainly includes loops with a smuon and a neutralino, as well as those with a muon-type sneutrino and a chargino [29][30][31][32]. The one-loop contributions to a SUSY µ in GNMSSM are given by [30,67] as equations (2.14): (2.14)

JHEP03(2022)203
where i = 1, · · · , 5, j = 1, 2 and l = 1, 2 denote the neutralino, chargino and smuon index, respectively.This gives us equations (2.15): 15) where N is the neutralino mass rotation matrix, X the smuon mass rotation matrix, and U c and V c the chargino mass rotation matrices defined by U c * M C V c † = m diag χ ± . F (x)s are the loop functions of the kinematic variables defined as x il ≡ m 2 /m 2 νµ , and take the following form given by equations (2.16)- (2.19): for the mass-degenerate sparticle case. In practice, it is instructive to understand the features of a SUSY µ through the mass insertion approximation [31]. Specifically, for the lowest order of the approximation, the contributions to a SUSY µ can be classified into four types: "WHL", "BHL", "BHR", and "BLR", where W , B, H, L, and R stands for Wino, Bino, Higgsino, and left-handed and right-handed Smuon fields, respectively. These are from the Feynman diagrams involving respectively, and take the following form [31,33,34] given by equations (2.20)-(2.23): where the loop functions are given by

JHEP03(2022)203
and they satisfy f C (1, 1) = 1/2 and f N (1, 1) = 1/6. Note that the Singlino fieldS can also enter the insertions. Because both theW −S andB 0 −S transitions and theμSμ L,R couplings vanish, the Singlino field only appears in the "WHL", "BHL" and "BHR" loops by two more insertions at the lowest order, which correspond to theH 0 d −S andS −H 0 d transitions in the neutralino mass matrix in eq. (2.11). Since a small λ is preferred by the DM physics, the Singlino-induced contributions are never significant [67]. Although in this case the GNMSSM prediction of a SUSY µ is roughly the same as that of the MSSM, except that the µ parameter of the MSSM should be replaced by µ tot , the two models predict different DM physics and different sparticle signals at the LHC. Thus, they are subject to different theoretical and experimental constraints. It should also be noted that, although there is a prefactor of the Higgsino mass µ in the expression of the "WHL", "BHL", and "BHR" contributions, the involved loop functions approach zero with the increase of |µ|, and consequently these contributions depend on µ in a complex way. By focusing on several typical patterns of sparticle mass spectrum with a positive µ, we found that the "WHL" contribution decreases monotonously as µ increases, while the magnitude of the "BHL" and "BHR" contributions increases when µ is significantly smaller than the slepton mass and decreases when µ is larger than the slepton mass. In addition, the "WHL" contribution is usually much larger than the other contributions ifμ L is not significantly heavier thanμ R .

Research strategy
This sector focuses on the h 2 scenario of µNMSSM. In order to analyze its characteristics in explaining the discrepancy, the following parameter space was scanned by MultiNest algorithm [135]: where the flat prior distribution was chosen for all input parameters and the active point number, n live , was set to be 6000. 5 Other dimensional parameters that are unimportant to this study were fixed at 2 TeV, and include SUSY parameters for the first and third generation sleptons, three generation squarks, and gluinos. In numerical calculations, the model file of GNMSSM was constructed using the package SARAH-4.14.3 [136][137][138][139].

JHEP03(2022)203
contains the value of a SUSY µ [67], and it takes the following form, given by equation (3.2): if restrictions unsatisfied, where the restrictions on each sample include: 1. DM relic abundance, 0.096 ≤ Ωh 2 ≤ 0.144. In implementing this constraint, the central value of the Planck-2018 data, Ωh 2 = 0.120 [149], was used, and a theoretical uncertainty of 10% in abundance calculation was assumed.
2. DM direct and indirect detections. Specifically, the SI and SD DM-nucleon scattering cross-sections should be lower than the bounds from the XENON-1T experiments [150,151], and the DM annihilation rate at present time should be consistent with dwarf galaxies observations from Fermi-LAT collaboration [152]. The method suggested in [153] was adopted in studying the latter constraint.
3. Higgs data fit. The properties of the next lightest CP-even Higgs boson h 2 (also denoted by h throughout this work) should be consistent at the 95% confidence level with corresponding data obtained by ATLAS and CMS collaborations. This condition was checked with the code HiggsSignal-2.2.3 [154] by requiring the sample's p value to be larger than 0.05.
4. Direct searches for extra Higgs bosons at LEP, Tevatron and LHC. This requirement was examined by the code HiggsBounds-5.3.2 [155].
5. Some B-physics observations. Specifically, the branching ratios of B s → µ + µ − and B → X s γ should agree with their experimental measurements, which were summarized in [156] at the 2σ level.
6. LHC searches for SUSY. In order to explain the discrepancy, the electoweakinos and sleptons in the GNMSSM can not be excessively heavy. Thus, they will be produced at the LHC to generate multi-lepton signals. The code SModelS-2.1.1 [157] was used to set limits on the signals in some simple topology cases. Sophisticated study of the constraints will be carried out in subsection 3.3, using the package 7. Vacuum stability for the scalar potential consisting of the Higgs fields and the last two generation slepton fields. This condition was checked by the Vevacious code [161,162], and its effect on the GNMSSM was recently discussed in [67].
In presenting the results, two-dimensional profile likelihood (PL) for the function L in eq.  where Θ i (i = 1, 2, . . .) denote the input parameters, Θ A,B are the variables of interest, and the maximization of L(Θ A , Θ B ) is achieved by scanning the parameters other than Θ A and Θ B . Related quantities includes 1σ and 2σ confidence intervals (CI), and the χ 2 function defined by χ 2 ≡ −2lnL(Θ A , Θ B ). These statistical measures were briefly introduced in [163], and they reflect the capability of the theory to explain the discrepancy.

Key features of the interpretation
All samples obtained in the scan were projected onto different parameter planes to show two-dimensional PLs, which could reveal the underlying physics of the h 2 scenario. Figure 1 illustrates that the scenario can interpret the discrepancy in a broad parameter space.  39.8 × 10 −10 , respectively. The upper right panel shows that the maximum reach of µ tot decreases monotonously with the increase of tan β, and it is about 500 GeV (260 GeV) for tan β = 10 (tan β = 60). The reason for such a behavior is that, in the case of a relatively small tan β, the second term in M 2 S,23 of eq. (2.4) is sizable, and can cancel the first term to suppress V S h , which is preferred by LHC Higgs data. As tan β increases, the cancellation effect becomes weak since the second term is suppressed by sin 2β, and tighter constraints are set on µ tot . 6 Moreover, analyzing the posterior probability of the scan results indicates that the scenario prefers small tan β region. Thus, most samples obtained in the scan predict tan β 30. 6 Throughout this work, A λ is fixed at 2 TeV. If a larger A λ , e.g., A λ = 10 TeV, was taken, it was found that tan β tended to become larger, while |µtot| and λ tended to be smaller [164]. This tendency is needed to suppress M 2 S, 23 and M 2 S,33 in eq. (2.4) simultaneously. In addition, the Bayesian evidence of the scenario decreases significantly as A λ increases [164], which means that setting a large A λ will cause a more subtle parameter tuning to obtain m h 1 125 GeV and correct electroweak symmetry breaking. In brief, even when A λ is treated as a variable in studying the parameter space, the natural realization of the h2 scenario to interpret the anomaly in the GNMSSM, as suggested by this work, has been tightly limited. This conclusion was verified by our alternative scans.

JHEP03(2022)203
The lower left and right panels of figure 1 depict the ranges of M 2 ,μ L , andμ R , which are determined by a SUSY µ in eqs. (2.20)-(2.23). They show that M 2 may be as large as 1.5 TeV, andμ L andμ R may be as large as 1 TeV. The lower left panel also exhibits that the mass of charginoχ ± 1 is less than about 350 GeV. It should be noted that the ranges of M 2 andμ L depend strongly on the value of tan β. For example, assuming that the theory explains the discrepancy of ∆a µ at 1σ level, it was found that M 2 andμ L must be less than about 400 GeV and 350 GeV, respectively, for tan β = 10. The upper bounds become 1.2 TeV and 700 GeV for tan β = 20, and 1.4 TeV and 1 TeV for tan β = 27. By contrast, mχ0 1 andμ R are not sensitive to tan β, e.g., mχ0 1 andμ R may vary in the range of 50 GeV mχ0 1 250 GeV and 100 GeV μ R 1 TeV for any value of tan β. The basic reason for the phenomenon is that the WHL contribution to a SUSY µ is usually the dominant one. It depends on M 2 , µ tot , andμ L , and is in particular proportional to tan β. Therefore, when tan β is relatively small, the invovled SUSY particles must be moderately light to predict a sizable a SUSY µ . As a result, the left sides of the lower panels usually correspond to a relatively small tan β, and the right sides correspond to a large tan β. Figure 2 focuses on the DM physics of the h 2 scenario, which involves the parameters λ, κ, µ tot , and the masses of singlet-dominated particles, i.e., mχ0 1 , m hs and m As . It reveals the following features: • 2mχ0 1 > m hs + m As for most of the parameter areas (see the upper right panel), which implies that in the early universe, the Singlino-dominated DM might annihilate into the singlet-dominated Higgs bosons h s and A s . As pointed out in [116], this annihilation proceeded by the s-channel exchange of Z boson and CP-odd Higgs bosons and the t-channel exchange of neutralinos. If the t-channel contribution to the annihilation rate was much larger than the s-channel contribution, |κ| 0.15 × (mχ0 1 /300 GeV) 1/2 could predict the measured abundance, while if the interference of the two contributions was significantly constructive/deconstructive, a smaller/larger |κ| could be fully responsible for the abundance. Given that 0.05 |κ| 0.25 on the upper left panel, we infer that the process played an important role in determining the abundance. In addition, it was verified in fewer cases that 2mχ0 1 < m hs + m As and/or |κ| 0.1, so that the DM obtained the measured abundance mainly by coannihilating with the Higgsino-dominated electroweakinos or µ-type sleptons.
• The SI and SD cross-sections of DM-nucleon scattering may be as low as 10 −50 cm and 10 −45 cm, respectively (see the lower left and right panels). The SD scattering proceeds only through the Z-mediated Feynman diagram, and the rate is proportional to (λv/µ tot ) 4 [116]. Thus, it is a small λ, e.g., λ ∼ O(0.01), that is responsible for the low SD cross-section (see the upper left panel). By contrast, the SI scattering is induced by three CP-even Higgs bosons, and it is the cancellation of h-and h smediated contributions that mainly accounts for the small SI cross section [115].
The DM physics in the h 2 scenario differs from those of the h 1 scenario, which were presented in figure 2 of [67], in three aspects. The first is that the DM is relatively light, i.e., 50 GeV mχ0  scenario. Two reasons may explain this phenomenon. One is that |mχ0 1 | must be less than µ tot , and as shown in figure 1, a moderately small µ tot is experimentally preferred for the h 2 scenario. The other reason is that |mχ0 1 | must be larger than (m hs + m As )/2 for most cases to proceed to the annihilationχ 0 1χ 0 1 → h s A s . A relatively lightχ 0 1 can meet this condition in the h 2 scenario (see the previous discussions). The second one is that |κ| is less than 0.25 in the h 2 scenario, while it is less than 0.4 in the h 1 scenario. The underlying reason for this is that a smaller |κ| can be fully responsible for the measured abundance in the h 2 scenario. The last aspect is that λ in the h 2 scenario may reach about 0.2, while it is at most 0.1 in the h 1 scenario. This is because the cancellation effect in the SI scattering is usually significant in the h 2 scenario, and consequently, a larger λ is still allowed by DM DD experiments.

JHEP03(2022)203
In figure 3, the mass distributions of the singlet-dominated Higgs states and the SUSY particles relevant to a SUSY µ are shown by a series of violin plots, which combines the advantages of the box plot and probability density distribution plot [165]. This figure shows that all SUSY particles except forχ 0 5 tend to be lighter than 500 GeV, and in particular,χ 0 2 ,χ 0 3 , andχ ± 1 are lighter than 500 GeV for nearly all samples obtained in the scan. The fundamental reason for the phenomenon, besides the explanation presented before, arises from the fact that a low tan β is preferred to predict the h 2 scenario. This tendency, once combined with the requirement of a sizable a SUSY µ , will necessitate light SUSY particles. 7 Given that these electroweakinos can be richly produced at the LHC, they have been restricted by searching for multi-lepton signals. This issue will be intensively studied in the following. 7 Without the a SUSY µ requirement,χ 0 1 may be very massive (e.g., |mχ0 1 | > 300 GeV [116]). In this case, the LHC constraints are significantly weakened.

LHC constraints
To comprehensively study the constraints from the LHC search for sparticles on the obtained parameter points, the following processes were analyzed in the Monte Carlo (MC) event simulation as given by equations (3.4): 8 In the calculation, the cross-sections of √ s = 13 TeV were obtained at the next-to-leading order (NLO) by the package Prospino2 [166]. The MC events were generated by the package MadGraph_aMC@NLO [167,168]  For each point, 10 6 MC events were generated in the simulations, and the LHC analyses listed in table 1 were used to test it. In particular, the following LHC analyses were included in our study, which played a crucial role in constraining the scenario:    The quantity R was used to describe the LHC's limitation on the samples in the discussion. It is defined by R ≡ max{S i /S 95 obs,i }, where S i stands for the number of simulated events in the i-th signal region (SR) of all the included analyses, and S 95 obs,i represents corresponding observed 95% confidence level upper limit. Accordingly, without considering the involved uncertainties, R > 1 represents that the considered parameter point is excluded due to the inconsistency with the LHC limit. Otherwise, it is allowed by the LHC searches [108].
The collider simulation results given by CheckMATE implied that the LHC searches for SUSY set strong restrictions on the h 2 scenario of the µNMSSM. In order to display the analysis results clearly, the points were classified by NLSP's dominated component, which may beμ R ,ν µ ,B,W , orH. The Histograms of R-value distribution for the different NLSP  is expressed in terms of its deviation from the FNAL measurement center value and is shown by different colors: turquoise for (−3σ, −2σ), orange for (−2σ, −1σ), pink for (−1σ, 0σ), green for (0σ, 1σ), blue for (1σ, 2σ), and violet for (2σ, 3σ). types were displayed in figure 4 from left to right and top to bottom, respectively. Points colored by turquoise, orange, pink, green, blue, and violet correspond respectively to the cases that a SUSY µ are in the range of (−3σ, −2σ), (−2σ, −1σ), (−1σ, 0σ), (0σ, 1σ), (1σ, 2σ) and (2σ, 3σ). The results were also summarized in table 2, which includes the number of samples obtained by the scan (denoted N tot ), that satisfying R < 1 (denoted by N pass ), and the more detailed classification of N tot and N pass by the ranges of a SUSY µ . According to table 2 and figures 4, the following conclusions are inferred: • Among the five types of NLSP, theH-dominated NLSP is the easiest one for explaining the discrepancy in the h 2 scenario, and by contrast, theB-dominated NLSP  is the least preferred one (see N tot in table 2). The LHC restrictions are extremely strong in excluding parameter points for any type of NLSP (see N pass in the table). In particular, they are strengthened significantly once the scenario is required to explain the discrepancy at 3σ level (see N a SUSY µ in the table). Specifically, it takes dozens of parameter points withH-dominated NLSP and only few points withν µ -dominated NLSP to interpret the discrepancy at the 3σ level. This situation reflects the difficulty of the h 2 scenario in explaining the discrepancy. One fundamental reason comes from the fact that the scenario prefers a relatively small tanβ and µ tot , and hence moderately light sparticles are predicted to obtain a sizable a SUSY µ .
It was verified that, among the experimental analyses, the analysis 4 usually sets the tightest constraints. In the case that the parameter points predict sizable signals with four or more leptons, analysis 5 could also impose the strongest restriction.
• In the case ofμ R -orν µ -dominated NLSP, Wino-and Higgsino-dominated electroweakinos will decay mainly into leptonic final states via slepton and/or sneutrino, which proliferates the lepton signals. As a result, R can reach 100 for lots of points, which is shown on the top left and right panels in figure 4. In addition, the LHC constraints on theμ R -dominated NLSP point are usually tighter than those on thẽ ν µ -dominated NLSP point because neutralinos will decay byχ 0 i →μ ± µ ∓ →χ 0 1 µ + µ − for the former case and byχ 0 i →μ ± µ ∓ ,ν µ ν →χ 0 1 µ + µ − ,χ 0 1 νν for the latter case. The former case can produce more µ leptons.
• In the case ofW -dominated NLSP, most points predict R < 20, but in very rare case R may reach 70. Detailed study indicated thatχ 0 2 decays mainly byχ 0 2 →χ 0 1 Z ( * ) ,χ 0 1 h 1,2 for most points, andχ 0 2 →χ 0 1 µ + µ − ,χ 0 1 νν are the dominant decay only for a small portion of the points. It also indicated thatχ ± 1 decays mainly byχ ± 1 →χ 0 1 W ( * ) for nearly all points, and smuons decay in a complex way, e.g., any of the channels µ →χ 0 i µ (i = 1, · · · , 5),χ − 1,2 ν may be the dominant decay. It is notable that R in theW -dominated NLSP case can not be exceedingly large. This conclusion comes from the fact that theW -dominated electroweakinos are forbidden to decay into sleptons directly, and thus, even in the optimum case, the lepton signal from the decayχ 0 2 →μ ± * µ ∓ →χ 0 1 µ + µ − is not much larger than the other final states. Consequently, pp →χ 0 2χ ± 1 , which is the largest sparticle production process, can not generate tri-lepton signal events efficiently. This feature results in a smaller signal rate than theB-dominated NLSP case, where the Wino-dominated electroweakinos may decay far dominantly into leptons. • By considering theH-dominated NLSP case it was found that the decay modes ofχ 0 2 ,χ ± 1 andμ 1,2 are similar to those of theW -dominated NLSP case, and the LHC constraints tend to be weaker than the other cases. This observation may be understood from four aspects [67]. First, since theH-dominatedχ 0 2,3 andχ ± 1 can not decay into sleptons, the leptonic signal rate is usually much smaller than the case whereμ R orν µ acts as NLSP. Second, the collider sensitive signal events are often diluted by the complicated decay chains of sparticles, given that heavy sparticles prefer to decay into the NLSP or other non-singlet-dominated sparticles first. They are diluted also by the decaysχ 0 2,3 →χ 0 1 h s ,χ 0 1 h, given that Br(h s /h → ± ∓ ) is much smaller than Br(Z → ± ∓ ). Third, the interpretation of ∆a µ requires that all crucial sparticles are usually in several hundred GeVs for a not too large tan β. Thus, for the parameter points surviving the LHC constraints, the mass splitting between sparticles is not large enough to produce high-p T signal objects, which can be significantly distinguished from the background in the collider. Last, in some rare cases, the leptonic signal of SUSY may mainly come from the NLSP. For this situation, the discussion of the LHC constraints can be simplified by considering the system that only contains NLSP and LSP. From this, it is evident that the constraints on thẽ H-dominated NLSP case are significantly weaker than those on theW -dominated NLSP case. In fact, we once scrutinized the property of all the samples surviving the LHC constraints. It was found that the above factors applied to these samples.

JHEP03(2022)203
In order to emphasize the characteristics of the parameter point withH-dominated NLSP, two benchmark points, P1 and P2, are chosen to present their detailed information in table 3. Both points satisfy the DM constraints and can explain the a µ discrepancy at 2σ level. The P1 point survives the LHC constraints, while the P2 point has been excluded by the LHC search for SUSY. These two benchmark points verify part of the discussions in this work.  Finally, it should be noted that the LHC search for τ -leptons plus missing momentum signal, such as the ATLAS analyses in [178] and [179], was not considered because of the massiveτ assumption in this study. Specifically, the assumption implies that the τ -leptons mainly come from the decay of the W/Z or Higgs bosons, which are the decay products of parent sparticles. For the former case, the final states containing e and/or µ are more efficient than the τ final state in restricting SUSY mass spectrum since all the lepton signal rates are roughly equal. For the latter case, the τ signal is usually less crucial in SUSY search because the branching ratio of the Higgs decay into ττ is significantly small (in comparison withτ decay). With the codes for the analyses in [178] and [179], which were implemented in our previous work [164] and CheckMATE-2.0.29, respectively, we studied JHEP03(2022)203 their prediction of R for the two benchmark points. We found that the analyses do not affect the results in table 3. As an alternative, ifτ or all sleptons are assumed to the NLSP (see, e.g., [41]), the production rate of the e/µ final states will be affected. In this case, R should be recalculated. In particular, the experimental analysis of the τ final state must be included in the study. It is expected that the LHC constraints are still strong because the h 2 scenario is featured by moderately light Higgsinos.

Explaining ∆a µ in the h 2 scenario of GNMSSM
The impact of the muon g-2 anomaly on the h 2 scenario of the GNMSSM is studied in this section. For this purpose, the parameter space including |µ | ≤ 1 TeV, −10 6 TeV 2 ≤ m 2 S ≤ 10 6 TeV 2 , and that in eq. (3.1), were scanned in a way similar to what we did in section 3. It was found that the DM was Singlino-dominated for all the obtained samples, and it annihilated mainly by a resonant Z, h, or h s /A s to obtain the measured abundance. These channels contributed to the total Bayesian evidence by about 43%, 19.6%, and 37%, respectively, before the MC simulations were implemented. The basic reason for such a behavior is thatχ 0 1 , m hs and m As in the GNMSSM can be changed freely by tuning µ , A κ , and m 2 S , respectively. Thus, the annihilations could easily happen. Given that the GMSSM might have different key features from the µNMSSM, various PL maps of the GNMSSM were surveyed in this study.
In figure 5, the two-dimensional profile likelihood function was projected onto tan β − µ tot and M 2 − µ tot planes. They show that the GNMSSM results are quite similar to the µNMSSM predictions. In particular, µ tot and mχ± 1 have upper bounds of about 310 GeV and 300 GeV, respectively. 9 The fundamental reason for this, as was emphasized before, is that the h 2 scenario prefers moderately small µ tot and tan β to predict m h 1 125 GeV and h 2 to be SM-like, which was verified by the posterior probability distribution function of the samples obtained from the scan. This characteristic, once combined with the requirement to explain the anomaly, will entail certain moderately light sparticles. In figure 6, the violin diagrams for the mass spectrum of the singlet-dominated Higgs bosons, the electroweakinos and µ-type sleptons are shown. The profiles for the sparticles are quite similar to those in figure 3 for the µNMSSM results, except that |mχ0 1 | can be as low as several GeV. This difference mainly comes from the DM annihilation mechanisms, and it usually makes the LHC's constraints much stronger. 9 We inferred that the 2σ CI of the µNMSSM covers a broader region on the parameter planes than that of the GNMSSM by comparing figure 5 with figure 1. This is contrary to the common sense that the former should be narrower than the latter since the parameter space of the µNMSSM is only a subset of the GNMSSM's parameter space. This phenomenon originates from the MultiNest algorithm utilized in the scan, which mainly collects the samples contributing significantly to the Bayesian evidence [135]. The parameter points of the µNMSSM correspond to µ = 0 and m 2 S = 0, and are relatively unimportant for the evidence, which was verified by the study of the two-dimensional posterior probability distribution function, P (µ , m 2 S ) [163]. Thus, only a few of them were considered in the sampling. It is expected that, with the increase of the setting n live , more samples of the GNMSSM will be collected, which will broaden the CI regions [164]. This process, however, is very computationally expensive, since a high-dimensional parameter space is surveyed.

JHEP03(2022)203
In table 4, the numbers of the samples surveyed by MC simulations and those passing the LHC constraints were presented in a way similar to table 2. This table shows that the LHC analyses have strongly constrained the parameter space of the GNMSSM.

Conclusion
The recent measurement of a µ by the FNAL corroborates further the long-standing discrepancy of a Exp µ from a SM µ . It can not only reveal useful information of the physics beyond the SM, but also place strong restrictions on certain theories. Recently, implications of the discrepancy were comprehensively discussed with respect to the GNMSSM, which is a theory that has the following attractive features: it is free from the tadpole problem and the domain-wall problem of the Z 3 -NMSSM, and it is capable of forming an economic secluded DM sector to naturally yield the DM experimental results [116]. It was found that the h 1 scenario of the GNMSSM could easily and significantly weaken the constraints from the LHC search for SUSY. It also predicted more stable vacuums than the Z 3 -NMSSM. As a result, the scenario can explain the discrepancy in a broad parameter space that is consistent with all experimental results, and at same time keeps the electroweak symmetry breaking natural [67]. By contrast, it is difficult for the popular MSSM and Z 3 -NMSSM to do this.
These theoretical advantages inspired us to consider the h 2 scenario of the GNMSSM, which is another well-known realization of the theory. It was shown by analytic formulae that, in order to obtain m h 1 125 GeV and an SM-like h 2 without significant tunings of relevant parameters, the scenario prefers a moderately light µ tot and tan β 30. This characteristic, if combined with the requirement to account for the anomaly, will entail some light sparticles, and sequentially make the LHC constraints rather tight. In this work, this speculation was tested using numerical results. Specifically, a special case of the GNMSSM called µNMSSM was first studied by scanning its parameter space with the MultiNest algorithm and considering the constraints from the LHC Higgs data, the DM experimental results, the B-physics observations, and the vacuum stability. Then, the samples obtained from the scan were surveyed by the LHC analyses in sparticle searches. Through sophisticated MC simulations, it was found that only a dozen of the samples, among about twenty thousand, passed the constraints, which corresponded to about 0.04% of the total Bayesian evidence. Given that the scan results have statistical significance, we conclude that the h 2 scenario of the µNMSSM is tightly constrained if it is intended to explain the anomaly. A similar study was carried out for the GNMSSM, and it was found that a smaller portion of the samples (about 0.008% of the total evidence) satisfied the LHC constraints. This difference arises from DM annihilation mechanisms: for the former case, the Singlino-dominated DM achieved the measured abundance mainly through the processχ 0 1χ 0 1 → h s A s , while for the latter case, it was through a resonant Z, h, or h s /A s annihilation to obtain the abundance. Since the latter case usually predicts a relatively light DM, the LHC constraints are stronger.
This work extends the research in [99] by considering a more general theoretical framework with more advanced and sophisticated research strategies. As a result, the conclusions JHEP03(2022)203 obtained in this work are more robust than those of the previous work, and apply to any realizations of the NMSSM.

A DM-nucleon scattering in the MSSM
In this section, we use analytic formulae to focus on Bino-dominated DM and study DMnucleon scatterings for three typical cases.
We begin with the neutralino mass matrix given by [180]: where g 1 = 2M Z s W /v, s β ≡ sin β and c β ≡ cos β. In terms of neutralino mass, mχ0 i , the eigenvectors are then exactly formulated by Parameterizing the couplings of the DM to the SM-like Higgs boson, h, and Z boson as the following form [26,93] L MSSM Cχ0

JHEP03(2022)203
we obtain [112] Cχ0 where α is the mixing angle of CP-even Higgs fields in forming mass eigenstates [94]. In the decoupling limit of the Higgs sector, i.e., m A v, the SI and SD cross-sections of the DM with nucleons are approximated by [112]: with C p 2.9 × 10 −41 cm 2 for protons and C n 2.3 × 10 −41 cm 2 for neutrons. In the following, we assume tan β 1 and m A v so that α β − π/2 [94], and investigate the dependence of σ SĨ to obtain the measured DM abundance, and C 1 is approximated by: Consequently, C hχ 0 where µ is parameterized by |µ| = (1 + δ 1 )|mχ0 1 |, with δ 1 denoting a small positive dimensionless number. These approximations indicate that the DM-nucleon scattering rates increase monotonously as the DM becomes heavier and/or |µ| departures from |mχ0 with the XENON-1T data on SI cross-section. For mχ0 1 = 300 GeV, it must be less than about 0.12. The Higgsinos in this narrow mass region contribute significantly to lepton signals at the LHC and the DM relic abundance. They also affect a µ . Consequently, such a situation needs tuning to satisfy all experimental constraints.
• Case II: the DM co-annihilated with the Wino-dominated electroweakinos to obtain the measured abundance, and |µ| is much larger than |mχ0 1 |. This case predicts [112]: and consequently the scattering rates decrease monotonously with the increase of |µ|.
• Case III: the DM co-annihilated with the Higgsino-dominated electroweakinos to obtain the measured abundance, and |M 2 | is much larger than |mχ0 1 |. For this case: where the second approximation is obtained by assuming |mχ0