Improved (g − 2)μ measurement and singlino dark matter in μ-term extended ℤ3-NMSSM

Very recently, a Fermilab report of muon g− 2 showed a 4.2σ discrepancy between it and the standard model (SM) prediction. Motivated by this inspiring result and the increasing tension in supersymmetric interpretation of the anomalous magnetic moment, it is argued that in the general next-to-minimal supersymmetric standard model (GNMSSM), a singlino-dominated neutralino can act as a feasible dark matter (DM) candidate in explaining the discrepancy naturally. In this case, the singlino-dominated DM and singlet-dominated Higgs bosons can form a secluded DM sector with χ~10χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\overset{\sim }{\chi}}_1^0{\overset{\sim }{\chi}}_1^0 $$\end{document}→ hsAs responsible for the measured DM relic abundance when mχ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_{{\overset{\sim }{\chi}}_1^0} $$\end{document} ≳ 150 GeV and the Yukawa coupling κ is around 0.2. This sector communicates with the SM sector by weak singlet-doublet Higgs mixing, so the scatterings of the singlino-dominated DM with nucleons are suppressed. Furthermore, due to the singlet nature of the DM and the complex mass hierarchy, sparticle decay chains in the GNMSSM are lengthened in comparison with the prediction of the minimal supersymmetric standard model. These characteristics lead to sparticle detection at the Large Hadron Collider (LHC) being rather tricky. This study surveys a specific scenario of the GNMSSM, which extends the ℤ3-NMSSM by adding an explicit μ-term, to reveal the features. It indicates that the theory can readily explain the discrepancy of the muon anomalous magnetic moment without conflicting with the experimental results in DM and Higgs physics, and the LHC searches for sparticles.


Introduction
Recently, experimentalists at Fermilab National Accelerator Laboratory (FNAL) reported the most accurate measurement of the muon anomalous magnetic moment a exp µ (FNAL) [1]. Combined with the previous Brookhaven National Laboratory E821 result a exp µ (BNL) [2], the statistical average a exp which reveals a 4.2σ discrepancy with the SM prediction a SM µ = 116591810(43)×10 −11 : Besides, the Run-I result indicates that the future complete results of Fermilab and/or Japan Proton Accelerator Research Complex (J-PARC) experiments are very likely to confirm the excess of a µ at 5σ discovery level. This expectation implies that the longstanding discrepancy of the muon anomalous magnetic moment between the SM prediction and experimental measurements, ∆a µ , may be the most promising hint of the new physics beyond the SM. As the best candidate theory for new physics, supersymmetric models predict the scalar partners of muon and µ-type neutrino. It was speculated about twenty years ago that the source of the observed ∆a µ might be the quantum effect contributed by JHEP09(2021)175 these supersymmetric particles (sparticles) [24][25][26]. Along this direction, numerous studies have been carried out in the minimal supersymmetric standard model (MSSM) and its extensions (see, e.g., [27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43]). At present, low-energy supersymmetric theories are subjected to the increasingly tight constraints from LHC experiments and DM search experiments. In the MSSM, the lightest neutralinoχ 0 1 , if it is the lightest supersymmetric particle (LSP), can act as a DM candidate accounting for the Planck measured relic density [44]. The likelihood analysis of the phenomenological MSSM in 11 free-parameter space showed that the DMχ 0 1 must be bino-dominated within the 1σ confidence level (see figure 12 and table 6 of ref. [45]). Aiming at the proper recasting of the LHC Run-II data for sparticle search and assuming thatχ 0 1 provides the full DM relic density, recent studies investigated the phenomenology of the bino DM co-annihilating with wino-dominatedχ ± 1 or with sleptons˜ L/R [46,47]. The result revealed that it is difficult to obtain the correct DM relic density and experimentally compatible DM-nucleon scatterings in the natural parameter space to explain ∆a µ . 1 The more recent research opened up the possibility of wino and higgsino DM by giving upχ 0 1 to provide the full DM relic density [48]. In the next-to minimal supersymmetric standard model with a Z 3 symmetry (Z 3 -NMSSM), which is another economic realization of supersymmetry, the DMχ 0 1 may be either bino-dominated or singlino-dominated [49][50][51][52][53][54][55][56][57][58][59]. The situation of the theory is similar to that of the MSSM, i.e., the parameter space to explain ∆a µ naturally have been constrained tightly by the LHC and DM experiments [60][61][62].
Considering the increasing tension between natural interpretations of a µ and the experimental constraints, the present paper aims to study the combined constraints on the singlino-dominated DM scenario in the GNMSSM from the DM relic density, spindependent (SD) and spin-independent (SI) direct detection experiments, sparticle direct searches at the LHC, and the existing muon (g − 2) measurement in eq. (1.2). The rest of this paper is organized as follows. In the next section, the basics of the GNMSSM and the annihilation mechanism of singlino-dominated DM is reviewed. In section 3, the relevant parameter space is scanned in a specific scenario of the GNMSSM, which extends the Z 3 -NMSSM by adding an explicit µ-term and is called µNMSSM hereafter. It is found that (g − 2) µ can be properly interpreted without conflicting with any experimental observations. In section 4, the constraints from the LHC searches for supersymmetry on the interpretation are studied comprehensively by specific Monte Carlo simulations to reveal their features. The theory's capabilities to relax the experimental constraints are emphasized. Conclusions about the µNMSSM interpretation of ∆a µ are drawn in section 5.

