Current status of the light neutralino thermal dark matter in the phenomenological MSSM

Abstract

In a previous publication, we studied the parameter space of the phenomenological Minimal Supersymmetric Standard Model (pMSSM) with a light neutralino thermal dark matter ($M_{\tilde{\chi}_{1}^{0}}\leq M_{h}/2$) and observed that the recent results from the dark matter and collider experiments put strong constraints on this scenario. In this work, we present in detail the arguments behind the robustness of this result against scanning over the large number of parameters in pMSSM. The Run-3 of LHC will be crucial in probing the surviving regions of the parameter space. We further investigate the impact of light staus on our parameter space and also provide benchmarks which can be interesting for Run-3 of LHC. We analyse these benchmarks at the LHC using the machine learning framework of XGBOOST. Finally, we also discuss the effect of non-standard cosmology on the parameter space.

I Introduction

Since the discovery of the 125 GeV Higgs boson ($h$) by the ATLAS and CMS collaborations in 2012, the precise quantification of its non-standard couplings has been a major cornerstone of the new physics search program at the LHC. Until now, clear evidence of physics beyond the Standard Model (SM) is yet to be observed at the LHC. Within experimental uncertainties, measurements indicate that the observed properties of $h$ are consistent with the expectations for the SM Higgs boson. While the couplings of the Higgs boson with the third-generation charged fermions and gauge bosons have been measured with considerable precision, uncertainties in the second-generation Yukawa coupling measurements are gradually reducing with improved statistics. Nonetheless, the current data still allow non-standard decays for the discovered Higgs boson. An exciting aspect of non-standard interactions of Higgs bosons is their decay into an invisible final state. Recent analyses of the LHC Run-II data by the ATLAS and CMS collaborations have constrained the branching fraction for Higgs invisible decay to 11% [1] and 15% [2] at $95\%$ CL, respectively. The invisible final states can be dark matter (DM) candidates, given that they are stable or have a lifetime larger than the age of the Universe, thus transmuting Higgs invisible searches at the LHC to potential probes for the elusive Dark Matter.

The R-parity conserving (RPC) Minimal Supersymmetric Standard Model (MSSM), with no explicit terms for baryon number and lepton number violation, provides a stable lightest supersymmetric particle (LSP), typically the neutralino $\tilde{\chi}_{1}^{0}$, which can be a WIMP (weakly interacting massive particle) DM candidate. Charged under the $SU(2)_{L}\times U(1)_{Y}$ gauge group, the LSP neutralino can typically generate the correct relic abundances at a mass of $\mathcal{O}(100)~{}$GeV, making it one of the most favorable and widely studied cold DM candidates. It is also worth noting that the MSSM can address the “naturalness” problem [3, 4, 5] while also accommodating a scalar boson consistent with the $h$ measurements, thus remaining one of the most attractive frameworks to pursue new physics beyond the scope of the Standard Model.

In this paper, we focus on the case of a light neutralino DM with mass $m_{\tilde{\chi}_{1}^{0}}\leq M_{h}/2$ such that the Higgs boson can decay invisibly through $h\to\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$. Since charged Higgsinos and Winos are constrained to have masses above $\gtrsim 100~{}$GeV by LEP searches [6, 7], the $\tilde{\chi}_{1}^{0}$ with mass below $\lesssim M_{h}/2$ is left with the sole possibility of a dominant Bino admixture. The observed DM abundance in our Universe is $\Omega_{DM}^{obs}h^{2}=0.120\pm 0.001$ as measured by the PLANCK collaboration [8]. However, the annihilation cross-section for Bino-dominated neutralinos is too small and leads to overclosure of the Universe ($\Omega_{LSP}h^{2}>\Omega_{DM}^{obs}h^{2}$), except under special circumstances, such as annihilation through an $s$-channel resonance with a mass of $\sim 2m_{\tilde{\chi}_{1}^{0}}$ and sfermion exchange. Our primary focus in this work will be the first scenario, however, we will also study the implications of light sfermions on our results. When the LSP contributes to the invisible decay of the SM-like Higgs boson, the available resonances through which it can annihilate are the $Z$ and $h$ bosons. The regions of parameter space where the LSP mass lies within a window of 3-5 GeV around half the $Z$ boson mass (45 GeV) or half the $h$ boson mass (62.5 GeV), are referred to as the $Z$ funnel region or the $h$ funnel region, respectively.

Several studies have explored the prospect of a light neutralino DM, in the constrained MSSM (cMSSM) and the phenomenological MSSM (pMSSM) considering the various experimental constraints at the time [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. The ATLAS and CMS Collaborations at the LHC have made available new results from searches of heavy Higgs bosons [32], direct searches of charginos and neutralinos [33, 34, 35, 36], as well as the invisible decay of the SM Higgs boson [37] at Run2. The XENON-1T, XENON-nT, PICO-60, PandaX-4T, and LUX-ZEPLIN (LZ) collaborations have also published limits on the DM direct detection (DD) cross-sections $-$ both spin-dependent (SD) and spin-independent (SI) [38, 39, 40, 41, 42, 43, 44]. Among these, the results from the LZ collaboration provide the most stringent bounds on the SI DD cross-sections [42] for DM masses in the 10 GeV - 1 TeV range. With these new and improved results, revisiting the MSSM parameter space containing light neutralino DM, which can also contribute to the invisible decay of the Higgs boson, becomes crucial.

In Ref. [31], we investigated the current status of the pMSSM parameter space with ten free parameters, that accommodates a light neutralino DM satisfying the upper bound on the relic density. We considered both positive and negative values for the Higgsino mass parameter $\mu$. The implications from the latest direct detection experiments in both spin-independent and spin-dependent interactions were studied in conjunction with the current bounds from Higgs invisible measurements, heavy Higgs searches as well as electroweakino searches at the LHC. We found that the latest direct detection limit from the LZ collaboration puts the $\mu>0$ scenario under severe tension. For the $\mu<0$ scenario, the LZ bound along with the constraints from electroweakino searches at colliders excludes most of the parameter space, except for very light Higgsinos, having masses $125-145$  GeV and $145-160$ GeV in the $Z$ and $h$ funnel regions, respectively.

In this paper, we present comprehensive and exhaustive arguments for the results found in Ref. [31]. To further concretise our findings, we study the interplay of direct detection and collider limits in a simplified scenario consisting of the SM extended by a spin-$1/2$ Majorana fermion DM in analogy with the neutralino $\tilde{\chi}_{1}^{0}$ DM in the pMSSM framework. We then investigate the light Higgsinos surviving in the negative $\mu$ scenario and the current analyses sensitive to them. We choose benchmarks from the different allowed regions of our parameter space and perform dedicated analysis of the $3l+\rm E{\!\!\!/}_{T}$ channel for light Higgsinos in the $\mu<0$ scenario using a machine learning algorithm, XGBOOST to explore its potential sensitivity. We further extend our previous work to study the impact of light staus. We present benchmarks where the Higgsino can decay to staus and perform an analysis with tau leptons in the final state to explore the sensitivity for Run-3. Finally, we discuss how the status changes for a thermal neutralino in non-standard cosmological scenarios.

We organise this paper as follows: In Sec. II, we summarize the Higgs and electroweakino sectors of the MSSM, most relevant to the present analysis. The parameter space of interest, scanning technique and the scan ranges, are also discussed in the same section. The impact of constraints from LEP measurements, flavor observables, and Higgs measurements at the LHC, are examined in Sec. III. We devote Sec. IV to exploring the implications from relic density bounds. The effect of constraints from direct detection measurements on the $Z$ and $h$ funnel regions are scrutinized in Sec. IV.1 and IV.2. Sec. IV.3 examines a simplified scenario with the SM extended by a stable Majorana fermion. We discuss the effect of electroweakino constraints in Sec. V. Details and results from our collider analysis targeted on smaller Higgsino mass $\mu\lesssim 160~{}$GeV are presented in Sec. VI. In Sec. VII and Sec. VIII, we respectively discuss how our results change when we have light staus and in scenarios of non-standard cosmology. We conclude in Sec. IX.

II The ${\rm p}$MSSM framework and the parameter space

In the pMSSM framework, the lightest neutralino $\widetilde{\chi}_{1}^{0}$ is a well-motivated DM candidate, provided it is the LSP. The $\widetilde{\chi}_{1}^{0}$ eigenstate can be written in terms of the Bino ($\tilde{B}$), neutral Wino ($\tilde{W}_{3}$) and neutral Higgsinos ($\tilde{H}_{1}^{0}$ and $\tilde{H}_{2}^{0}$),

$$\tilde{\chi}_{1}^{0}=N_{11}\tilde{B}+N_{12}\tilde{W}^{3}+N_{13}\tilde{H}_{1}^{0}+N_{14}\tilde{H}_{2}^{0},$$ (1)

where, $N_{1i}$ represents the amount of Bino ($i=1$), Wino ($i=2$) and the Higgsino ($i=3,4$) admixtures. In the present work, we are interested in the region of the pMSSM parameter space where $\widetilde{\chi}_{1}^{0}$ is the LSP and ‘light’, i.e., $M_{\widetilde{\chi}_{1}^{0}}\leq M_{h}/2$, such that it is kinematically feasible for the SM-like Higgs boson $h$ to decay into a pair of $\widetilde{\chi}_{1}^{0}$’s, thus contributing to its invisible decay. $\widetilde{\chi}_{1}^{0}$ interacts with other electroweakinos via the exchange of SM $Z/W^{\pm}$ bosons and the pMSSM Higgs bosons, with the coupling strengths determined by their electroweakino composition. These various interactions of the LSP become important factors in the calculation of its relic density and direct detection cross-sections. The $\widetilde{\chi}_{1}^{0}\widetilde{\chi}_{1}^{0}Z$ coupling can be expressed as follows [45],

$$g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}Z}^{L}=-\frac{g}{2{\rm cos}\theta_{W}}\left(|N_{13}|^{2}-|N_{14}|^{2}\right),~{}~{}~{}g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}Z}^{R}=-g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}Z}^{L},$$ (2)

where $g$ is the SU(2) gauge coupling and $\theta_{W}$ is the Weinberg angle. Similarly, the coupling of $\widetilde{\chi}_{1}^{0}$ with the three neutral Higgs bosons of the MSSM Higgs sector can be expressed as [45]:

$$g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}(h/H/A)}^{L}=g\left(N_{12}-{\rm tan}\theta_{W}N_{11}\right)\left(e_{h/H/A}N_{13}+d_{h/H/A}N_{14}\right),$$

$$g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}(h/H)}^{R}=g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}(h/H)}^{L},~{}~{}~{}g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}A}^{R}=-g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}A}^{L},$$ (3)

where,

$$e_{h}=-{\rm sin}\,\alpha,~{}~{}e_{H}={\rm cos}\,\alpha,~{}~{}e_{A}=-{\rm sin}\,\beta,$$

$$d_{h}=-{\rm cos}\,\alpha,~{}~{}d_{H}=-{\rm sin}\,\alpha,~{}~{}d_{A}={\rm cos}\,\beta,$$ (4)

where, $\alpha$ is the mixing angle in the CP–even neutral Higgs sector and ${\rm tan}\beta$ is the ratio of the vevs of the two Higgs doublets.

We require the light CP-even neutral Higgs boson $h$ to be consistent with the properties of the observed SM-like Higgs boson. $h$ can decay via $h\to\widetilde{\chi}_{1}^{0}\widetilde{\chi}_{1}^{0}$, provided $M_{\widetilde{\chi}_{1}^{0}}\leq M_{h}/2$ and $|g_{\widetilde{\chi}_{1}^{0}\widetilde{\chi}_{1}^{0}h}|>0$. As previously discussed, the former condition requires a Bino-dominated $\widetilde{\chi}_{1}^{0}$ in order to evade the lower bounds on $M_{{\widetilde{\chi}_{1}}^{\pm}}$ from LEP measurements [7]. We also notice from Eqn. 3 that $g_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}h}$ becomes zero when $\widetilde{\chi}_{1}^{0}$ is a pure Gaugino, referred to Bino and Wino collectively, or a pure Higgsino. Therefore, a Gaugino-dominated $\widetilde{\chi}_{1}^{0}$ must have some Higgsino admixture to couple with the SM-like Higgs boson, which typically entails a dominant Higgsino-like next-to-lightest supersymmetric particle (NLSP). Furthermore, the lower bounds on the masses of Higgsinos are less stringent than that of the Winos. Accordingly, in the present study, we are led to the parameter space where $\widetilde{\chi}_{2}^{0},~{}\widetilde{\chi}_{3}^{0}$ and ${\widetilde{\chi}_{1}}^{\pm}$ have a predominant Higgsino composition, while $\widetilde{\chi}_{4}^{0}$ and ${\widetilde{\chi}_{2}}^{\pm}$ have a dominant Wino admixture.

Additionally, $\tilde{\chi}_{1}^{0}$ also interacts with the SM fermions and sfermions. Searches at the LHC have derived strong lower bounds on the first two generations of squarks. For example, a single non-degenerate squark is constrained above $\gtrsim 1200~{}$GeV for $m_{\widetilde{\chi}_{1}^{0}}\sim 60~{}$GeV [46]. Therefore, the effect of their interactions on the observables related to the LSP will be negligible and motivates fixing their mass parameters to high values, say 5 TeV.

The sleptons are relatively weakly constrained from collider searches and their couplings with $\tilde{\chi}_{1}^{0}$ can be expressed as [45],

$$g_{\tilde{\chi}_{1}^{0}l\tilde{l}}=\sqrt{2}g{\rm sin}\theta_{W}\left[Q_{l}\left(N_{11}{\rm cos}\theta_{W}+N_{12}{\rm sin}\theta_{W}\right)+\left(I_{l}^{3}-Q_{l}{\rm sin}^{2}\theta_{W}\right)\frac{N_{12}{\rm cos}\theta_{W}-N_{11}{\rm sin}\theta_{W}}{{\rm cos}\theta_{W}{\rm sin}\theta_{W}}\right],$$ (5)

where $Q_{l}$ and $I_{l}^{3j}$ are the charges of the lepton and third component of isospin of the lepton, respectively. Among the three generations of sleptons, the staus have the weakest limits from collider searches. Therefore, the presence of light staus can impact the parameter space of light neutralino thermal DM. We perform our scan in two parts $-$ without and with light staus.

The relevant input parameters are: the Gaugino masses, $M_{1}$ (Bino) and $M_{2}$ (Wino), the Higgsino mass parameter $\mu$, the ratio of the Higgs vacuum expectation values $\tan\beta$, the pseudoscalar mass $M_{A}$, the mass of the third generation squarks $\{M_{\tilde{Q}_{3l}}$, $M_{\tilde{t}_{R}}$, $M_{\tilde{b}_{R}}\}$, the trilinear coupling of the stop, $A_{t}$, and the gluino mass parameter, $M_{3}$. We perform a random scan over these ten input parameters in the ranges specified below:

$$30~{}{\rm GeV}<M_{1}<100~{}{\rm GeV},~{}1~{}{\rm TeV}<M_{2}<3~{}{\rm TeV},$$

$$100~{}{\rm GeV}<|\mu|<~{}2~{}{\rm TeV},$$

$$2<\tan{\beta}<50,~{}100~{}{\rm GeV}<M_{A}<5~{}{\rm TeV},$$

$$3~{}{\rm TeV}<M_{\tilde{Q}_{3L}}<20~{}{\rm TeV},~{}3~{}{\rm TeV}<M_{\tilde{t}_{R}}<20~{}{\rm TeV},$$

$$3~{}{\rm TeV}<M_{\tilde{b}_{R}}<20~{}{\rm TeV},~{}-20~{}{\rm TeV}<A_{t}<20~{}{\rm TeV},$$

$$2~{}{\rm TeV}<M_{3}<5~{}{\rm TeV}.$$ (6)

For the first part of the scan, i.e., without the light staus, we decouple the first two generations of squarks and all the three generations of sleptons from the spectrum and set the following pMSSM input parameters to a fixed value,

$$M_{\tilde{Q}_{1,2L}}=M_{\tilde{u}_{1,2R}}=M_{\tilde{d}_{1,2R}}=5\,{\rm TeV},~{}A_{u/d/c/s/b}=0,$$

$$M_{\tilde{L}_{1,2,3L}}=M_{\tilde{e}_{1,2,3R}}=2\,{\rm TeV},~{}A_{e/\mu/\tau}=0.$$ (7)

Here, $\{M_{\tilde{Q}_{1,2L}},M_{\tilde{u}_{1,2R}},M_{\tilde{d}_{1,2R}}\}$ are the first and second generation squark mass parameters, and $\{M_{\tilde{L}_{1,2,3L}},M_{\tilde{e}_{1,2,3R}}\}$ are the left and right-handed slepton mass parameters.

We perform separate scans for the positive and negative values of $\mu$ to examine the role of $sgn(\mu)$ on the results. It is also worth noting that within the parameter space of our interest, DM relic density and direct detection constraints restrict $M_{\widetilde{\chi}_{1}^{0}}$ to the $Z$ and $h$ funnel regions only [29, 30, 47]. We begin by performing a random scan over the specified parameter space. Consequently, in order to sufficiently populate the funnel regions, we perform a dedicated scan where we dynamically tune $M_{1}$ such that $m_{\tilde{\chi}_{1}^{0}}$ is within $m_{Z}/2\pm 5~{}$GeV and $M_{h}/2\pm 3~{}$GeV. Additionally, we extract the pole mass of the top quark, $M_{t}$, randomly from a gaussian distribution with a central value of 173.21 GeV and a standard deviation of 0.55 GeV [48]. Until this point, we scan over $2\times 10^{8}$ points in total.

To study the effect of light staus, we perform a second scan where we vary the stau mass parameters, in a way that the staus become the NLSPs. We can accommodate single left-handed (LH) or right-handed (RH) staus with masses around 100-150 GeV according to Ref. [49]. According to Eqn. 5, the coupling of sleptons with the LSP neutralino depends on the third component of isospin of the slepton, which is higher for RH staus as compared to the LH ones. Therefore, RH staus will have greater impact on the relic density. To avoid constraints from additional light LH staus and light sneutrino, we consider only light RH staus. This is a minimal extension to our earlier scan [31], which is achieved by lowering the parameter $M_{\tilde{e}_{3R}}$, and varying it between 85 GeV and 500 GeV, while keeping $M_{\tilde{L}_{3L}}$ and $A_{\tau}$ fixed at 2 TeV and 0, respectively.

We use FeynHiggs 2.18.1 [50, 51, 52, 53, 54, 55, 56, 57] to generate the particle spectrum of the SUSY particles and of the Higgs bosons for each set of input parameters  and branching fractions for the decay of the Higgs bosons. MicrOMEGAS 5.2.13 [58, 59, 60, 61, 62, 63] is used to compute the LEP, flavor physics, and dark matter observables, as further discussed in Sec. III, where we also describe the various relevant constraints and their impact on the scanned parameter space.

III Constraints from LEP, flavor observables and the Higgs sector

As previously discussed, we associate the lightest CP-even Higgs boson $h$ with the discovered Higgs boson at the LHC, and require that the masses and branching of the two match within experimental and theoretical uncertainties. Measurements at the LHC have put the mass of the Higgs boson at $125.38\pm 0.14~{}$GeV [64]. We require that the mass of $h$, as computed by FeynHiggs 2.18.1, must be in the range $122~{}\rm{GeV}-128~{}$GeV. Considering theoretical uncertainties stemming from the dependence on the renormalisation scheme and scale, from the assumption of zero external momentum in two-loop corrections, and also absence of higher order corrections, we allow for a conservative 3 GeV window around the experimentally measured value [65, 66, 67] . At low values of tan $\beta$, consistency with the Higgs boson mass constraint requires large $A_{t}$ and stop masses. However, large $A_{t}$ can give rise to color and charge-breaking minima (CCB) [68, 69, 70], where the scalar partners of top quarks having color and electric charges develop a non-zero vev and the corresponding minima is lower than the minima of the Higgs field. This can be evaded, given, $|X_{t}|<\sqrt{6M_{t_{1}}M_{t_{2}}}$ [69], where $X_{t}=A_{t}-\mu/\rm{tan}\,\beta$ and $M_{t_{1,2}}$ represents the stop masses. In the present study, we see that the CCB condition has no significant effect on the allowed parameter space.

We also apply limits on the invisible decay width of $Z$-boson $\Gamma_{Z}\leq 2~{}$MeV [71], chargino mass $m_{\widetilde{\chi}_{1}^{0}}\geq 103.5~{}$GeV [7], and cross-section of neutralino pair production $\sigma(\widetilde{\chi}_{2}^{0}\widetilde{\chi}_{1}^{0})\leq 0.1~{}$pb in final states with jets $+\rm E{\!\!\!/}_{T}$ [7], as obtained from LEP. We impose constraints on various flavor physics observables, such as, $Br(b\rightarrow s\gamma)=(3.32\pm 0.16)\times 10^{-4}$ [72], $Br(B_{s}\rightarrow\mu^{+}\mu^{-})=3.0^{+0.67}_{-0.63}\times 10^{-9}$ [73], and $Br(B\rightarrow\tau\nu)=1.28\pm 0.25$ [74], allowing $2\sigma$ uncertainty around the best-fit values. We have used MicrOMEGAS 5.2.13 [58, 59, 60, 61, 62, 63] to calculate both the LEP and flavor physics observables.

We also impose the Higgs signal strength constraints on the parameter space using the HiggsSignal 2.6.2 [75, 76, 77] package, while limits from the heavy Higgs boson searches at the LHC are imposed using the HiggsBounds 2.10.0 [78, 79, 80, 81, 82] package. The parameter space is also required to satisfy the most stringent upper bound on the Higgs invisible branching ratio $Br(h\to\rm{inv})\leq 11\%$, as measured by the ATLAS collaboration [1], which is stronger than the current CMS bound (15% [2]). Hereafter, constraints on the mass, signal strength (imposed through HiggsSignal 2.6.2), and the invisible branching ratio of the Higgs boson are combinedly referred to as the constraints on “Higgs properties”. We summarize the constraints in Table 1.

Sr. No.	Observable	Calculated by	Constraint
(1)	Light Higgs boson mass	FeynHiggs 2.18.1	122 GeV$<M_{h}<$128 GeV
(2)	Higgs signal strength	HiggsSignal 2.6.2	111 channels, $p$-value$>0.05$
(3)	Heavy Higgs bosons	HiggsBounds 2.10.0	Constraints from collider searches of heavy
			Higgs bosons implemented in HiggsBounds
(4)	Invisible decay of Higgs boson	FeynHiggs 2.18.1	Br$(h\rightarrow\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0})<0.11$
(5)	Invisible decay of $Z$ boson from LEP	MicrOMEGAS 5.2.13	$\Gamma(Z\rightarrow\text{invisible})<2$ MeV
(6)	Chargino mass limit from LEP	FeynHiggs 2.18.1	$m_{\chi_{1}^{\pm}}>103$ GeV
(7)	LEP limits on neutralino	MicrOMEGAS 5.2.13	$\sigma(e^{+}e^{-}\rightarrow\tilde{\chi}_{1}^{0}\tilde{\chi}_{2}^{0})\times{\rm Br}(\tilde{\chi}_{2}^{0}\rightarrow\tilde{\chi}_{1}^{0}+{\rm jets})+$
	in dijet + MET final states		$\sigma(e^{+}e^{-}\rightarrow\tilde{\chi}_{1}^{0}\tilde{\chi}_{3}^{0})\times{\rm Br}(\tilde{\chi}_{3}^{0}\rightarrow\tilde{\chi}_{1}^{0}+{\rm jets})<0.1$ pb
(8)	Flavour observables	MicrOMEGAS 5.2.13	$3.00\times 10^{-4}<{\rm Br}(b\rightarrow s\gamma)<3.64\times 10^{-4}$
			$1.66\times 10^{-9}<{\rm Br}(B_{s}\rightarrow\mu^{+}\mu^{-})<4.34\times 10^{-9}$
			$0.78<\frac{({\rm Br}(B\rightarrow\tau\nu))_{obs}}{({\rm Br}(B\rightarrow\tau\nu))_{SM}}<1.78$

Table 1: Summary of constraints from Higgs sector, LEP, and flavor observables.

Figure 1: Left: Parameter space in the $M_{A}$-${\rm tan}\beta$ plane satisfying the Higgs properties, LEP constraints, flavor constraints, and results of heavy Higgs boson searches implemented in HiggsBounds; Right: Parameter space in the $M_{A}$-${\rm tan}\beta$ plane zoomed-in to show the effect of the different flavor observables for $M_{A}<1$ TeV. The top and bottom panels are for $\mu>0$ and $\mu<0$ respectively.