The basics of GNMSSM
Compared with the MSSM, the NMSSM introduces a singlet Higgs superfieldŜ. Given the superfield composition, the general form of the NMSSM superpotential is [63,64]

1)
1 When we mention the term "natural" in this work, it means that higgsinos lighter than about 500 GeV are preferred to predict Z boson mass without causing serious offsets between the different contributions of mZ .

JHEP09(2021)175
where the Yukawa terms W Yukawa are the same as those in MSSM. The scenario of µ = ν = 0 possess an accidental Z 3 -symmetry, and the theory is defined in a scaleinvariant form. Clearly, if the dimensional parameters µ and ν are non-vanishing, they should be at the weak or supersymmetry-breaking scale to break the electroweak symmetry without fine-tuning. These Z 3 broken terms µ and ν are introduced in some works to solve the tadpole problem [64,65] and the cosmological domain-wall problem [66][67][68].
Past studies [66,[69][70][71][72] demonstrated that the electroweak scale µ and ν may come from the discrete Z R 4 or Z R 8 symmetry breaking at high-energy scale. Furthermore, the scaleinvariant Z 3 -NMSSM allows it to be embedded into canonical superconformal supergravity in the Jordan frame. The superconformal symmetry of matter multiplets can be broken via a non-minimal interaction χĤ u ·Ĥ d R, where R is the supersymmetric version of the Ricci scalar. As shown in ref. [73], the dimensionless coupling χ can drive inflation in the early Universe, and can also provide a µ term correction to the Z 3 -NMSSM superpotential, where µ = 3 2 m 3/2 χ, with m 3/2 denoting the gravitino mass. This work treats µ and ν as free theoretical input parameters, irrespective of their physical origin. Particularly for the sake of brevity without loss of generality, a specific scenario of the GNMSSM characterized by ν ≡ 0 is investigated [74]. This so-called µ-term extended Z 3 -NMSSM, i.e., µNMSSM, is defined by its superpotential and the soft breaking Lagrangian as follows [75,76] where H u , H d , and S are the scalar parts of superfieldsĤ u ,Ĥ d , andŜ, respectively. After the electroweak symmetry breaking, the neutral Higgs fields acquire non-zero vacuum expectation values (vevs), GeV. In practice, the free input parameters of the Higgs sector can be taken as follows 2 In the field convention that [77,78], the elements of the CP -even

JHEP09(2021)175
Higgs boson mass matrix M 2 S in the bases (H NSM , H SM , Re(S)) are read as (2.5) Dropping the Goldstone mode, those elements for CP -odd Higgs fields in the bases (A NSM , Im(S)) are given by (2.6) Three The fermion parts of Higgs superfields H u ,H d ,S and gauginos B ,W form five neutralino states and two chargino states, and they are referred as electroweakinos (EWinos) in general. The symmetric neutralino mass matrix in the gauge eigenstate bases of where the abbreviations s W = sin θ W and c W = cos θ W are used, with θ W being the weak mixing angle, and s β = sin β and c β = cos β. Similarly, the chargino mass matrix in the JHEP09(2021)175 After diagonalization, one arrives at the neutralinoχ 0 i and charginoχ ± i as mass eigenstates, with increasing mass for a higher label i.
The smuon mass matrix in the gauge eigenstate bases (μ L ,μ R ) is given as where A µ , mμ L , and mμ R are muon-type soft breaking parameters. Eq. (2.10) indicates that the left-right mixing term is dominated by (µ + µ eff ) tan β, so A µ is fixed as zero in the following. The muon-type sneutrino mass is (2.11)

Muon g − 2 in µNMSSM
The SUSY contribution a SUSY µ , in which the muon lepton number is carried byμ orν µ in the loops, 3 can be the source of ∆a µ [24,93,94]. The expression of a SUSY µ in the µNMSSM is similar to that in the MSSM [79], which is given by [24]: where i = 1, · · · , 5, k = 1, 2, and l = 1, 2 denote the neutralino, chargino, and smuon index, respectively, and In the NMSSM, besides the contributions from the SM particles in the loops, heavy doublet-dominated Higgs bosons can also mediate the contribution to aµ. However, after considering the constraints from the LHC searches for extra Higgs bosons and the measurements of the branching ratios for B → Xsγ and Bs → µ + µ − , this contribution is negligibly small because the Higgs bosons should be very massive for a large tan β. In addition, although the contribution from light singlet-dominated Higgs bosons might reach O(10 −10 ) as pointed out in [79], it is negligible in this study. The reason is that the constraints from the DM direct detection experiments strongly favor a small λ for a singlino-dominated DM, and consequently, the singlet-doublet Higgs mixings and their relatedμµAs andμµhs couplings are suppressed significantly. It was testified numerically that the total Higgs-mediated contributions are less than 10 −10 for the samples obtained in this study. Besides, a SUSY µ has two-loop contribution [80][81][82][83][84][85][86][87][88][89][90][91]. A recent analysis revealed that the correction is less than 4 × 10 −10 [92]. We anticipate that it can be further suppressed if the restrictions from the LHC search for SUSY and DM physics are considered.

JHEP09(2021)175
Here, 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 * X 2×2 V c † = m diag χ ± . The kinematic loop functions F (x)s depend on the variables x il ≡ m 2 /m 2 νµ , and are given by (2.14) (2.17) for mass-degenerate sparticle case. It is instructive to point out that, although the µNMSSM predicts five neutralinos, the singlino-induced contribution never makes sense, and the µNMSSM prediction of a SUSY µ is roughly the same as that of the MSSM except that the µ parameter of the MSSM should be replaced by µ + µ eff . This feature can be understood by noting the fact that the field operator for a µ involves chiral flipped Muon leptons and adopting the mass insertion approximation in the calculation of a µ [93]. Specifically, the contributions to a SUSY µ in the MSSM can be classified into four types: "WHL", "BHL", "BHR", and "BLR", where W , B, H, L, and R stands for wino, bino, higgsino, left-handed and right-handed smuon field, respectively. They arise from the Feynman where the loop functions are given by and they satisfy f C (1, 1) = 1/2 and f N (1, 1) = 1/6. In the µNMSSM, the singlino fieldS can also enter the insertions, but because both theW −S andB 0 −S transitions and thē µSμ L,R couplings vanish, it only appears in the "WHL", "BHL" and "BHR" loops by two JHEP09(2021)175 more insertions at the lowest order, which corresponds to theH 0 d −S andS −H 0 d transitions in the neutralino mass matrix in eq. (2.8), respectively. Since a massive singlino and a small λ are preferred by DM physics (see discussion below), the singlino-induced contribution can not be significantly large. 4 It should be emphasized that, although a SUSY µ has roughly the same properties in the µNMSSM and MSSM, it is subject to significantly relaxed experimental and theoretical limitations in the µNMSSM (see following discussions). Thus, the µNMSSM is more readily to explain the a µ discrepancy, which is the main motivation of this work. In addition, the WHL contribution is usually much larger than the other contributions ifμ L is not significantly heavier thanμ R .

Singlino-dominated DM
This work aims to reveal DM physics in the natural parameter space of interpreting ∆a µ . In the Z 3 -NMSSM, the properties of the singlino-dominated DM are mainly determined by three parameters: λ, µ, and mχ0 1 [62]. The Yukawa coupling κ is related with mχ0 1 and satisfies 2|κ| < λ to ensure thatχ 0 1 is singlino-dominated. The DM obtained the correct relic abundance mainly through co-annihilation with higgsinos in the early Universe, and λ should be less than 0.1 to suppress the DM direct detection rate [61,62]. As a result, the parameters in the Z 3 -NMSSM are highly constrained. In contrast, due to the introduction of the additional µ term in the µNMSSM, the DM properties are described by four Higgs parameters: λ, κ, µ, and mχ0 1 2κ λ µ eff , and a singlino-dominated DM does not require 2|κ| < λ [74]. This causes the singlino-dominated DM properties in the µNMSSM to be significantly different from those in the Z 3 -NMSSM.
Assuming a standard thermal history of the Universe with a singlino-dominated DM candidate in thermal equilibrium until it freezed out, DM annihilation rate must be sufficiently large to be consistent with the Planck observation Ω DM h 2 = 0.120 ± 0.001. Various processes, as follows, may provide a sufficiently large annihilation cross-section [74].
This process is mainly carried out through s-channel exchange of Z and CP -odd Higgs A s and t-channel exchange of neutralinos [52,97]. As shown in eq. (2.19) in [74], the t-channel annihilation cross-section is roughly proportional to κ 4 . In the t-channel dominated case, the parameters should satisfy the relationship [52] |κ| ∼ 0.15 × to obtain the measured density. For the case of m As 2 × mχ0 1 , the cross-section is enhanced by the s-channel pole. However, the A s pole enhancement may lead to a very light singlet Higgs boson h s . This h s must satisfy the constraints from LEP Higgs searches and predict a small BR(h → h s h s ) to fit SM-like Higgs data.

JHEP09(2021)175
Besides, the SI direct detection rate mediated by h s must satisfy the current bound from Xenon-1T [98] and PandaX-II [99]. A comparative study of light h s scenario (m hs < m h ) and heavy h s scenario (m hs > m h ) indicated that the latter scenario is preferred by a Bayes factor 2.42 [74]. 5 For the above reasons, only the m hs > m h scenario will be considered (i.e., h 1 = h and h 2 = h s ) in the following numerical study.
This process is mediated by the s-channel exchange of Z and Higgs bosons. It is significant only when λ is sizable, but this usually leads to a sizable DM-nucleon scattering rate [62].
• Co-annihilation with EWinos. In principle, this channel affects the abundance when the mass splitting betweenχ 0 1 and a co-annihilation particle is less than approximately 10% [62].
• Co-annihilation with smuonsμ L /μ R or µ-type sneutrinoν µ . Within the interpretation of ∆a µ , smuons should not be too heavy, so this co-annihilation channel can be opened when mμ mχ0 1 .
In the heavy squark limits, SI DM-nucleon scattering is mainly from t-channel exchange of CP -even Higgs bosons, and the cross-section is given as [102,103] is the reduced mass of the DM-nucleon system, and C N N h i is the coupling of a Higgs boson with a nucleon, In the heavy H case, the SI cross-section is dominated by two light Higgs contributions, and the h-mediated contribution is usually significantly larger than the h s -mediated contribution. The couplings Cχ0 It can be interpreted as a measure of the strength of evidence in favor of one theory among two competing theories [100,101]. A factor of 2.42 is generally considered as a decisive result to indicate the preference of one theory over another.

JHEP09(2021)175
are given by [74] Cχ0 In contrast, the SD scattering cross-section takes the following simple from [104,105]: with C p 4.0 for the proton and C n 3.1 for the neutron. The above formulae indicate that the DM direct detection rate is positively related to the Higgs coupling λ, roughly in terms of λv/( √ 2(µ + µ eff )).

Research strategy
The following relevant parameter space was scanned using the MultiNest technique with the setting n live = 30000 [106,107] 6 to explore the features of the µNMSSM interpretation JHEP09(2021)175 of ∆a µ : All of the input parameters are flatly distributed beforehand. Other SUSY parameters, such as A λ , the parameters for the first and third generation sleptons, three generation squarks, and gluino, are fixed at 2 TeV. In the numerical calculation, the model file of the µNMSSM is constructed through the package SARAH-4.14.3 [108][109][110][111]. The particle mass spectra and low-energy observables, such as a SUSY µ , are generated by the codes SPheno-4.0.3 [112,113] and FlavorKit [114]. The DM relic density and direct detection cross-sections are calculated using package MicrOMEGAs-5.0.4 [115][116][117][118][119][120]. The following likelihood function if restrictions unsatisfied.
was constructed to guide the scan, where the restrictions on each sample include: • Higgs data fit. As mentioned above, the lightest CP -even Higgs boson h 1 corresponds to the SM-like state, and this state must satisfy the constraints from the LHC data using the code HiggsSignal-2.2.3 [121]. The p value in the fit is required to be larger than 0.05, which implies that the Higgs property coincides with the data at the 95% confidence level. The extra Higgs states must pass the constraints from the direct searches at the LEP, Tevatron, and LHC, which is implemented by the code HiggsBounds-5.3.2 [122].
• Constraints from B-physics observation BR(B s → µ + µ − ) and BR(B s → X s γ) are taken into consideration [100]. These two B-physics observables must fall into the 2σ bounds.
• DM relic density constraints. The samples are required to have a neutralino LSP, and the predicted relic density of samples must agree with the Planck measurement [44], 7 i.e., 0.096 ≤ Ωh 2 ≤ 0.144.
• Direct detection limits on DM. The SI DM scattering cross-section σ SI p is required to be below the constraint from the Xenon-1T experiment [98]. The SD cross-section σ SD n is needed to pass the limits of the Xenon-1T report [123].
• Constraints from LHC sparticles direct searches. In the µNMSSM explanation of ∆a µ , the EWinos and sleptons can be produced at the LHC, and thus restricted JHEP09(2021)175 by the searches for multi-lepton signals. The constraints implemented in this step are produced using the code SModelS-1.2.3 [124], which contained the experimental analyses in simplified models that are summarized in appendix A.
For each sample obtained in the scan, the stability of its vacuum for the scalar potential consisting of Higgs and the last two generation slepton fields was finally checked by the code Vevacious [125,126]. Compared with the MSSM, the vacuum in the µNMSSM is more stable due to the addition of the singlet Higgs field as dynamical degree of freedom, especially in the case of a small λ and a large µ eff [75], but tremendously large softbreaking trilinear coefficients may still cause its destabilization [127]. Generally speaking, the vacuum destabilization occurs in the following situations: • One or more tachyonic Higgs masses are predicted. Tachyonic masses are related to the fact that the electroweak point, around which the potential is expanded, is not a local minimum in the scalar potential, but rather resembles a saddle point or even a local maximum [75]. In this case, the true vacuum lies at a deeper point along this tachyonic direction. Consequently, the true vacuum has vevs different from the input values, and the electroweak breaking condition does not select a minimum.
In this study, it was found that more than half of the samples encountered in the scan correspond to the tachyonic mass case. They are abandoned directly in the calculation since they can not predict physical mass spectra.
• The formation of non-standard minima which break the electric and/or color charges, known as charge-and color-breaking (CCB) minima. In the MSSM, the following condition should be satisfied to avoid the CCB vacuum [128,129]: where y ,eff with = τ, µ are the lepton Yukawa couplings including radiative corrections [130], η τ 1 and η µ 0.88. If the singlet-field direction were neglected, the formula could be directly applied to the µNMSSM, keeping v s = 0 GeV and replacing µ → µ + µ eff . However, with the singlet as dynamical degree of freedom, the stability of the electroweak vacuum is improved as the only singlet−slepton contribution is actually a quadrilinear term λY SH 0 u˜ * L˜ R , and the occurrence of a vacuum with ˜ L = 0 and ˜ R = 0 can enhance the scalar potential [75].
It should be pointed out that, for the parameter space in eq. (3.1), samples obtained in the scan satisfy the inequality automatically. The reason is the ATLAS measurement of the properties of the discovered Higgs particle has required 0.77 ≤ Y τ,eff /Y τ,SM ≤ 1.37 [131] and Y µ,eff /Y µ,SM ≤ 2.4 [132] at 2σ confidence level. This
• The electroweak vacuum of the sample, which was called the desired symmetry breaking (DSB) vacuum in literature [75,127], corresponds to a local minimum of the scalar potential instead of a global minimum. In this case, the DSB vacuum could undergo quantum tunneling to the true vacuum. If the tunneling time is short enough in reference to the age of Universe, the DSB vacuum would decay completely [125,126]. Such vacuum was called metastable but short-lived. Evidently, the occurrence of the metastable vacuum depends on the contour of the scalar potential, which is mainly decided by the parameters λ, κ, µ eff , and A κ [75]. The calculation of the code Vevacious indicated that only about 1% of the scanned samples predict short-lived vacuum.

Numerical results
The discussion begins with the normalized two-dimensional profile likelihood (PL), L(θ 1 , θ 2 ), for the likelihood function in eq. . Within the 1σ level in ∆a µ , the contours in the (µ + µ eff ) − M 2 plane imply that the mass of charginoχ ± 1 cannot be larger than approximately 700 GeV (in particular, the higgsino mass may be less than 500 GeV to predict m Z naturally), while both mμ L and mμ R can be larger than 800 GeV.
As shown in figure 2, λ is less than 0.1, the absolute values of κ vary from 0.1 to 0.4, and all of the samples meet the condition mχ0 1 (m hs + m As )/2. This feature indicates that the strengths of DM coupled to other non-singlet fields are relatively weak, and the main DM annihilation channel isχ 0 1χ 0 1 → h s A s . Thus, the singlet-dominated particles form a secluded DM sector where the singlet Higgs states h s and A s act as the mediators between DM and SM particles [133], which was pointed out in our previous work [74]. Owing to the smallness of λ, the DM direct detection rates σ SI p and σ SD n in the ∆a µ favored parameter space can be far below the current detection limits. The above discussion shows that the DM phenomenology and natural interpretations of ∆a µ are very weakly connected JHEP09(2021)175 in the µNMSSM . This situation is significantly different from the MSSM case, where the bino-dominatedχ 0 1 mainly co-annihilated with the other sparticles to obtain the measured relic density, and consequently the parameter space to explain ∆a µ is limited [46,47].
In figure 3, the mass distributions of the singlet Higgs states and SUSY particles are shown with violin plots. 8 From figure 3, one can find that the masses of sleptons and EWinos can be lower than 500 GeV. They can be produced at the LHC, and they are thus restricted by searching for multi-lepton signals. Compared with the MSSM prediction, the sparticles in the µNMSSM have the following distinct features.
• As shown in figures 1 and 2,χ 0 1 is moderately massive mχ0 1 > 150 GeV from the current results because it must be heavier than (m hs + m A S )/2 to proceed with χ 0 1χ 0 1 → h s A s . 8 A violin plot combines the advantages of the box plot and probability density distribution plot [134]. • Sinceχ 0 1 is singlino-dominated, and thus coupled very weakly to the other sparticles, heavy sparticles prefer to decay into next-to-LSP (NLSP) or next-next-to-LSP (NNLSP) first. Consequently, their decay chain is lengthened and their signals become complicated.

JHEP09(2021)175
• Since singlet-dominated Higgs bosons are preferred light, they may serve as sparticle decay products and enrich the decay channels of sparticles.
These characteristics make sparticle detection at the LHC rather tricky, and the corresponding constraints should be weaker than those in the simplified models and in the MSSM. This issue will be intensively studied in the following.

LHC constraints
Given that the capability of the SModelS package in implementing the LHC constraints is limited by its database and the strict prerequisites to use it, the LHC detection of JHEP09(2021)175 sparticles is further studied by Monte Carlo event simulations. 9 The following processes are considered in the simulation: The cross-sections at √ s = 13 TeV were obtained at the NLO using the package Prospino2 [135]. The signal events were generated by MadGraph_aMC@NLO [136,137] with the package PYTHIA8 [138] for parton showers, hadronization, and sparticle decay. Finally, the event files were put into CheckMATE-2.0.29 [139][140][141] with the embedded package Delphes [142] for detector simulation.
In the simulations, 10 5 signal events were generated for each µNMSSM parameter point, and the following latest LHC searches among the analyses mentioned before were found crucial in testing the samples: 9 More than 88000 samples were obtained in the previous scan. In order to save time and at the same time make our conclusions as general as possible, a smaller n live (n live = 3000) was chosen to repeat the scan and obtained about 7500 samples for studying the LHC constraints. These samples satisfy the vacuum stability requirement and can interpret ∆aµ within 2σ level. Besides, the analyses of √ s = 8 TeV pp collisions were not considered in this work because the LSP is relatively heavy so that their constraints on the theory are weak.
• Search for chargino-neutralino production with involved mass splittings near the electroweak scale in three-lepton final states in √ s = 13 TeV pp collisions with the ATLAS detector (see Report No. CERN-EP-2019-263) [144].
• Search for direct production of electroweakinos in final states with one lepton, missing transverse momentum and a Higgs boson decaying into two bb jets in pp collisions at √ s = 13 TeV with the ATLAS detector (see Report No. CERN-EP-2019-188) [145].
The R values obtained from the code CheckMATE were applied to implement the LHC constraints. Here, R ≡ max{S i /S 95 obs,i } for all the involved analyses, where S i represents the simulated event number of the ith signal region (SR), and S 95 obs,i is the corresponding 95% confidence level upper limit. Therefore, R > 1 indicates that the sample is excluded by the LHC searches if the involved uncertainties are not considered (see the discussion in footnote 9 of this work).
The collider simulation results via the CheckMATE show that the LHC searches for supersymmetry set clearly limits on the mass spectrum. In order to show details of the results, the samples are classified by their NLSP's dominant component, which may beν µ , µ R ,B,W , orH. In figures 4, 5, and 6, samples featured byν µ -,μ R -, andB-dominated NLSP, respectively, are projected on mχ0 LHC experiments with 0.67 ≤ R < 1 and R < 0.67, respectively. 10 In figures 7 and 8, similar diagrams are plotted for the samples withW -andH-dominated NLSP, respectively. The total number of each type of samples surveyed by specific Monte Carlo simulations and the corresponding number satisfying R < 1 are summarized in table 1, which are denoted by N tot and N pass , respectively. The lower limits of the representative parameters (µ + µ eff ), M 2 , mχ0 1 , mμ L and mμ R after considering the LHC constraints are also presented in the table.
From figures 4-8 and table 1, the following conclusions are inferred: • Due to the singlet nature of the LSP, sparticle prefers to decay first into lighter sparticles other than the LSP. Whenν µ orμ R acts as NLSP, wino-dominated and 10 0.67 ≤ R < 1 means that the sample's signal is close to exclusion, but a full accounting of uncertainties (originating from e.g., parton distribution function sets, the choice of renormalization and factorisation scale, the details of parton showering or the finite Monte Carlo statistics) would certainly place it within error bars [54]. On the other hand, if R < 0.67, the sample appears to be essentially compatible with the experimental results. It corresponds to the number of signal events which is below the 95% C.L. upper bound divided by 1.5. higgsino-dominated EWinos decay mostly via slepton and/or sneutrino into leptonic final states, which can enhance the multi-lepton signals. The sparticle's signal is similar to the prediction of the simplified models adopted by ATLAS and CMS collaborations in analyzing experimental data, and it is highly restricted by the current LHC searches.
• When the bino-dominatedχ 0 2 is NLSP, EWinos and sleptons prefer to decay into it with significant branching ratios becauseχ 0 2 couples to these particles with unsuppressed gauge interactions. For samples featured by ∆ ≡ mχ0 2 − mχ0 1 less than dozens of GeV,χ 0 2 appears as a missing track at the collider detector. The signal is similar to the MSSM prediction with a bino-dominated LSP, and it leads to strong constraints on the samples from the LHC searches for supersymmetry. However,χ 0 2 will decay byχ 0 2 →χ 0 1 h,χ 0 1 Z if ∆ > m h , and theχ 0 1 h channel is usually dominant. In this case, it is hard to detect the sparticles partially due to the complexity of the decay chain.
• For samples with wino-or higgsino-dominatedχ 0 2 acting as NLSP, the collider constraints are relatively weak. This observation comes from at least three facts. First, since the wino-or higgsino-dominatedχ 0 2 can not decay into sleptons, the leptonic sig- nal rate is usually much smaller than the case whereν µ orμ R 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 first. Third, the interpretation of ∆a µ requires that EWinos and smuons are in several hundred GeVs. So for most of the samples 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. For example, it was found that the samples featured by wino-dominated NLSP receive the weakest restriction, and most of the surviving samples roughly satisfy mχ0 In addition, it is noticeable that most of the surviving samples are characterized by µ + µ tot < 600 GeV, which can predict m Z naturally. It is a distinct difference between the µNMSSM and MSSM.
• The a SUSY µ enhancement factor tan β is apparently restricted by the current LHC • It should be noted that some of the latest LHC analyses, e.g., the analysis in [148], are not included in this work. However, one can make some rough estimations by thinking that the red points in figures 4-8 will be excluded by the latest or the near future LHC analyses. There is still a relatively large parameter space of the µNMSSM that can interpret ∆a µ without conflicting with the newest LHC analyses. So the LHC constraints considered in this article are of a reference value.
In order to emphasize the properties of the samples with higgsino-dominatedχ 0 2 and wino-dominatedχ 0 2 , two benchmark points, P1 and P2, are chosen to present their detailed information in table 2. Both points predict a SUSY µ values approximately equal to 2.51×10 −9 and pass all the experimental constraints. These two benchmark points verify our previous discussions.
Finally, the mass spectra of the sparticles surviving the LHC constraints are presented in figure 9. Comparing it with figure 3, it was found that the LHC constraints are very effective in excluding relatively light sparticles, especially light sleptons. In fact, this conclusion is also reflected in the last three panels of figures 4-8. Another important conclusion is that theW orH-dominatedχ ± 1 is always lighter than about 700 GeV. Taking into account its production cross-section at linear e + e − colliders [149], one can infer thatχ ± 1 is very likely to be discovered at the CLIC, whose collision energy can reach 3 TeV [150][151][152][153]. This point was recently emphasized by the authors of [46,47]. Moreover, as mentioned in [96,154], future high luminosity LHC can significantly extend the LHC Run-II's capability in sparticle detection. The preliminary analyses carried out in, e.g., [155][156][157], have proven this point. These machines provide an opportunity to test the µNMSSM interpretation of ∆a µ once the deviation of a exp µ from the SM prediction is confirmed. This issue will be studied in detail in our future work.

Summary
In this work, the phenomenology of the new Fermilab result of ∆a µ interpreted in the µNMSSM was investigated. The obtained results show the following features: • Compared with the MSSM or the Z 3 -NMSSM, the strong exclusivity from DM physics and natural interpretations of ∆a µ is weak in the µNMSSM.

JHEP09(2021)175
Mass (GeV) Test Test Figure 9. Same as figure 3, but for the samples surviving the LHC constraints.
• A singlino-dominated DM candidate is preferred in the interpretation. Owing to the smallness of the singlet-doublet Higgs coupling λ, the singlino-dominated neutralino and singlet-dominated Higgs bosons may form a secluded DM sector in which the annihilation channelχ 0 1χ 0 1 → h s A s is responsible for the measured abundance by adopting an appropriate singlet Yukawa coupling κ.
• The secluded sector communicates with the SM sector by weak singlet-doublet Higgs mixing, so the scatterings of singlino-dominated DM with nucleons are suppressed below the current experimental limits.
• The mass of DM mχ0 1 must be heavier than about 150 GeV to proceed with the annihilation, and must be lighter than about 550 GeV to explain ∆a µ at 1σ level.
• Owing to the singlet nature of DM and the complex mass hierarchy of sparticles, the decay chains of EWinos and sleptons are lengthened in comparison with the MSSM prediction. Moreover, the singlet Higgs bosons h s and A s in the final states of the searching channels weaken the LHC detection capability. These characteristics make sparticle detection at the LHC rather tricky.
This study shows that the proposed theory can readily explain the discrepancy of the muon anomalous magnetic moment between its SM prediction and experimentally measured value, without conflicting with DM and Higgs experimental results and the LHC searches for supersymmetry. Among the interpretations, it is remarkable that the higgsino mass is less than 500 GeV in most cases so that the Z boson mass can be naturally predicted.    Table 2. Detailed information of two benchmark points that agree well with all of the DM and Higgs experiments and predict a SUSY µ 2.51 × 10 −9 . Numbers after annihilation processes represent their fractions in contributing to total DM annihilation cross-section at freeze-out temperature. Numbers after sparticle decay channels denote their branching ratios.

A Fast simulation via SModelS
SModelS [124,[166][167][168][169][170] enables the fast interpretation of the LHC data via simplified model results from ATLAS and CMS searches for SUSY particles. It decomposes all of the signatures occurring in a given SUSY model into simplified model topologies 11 via a generic procedure. In practice, the cross-section upper limits and efficiency maps are used in re-interpreting each topology result. Compared with Monte Carlo simulation, SModelS is much easier and faster. It not only allows for re-interpreting searches of the cut-andcount methodology, but it also allows for other searches, such as those relying on boosted decision tree (BDT) variables. The power of SModelS comes from its superfast speed and 11 Simplified model topologies are also referred to as simplified model spectra. An event topology is defined by the vertex structure as well as the SM and BSM final states. In each topology, the intermediate Z2 odd BSM particles are characterized only by their masses, production rates, and decay modes.   its large and continuously updated database. The applicability of SModelS is limited by the simplified model results available in the database. Moreover, when the tested spectra split into many different channels, as is often the case in a complex model, the derived results are generally conservative.

JHEP09(2021)175
In this work, the samples are refined with the analyses at 13 TeV LHC in SModelS, which are summarized in table 3.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.