We show the scanned points in the $M_{A}-\tan\beta$ plane for $\mu>0$ (top) and $\mu<0$ (bottom) and the effect of various constraints in Fig. 1. In the left panels, we show the parameter points allowed by the Higgs boson measurements and LEP constraints (yellow), flavor physics bounds (light-green), and the constraints from heavy Higgs searches included via HiggsBounds (dark-green), applied consecutively. In the right panels, we zoom in on the $M_{A}\leq 1~{}$TeV region, showing the parameter space excluded by different flavor-changing processes. The flavor physics observables combinedly exclude points with low $M_{A}~{}$($\lesssim 700~{}$GeV). The low $\tan\beta$ region is excluded by $Br(b\to s\gamma)$ while the other two flavor observables are more impacting in the high $\tan\beta$ region. Points in the high $\tan\beta$ region for $700~{}\textrm{GeV}\lesssim M_{A}\lesssim 2~{}$TeV are excluded by constraints from heavy Higgs boson searches at the LHC applied via the HiggsBounds package. Among the various heavy Higgs searches considered, results from the search for heavy Higgs bosons decaying into a pair of tau leptons at $\sqrt{s}=13~{}$TeV with $\mathcal{L}=139~{}\textrm{fb}^{-1}$ data [32] by the ATLAS collaboration resulted in the most stringent constraints on the parameter space of our interest. We further observe that the LEP constraints (Table 1 (2), (3), and (4) applied together), HiggsSignal and the upper limit on the invisible branching of $h$ to invisible final states do not show any specific trend in excluding the parameter space in the $M_{A}$-${\rm tan}\beta$ plane. All these observations are irrespective of the sign of the $\mu$ parameter.

The observed discrepancy of the muon $g-2$ measurement with the SM prediction requires an additional contribution of $(24.9\pm 4.8)\times 10^{-10}$ [83, 84] from new physics. In the MSSM, the sign of the contribution to the muon $g-2$ depends on the sign of $\mu$. Therefore, to resolve the muon $g-2$ anomaly within the MSSM, one prefers the positive sign of $\mu$. However, in our analysis the sleptons are fixed to have a high mass around 2 TeV, therefore their contribution to the muon $g-2$ is negligible. For the positive $\mu$ benchmark that we obtain later in Sec. VI, the MSSM contribution comes to be around $\sim 2\times 10^{-11}$, which is two orders of magnitude away from the required value.

IV Dark Matter constraints

As discussed previously, the lightest supersymmetric particle in the pMSSM, here $\tilde{\chi}_{1}^{0}$, is a viable DM candidate. It can have a thermal production in the early Universe, which freezes-out. In the standard cosmology, we require the relic density of the LSP $\Omega_{\widetilde{\chi}_{1}^{0}}$ to be equal to the observed DM relic density as measured by the PLANCK collaboration $\Omega^{obs}_{DM}h^{2}=0.120\pm 0.001$ [85], which assuming a $2\sigma$ interval can vary from 0.118 $\lesssim\Omega^{obs}_{DM}h^{2}\lesssim$ 0.122. Lifting up the requisite that the neutralino LSP forms 100% of the observed DM relic owing to the possibility of multicomponent DM, we can modify the relic density constraint to $\Omega_{\widetilde{\chi}_{1}^{0}}\lesssim 0.122$. We use the MicrOMEGAS 5.2.13 [58, 59, 60, 61, 62, 86, 63] package to compute the relic density of $\tilde{\chi}_{1}^{0}$.

In addition to the relic density constraint, we need to take into consideration the limits from dark matter direct detection (DD) experiments which constrain the spin-dependent DM-neutron (SDn), DM-proton (SDp) and spin-independent (SI) DM-nucleon interaction cross-sections as a function of mass of the DM. We use MicrOMEGAS 5.2.13 to compute these cross-sections and then compare them with the 90% confidence level (CL) upper limits quoted by the PICO-60 (SDp [40]), PandaX-4T (SDn [43]), and LZ (SI [42]) experiments, since these are the strongest available bounds for each category in the DM mass range of 10 GeV to 1 TeV at present. The LZ collaboration sets an upper limit on the SI cross-section of a DM particle with mass in the $Z$-funnel to be $1.06\times 10^{-47}$ cm${}^{2}$, and in the $h$-funnel to be $1.54-1.64\times 10^{-47}$ cm${}^{2}$ for $M_{h}\in 122-128$ GeV. We tabulate the DM related constraints applied on the parameter space of our interest in Table 2.

Sr. No.	Observable	Calculated by	Constraint
(9)	Relic density	MicrOMEGAS 5.2.13	$\Omega<0.122$ (PLANCK)
(10)	Limits on direct detection	MicrOMEGAS 5.2.13	Spin-dependent proton: PICO-60
	cross-sections scaled with $\xi$ (Eqn. 8)		Spin-dependent neutron: PandaX-4T
			Spin-independent: LUX-ZEPLIN (LZ)

Table 2: Summary of DM relic density and direct detection constraints on the LSP neutralino.

The DD limits from the experimental collaborations are placed assuming that a single DM candidate constitutes the entire relic. Therefore, if the neutralino DM is underabundant, i.e., $\Omega_{\widetilde{\chi}_{1}^{0}}<0.118$, then the DD limits are applied on scaled cross-sections. The scaling factor $\xi$ is taken unity when the LSP relic is within the experimental uncertainty, i.e., $0.118<\Omega_{LSP}<0.122$. For $\Omega_{LSP}<0.118$, it is scaled by the ratio of the central value of the observed relic density to the computed relic density for $\widetilde{\chi}_{1}^{0}$, as follows:

$$\xi=\frac{\Omega_{\widetilde{\chi}_{1}^{0}}}{0.120}$$ (8)

The recent upper limit on the SI cross-section $\sigma_{\textrm{SI}}$ derived by the LZ collaboration is roughly $\sim 3-4$ times stronger than the previous most stringent limits from PandaX-4T in the region of DM masses considered in this work . To demonstrate the role of the LZ result, we divide the constraints on our scanned parameter space into “Before LZ” which includes constraints from LEP, flavor, Higgs properties, heavy Higgs searches using HiggsBounds, relic density, and the DD experiments XENON-1T, PICO-60, and PandaX-4T, and “After LZ” with the constraint from the LZ experiment.

Although our scan is over a ten-dimensional parameter space, not all parameters contribute to the individual observables. The relic density is determined by the annihilation channels of the DM, which in the present scenario will dominantly proceed through the $s$-channel diagrams involving the $Z$ and the Higgs bosons as propagators. Scattering between the neutralino and the SM quarks and gluons, which forms the basis of the DD experiments, will involve the same propagators in the $t$-channel. It is worth noting that diagrams involving squark exchange do not play an important role due to strong lower limits on squark masses from searches at the LHC. Therefore, the most important parameters for the DM constraints from relic density and the DD experiments are $-$ $M_{1}$, $M_{2}$, $|\mu|$, ${\rm tan}\beta$, and $M_{A}$, since these affect the couplings of the DM with the $Z$ and Higgs bosons as shown in Eqns. 2 and 3. Large couplings are excluded by the DD experiments, whereas small values of couplings are excluded by the observed relic density constraint unless the mass of DM lies within a narrow window around half the mediator mass resulting in resonant enhancement of the DM annihilation cross-section.

In the $Z$ funnel, the coupling depends only on the Higgsino components in the LSP, which decreases as we move to higher values of $\mu$. Rewriting and simplifying Eqn. 3 for the lighter CP-even Higgs boson in the limit $M_{A}\gg M_{Z}$ (where cos $\alpha\,\sim\,$sin $\beta$), we have

$$g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}\approx-g\left(N_{12}-{\rm tan}\theta_{W}N_{11}\right)\left({\rm sin}\beta N_{14}-{\rm cos}\beta N_{13}\right)$$ (9)

Fig. 2 shows the coupling of the SM Higgs boson with the lightest neutralino as a function of the $\mu$ parameter (top panel) and tan$\beta$ (bottom panel) for both positive (left panel) and negative (right panel) values of $\mu$. For positive $\mu$, the coupling is always negative and can increase in magnitude with decreasing value of $\mu$ and tan$\beta$. For negative $\mu$, the maximal value of the coupling increases for small $|\mu|$. Moreover, it can have either sign depending on the value of tan $\beta$ $-$ the coupling is negative at large tan $\beta$ and it increases with decreasing tan$\beta$, eventually becoming positive at low tan$\beta$, around $\sim 10$. Therefore, we can have large magnitude of coupling for both high and low values of tan$\beta$.

Figure 2: Variation of the $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ coupling with $\mu$ (top) and tan$\beta$ (bottom) for $\mu>0$ (left) and $\mu<0$ (right) for points satisfying Higgs constraints, LEP and flavor constraints.

The SI DD cross-section will further receive contribution from the heavy CP-even Higgs boson in the pMSSM ($H$) present in the $t$-channel, which has the following coupling to $\tilde{\chi}_{1}^{0}$ (from Eqn. 3):

$$g_{H\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}\approx g\left(N_{12}-{\rm tan}\theta_{W}N_{11}\right)\left({\rm sin}\beta N_{13}+{\rm cos}\beta N_{14}\right) ,~{}~{}M_{A}\gg M_{Z},$$ (10)

and couples to the SM up-type and down-type quarks as follows:

$$g_{Huu}=i\frac{m_{u}}{v}\frac{{\rm sin}\alpha}{{\rm sin}\beta}\propto-{\rm cot}\beta ,~{}~{}M_{A}\gg M_{Z}$$

$$g_{Hdd}=i\frac{m_{d}}{v}\frac{{\rm cos}\alpha}{{\rm cos}\beta}\propto{\rm tan}\beta ,~{}~{}M_{A}\gg M_{Z}$$ (11)

The DM-quark scattering cross-section involves the product of the coupling $g_{H\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ with $g_{Huu}$ or $g_{Hdd}$. For $\mu>0$, the $h$ and $H$ contributions have opposite signs for up-type quarks, whereas for down-type quarks they add up, moreover, $g_{Hdd}$ is tan$\beta$ enhanced. The constructive interference between the contributions for down-type quarks is more effective than the destructive interference between the up-type quark contributions from $h$ and $H$. For $\mu<0$, when tan$\beta$ is large and $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ is negative, contributions from the two Higgs bosons destructively interfere for down-type quarks, and add up for up-type quarks. Since the cancellation is for coupling with down-type quarks, it is more effective for larger values of tan$\beta$. For small tan$\beta$, $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ turns positive, and follows the same trend as for $\mu>0$.

Having discussed the trends of the various couplings of the LSP DM, let us have a closer look at the $Z$ and $h$ funnels of both positive and negative $\mu$ and study how the recent LZ limit has affected these scenarios.

IV.1 $Z$ funnel

In the $Z$ funnel, the observed relic density bound restricts $|\mu|$ to small values since the $g_{Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ coupling depends only on the Higgsino components of $\tilde{\chi}_{1}^{0}$. Fig. 3 shows the fraction of DM satisfied by the LSP, $\xi$, as a function of the NLSP neutralino mass ($M_{\tilde{\chi}_{2}^{0}}$) for both $\mu>0$ and $\mu<0$ with the “Before LZ” set of cuts. In both cases, beyond $M_{\tilde{\chi}_{2}^{0}}\sim 450$ GeV, DM becomes overabundant assuming the standard cosmological model.

For interpreting the constraints from the LZ experiment, we need to consider the $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ coupling and the heavy Higgs contribution as well, which we discussed earlier. For $\mu>0$, in the low tan$\beta$ region, $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ increases and in the high tan$\beta$ region, the $g_{Hdd}$ coupling gets enhanced, both increasing the SI DD cross-section. The left panel of Fig. 4 shows the variation of $\sigma_{SI}\times\xi$ with $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ for $\mu>0$, with tan$\beta$ in the colorbar. It shows that even for smaller magnitude of couplings, we get values of $\sigma_{SI}\times\xi$ as comparable with those at large $|g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}|$, and this usually happens for a large value of tan$\beta$. The DD cross-sections are scaled with $\xi$, which is determined by the $g_{Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ coupling. Since this coupling does not affect the SI cross-section, $\xi$ can attain very small values, and can help scale down large SI DD cross-sections. However, we have seen from Fig. 3 that $\xi$ attains small values only for small NLSP masses, where $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ is also large. This puts the $Z$ funnel of $\mu>0$ region under severe tension in the pMSSM, where we do not find any region of parameter space satisfying the relic density constraint and the LZ DD limit of $\sigma_{SI}\times\xi\lesssim 0.106\times 10^{-10}$ pb, simultaneously.

Figure 3: Fraction of DM satisfied by $\tilde{\chi}_{1}^{0}$ ($\xi=\frac{\Omega_{LSP}}{0.120}$) as a function of $M_{\tilde{\chi}_{2}^{0}}$ for both $\mu>0$ (left) and $\mu<0$ (right) with the “Before LZ” set of cuts.

Figure 4: Variation of the scaled SI DM-nucleon cross-section ($\sigma_{SI}\times\xi$) with $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ for both $\mu>0$ (left) and $\mu<0$ (right) with the “Before LZ” set of cuts.

The right panel of Fig. 4 shows a plot similar to the one in the left panel but for $\mu<0$. For $\mu<0$, we have seen that $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ can attain very small values as it crosses zero coupling. Since we are in the $Z$ funnel, the relic density is determined by $g_{Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ and smaller values of $|g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}|$ are allowed. These lead to very small $\sigma_{SI}\times\xi$, well below the present LZ limit. Moreover, for large tan$\beta$, we have negative values of $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ which leads to tan$\beta$ enhanced cancellations from the heavy Higgs contribution. Therefore, the negative couplings have relatively smaller values of $\sigma_{SI}\times\xi$, as compared to the positive couplings, even when the magnitude of the coupling is the same. The region of parameter space with negative couplings, where there is interference between the $H$ and $h$ contributions, all satisfy the LZ limit. The future SDn DD experiments, which constrain the $g_{Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ coupling, play a much more crucial role in probing the parameter space in the $Z$ funnel of $\mu<0$, as we will later see in Fig. 11.

IV.2 $h$ funnel

As we move to the Higgs funnel, the LZ limit becomes slightly weaker than in the $Z$ funnel. Moreover the relic density bound can be satisfied for Higgsinos heavier than $\sim 450$ GeV, unlike the $Z$ funnel. This is because at high $\mu$, the coupling $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ can be large provided tan $\beta$ is small (see Fig. 2). The left panel of Fig. 5 shows the parameter space points in the $h$ funnel of $\mu>0$ in the $\mu$-tan$\beta$ plane with the colorbar showing the value of $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$. The “Before LZ” set of constraints, especially the bound on relic density restricts large $\mu$ to have only small tan$\beta$ values. The points shown in green in the left panel of Fig. 5 are allowed by the LZ limit. These points predominantly occupy regions characterised by smaller $g_{h\widetilde{\chi}_{1}^{0}\widetilde{\chi}_{1}^{0}}$, $\mu\gtrsim 400$ GeV and $\tan\beta\lesssim 35$. Note that even though the points at higher values of tan$\beta$ have very small couplings, they do not survive the LZ constraint due to the added tan$\beta$ enhanced heavy Higgs contribution.

Figure 5: Left: The $\mu>0$ $h$ funnel parameter space in the $\mu$-tan$\beta$ plane before (yellow) and after (green) applying the LZ limit [42] with $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ in the colorbar for “Before LZ”. Right: The full $\mu>0$ parameter space in the $M_{A}$-tan$\beta$ plane that satisfy the relic density constraint (grey) and that satisfy all the constraints including LZ (green) with $M_{\tilde{\chi}_{2}^{0}}$ in the colorbar.

To further demonstrate the significant impact of the LZ result on the positive $\mu$ scenario, the right panel of Fig. 5 shows the $\mu>0$ parameter space in the $M_{A}$-tan$\beta$ plane satisfying the relic density constraint in grey. The parameter space obtained after imposing both the relic density and LZ bounds survives all the other set of cuts applied till this point, i.e., the constraints from Higgs properties, LEP, flavor and DD bounds on SDn and SDp cross-sections of the DM. These points are shown in shades of green. The colorbar in the right panel of Fig. 5 shows that the minimum allowed $M_{\tilde{\chi}_{2}^{0}}$ value is around 350 GeV. The SI DD constrain large values of $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$, hence rules out lighter Higgsinos (see Fig. 2). Moreover, Higgsinos below 500 GeV are found at large values of tan $\beta$ where the coupling $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ is reduced and large $M_{A}$ as can be seen in Fig. 5. Higher values of tan$\beta$ are allowed as we go to large $M_{A}$ where the heavy Higgs contribution decreases.

Figure 6: Parameter space points for the $\mu>0$ scenario in the $M_{A}-\tan\beta$ plane obtained from a dedicated scan for large $M_{A}$ ranging from 10-50 TeV. The points satisfying the relic density constraints are shown in grey, while the green colored points are allowed by the other current constraints with $M_{\widetilde{\chi}_{2}^{0}}$ represented in the colorbar.

We have scanned $M_{A}$ up to 6 TeV, and a natural question which arises is whether even larger $M_{A}$ can decouple the effect of $H$ and extend the allowed parameter space to include larger values of tan $\beta$. This motivates us to perform a dedicated scan in the 10-50 TeV $M_{A}$ region. The result in Fig. 6 shows that we do reach higher tan$\beta$ values, which are allowed by both the observed relic density and LZ bound, and survive all the other constraints described in the previous section. In Fig. 6, the colorbar shows that it is really difficult to get Higgsinos having masses less than $\sim 250$ GeV, even when we go to $M_{A}$ as large as 50 TeV. In the next section, we study the mass ranges of Higgsinos allowed by the present electroweakino searches at the LHC.

Figure 7: Left: Parameter space points in the $h$ funnel of $\mu<0$ in the $\mu$-tan$\beta$ plane when all the constraints except the DD constraints are applied (in yellow), and in addition when the strongest available SDn DD constraint from the PandaX-4T experiment is applied (in green). Right: $\xi$ values of the points passing the PandaX-4T SDn bound with tan$\beta>10$ as a function of tan$\beta$ with the colorbar showing $\mu$.

In the left panel of Fig. 7, we show the parameter points in the $h$ funnel region of the $\mu<0$ scenario in the $\mu-\tan\beta$ plane. The points are depicted under two conditions: firstly, all the constraints except the DD constraints are applied (yellow), and secondly, the most stringent SDn DD constraint from PandaX-4T is additionally applied (green). Note that the LZ limit is not applied on this parameter space yet. We observe that $\tan\beta$ is restricted to smaller values for large $|\mu|$, similar to the $\mu>0$ scenario. The SDn DD cross-section depends on the coupling of the LSP with the $Z$ boson, which grows with decreasing magnitude of $\mu$, and therefore, the PandaX-4T SDn bound excludes low values of $|\mu|$. Regions of large tan$\beta$ evade this bound, although $g_{Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ does not have any significant dependence on tan$\beta$. The tan$\beta$ dependence comes from the scaling factor $\xi$, since the tan$\beta$ dependent $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ coupling determines the relic density in the Higgs funnel.

We show values of $\xi$ for the points passing the PandaX-4T SDn bound with tan$\beta>10$ as a function of tan$\beta$ in the right panel of Fig. 7 with the colorbar showing $\mu$. The purpose of this is to show how points with low $\mu$ and high tan $\beta$ survive the PandaX-4T bound on SDn cross-sections. It was shown in Fig. 2 that for $\mu<0$, large magnitudes of $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ are possible for large tan$\beta$ and small $\mu$. This makes $\xi$ small and reduces the scaled SDn direct detection cross-section, even for low values of $|\mu|$, where otherwise the SDn cross-sections are high due to large $g_{Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ coupling. These high coupling values at high tan$\beta$ also survive the SI DD bounds due to destructive contribution from the heavy Higgs boson. This region is particularly interesting since it involves light Higgsinos, which provide important benchmarks for Run-3 of LHC. For larger $|\mu|$, the $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ coupling is smaller and larger tan$\beta$ is required to increase $\xi$ in order to satisfy the bound on SDn DD cross-section.

Figure 8: Left: The $h$ funnel of the $\mu<0$ in the $M_{\tilde{\chi}_{1}^{0}}$-$\sigma_{SI}\times\xi$ plane before and after the LZ limit. Right: The $h$ funnel of the $\mu<0$ in the $M_{A}$-tan$\beta$ plane satisfying various cuts, and allowed points with $M_{\tilde{\chi}_{2}^{0}}$ in the colorbar.

To understand the effect of the LZ result, the left panel of Fig. 8 shows the parameter space in the $M_{\tilde{\chi}_{2}^{0}}$-$\sigma_{SI}\times\xi$ plane for the $h$ funnel for $\mu<0$ before and after the LZ result. This includes the SDn constraint from PandaX-4T, which we discussed in Fig 7. The colorbar in the plot indicates the corresponding $\tan\beta$ values. LZ excludes a significant part of the parameter space, especially for high $\mu$ values corresponding to large $M_{\tilde{\chi}_{2}^{0}}$, where only a narrow strip of allowed region with very small tan $\beta$ remains. With only a 20% improvement in the LZ limit in the future, this region can be fully probed. The very low SI DD cross-sections found for Higgsinos in the mass range 200-400 GeV is due to the destructive interference between the $h$ and $H$ contributions which is more important at large tan $\beta$. The right panel of the same figure shows the allowed parameter space in the $M_{A}$-tan$\beta$ plane, with $M_{\tilde{\chi}_{2}^{0}}$ in the colorbar. The gap in the allowed region around tan $\beta$ values of 6-18 is due to the recent bound on SDn cross-section by the PandaX-4T collaboration, as we observed previously in Fig. 7.

The upshot for $\mu>0$ scenario is that after the DM constraints, the $Z$-funnel is excluded by the LZ result, however, a region of the $h$-funnel survives these constraints. In this region, the lightest allowed Higgsino is around 350 GeV. For the $\mu<0$ scenario, we found allowed points both in the $Z$ and $h$ funnel, and Higgsinos as light as 125 GeV survive DM constraints. The major factor that creates a difference between the results of positive and negative $\mu$ is the effect of the heavy Higgs bosons in the DD cross-sections.

IV.3 Comparison with a simplified model

The previous section has shown the importance of the heavy Higgs boson contribution to the SI DD cross-section and its dependence on tan $\beta$ in pMSSM. As we move towards heavier Higgsinos, the relic density upper bound can only be satisfied at very low values of tan $\beta$ for both $\mu>0$ and $\mu<0$. In this region, the effect of $H$ is decoupled, and the lighter Higgs boson plays the dominant role. It is worth studying whether our results for high $\mu$ in pMSSM with the DM constraints, especially the interplay of relic density and the recent LZ upper limit on the SI DD cross-section, generalise to any BSM theory consisting of a Majorana fermion coupling with only the light Higgs boson, which resembles the discovered Higgs boson at LHC.

We consider a simplified model where the SM is extended by a single Majorana fermion, $\chi$, which is the DM candidate. It has coupling with the SM Higgs boson ($g_{\chi h}$) and has the following Lagrangian:

$$\mathcal{L}=\mathcal{L}_{SM}+g_{\chi h}\bar{\chi^{C}}\chi h.$$ (12)

We scan over the mass of the DM and its coupling with $h$, i.e., in the $M_{\chi}-g_{\chi h}$ plane. The Higgs boson mass is fixed at 125 GeV, and the total width of $h$ is calculated from the model. Since this minimal model has only two input parameters, $M_{\chi}$ and $g_{\chi h}$, the total width of the Higgs boson has negligible variation and attains a value around $4.02$ MeV when the model is implemented in MicrOMEGAS 5.2.13.

Figure 9: Result from the scan of a simplified model where we extend the SM by a Majorana fermion coupled to the $h$ boson in the $g_{\chi h}$-$M_{h}-2\times M_{\chi}$ plane. The blue points satisfy the relic density upper bound with $\xi$ shown in the colorbar, and the green points also satisfy the LZ limit.

The DM relic density depends on the couplings involved in the annihilation process, the difference between twice the DM mass and the resonance in the funnel region, and the width of the mediator particle. The coupling of the Higgs boson to the SM particles and the width of the Higgs boson do not vary much in this set-up, the only variables are therefore $M_{h}-2M_{\chi}$ and $g_{\chi h}$. The contribution of DM to the total width of $h$ is negligible. Moreover, the scattering cross-section of $\tilde{\chi}$ on nucleons relevant for DD depends on the same coupling. We present our results of the scan in the plane of $M_{h}-2M_{\chi}$ versus $g_{\chi h}$. The left panel of Fig. 9 shows the region of parameter space surviving the relic density constraint in blue with the colorbar showing the fraction of DM constituted by this additional fermion, and the allowed region of parameter space after the LZ DD bound in green. We find that the coupling allowed by both relic density and DD is minimum when $M_{h}-2M_{\chi}\sim 2$. Furthermore, $g_{h\chi}$ cannot be much smaller than $\sim~{}0.005$ in order to satisfy these two constraints. The region of parameter space that satisfies current constraint can be probed by improving the DD limit by $\sim$11%, which should be achieved in a few days of running of LZ.

The result for the simplified scenario makes it clear that the relic density and the recent LZ limit strongly constrain the coupling and mass plane of the DM when it connects to the SM particles through a Higgs portal. Therefore, even in a BSM theory like pMSSM, where we have a wide range of input parameters, the parameter space is equally constrained. The right panel of Fig. 9 shows an analogous plot for the $\mu>0$ scenario of pMSSM for heavier Higgsinos ($\mu\gtrsim 850$ GeV). For $\mu<0$ scenario of pMSSM where $|\mu|$ is kept higher than $\sim 850$ GeV, we get a similar result. The key difference that we observe in the pMSSM for large $|\mu|$ is that few points with larger coupling ($g_{\chi h}\sim 0.01-0.025$) are still allowed by the LZ limit, unlike the simplified scenario. This is due to a larger variation of $\Gamma_{h}$ in the pMSSM than in the simplified model as a result of more input parameters in the former. These points marginally pass the LZ limit in pMSSM due to a $\sim 5$% change in the width which lowers the relic density by the same amount, thereby reducing $\xi$. However, the constraints on $\Gamma_{h}$ from Higgs signal strength measurements do not allow it to vary from the SM prediction by a large fraction, which limits its effect on our results.

V Electroweakino constraints and the allowed parameter space

The electroweakinos can be directly produced at the colliders, where the NLSP Higgsinos ($\tilde{\chi}_{1}^{\pm}$/$\tilde{\chi}_{2}^{0}$/$\tilde{\chi}_{3}^{0}$) can decay to final states involving the LSP neutralino ($\tilde{\chi}_{1}^{0}$) along with $W$, $Z$ or $h$ bosons, which can have both leptonic and hadronic decays. We use the SModelS 2.2.1 [87, 88, 89, 90, 91, 92, 93, 94] package to implement the electroweakino search constraints on our scanned parameter space. Recently many analyses have updated their results with the full Run-2 data and this version of SModelS includes results from the recent search for electroweakinos in the leptonic final states at CMS [33] and ATLAS [34] and in the hadronic final states at ATLAS [35]. The constraints from the recent searches play a significant role in excluding a large range of $m_{\tilde{\chi}_{1}^{\pm}}$, $m_{\tilde{\chi}_{2}^{0}}$ and $m_{\tilde{\chi}_{3}^{0}}$, extending the sensitivity to higher masses, especially with the ATLAS analysis of the hadronic final states.

Figure 10: Left: Scaled SI DM-nucleon cross-section ($\sigma_{SI}\times\xi$) for $\mu>0$ as a function of the mass of the LSP neutralino DM in the region of parameter space satisfying LEP, flavor, Higgs constraints, relic density, DM DD constraints from the XENON-1T, PICO-60, PandaX-4T, and LZ experiments, as well as the electroweakino limits implemented in SModelS; Right: Parameter space satisfying all the constraints defined as “After LZ” and the region surviving the electroweakino searches in the $M_{\tilde{\chi}_{1}^{0}}$-$M_{\tilde{\chi}_{2}^{0}}$ plane for $\mu>0$, where the colorbar shows tan $\beta$.

Figure 11: Top: Scaled SI DM-nucleon cross-section ($\sigma_{SI}\times\xi$, left) and scaled SD DM-neutron cross-section ($\sigma_{SDn}\times\xi$, right) for $\mu<0$ as a function of the mass of the LSP neutralino DM in the region of parameter space satisfying LEP, flavor, Higgs constraints, relic density, and DM DD constraints from the XENON-1T, PICO-60, PandaX-4T, and LZ experiments (shown in yellow and marked “After LZ”) and after applying the electroweakino limits implemented in SModelS; Bottom: Parameter space satisfying all the constraints defined as “After LZ” and the regions surviving the electroweakino searches in the $M_{\tilde{\chi}_{1}^{0}}$-$M_{\tilde{\chi}_{2}^{0}}$ plane for $\mu<0$, where the colorbar shows tan $\beta$.

We apply the limits from electroweakino searches on the parameter space surviving all the constraints discussed previously in Sections 3 and 4, which are defined as “After LZ”. We identify regions of the parameter space surviving all the constraints: $A)$ high values of $|\mu|$ ($\gtrsim 800$ GeV) with low tan $\beta$ ($\lesssim 10$) in the $h$ funnel of both positive and negative $\mu$, and B) low values of $|\mu|$ ($\lesssim 200$ GeV) in both the $Z$ and $h$ funnels of negative $\mu$. We perform dedicated scans over these regions again with an additional sample of size $\sim 10^{8}$, which makes the total number of points scanned for our analysis to be $\sim 3\times 10^{8}$.

Let us first discuss the results in the positive $\mu$ scenario. Fig. 10 (left panel) shows the scaled SI DM-nucleon cross-section ($\sigma_{SI}\times\xi$) with the mass of the LSP neutralino DM for the allowed parameter space for $\mu>0$. In the right panel of Fig. 10, we show these points in the $M_{\tilde{\chi}_{1}^{0}}$-$M_{\tilde{\chi}_{2}^{0}}$ plane with $\tan\beta$ in the colorbar. In Section IV, we have discussed that the LZ limit excludes the $Z$ funnel region for $\mu>0$. As a result, both the panels of Fig. 10 have allowed points only in the $h$ funnel region.

The right panel of Fig. 10 reveals that on applying the electroweakino constraints implemented in SModelS, $M_{\tilde{\chi}_{2}^{0}}\lesssim 850$ GeV are excluded by various collider searches. These heavy Higgsinos populate very small values of tan $\beta$, which results from the DM constraints as we have discussed earlier in Section IV.2. The region of parameter space surviving the constraints from electroweakino searches can be probed by improving the LZ limit by 20% which can be achieved with just a few more days of running, as can be seen from the left panel of Fig. 10. For the simplified scenario, note that only 11% improvement in the LZ bound is required to probe the allowed parameter region, which is smaller than that required for the pMSSM parameter space. The difference in the two scenarios is again due to the variation in the Higgs boson total decay width. The coupling $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ for heavier Higgsinos becomes constant, and since the surviving regions are restricted to low tan$\beta$ values, there is no significant contribution from the heavy Higgs. This parameter space cannot have any smaller values of $\sigma_{SI}\times\xi$, as the relic density constraint does not allow for lower values of the coupling.

For the negative $\mu$ scenario, the top panels of Fig. 11 show the scaled SI DM-nucleon cross-section ($\sigma_{SI}\times\xi$, left) and scaled SD DM-neutron cross-section ($\sigma_{SDn}\times\xi$, right) as a function of the mass of the LSP neutralino DM. We observe that we have regions of parameter space which survive all the constraints in both the $Z$ and the $h$ funnels. The top left panel of Fig. 11 shows that the allowed region in the $h$ funnel is well within the reach of the next few days of LZ data which can improve the limit on SI DD cross-section by 80%, and from the top right panel, we infer that the allowed parameter space in the $Z$ funnel can be probed by the SDn result projected by the XENON-nT collaboration.

The bottom panel of Fig. 11 shows the parameter space in the $M_{\tilde{\chi}_{1}^{0}}$-$M_{\tilde{\chi}_{2}^{0}}$ plane for $\mu<0$, with the colorbar representing tan $\beta$, for the points allowed by all the constraints. The electroweakino searches restrict the allowed parameter space to either heavy Higgsinos, having masses $\gtrsim 850$ GeV, or to a narrow region of parameter space with light Higgsinos with masses in the range of 125-145 GeV in the $Z$-funnel and 145-160 GeV in the $h$-funnel, many of which have very small $R$-values . We further investigate the allowed region of such light Higgsinos in the following section. Fig. 10 and Fig. 11 are similar to Figs. 1 and 2 of Ref. [31], they are shown here for completeness. The right panel of Fig. 10 and the bottom panel of Fig. 11 also show the tan $\beta$ range of the allowed regions in the colorbar.

The future lepton colliders like ILC and CEPC will be crucial for precision measurements of Higgs boson. The projected upper limit on the invisible branching of the Higgs boson is 0.4% at ILC [95] and 0.3% at CEPC [96]. Although these can probe a significant part of the allowed parameter space in the $\mu<0$ case, as shown in Fig. 21 in Appendix A, we still have regions with Br($h\rightarrow$invisible)$<0.003$ in both the $Z$ and $h$-funnels. In the $\mu<0$ case, the partial decay width of the $Z$ boson to $\tilde{\chi}_{1}^{0}$ ($\Gamma_{\rm inv}^{\rm new}$) is always less than 0.1 MeV for the allowed parameter region that we obtain. Therefore, we do not expect the Giga-$Z$ option of ILC, which is expected to have a modest improvement over LEP [97], to be sensitive to this region.

The relic density can have theoretical uncertainties in its calculation [98, 99, 100]. If we overestimate the relic density by few percent as compared to its actual value due to the theoretical uncertainties, the scaling factor, $\xi$, reduces, resulting in a lower SI DD cross-section. Fig. 12 shows the scaled SI DD cross-section ($\sigma_{SI}\times\xi$) with the DM mass of the parameter space for $\mu>0$ surviving the “After LZ” set of constraints assuming that the relic density is overestimated by 20%, with the colorbar showing the mass of $\tilde{\chi}_{2}^{0}$. We observe that a small allowed region has opened up in the $Z$ funnel which survives the present LZ result. However, the Higgsinos have masses above $\sim 200$ GeV, and get excluded by the collider bounds from electroweakino searches, as we have seen from the bottom panel of Fig. 11. In the $h$ funnel, the allowed region extends further down to $\sigma_{SI}\times\xi\sim 0.112\times 10^{-10}$ pb as compared to our result with no theoretical uncertainty on the relic density, where the allowed region in the $h$ funnel had the lowest scaled SI DD cross-section of around 0.137$\times 10^{-10}$ pb. Still, this region is well within the projected limit from the full 1000 days of the LZ experiment.

Figure 12: Scaled SI DM-nucleon cross-section ($\sigma_{SI}\times\xi$) for $\mu>0$, with 20% overestimation of relic density due to theoretical uncertainty, as a function of the mass of the LSP neutralino DM in the region of parameter space satisfying LEP, flavor, Higgs constraints, relic density, and DM DD constraints from the XENON-1T, PICO-60, PandaX-4T, and LZ experiments, defined as “After LZ”.

Let us now focus on the regions of light Higgsinos allowed by the present electroweakino searches.

VI Collider analyses for probing the light Higgsinos

In this section, we study the surviving regions of parameter space with lighter Higgsinos in the standard cosmological scenario. In Fig. 13, we present the $R$-value as given by the SModelS package for the most sensitive analyses on the allowed parameter space in the $Z$ (left) and $h$ funnel (right) of the negative $\mu$ scenario. This includes ATLAS analyses for the final states 2 leptons ($e$,$\mu$) + MET and 3 leptons ($e$,$\mu$,$\tau$) + MET with 20.3 fb${}^{-1}$ of data, along with final states having jets + MET and 3 leptons + MET with an integrated luminosity of 139 fb${}^{-1}$. The relevant CMS analyses include searches for electroweakinos with decays to leptons, $W$, $Z$, and Higgs bosons with a luminosity of 19.5 fb${}^{-1}$ and multilepton final states with 35.9 fb${}^{-1}$. We observe that many of these points have very low $R$-values which indicate low sensitivity of the collider searches. The sensitivity drops around the region where the mass difference between the $\tilde{\chi}_{1}^{\pm}$/$\tilde{\chi}_{2}^{0}$ and the $\tilde{\chi}_{1}^{0}$ is close to the mass of the $Z$ boson. Since in this region the decay products of the NLSP electroweakinos are produced at rest, the $\tilde{\chi}_{1}^{0}$ does not carry significant momentum, thereby, decreasing the effectiveness of the $p_{T}^{miss}$ variable to differentiate the signal from SM backgrounds. Therefore, the experimental results suffer from low sensitivity near the $Z$ boson mass threshold.

Representative benchmarks from each of the allowed regions of the parameter space are presented in Table 3. These benchmarks have very small uncertainty in the Higgs boson mass as estimated by FeynHiggs ($\Delta_{M_{h}}^{FH}\lesssim\mathcal{O}(1)~{}$GeV), and have SModelS $R$-values below 0.5 . They are also allowed when tested with CheckMATE 2 [101], another package that implements the constraints from electroweakino searches. We find that the Tevatron searches for light charginos [102] are also not sensitive to these benchmarks.

Figure 13: $R$-values from SModelS-2.2.1 of the most sensitive analyses versus the $M_{\tilde{\chi}_{2}^{0}}$ for the light Higgsino scenarios allowed for $\mu<0$ in the $Z$ funnel (left) and $h$ funnel (right). The different colors depict the most sensitive analysis for each point.

Benchmarks (mass parameters in GeV)				$M_{h}[\Delta_{M_{h}}^{FH}]~{}$[GeV]	$\sigma_{SI}\times\xi\times 10^{-10}$ [pb]
$\mu>0$	$h$-funnel	BP1	$M_{t}=173.21$, $M_{1}=62.5$, $M_{2}=2000$, $\mu=1000$, tan$\beta=5$, $M_{A}=3000$,	125.38 [$\pm$0.97]	0.151
			$M_{\tilde{Q}_{3L}}=10000$, $M_{\tilde{t}_{R}}=10000$, $M_{\tilde{b}_{R}}=10000$, $A_{t}=10000$, $M_{3}=3000$
$\mu<0$	$Z$-funnel	BP2	$M_{t}=173.21$, $M_{1}=44$, $M_{2}=2000$, $\mu=-124$, tan$\beta=5$, $M_{A}=3000$,	125.88 [$\pm$0.96]	$7.46\times 10^{-4}$
			$M_{\tilde{Q}_{3L}}=10000$, $M_{\tilde{t}_{R}}=10000$, $M_{\tilde{b}_{R}}=10000$, $A_{t}=10000$, $M_{3}=3000$
	$h$-funnel	BP3	$M_{t}=173.21$, $M_{1}=68$, $M_{2}=2000$, $\mu=-150$, tan$\beta=50$, $M_{A}=3000$,	125.67 [$\pm$0.63]	0.143
			$M_{\tilde{Q}_{3L}}=5000$, $M_{\tilde{t}_{R}}=5000$, $M_{\tilde{b}_{R}}=5000$, $A_{t}=-5000$, $M_{3}=3000$
		BP4	$M_{t}=173.21$, $M_{1}=$, $M_{2}=2000$, $\mu=-1000$, tan$\beta=4.5$, $M_{A}=3000$,	125.15 [$\pm$0.99]	0.150
			$M_{\tilde{Q}_{3L}}=10000$, $M_{\tilde{t}_{R}}=10000$, $M_{\tilde{b}_{R}}=10000$, $A_{t}=10000$, $M_{3}=3000$

Table 3: Parameters corresponding to four benchmark points satisfying all the present constraints from the $\mu>0$ and $\mu<0$ scenarios along with their scaled SI DD cross-sections. The mass of the Higgs boson $M_{h}$ and the uncertainty in $M_{h}$ computed by FeynHiggs ($\Delta_{M_{h}}^{FH}$) are also shown.

For probing the benchmarks with heavy Higgsinos, BP1 and BP4, the hadronic decay channels of the $W$ and $Z$ bosons are more sensitive than the leptonic ones. The current ATLAS result for the hadronic final state excludes Higgsinos below 850 GeV (right panel of Fig. 10). Assuming that the upper limit on the cross-section improves by a factor of $\sqrt{\mathcal{L}}$ with increasing luminosity, $\mathcal{L}$, Run-3 will be able to probe Higgsinos up to a mass of 900-925 GeV and HL-LHC will further increase the sensitivity to $\sim 1100$ GeV. Therefore, BP1 and BP4 can be interesting benchmarks to be probed in the HL-LHC runs.

Our benchmarks for light Higgsinos survive the electroweakino constraints as implemented in recasting frameworks like SModelS and CheckMATE. However, it is not guaranteed that all the recent relevant public results from the experimental collaborations have been added in the database of these packages. This motivates a detailed analysis for such light Higgsinos, which we present in the subsequent section.

To estimate the prospects for probing the region with light charginos and neutralinos at the LHC, we perform an analysis of the low mass Higgsino-like electroweakinos in the leptonic $3l+\rm E{\!\!\!/}_{T}$ final state at $\sqrt{s}=14$ TeV using the XGBOOST [103] framework. We study the process $pp\rightarrow\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{2}^{0}/\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{3}^{0},~{}\tilde{\chi}_{1}^{\pm}\rightarrow ff^{\prime}\tilde{\chi}_{1}^{0},~{}\tilde{\chi}_{2}^{0}/\tilde{\chi}_{3}^{0}\rightarrow f\bar{f}\tilde{\chi}_{1}^{0}$ with $m_{\tilde{\chi}_{1}^{\pm}}=125.1$ GeV, $m_{\tilde{\chi}_{2}^{0}}=129.9$ GeV, $m_{\tilde{\chi}_{3}^{0}}=133.5$ GeV, and $m_{\tilde{\chi}_{1}^{0}}=44.6$ GeV (benchmark 2 from Table 3) where $f,f^{\prime}$ are SM fermions. We restrict to the leptonic final state which is cleaner for a lighter benchmark, such as ours. The SM background processes studied in the analyses are summarised in Table 4 with their respective cross-sections and simulation details.

Background	Cross section [pb]	Generated using	Total generated
$lll\nu$	0.4684$\times 1.2$	MadGraph 2.7.3	$9.98\times 10^{6}$
$WZ$, leptonic, $2j$ matched	1.253$\times 1.2$	MadGraph 2.7.3	$4.97\times 10^{6}$
$ZZ$, leptonic, $2j$ matched	0.1186$\times 1.2$	MadGraph 2.7.3	$1.25\times 10^{6}$
$t\bar{t}$, leptonic	55.36$\times 1.74$	MadGraph 2.7.3	$6\times 10^{7}$
$VVV$, inclusive	0.2678$\times 1.2$	MadGraph 2.7.3	$2.5\times 10^{6}$
$Wh$, inclusive	1.504 [104]	Pythia 8.306	$5\times 10^{6}$
$Zh$, inclusive	0.883 [104]	Pythia 8.306	$5\times 10^{6}$
ggF $h\rightarrow ZZ$, leptonic	0.0137	Pythia 8.306	$5\times 10^{6}$
VBF $h\rightarrow ZZ$, leptonic	0.00115	Pythia 8.306	$5\times 10^{6}$
$t\bar{t}h$, inclusive	0.6113 [104]	Pythia 8.306	$5\times 10^{6}$
$t\bar{t}W$, leptonic	0.01387$\times$1.22	MadGraph 2.7.3	$2.5\times 10^{6}$
$t\bar{t}Z$, leptonic	0.00644$\times$1.23	MadGraph 2.7.3	$2.5\times 10^{6}$

Table 4: Details of the background simulation and cross-sections

We perform an analysis of the $3l+\rm E{\!\!\!/}_{T}$ final state where we require exactly three leptons satisfying $p_{T}>25,25,20$ GeV and $|\eta|<2.4$, and we have put a veto on $b$-jets with $p_{T}>30$ GeV and $|\eta|<2.5$. In our signal benchmark, BP2, since we do not have any on-shell $Z$-boson, we also veto events where the invariant mass of a pair of same flavor opposite sign (SFOS) leptons lie within 10 GeV of the $Z$ mass. After these preselections, we train our signal and background samples using XGBOOST with a set of the following variables:

• Transverse momenta ($p_{T}$) of the three leptons

• Transverse mass ($M_{T}$) and contransverse mass ($M_{CT}$) of each of the three leptons with the $\rm E{\!\!\!/}_{T}$

• Minimum and maximum values of $\Delta R$ between opposite sign lepton pairs along with their $\Delta\eta$ values

• Invariant mass of the opposite sign lepton pairs with minimum and maximum $\Delta R$

• Missing transverse momentum

• Number of jets in the event with the $p_{T}$ of the two leading jets

• Scalar sum of $p_{T}$ of all the jets in the event ($H_{T}$)

• Invariant mass of the three leptons

In benchmark BP3, we study the process $pp\rightarrow\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{2}^{0}/\tilde{\chi}_{1}^{\pm}\tilde{\chi}_{3}^{0},~{}\tilde{\chi}_{1}^{\pm}\rightarrow W^{\pm}\tilde{\chi}_{1}^{0},~{}\tilde{\chi}_{2}^{0}/\tilde{\chi}_{3}^{0}\rightarrow Z\tilde{\chi}_{1}^{0}$ with $m_{\tilde{\chi}_{1}^{\pm}}=125.9$ GeV, $m_{\tilde{\chi}_{2}^{0}}=155.3$ GeV, $m_{\tilde{\chi}_{3}^{0}}=157.4$ GeV, and $m_{\tilde{\chi}_{1}^{0}}=62.2$ GeV (benchmark 3 from Table 3). For this benchmark, we apply the preselections described for BP2 above. In this benchmark, we have an on-shell $Z$ boson in the final state. We, therefore, select those events where the invariant mass of a pair of same flavor opposite sign (SFOS) leptons lie within 10 GeV window of the $Z$ mass, and we define these two leptons as the SFOS pair of leptons. We use the following variables for training the XGBOOST framework:

• Transverse momenta ($p_{T}$) of the three leptons

• Transverse mass ($M_{T}$) and contransverse mass ($M_{CT}$) of the lepton, which is not part of the SFOS pair of leptons, with the $\rm E{\!\!\!/}_{T}$

• $\Delta R$ and $\Delta\eta$ between the SFOS lepton pair

• $\Delta\phi$ and $\Delta\eta$ between the SFOS lepton pair system and the unpaired lepton

• $\Delta\phi$ between the SFOS lepton pair system and $\rm E{\!\!\!/}_{T}$

• $\Delta\phi$ between the unpaired lepton and $\rm E{\!\!\!/}_{T}$

• Missing transverse momentum

• Number of jets in the event with the $p_{T}$ of the two leading jets

• Scalar sum of $p_{T}$ of all the jets in the event ($H_{T}$)

• Invariant mass of the three leptons

We train our XGBOOST model using the following hyperparameters:

‘objective’:‘multi:softprob’, ‘colsample_bytree’:0.3, ‘learning_rate’:0.1,

‘num_class’:12, ‘max_depth’:7, ‘alpha’:5, ‘eval_metric’:‘mlogloss’,

‘num_round’:1000, ‘early_stopping_rounds’:3

We divide our total sample in two parts $-$ one for training and one for validation. The background events are merged with a weight factor calculated using the fraction of the number of events expected at the LHC for a particular luminosity and the number of events generated for each background process. The weights are then normalised such that the sum of the weights of all the background processes becomes unity. For each epoch, we train on the training data and test the training on the validation sample. The model minimises its loss function unless the loss on the validation sample does not decrease in three consecutive iterations. The XGBOOST models are separately trained with 21 kinematic variables for BP2 and 18 kinematic variables for BP3. These trained models are then used to discriminate the signal benchmarks from each background class by computing the significance of observing the signal over the background events. At the $\sqrt{s}=14$ TeV LHC with 137 fb${}^{-1}$ of integrated luminosity ($\mathcal{L}$), Table 5 shows the expected number of our two signal benchmark points and background events for a threshold of 0.9 on our XGBOOST output.

Number of events for $\mathcal{L}=137$ fb${}^{-1}$		BP2	BP3
Backgrounds	$lll\nu$	205.6	$-$
	$WZ$, leptonic, $2j$ matched	$-$	46.7
	$ZZ$, leptonic, $2j$ matched	14.7	5.8
	$t\bar{t}$, leptonic	677.6	21.8
	$VVV$, inclusive	13.0	2.3
	$Wh$, inclusive	46.5	1.4
	$Zh$, inclusive	7.4	1.4
	ggF $h\rightarrow ZZ$, leptonic	2.2	0.002
	VBF $h\rightarrow ZZ$, leptonic	0.2	6.0$\times 10^{-4}$
	$t\bar{t}h$, inclusive	8.2	0.3
	$t\bar{t}W$, leptonic	9.2	0.5
	$t\bar{t}Z$, leptonic	2.5	1.0
	Total	987.1	81.2
Signal		763.4	112.1
Significance with 20% systematic uncertainty		3.1	4.5
Significance with 50% systematic uncertainty		1.3	1.98

Table 5: Number of events from individual background processes and the signal surviving a threshold of 0.9 on the XGBOOST output from two models trained on benchmarks BP2 and BP3 respectively, along with the signal significance for $\mathcal{L}=137$ fb${}^{-1}$.

We quote our results by assuming a 20% (50%) systematic uncertainty, where the signal significance is estimated using the formula in Ref. [105]. We present our results for $\sqrt{s}=14$ TeV to make it easier to translate to the case of Run-3 ($\sqrt{s}=13.6$ TeV) and HL-LHC ($\sqrt{s}=14$ TeV) as the cross-sections for direct electroweakino production are not expected to change much. We find that the result sensitively depends on the systematic uncertainty, which can have a significant impact for light electroweakinos. Our result shows that these light Higgsinos are within the reach of LHC and could be probed with upcoming analyses of the Run-2 data or at Run-3 of the LHC, provided the systematic uncertainties can be controlled.

VII Impact of light staus on the spectrum

In our previous scan, we had fixed the soft parameters related to the first and second generation squarks and all the three generations of sleptons. The former are fixed at masses around 5 TeV and the latter at masses around 2 TeV, with all the trilinear couplings associated with these squarks and sleptons set to zero. Lighter squarks and sleptons can enter various processes of the neutralino DM and affect its relic density and in turn, impact of DD experimental constraints due to the scaling factor. For squarks, as discussed previously, the strong limits from the collider searches reduce their effect on these observables. According to Ref. [46], a single non-degenerate squark has to be heavier than $\sim 1200$ GeV for $M_{\tilde{\chi}_{1}^{0}}\sim 60$ GeV. We have found that the presence of a 1200 GeV squark has very little effect on the relic density of DM (less than 2%) and negligible effect on the DD cross-section.

Among the three generations of sleptons, staus have the weakest limits. As motivated earlier in Section II, for light staus, we are more interested in studying the effect of RH light staus $-$ for this we vary the parameter $M_{\tilde{e}_{3R}}$ from 85-500 GeV. We find that the present searches of stau leptons at the LHC, which are already recasted in the SModelS package, does not constrain the scenario of RH staus as the NLSP for the DM mass range under consideration, as shown in Fig. 22 of Appendix B.

Figure 14: Fractional change in the relic density of the LSP with and without light RH staus as a function of the mass of the lightest stau for $\mu=+1$ TeV (left) and $\mu=-1$ TeV (right). These parameter points survive all the constraints described earlier and the colorbar shows the tan $\beta$ of these points.

Figure 15: Scaled SI DM-nucleon cross-section ($\sigma_{SI}\times\xi$) for $\mu>0$, with light RH staus ($M_{\tilde{e}_{3R}}=85$ GeV), as a function of the mass of the LSP neutralino DM in the region of parameter space satisfying LEP, flavor, Higgs constraints, relic density, and DM DD constraints from the XENON-1T, PICO-60, PandaX-4T, and LZ experiments, defined as “After LZ”.

Since light RH staus are still allowed by the collider constraints, we study their impact on the relic density of the Bino-like LSP that we have studied so far. We observe that for lighter Higgsinos, the couplings $g_{Z\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ and $g_{h\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}}$ are larger, therefore, we do not expect a large effect from lighter RH staus. For heavier Higgsinos, say having $|\mu|=1$ TeV, these coupling have small values due to reduced Higgsino components in the lightest neutralino, and the effect of light RH stau becomes important. To demonstrate this, Fig. 14 shows the fraction of change in the relic density of the LSP DM with and without light RH staus for varying masses of the lightest stau for both positive (left) and negative (right) $\mu$, with the Higgsino mass parameter having a value of 1 TeV in the $h$ funnel. We observe that the sign of $\mu$ does not play a significant role, and light staus of 100 GeV can reduce the relic density by 30-40%.

To maximize the effect of a light stau, we fix $M_{\tilde{e}_{3R}}=85$ GeV corresponding to $M_{\tilde{\tau}_{1}}=90$-$95$ GeV, a value above the LEP bound [106, 107, 108, 109, 110, 111]. We then redo the scan to examine the impact of adding the light stau on the parameter space of the light neutralino thermal dark matter. Fig. 15 shows the allowed parameter space with the “After LZ” set of cuts in the $\sigma_{SI}\times\xi-M_{\tilde{\chi}_{1}^{0}}$ plane for $\mu>0$ We find that the light staus reduce the relic density and therefore, brings down the scaled SI DD cross-sections, thereby allowing a small region of parameter space in the $Z$ funnel and extending the allowed region in the $h$ funnel further down, both of which are within the reach of the LZ 1000 days projection.

Our analysis is for a 100% branching fraction of Higgsinos to the $WZ$ final state. Experimental collaborations also quote their exclusion boundaries assuming 100% branching to a specific final state and a particular mass hierarchy. The prospects of the presence of other light SUSY particles might also affect the collider constraints on Higgsinos, if the latter decay into the former with significant branching fractions. One such possibility which we mention in our paper is the presence of light staus having masses between the Bino-like LSP and Higgsinos, which are still allowed by the searches at LHC [49] (also see Fig. 22 in Appendix B).

Figure 16: SModelS $R$-values in the plane of mass of the Higgsino-like neutralino and the lightest stau, which is RH, for a range of masses where the Higgsino can kinematically decay to final states involving stau leptons for tan $\beta=5$ (left) and tan $\beta=20$ (right). The points marked in red are already excluded.

In order to study the impact of light staus on the exclusion limit of Higgsinos, we perform a scan over the Higgsino-RH stau parameter space, where we fix the Bino mass parameter at 60 GeV and vary $\mu$ in the range [$100$-$1000$] GeV and $M_{\tilde{e}_{3R}}$ in the range [85-500] GeV. Fig. 16 shows the $R$-values of the present anaylses as implemented in SmodelS in the plane of mass of the Higgsino-like neutralino and the lightest stau. We present the result for a range of masses where the Higgsino can kinematically decay to final states involving stau leptons for tan $\beta=5$ (left) and tan $\beta=20$ (right). The points marked in red are already excluded, i.e., have a $R$-value greater than unity. We observe that the presence of such light staus affect the exclusion limit of Higgsino NLSPs. The impact of light RH staus become stronger for higher tan $\beta$ values where the branching of Higgsinos to stau leptons increase, thereby, weakening the result.

Benchmarks with light staus (mass parameters in GeV)				$M_{h}[\Delta_{M_{h}}^{FH}]~{}$[GeV]	$\sigma_{SI}\times\xi\times 10^{-10}$ [pb]
$\mu<0$	$Z$-funnel	BP5	$M_{t}=173.21$, $M_{1}=44$, $M_{2}=2000$, $\mu=-124$, tan$\beta=5$, $M_{A}=3000$,	125.86[$\pm$0.96]	$7.45\times 10^{-4}$
			$M_{\tilde{Q}_{3L}}=M_{\tilde{t}_{R}}=M_{\tilde{b}_{R}}=A_{t}=10000$, $M_{3}=3000$, $M_{\tilde{e}_{3R}}=85$
	$h$-funnel	BP6	$M_{t}=173.21$, $M_{1}=68$, $M_{2}=2000$, $\mu=-150$, tan$\beta=50$, $M_{A}=3000$,	125.65[$\pm$0.63]	0.137
			$M_{\tilde{Q}_{3L}}=M_{\tilde{t}_{R}}=M_{\tilde{b}_{R}}=5000$, $A_{t}=-5000$, $M_{3}=3000$, $M_{\tilde{e}_{3R}}=85$
$\mu>0$	$Z$-funnel	BP7	$M_{t}=173.21$, $M_{1}=44$, $M_{2}=2000$, $\mu=500$, tan$\beta=50$, $M_{A}=6000$,	125.11[$\pm$0.99]	0.095
			$M_{\tilde{Q}_{3L}}=M_{\tilde{t}_{R}}=M_{\tilde{b}_{R}}=4500$, $A_{t}=4000$, $M_{3}=5000$, $M_{\tilde{e}_{3R}}=85$
	$h$-funnel	BP8	$M_{t}=173.21$, $M_{1}=62$, $M_{2}=2000$, $\mu=500$, tan$\beta=20$, $M_{A}=6000$,	124.77[$\pm$0.97]	0.152
			$M_{\tilde{Q}_{3L}}=M_{\tilde{t}_{R}}=M_{\tilde{b}_{R}}=4500$, $A_{t}=4000$, $M_{3}=5000$, $M_{\tilde{e}_{3R}}=150$

Table 6: Parameters corresponding to four benchmark points satisfying all the present constraints from the $\mu>0$ and $\mu<0$ scenarios along with their scaled SI DD cross-sections. The mass of the Higgs boson $M_{h}$ and the uncertainty in $M_{h}$ computed by FeynHiggs ($\Delta_{M_{h}}^{FH}$) are also shown.

Subsequently, we study the prospect of a XGBOOST based analysis for benchmark points where we have a light stau, which are still allowed by the SUSY searches implemented in SModelS-2.3.0. We select four benchmarks $-$ from the $Z$ and $h$ funnel regions of $\mu>0$ and $\mu<0$, each. These are listed in Table 6 with the relevant soft parameters. BP5 and BP6 correspond to the light Higgsino benchmarks studied in the previous section (BP2 and BP3 in Table 3) to which we add a light stau with a physical mass around 90 GeV. Due to light RH staus, a region of parameter space in the previously excluded $Z$ funnel for positive $\mu$ opens up, and survives both LZ and electroweakino direct search bounds. We select a benchmark from this region, called BP7, with moderate Higgsino mass around 500 GeV and a very light stau with $M_{\tilde{e}_{3R}}=85$ GeV. For the $h$ funnel of $\mu>0$, Higgsinos up to 850 GeV masses were excluded for a 100% branching to the LSP. If RH staus are brought below the Higgsino, then Higgsinos around 500 GeV masses can still satisfy the collider limits for particular masses of the RH stau, like with $M_{\tilde{e}_{3R}}=150$ GeV, which we choose as BP8.

Number of events for $\mathcal{L}=300$ fb${}^{-1}$		BP5	BP6
Backgrounds	$lll\nu$	190.7	105.6
	$ZZ$, leptonic, $2j$ matched	39.4	26.5
	$t\bar{t}$, leptonic	3500.0	1520.5
	$VVV$, inclusive	14.9	7.1
	$Wh$, inclusive	61.0	27.7
	$Zh$, inclusive	22.8	13.5
	ggF $h\rightarrow ZZ$, leptonic	1.2	0.5
	VBF $h\rightarrow ZZ$, leptonic	0.2	0.05
	$t\bar{t}h$, inclusive	11.0	6.4
	$t\bar{t}W$, leptonic	4.7	2.1
	$t\bar{t}Z$, leptonic	2.4	1.4
	Total	3848.3	1711.4
Signal		5937.9	3513.9
Significance with 20% systematic uncertainty		5.51	6.81
Significance with 50% systematic uncertainty		2.21	2.74

Table 7: Number of events from individual background processes and the signal surviving a threshold of 0.9 on the XGBOOST output from two models trained on benchmarks BP5 and BP6, respectively, along with the signal significance for $\mathcal{L}=300$ fb${}^{-1}$.

Number of events for $\mathcal{L}=300$ fb${}^{-1}$		BP7	BP8
Backgrounds	$lll\nu$	40.1	35.7
	$ZZ$, leptonic, $2j$ matched	4.9	3.43
	$t\bar{t}$, leptonic	1860.1	1659.2
	$VVV$, inclusive	16.7	16.6
	$Wh$, inclusive	13.8	10.5
	$Zh$, inclusive	3.0	2.1
	ggF $h\rightarrow ZZ$, leptonic	0.02	0.02
	VBF $h\rightarrow ZZ$, leptonic	0.004	0.004
	$t\bar{t}h$, inclusive	12.6	10.6
	$t\bar{t}W$, leptonic	8.7	9.2
	$t\bar{t}Z$, leptonic	2.7	2.8
	Total	1962.6	1750.0
Signal		406.8	170.2
Significance with 5% systematic uncertainty		3.56	1.70
Significance with 10% systematic uncertainty		1.90	0.92

Table 8: Number of events from individual background processes and the signal surviving a threshold of 0.98 on the XGBOOST output from two models trained on benchmarks BP7 and BP8, respectively, along with the signal significance for $\mathcal{L}=300$ fb${}^{-1}$.

When the Higgsino decays to staus, we have final states enriched with tau leptons. They in turn decay to electrons, muons, or pions. In our analysis for these benchmarks, we perform a similar analysis like the $3l+$MET, including the hadronic decays of the tau leptons. Table 7 shows the expected number of our signal benchmark points BP5 and BP6 along with the background events for a threshold of 0.9 on our XGBOOST output. We also quote the significance by assuming a 20% (50%) systematic uncertainty. We find that both BP5 and BP6, belonging to the $\mu<0$ scenario, can be probed with our analysis at the Run-3 of the LHC using 300 fb${}^{-1}$ of data, with a signal significance $\gtrsim 2\sigma$, despite a large systematic uncertainty of $50\%$. For positive $\mu$, the two benchmarks have higher Higgsino masses (500 GeV), and therefore, lower production cross-sections. Hence, we put a stronger XGBOOST threshold to reduce the backgrounds further. Table 8 shows the expected number of our signal benchmark points BP7 and BP8 along with the background events for a threshold of 0.98 on our XGBOOST output. For heavier Higgsinos, the systematic uncertainties might be much smaller. We, therefore, quote the significance by assuming a 5% (10%) systematic uncertainty. For BP7, we find that if the uncertainty can be brought down to 5%, we can achieve more than 3$\sigma$ significance, while for BP8, we require the uncertainty to be around 2% to have 3$\sigma$ significance.

VIII The thermal neutralino in non-standard cosmology

Until now, we have worked within the framework of standard cosmological scenario. However, if the neutralino DM is produced thermally in a non-standard cosmology, then the relic density constraint can be relaxed. This can happen, for example, due to entropy injection in the Universe from the late decay of some particle after the DM freezes out. In this scenario, even if the relic density of the DM at freeze-out is much larger than the present observed relic, it can be diluted due to the increase in the entropy density of the Universe. Ref. [112, Fig. 8] shows that for a mass of WIMP DM in the range 40-60 GeV, the quantity $\langle\sigma v\rangle$, which is the annihilation cross-section multiplied by the DM velocity, can be reduced from the usual value of $10^{-9}$ GeV${}^{-2}$ in standard cosmological scenarios to a value below $10^{-15}$ GeV${}^{-2}$ in non-standard cosmologies. Therefore, in the non-standard cosmology, we can allow for very small DM annihilation cross-sections. The allowed parameter space is not restricted to the funnel regions and the Higgsinos can have masses as large as 2 TeV, or even heavier in both the $Z$ and $h$ funnels. Having very small couplings, these points will also satisfy the DD bounds and therefore, will provide interesting benchmarks for probing non-standard cosmologies. Fig. 17 shows the parameter space surviving the LEP, flavor, Higgs and DM DD constraints in the mass and relic density of the LSP plane, with the colorbar showing the mass of $\tilde{\chi}_{2}^{0}$.

Figure 17: The relic density at freeze-out of the LSP neutralino DM as a function of its mass ($x$-axis) and the mass of the Higgsino-like $\tilde{\chi}_{2}^{0}$ (colorbar), in the region of parameter space satisfying LEP, flavor, Higgs constraints, and DM DD constraints from the XENON-1T, PICO-60, PandaX-4T, and LZ experiments. It is assumed that entropy injection in non-standard cosmology can reduce the relic density of overabundant DM.

Figure 18: Scaled SI DM-nucleon cross-section ($\sigma_{SI}\times\xi$) for $\mu>0$ in the non-standard cosmological scenario, where the entropy dilution can provide the correct observed DM relic density upto a freeze-out density of around 10, as a function of the mass of the LSP neutralino DM in the region of parameter space satisfying LEP, flavor, Higgs constraints, relic density ($0.122<\Omega<10$), and DM DD constraints from the XENON-1T, PICO-60, PandaX-4T, and LZ experiments.

Fig. 18 shows the scaled SI DD cross-section as a function of the DM mass for scenarios where DM is found to be overabundant in the standard cosmological scenario, however, can satisfy the observed relic density due to non-standard cosmology. We apply all the constraints from the set “After LZ”, where only the relic density constraint is modified to $0.122<\Omega h^{2}<10$, where we assume that the non-standard cosmology can dilute the DM relic density by a factor of about 100. We observe that new regions of the parameter space now survive the experimental constraints, in both the $Z$ and $h$ funnel regions. This would allow us to identify scenarios of non-standard cosmology, depending on the nature of the observed signals. For instance, if one observes a DM in the $Z$-funnel region and simultaneously a LHC signal for Higgsinos heavier than 500 GeV, it might indicate non-standard cosmology with thermal production of the neutralino DM. Even in the $h$-funnel region, we have an idea of the minimum DD cross-section values that can survive in standard cosmology. Observing a signal with the future LZ data might hint towards a non-standard cosmological picture.

Figure 19: Branching fractions of the two Higgsino-like neutralinos with $\mu=1$ TeV to $\tilde{\chi}_{1}^{0}Z$ (left) and $\tilde{\chi}_{1}^{0}h$ (right) when the staus are heavy.

Figure 20: Branching fractions of the Higgsino-like chargino (left) and the two neutralinos (right) to final states involving various masses of RH staus.

Additionally, for a DM signal from DD experiments in the $h$ funnel, it is interesting to explore whether collider experiments can provide any hint of non-standard cosmology. If we observe a signal for heavy Higgsinos, having masses around a TeV for high tan $\beta$, then it might suggest a non-standard cosmological scenario. To answer the question whether we can get an idea of the tan $\beta$ value from the signal, we studied the variation in branching of the Higgsinos to various final states as a function of tan $\beta$. Fig. 19 shows the branching of the neutral Higgsinos to the LSP and Z/h. We find that although there is a tan $\beta$ dependence, the sum of Br$(\tilde{\chi}_{2}^{0}\rightarrow\tilde{\chi}_{1}^{0}Z)+$Br$(\tilde{\chi}_{3}^{0}\rightarrow\tilde{\chi}_{1}^{0}Z)$ (or with Higgs boson in the final state), has no variation with tan $\beta$. Therefore, it won’t be possible to estimate the tan $\beta$ value from the branching fraction of Higgsinos to these final states. We then turn to the case of light staus, such that the staus have masses between the Bino and the Higgsinos. Fig. 20 shows the branching of the chargino and the two neutralinos decaying to stau and tau for three different stau masses. We observe a very clear tan $\beta$ dependence in this case, especially for lighter staus. If we observe a signal of Higgsinos and could identify the final state with staus, we can get an idea of the tan $\beta$ from the branching fraction of Higgsinos decaying to staus. The branching to staus is higher for high tan $\beta$, where relic density cannot be satisfied for a TeV scale Higgsino within standard cosmology, and therefore, points toward a non-standard cosmological scenario.

IX Conclusion

In summary, this study shows that the current experiments, especially the recent results from electroweakino searches at the LHC and dark matter DD measurements from LZ, have severely constrained the $\mu>0$ scenario for a light neutralino thermal DM in the pMSSM with 10 free parameters. The DD result of the LZ collaboration is a strong constraint and it affects different regions of the parameter space in pMSSM depending on the constructive and destructive interferences between the light and the heavy CP-even neutral Higgs bosons. For heavy staus, the $Z$-funnel is completely excluded, and only heavy Higgsinos ($M_{\tilde{\chi}_{1}^{0}}\gtrsim 850$ GeV) are allowed in the $h$-funnel region. In the $\mu<0$ scenario of the same model, the allowed parameter space consists of either Higgsinos heavier than $\sim 850$ GeV in the $h$-funnel or restricted to a narrow region of light Higgsinos having mass of 125-160 GeV in the $Z$ and $h$-funnels, unlike the $\mu>0$ scenario. For light Higgsinos, there is a constructive interference between the $h$ and $H$ contributions to the DD cross-section for positive $\mu$. Therefore, it is not possible to evade the LZ bounds. Moreover, the DM constraints on relic density and the DD SI cross-sections severely constrain any simplified extension of SM with a Majorana fermionic DM coupling only to the discovered Higgs boson. There is a destructive interference between the contributions from the two CP-even neutral Higgs bosons for negative $\mu$, thereby, opening up allowed parameter space with light Higgsinos. Our XGBOOST analysis with the $3l+$MET signature shows that these light Higgsinos could be probed with the Run-2 data of 137 fb${}^{-1}$, if the systematic uncertainties lie within 20-30%. These benchmarks are still allowed by the available recasting frameworks, like SModelS and CheckMATE, which can translate the result of Wino-like NLSPs, provided by experimental collaborations, for Higgsino-like NLSPs. Thus, they form an important target for Run-3 searches.

The situation changes when we have light RH staus as NLSPs. They provide an additional annihilation channel for the DM and therefore, reduces the relic density by a factor depending on the mass of the stau. Furthermore, they also affect the exclusion limits of Higgsinos in collider searches, when the Higgsino can decay to staus with significant branching fractions. Due to both these effects, we get allowed region of parameter space even in the $Z$ funnel of $\mu>0$ with lighter staus. In the $h$ funnel of positive $\mu$, the presence of light staus relax the lower limit on Higgsino masses, and 500 GeV Higgsinos are still allowed in such a scenario. Our preliminary analysis of these benchmarks show that they are accessible at Run-3 of LHC, if the systematic uncertainties can be controlled. The future direct detection experiments, e.g. 1000 days of LZ, can probe both the funnels with light staus, as can be seen from Fig. 15. The status further changes when we go to non-standard cosmological scenarios, where the relic density can be satisfied by late injection of entropy in the Universe. Therefore, a large region of the parameter space becomes available. We discuss some ways to identify a non-standard history of our Universe from the combination of signals in future DM DD and collider experiments, which could not be realised assuming the standard cosmological model. The future DM DD experiments and the LHC Run-3 have promising prospects in exploring the remaining corners of the pMSSM parameter space with a light neutralino thermal dark matter.

Appendix A Impact of the constraints on the invisible branching of the Higgs boson

Figure 21: Allowed region of parameter space in the $M_{\tilde{\chi}_{1}^{0}}$-Br($h\rightarrow\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$) plane for $\mu>0$ (left) and $\mu<0$ (right) with all the constraints ’Before LZ’ (yellow), ’After LZ’ (light green), and from electroweakino searches (dark green). Dashed and solid black lines show the projected sensitivity for Br($h\rightarrow\tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$) from the ILC and CEPC experiments, respectively.

Appendix B Collider limit on right-handed staus

Fig. 22 shows the $R$-values of present analyses searching for stau leptons at the LHC, as recasted by the SModelS package in the plane of mass of the lightest stau which is RH, $M_{\tilde{\widetilde{\tau}}_{1}}$ and the mass of the lightest neutralino which is Bino-like, $M_{\tilde{\chi}_{1}^{0}}$ for a Higgsino mass parameter of 1 TeV. We observe that the $R$-values are all less than unity, implying that the current analyses are not sensitive to this region of parameter space. As a result, light staus in the mass range of 90-400 are still allowed, the lower limit coming from the LEP experiment [106]. We translate these $R$-values with the square root of luminosity for the future runs of LHC. We find that Run-3 with 300 fb${}^{-1}$ integrated luminosity will be sensitive to a small region of parameter space for $M_{\tilde{\chi}_{1}^{0}}\lesssim 20$ GeV and $M_{\tilde{\widetilde{\tau}}_{1}}$ between 200-250 GeV, and the HL-LHC run with 3 ab${}^{-1}$ of luminosity will be able to probe this whole region of parameter space.

Figure 22: SmodelS $R$-values of present analyses searching for stau leptons at the LHC in the $M_{\tilde{\widetilde{\tau}}_{1}}$-$M_{\tilde{\chi}_{1}^{0}}$ plane, where the lightest stau is RH and $M_{\tilde{\chi}_{1}^{0}}$ is Bino-like for $\mu=1$ TeV.