Supersymmetric explanation of the muon g-2 anomaly with and without stable neutralino

In this paper we explore the possibility of explaining the muon $g-2$ anomaly in various types of supersymmetric extensions of the Standard Model. In particular, we investigate and compare the phenomenological constraints in the MSSM with stable neutralino and the other types of scenarios where the neutralino is unstable. For the latter case we study the Gauge Mediated SUSY Breaking (GMSB) scenario with very light gravitino and the $UDD$-type R-Parity Violating (RPV) scenario. In the MSSM with stable neutralino, the parameter region favoured by the $(g-2)_\mu$ is strongly constrained by the neutralino relic abundance and the dark matter direct detection experiments, as well as by the LHC searches in the lepton plus missing transverse energy channel. On the other hand, the scenarios without stable neutralino are free from the dark matter constraints, while the LHC constraints depend strongly on the decay of the neutralino. We find that in GMSB the entire parameter region favoured by the muon $g-2$ is already excluded if the Next Lightest SUSY Particle (NLSP) is the neutralino, while some regions are still allowed if the NLSP is stau. In the RPV scenario, the LHC constraints are much weaker than the other scenarios and a wide region of the parameter space is still open for the muon $g-2$.


Introduction
The 4σ-level anomaly on the muon anomalous magnetic moment (muon g − 2), established by the measurements at the Brookhaven National Laboratory [1][2][3] and more recently at the Fermilab [4], is arguably one of the most interesting hints for physics beyond the Standard Model (BSM). The anomaly, expressed in terms of a µ = (g µ − 2)/2, amounts to ∆a µ = a BNL+FNAL µ − a SM µ = (25.1 ± 5.9) × 10 −10 , (1.1) where we use the averaged value of the two measurements as the experimental value and the expected value based on the Standard Model (SM) is taken from the white paper . 1 Taking the 4σ-level anomaly as a face value, it may be attributed to the effect of physics beyond the Standard Model. There is a wealth of literature pursuing this possibility. The list of BSM scenarios that can explain the (g−2) µ anomaly includes, for instance, models with leptoquarks [31][32][33][34], vector-like leptons [35][36][37], axion-like particles [38], and dark matter (DM) [39]. Among these scenarios, supersymmetry (SUSY) [40], in particular the Minimal Supersymmetric extension of the Standard Model (MSSM) [41], is especially attractive since it also offers other benefits, such as providing candidates for DM [42,43], gauge coupling unification [44], and radiative breaking of the electroweak symmetry [45].
The SUSY explanation of the (g−2) µ anomaly has been widely studied in the literature [46][47][48][49][50][51][52][53]. 2 The majority of these analyses assumes that the lightest supersymmetric particle (LSP) is the lightest neutralino,χ 0 1 , and it is stable due to the conservation of a discrete Z 2 symmetry called the R-parity. Moreover, many of the studies try to identify the lightest neutralino as the DM particle [58][59][60][61][62], selecting the parameter region in which the relic abundance of the neutralino, Ωχ0 1 , agrees with the observed density of the DM in the present Universe, Ω DM . In this setup, the parameter region favoured by the (g − 2) µ anomaly is severely constrained by the requirement Ωχ0 1 Ω DM and the upper bounds on the spin-independent DM-nucleon scattering cross-section set by the DM direct detection experiments [60][61][62]. Even if one gives up the possibility to identify the neutralino as the DM candidate, the requirement that the relic neutralinos do not overclose the Universe, Ωχ0 1 < Ω DM , and the constraint from the DM direct detection experiments apply and these exclude a large part of the parameter region preferred by the (g − 2) µ anomaly. The tension between the SUSY explanation of the (g − 2) µ anomaly and the DM constraint is not a coincidence. As we will discuss in the next section, the requirement on the SUSY parameters to give large contributions to (g − 2) µ is directly related to the enhancement of the spin-independent DM-nucleon cross-section and suppression of theχ 0 1 -χ 0 1 annihilation cross-section, leading to the overproduction of the relic neutralinos.
Another important constraint comes from the direct BSM searches by the ATLAS and CMS experiments at the LHC. The SUSY contribution to (g − 2) µ is large enough only when the masses of sleptons and electroweak (EW) gauginos are below O(500) GeV. Such light sleptons and electroweakinos (EWinos) would be produced at the LHC and detected in search channels looking for, e.g., multiple leptons and a sizable missing transverse energy (E miss T ). Providedχ 0 1 is the stable LSP, the slepton and electroweakino mass spectrum must be compressed in order to evade the LHC exclusion, so that decays of the produced SUSY particles do not produce high-p T leptons. where µ is the Higgsino mass parameter, tan β ≡ H 0 u / H 0 d is the ratio of the vacuum expectation values of the two Higgs doublets, and M 1 , M 2 ,m l L andm l R are the soft SUSY breaking masses for the Bino (U(1) Y gaugino), the Wino (SU(2) L gaugino), the left-handed slepton doublet and the right-handed slepton singlet, respectively. To reduce the number of free parameters we adopt the universal slepton mass assumption:m l 1 =m l 2 =m l 3 ≡m l L andm e R =m µ R =m τ R ≡m l R . With this assumption, the model can evade the strong constraints from lepton-flavour violating (LFV) processes such as µ → eγ. Strictly speaking, the universal slepton mass assumption is not necessary to avoid the LFV constraints. It is sufficient to assume that the charged slepton mass matrices, m l L and m l R , are diagonal in the same basis in which the charged lepton mass matrix is diagonal. However, in many known scenarios of SUSY breaking mediation, flavour universal sfermion masses are realised at the mediation scale. Although some of those scenarios predict particular gaugino mass relations, these ratios are relatively unconstrained when there are multiple sources of SUSY breaking mediation. 5 Our analysis is also applicable to such situations. We also assume for simplicity that the SUSY-breaking parameters and µ are real and do not contribute to the CP violating observables. 6 For large or moderate values of tan β (5 tan β 50), the one-loop SUSY contribution is approximated by [48] a SUSY µ a WHL µ + a BHL µ + a BHR µ + a BLR µ .

(2.2)
In the mass-insertion approximation, each term on the right-hand side is represented by the re- 5 For example, in the minimal gauge mediation, the messenger fields are a vector-like pair of SU (5), (5,5), and the gaugino masses are proportional to the square of the corresponding gauge couplings, Mi ∝ g 2 i . However, this relation can be easily modified by allowing mass splitting between the triplet and doublet parts of (5,5), while keeping the flavour universality of sfermion masses intact. The mixed modulus and anomaly mediation (or the mirage mediation) [77][78][79][80][81] and a non-singlet F-term breaking scenario (see e.g. [82]) are other examples to achieve non-standard gaugino mass relations. 6 Again, such an assumption is not required once we allow the possibility of cancellation among the imaginary parts of various parameters. However, we do not consider such a contrived case and simply adopt real parameters. spective diagram in Fig. 1, which is given by (2.3d) The functions f C and f N are given by and mμ L = (m 2 l L + m 2 µ + ∆ L ) 1 2 and mμ R = (m 2 l R + m 2 µ + ∆ R ) 1 2 with ∆ L ≡ m 2 Z (− 1 2 + sin 2 θ W ) cos 2β and ∆ R ≡ m 2 Z sin 2 θ W cos 2β are the diagonal (LL and RR) entries of the smuon mass matrix, and mν µ , m µ and θ W are the masses of the sneutrinos and the muon and the weak mixing angle respectively.
If all SUSY particles have the same mass, the a WHL µ is an order of magnitude larger than the other contributions. Although we define our parameter planes based on the four 1-loop diagrams described above, we do not use the mass insertion approximation in our numerical study. To compute a SUSY µ we use GM2Calc [83] instead, which includes the full 1-loop contribution as well as the leading 2-loop effect.

Parameter planes
Since there are six parameters in the problem, exploration and visualization of the whole parameter space are challenging. In this paper, we instead define eight two-dimensional parameter planes motivated by the above one-loop contributions and systematically study the interplay between a SUSY µ and the other phenomenological constraints. As a SUSY µ at the one-loop level is proportional to tan β, we take tan β = 50 as the default value in order to maximise the SUSY contributions and specify three of the remaining five parameters to define two-dimensional parameter spaces.
We set the masses of gluinos, squarks, and heavy Higgs bosons arbitrary large in our analysis because a SUSY µ at the one-loop level is independent of their mass and recent SUSY searches at the LHC have placed stringent lower bounds on the masses of the coloured SUSY particles [84][85][86].
The eight parameter planes are summarized in Table 1. In each plane, three mass parameters among the five (cf. Eq. (2.1)) are set as light as O(100) GeV and the remaining two are set heavier, so that one of the four contributions in Eq. (2.3) dominates. We prepare two planes, which result in different collider phenomenology, for each of the four scenarios. In these parameter planes the LSP among the MSSM spectrum is always given by the lightest neutralino,χ 0 1 . Beyond MSSM, the gravitino can be lighter thanχ 0 1 , which we study in detail section 7.

WHL planes
The WHL planes are characterized by Thanks to the first requirement, the WHL contribution a WHL µ yielded by the WHL diagram becomes large enough to explain the muon g − 2 anomaly, where the last requirement guarantees a WHL µ is positive (cf. Eq. (2.3a) and Fig. 1 As discussed in Sec. 1, if the LSP isχ 0 1 and stable, the parameter space under these requirements receives stringent constraints from SUSY searches at the LHC, such as searches for multi-lepton plus large missing transverse momentum (E miss T ) signature. In order to avoid emergence of high-p T leptons, we require two of the three small parameters are close to each other and define the following two planes: • WHL µ (M 2 vs µ) plane: m l L = min(M 2 , µ) + 20 GeV, tan β = 50, (M 1 ,m l R ) > a few TeV In both planes the MSSM-LSP isχ 0 1 . Bothχ 0 1 andχ ± 1 are composed either of the Wino (M 2 < µ) or the Higgsino (µ < M 2 ). Theχ 0 1 andχ ± 1 are almost mass degenerate and the scale of the mass splitting is given by v 2 EW /mχ0 1 . On the WHL µ plane,m l L is set close toχ 0 1 . By this assumption, the decayl L → l +χ 0 1 , which may occur in the cascade-decay chain from the SUSY particles produced at the LHC, does not yield a high-p T lepton. The models are therefore elusive at the LHC, even though the electroweakino pair-productions have sizable cross sections. At the same time, ifχ 0 1 is stable and mostly Wino-like, χ ± 1 might be long-lived enough so thatχ ± 1 can be detected by the disappearing track searches at the LHC.
On the WHL L plane,χ 0 1 is close in mass either to the other neutralinosχ 0 2,3 and charginoχ ± 2 (M 2 <m l L ), or tol L (m l L < M 2 ). The latter case is similar to the WHL µ plane and thus elusive at the LHC. In the former case with M 2 <m l L , high p T leptons may be produced from the slepton 7 Only the relative sign between µ and Mi (i = 1, 2, 3) is physical. We take Mi > 0 throughout the paper.
decays in the pp →l Ll * L process. However, the cross-section of the slepton pair production is much smaller compared to the electroweakino pair production at the same mass and the LHC constraint in this region can be relatively easily evaded.

BHL planes
The BHL diagram is obtained from the WHL diagram by the replacement ( W ± , H ± ,ν µ ) → ( B, H 0 ,μ L ), and thus a BHL µ dominates over the other three contributions if where the last condition guarantees a BHL µ > 0. Similarly to the WHL scenarios, we define two planes by Unlike the WHL planes,χ ± 1 is always Higgsino-like in the BHL planes and short-lived enough not to be detected by the disappearing track searches at the LHC. Also,χ 0 1 is composed of Bino or neutral Higgsino. Because the pure-Bino cross-section is vanishingly small when squarks are decoupled, SUSY-particle production is only through the Higgsino pair-production or slepton pair-production. Their cross sections are smaller than the Wino pair-production (compared at the same masses) and the parameter space tends to be less constrained by the current LHC searches.

BHR planes
The BHR diagram can be obtained from the BHL diagram simply by flipping the chirality of the muon and smuon. Hence a BHR µ becomes dominant if where, as can be seen in Eq. (2.3c), we take µ < 0 and M 1 > 0 to have a BHR µ > 0, which fits the observed (g − 2) µ anomaly.
The BHR planes are defined as follows: The nature of MSSM-LSP and the LHC phenomenology is similar to those in the corresponding BHL planes, but sneutrinos are always decoupled and do not participate in LHC signatures.

BLR planes
The BLR diagram in Fig. 1 is special because, unlike the other diagrams, the chirality flip occurs in the smuon sector. The size of the smuon chirality flip is given by the off-diagonal element of the smuon mass-squared matrix and thus proportional to the µ parameter (cf. Eq. (2.5)). The BLR contribution is therefore proportional to µ as well as to tan β, whereas all the other one-loop contributions are vanishing in the limit of |µ| → ∞. Therefore, a BLR µ becomes dominant if smuons and Bino are as light as O(100) GeV and µ tan β is large enough.
Under these conditions, we have to pay attention to two constraints from the stau sector, i.e., the LEP limit on the stau mass and the vacuum stability. These constraints are critical because the SUSY breaking mass parameters are common for all generations. The mass-squared matrices for the left-and right-handed charged sleptons are given by at the tree-level, where l = e, µ, τ . Here, m l is the charged lepton mass, θ W is the weak mixing angle, and A l is the SUSY-breaking trilinear coupling of the slepton, which we set to be zero. Becausem l L andm l R are common for all generations and the off-diagonal element is proportional to the lepton mass, the left-right mixing of staus is larger than that of smuons and the lighter stau (heavier stau) becomes lighter (heavier) than the lighter smuon (heavier smuon). The LEP bounds on the stau mass, mτ 1 > 81.9 GeV [87], is therefore critical in this scenario. In addition, with the large left-right mixing of the staus, the effective potential may develop a global minimum with non-zero vacuum expectation values of the stau fields. If the vacuum transition rate from our EW vacuum to this charge-breaking global minimum is significantly faster than the age of the Universe, the parameter point is excluded. The vacuum stability thus provides an upper bound on |µ tan β| [88,89]. Accordingly, the BLR planes are characterized by where (µ tan β) max is the maximal value of |µ tan β| under the vacuum-stability condition given in Refs. [88,89], which is dependent onm l L andm l R . The last requirement guarantees a BLR µ > 0 to explain the muon g − 2 anomaly and a WHL µ is subdominant because of the third condition, while a BLR µ is enhanced by the first and second ones. However, one should note that a BHL µ or a BHR µ may not be negligible because |µ| has an upper bound µ max ≡ (µ tan β) max / tan β.
We define two BLR planes, in whichχ 0 1 is the MSSM-LSP, as follows: • BLR 10 (m l L vsm l R ) plane: In both planes, we fix µ = µ max so that BLR contribution is maximised without violating the vacuum stability condition. SUSY-particle production at the LHC is through the slepton pairproduction or the Higgsino pair-production, but the latter is available only if µ is at O(100) GeV due to the vacuum stability condition. On the BLR 10 plane, the condition is relaxed due to the non-default value of tan β, and thus the LHC constraints originated from the Higgsino production is relaxed.

Phenomenological constraints 3.1 The neutralino relic density
In the cases with the lightest neutralinoχ 0 1 being stable, one must consider the phenomenological constraint that the relic neutrinos do not overclose the Universe. In the early Universe, neutralinos are coupled to the thermal bath and their number density is fixed by the temperature. While the Universe expands and the temperature drops, their number density decreases and at some point the lightest neutralino freezes out as the dark matter. Those relic neutralinos contribute to the energy density of the dark matter in the present Universe. Clearly, the energy density of the relic neutralino, Ωχ0 1 , cannot exceed the measured DM density, Ω DM . Namely, we require Ωχ0 1 ≤ Ω DM in our analysis for the case of stable neutralino. If the equality is satisfied in the above condition, the thermal relic neutralino explains the observed density of the dark matter, while if the inequality is unsaturated, one needs extra components (e.g. axions) to explain the dark matter abundance.
It is well known that if the neutralino is Higgsino-or Wino-like, they can annihilate rapidly into a pair of EW bosons. Due to the relatively large annihilation rate, they can stay in the thermal bath rather longer and the number density at the time of thermal decoupling tends to be sufficiently low. For the pure-Higgsino neutralino, the bound Ωχ0 1 ≤ Ω DM is translated to mχ0 1 ∼ |µ| 1 TeV [90]. For the case of pure-Wino, this condition corresponds to mχ0 1 ∼ |M 2 | 3 TeV [90]. These mass scales are significantly higher than those required to fit the (g − 2) µ anomaly. Therefore, the neutralino thermal relic abundance does not give strong constraints on the SUSY (g − 2) µ solution when the neutralino is Higgsino-or Wino-like Unlike Higgsino-and Wino-like neutralinos, the annihilation cross-section is suppressed for the Bino-like neutralino, since Bino is gauge-singlet and does not directly couple to EW bosons. The leading annihilation process comes from a t-channel sfermion exchange diagram, but the Majorana nature of Bino forces the s-wave cross-section to be chirally suppressed. In order to bring the Bino relic abundance into the acceptable level, one needs to resort the coannihilation mechanism [91], in which the Bino-number changing process, ξ SM B ↔ ξ SM η, is operative, where ξ SM and ξ SM are SM particles and η is a coannihilation partner. For this mechanism to work, the coannihilation partner must have a large annihilation rate and the mass gap between η and B must be smaller than 5 − 10 % of the Bino mass.
The Bino-like neutralino appears on the BLR plane and the BHL and BHR planes with M 1 < |µ|. These regions are tightly constrained by the neutralino relic abundance condition. In the BLR planes, the Bino mass is placed 20 GeV below mτ 1 , such that the Bino may efficiently coannihilate withτ 1 .

The dark matter direct detection
In the cases with the lightest neutralinoχ 0 1 being stable, the thermal relic neutralinos in the present Universe may be detected at the dark matter direct detection (DMDD) experiments. 8 These experiments are located in deep underground and look for a delicate signature from nucleons recoiled by dark matter particles. The DMDD experiments have not found any convincing evidence for such events so far, which in turn set a strong upper limit (UL) on the DM-nucleon scattering cross-section as a function of the DM mass. For the neutralino DM, the most stringent constraints are placed on the spin-independent (SI) cross-section and we have σχ We generally adopt this condition throughout our analysis for stable neutralinos. The spin-independentχ 0 1 -nucleon scattering is obtained by the s-channel squark exchange or the t-channel Higgs (h, H, A) exchange diagrams. Since our analysis adopts decoupled squarks and heavy Higgses, the contribution comes entirely from the h exchange process, which requires theχ 0 1 -χ 0 1 -h 3-point interaction. In the MSSM, however, the pure-Higgsino and pure-gaugino do not have this interaction, i.e. the coupling for the H-H-h and λ-λ-h are vanishing in the gauge eigenbasis, where λ is the EW gauginos ( B, W 0 ), while the couplings for H-λ-h are non-vanishing. This implies that the SI cross-section is enhanced when both Higgsino and one of the EW gauginos are light and largely mixed inχ 0 1 . It happens that this is precisely what is required to have large WHL, BHL and BHR contributions to (g − 2) µ as discussed in the previous section. As we will see in the next section, the current constraint from the DMDD experiments excludes large parts of the WHL, BHL and BHR planes if the neutralino is stable.

The LHC constraints
As discussed in the previous section, both sleptons and electroweakinos need to be as light as O(100) GeV. More precisely, at least one ofμ L andμ R needs to be light in any scenarios. Higgsinos and one of the EW gauginos must be light in the WHL, BHL, and BHR scenarios, whereas only the Bino needs to be light in the BLR scenario. Since these light particles should be produced at the LHC, the models are constrained by the direct BSM searches by the ATLAS and CMS collaborations at the LHC.
If the LSP is the lightest neutralinoχ 0 1 and is stable, production of these SUSY particles contributes to the lepton plus E miss T channel. There are ATLAS and CMS analyses targeting exactly the same kinematics and final states, and they set the upper limit on the cross-section for their simplified model as a function of the masses relevant to the kinematics. Whenever this is the case, we assess the exclusion of our scenario at a particular mass point using the methodology outlined in [93,94]. Namely, we confront the simplified model cross-section limit published by an experimental collaboration with the production cross-section times the branching ratios to the same final state of our model. This comparison must be performed with consistent mass assumptions between the simplified model and our scenario. The procedure is described in detail in section 5.1.
Ifχ 0 1 is not the LSP, or if it is the LSP but unstable, the LHC signature drastically changes. In the baryonic RPV scenario studied in section 6,χ 0 1 decays into three light-flavour (anti-)quarks, χ 0 1 → qqq orqqq. Events at the LHC then lose their characteristic E miss T signature, but can still be searched for in the multi-jet with small E miss T channel if E miss T is yielded by neutrinos produced in the upper stream of the SUSY-particle cascade decays. In the gravitino ( G) LSP scenarios, the signature largely depends on the type of the next-to-the-lightest SUSY particle (NLSP). If the gravitino is almost massless and NLSP isτ 1 , the large E miss T in the stable neutralino signature is replaced by the moderately large E miss T and energetic taus produced from the NLSP decays τ 1 → τ G.
These unconventional signatures cannot be analysed systematically in the simplified model framework by ATLAS and CMS, and the methodology mentioned above cannot be used to assess the exclusion. In order to confront these scenarios with the LHC data, we resort the Monte Carlo simulation. At each mass point, we generate inclusive signal events with Pythia 8 [95,96] for the

Analysis procedure
We use GM2Calc [83] for calculation of a SUSY µ , which includes the leading two-loop contributions. For consistency, we take physical masses and mixing matrices of non-colored SUSY particles from the GM2Calc output. These observables are calculated at tree-level but the corrections to the lepton Yukawa couplings that can be sizable for large µ tan β [89,101] are included. We also include the EW radiative correction to the chargino mass, mχ± 1 , using the fitting function provided in [102].
This correction is very small, mχ± 1 ∼ O(100) MeV, and barely affects a SUSY µ , whilst it may give significant effect on the chargino lifetime through the mass splitting, ∆m ±,0 = mχ± is Wino-like. The decay widths and branching ratios of sparticles are calculated with SDecay [103] within the SUSY-HIT package [104]. In the compressed mass region, ∆m ±,0 < 1 GeV, we calculate and include theχ + 1 → π +χ0 1 mode since SDecay does not support this decay mode. For the scenarios whereχ 0 1 is stable, the neutralino relic density and σχ

Stable neutralino
In this section, we discuss the eight parameter planes defined in section 2.1 (cf. Table 1) assuming χ 0 1 is LSP and stable due to the R-parity conservation. The constraints on the neutralino relic density (section 3.1) and the DMDD (section 3.2) are of importance, while the LHC experiments benefit from the large E miss T signature discussed in section 3.3.

Implementation of LHC constraints
In the stable neutralino case, the main LHC constraint comes from the lepton plus E miss T channel. In order to constrain our parameter planes we use the simplified model limits on the cross-section and the limit on the mass vs lifetime plane from the disappearing track searches for the Wino-like LSP scenario listed in Table 2.
CMS + − (search with an opposite charge same flavour lepton pair plus E miss T ) [107] interpreted their data for the˜ +˜ − → ( +χ0 1 )( −χ0 1 ) topology and provided the cross-section limit on the (m˜ vs mχ0 1 ) plane, assuming a 100% branching ratio. We impose this constraint not only on the crosssection of˜ +˜ − → ( +χ0 1 )( −χ0 1 ) process, but also on the cross-section times branching ratios 9 of This is because with this condition the decay products of η ( ) →χ 0 1 are very soft and do not alter the signal efficiency significantly. Whenχ 0 1 is Wino-or Higgsino-like, χ ± 1 (as well as χ 0 2 in the Higgsino-like case) is quasi mass degenerate with χ 0 1 and the above condition is generally satisfied. The corresponding cross-section limit is read from the original (m˜ vs mχ0 1 ) plane by mapping (mξ, mη) → (m˜ , mχ0 1 ). In the same philosophy, we also use this simplified model limit to constrain the cross-section times branching ratios ofηη The mass degeneracy condition here is often satisfied since the slepton mass is placed 20 GeV above the LSP mass parameter in some of our parameter planes.
When χ 0 1 is Wino-like, the mass difference ∆m ±,0 = mχ± 1 − mχ0 1 becomes very small and the decay width ofχ ± 1 is phase-space suppressed. In this region, the lifetime ofχ ± 1 can become as large as cτχ± 1 decays inside of the tracking system into unobservably soft particles, X soft , plus invisibleχ 0 1 , it can give a disappearing track (DT) signature. ATLAS and CMS have been looking for this type of signature [110,111]. To constrain our parameter space with their experimental results, we use the 95% CL limit on the (mχ±  by the latest a µ experimental results [4] is depicted with green and yellow bands, corresponding to one and two sigma agreement respectively. The bottom left region shaded by the dull yellow colour corresponds to a µ bigger than experimental value by more than 2 sigma. Hatched region is excluded by dark matter abundance criterion. Blue and red shaded regions are excluded by LHC constraints. Light grey shaded area is excluded by Xenon1T experiment. Expected exclusion range by future DDMD experiments is also shown with appropriately labelled contours.

Results
In this subsection we identify the parameter region that can explain the (g − 2) µ anomaly and confront it with relevant phenomenological constraints. We begin with the WHL µ (left) and WHL L (right) planes shown in Fig. 2. Here, and also in the other plots below, the area highlighted with green (yellow) corresponds to the region where the predicted value of a µ agrees with the measurement within the 1 (2) σ accuracy, while in the bottom left region shaded by the dull yellow colour, the SUSY contribution to a µ is too large and the model is disfavoured by more than 2 σ.
The hatched region surrounded by the black curve indicates the region excluded by the neutralino overproduction, Ωχ0 1 ≥ Ω DM . At the boundary of this region (i.e. on the black curve), the relic neutrino can explain the full amount of observed dark matter density. One can see that in the WHL planes the (g − 2) µ region is located outside of the black hatched region, implying that extra components, such as axions, are needed to explain the dark matter.
The region shaded by light grey corresponds to the excluded area due to the upper bound on the spin-independent DM-nucleon scattering cross-section, σ UL SI , set by the latest XENON1T experiment [112]. More precisely, the condition Eq. (3.1) is violated in this area. As discussed in the previous section, σχ 0 1 SI is enhanced whenχ 0 1 has both Higgsino and gaugino components. In the WHL scenario, this happens when µ ∼ M 2 . In the WHL µ plane (Fig. 2 left), we can explicitly see the region with µ ∼ M 2 is excluded by the XENON1T constraint. In the WHL L case (Fig. 2  right), a half of the plane withm l L M 2 is excluded since µ ∼ M 2 is realised with the condition µ = min(M 2 ,m l L ) − 20 GeV.
The regions that can be probed by the future DMDD experiments are also shown in Fig. 2. In the WHL µ plane, the area between two green (blue) dashed lines are sensitive to PandaX-4t [113] (LZ [114] and XENONnT [115]). Also, the area below the magenta dashed line can be explored by ARGO [116]. In the WHL L plane, the region above the green and blue dashed lines can be probed by the PandaX-4t and LZ/XENONnT, respectively, and the entire region in the plane will be explored by ARGO. One can see that almost entire (g − 2) µ region in the WHL L plane will be tested in the future DMDD experiments.
Constraints from direct BSM searches at the LHC are also visible in Fig. 2. The red shaded area corresponds to the 95 % CL excluded region estimated using the slepton simplified model limit published in the CMS + − analysis [107], while the blue shaded area indicates the 95 % CL exclusion obtained from the published limit for a compressed slepton scenario given in the ATLAS soft-analysis [109] (see Table 2). In the WHL µ plane we also see the exclusion from the ATLAS DT search [110], which is shaded by the orange colour. We find that only these three constraints provide 95 % CL exclusions on the WHL parameter planes, although other potentially relevant constraints listed in Table 2 are also included in our analysis.
In the WHL µ plane we see that the region with µ M 2 1 TeV is excluded by the CMS + − analysis. In this region, the main SUSY processes, which contribute to the dilepton plus E miss T signal region targeted in this analysis, are pp → W +,0 W −,0 followed by e.g. W ± → ±ν and W 0 → ±˜ ∓ . Note that sincem l L = µ + 20 GeV in this region, the decay products of the slepton/sneutrino into the Higgsino-like LSP are too soft and may be neglected for this analysis. On the other side of the diagonal dotted line, we do not see the corresponding excluded region with M 2 µ for the Higgsino production, pp → H +,0 H −,0 . This is because the Higgsino production cross-section is smaller than the Wino and the main decay is given by H ±,0 → h W ±,0 , producing leptons much less often.
Exclusion from the disappearing track signature from long-lived Winos is visible in the top left corner of the WHL µ plane with µ M 2 . In this region, ∆m ±,0 = mχ± We see that the current ATLAS DT analysis [110] excludes the Winolike LSP up to M 2 ∼ 760 GeV. Exclusion from the disappearing track signature does not show up in the WHL L plane, since theχ 0 1 is Higgsino-like across this plane. The exclusion from the ATLAS soft-analysis can be seen in the regions with (i) µ 250 GeV and (ii) M 2 220 GeV, excluding the M 2 ∼ µ region. In these areas, the relevant SUSY process is ξξ → ( +η )( −η ) whereξ ( ) = (˜ orν) is a left-handed slepton doublet andη ( ) = (χ ± 1 orχ 0 1/2 ) is dominantly Higgsino-and Wino-like in regions (i) and (ii), respectively. Although the typical mass difference betweenξ andη is small (∼ 20 GeV), the ATLAS soft-analysis can detect the signal by exploiting soft leptons.
Taking both the LHC and DM constraints into account, one can find two (g − 2) µ favoured regions that are phenomenologically allowed in the WHL µ plane; (a) 250 GeV M 2 ∼m l L 500 GeV, 500 GeV µ 800 GeV and (b) 1 TeV M 2 2 TeV, 250 GeV µ ∼m l L 500 GeV. Both of these regions will be probed in future by the next generation DMDD experiments.
In the WHL L plane two regions are excluded by the CMS + − analysis. One region with M 2 >m l L is identical to the corresponding excluded region in the WHL µ plane, which has already been discussed above. The second excluded region has M 2 <m l L 700 GeV and µ = M 2 −20 GeV. In this region, the main SUSY signature isξξ → ( +η )( −η ) whereξ ( ) = (˜ orν) andη ( ) = W ±,0 . The blue shaded region is excluded by the ATLAS soft-, which is identical to the corresponding region in the the WHL µ plane discussed above. In the WHL L plane the allowed (g − 2) µ region appear around M 2 ∼ 1.5 TeV, µ ∼ 350 GeV, which has been already found in the WHL µ plane, and the entire (g − 2) µ region will be probed in the early stage of next generation DMDD experiments.
We show the results of our analysis for the BHL and BHR scenarios in the top and bottom panels of Fig. 3, respectively. Note that we take the sign of µ to be negative for the BHR planes in order to fit the (g − 2) µ anomaly due to the minus sign in Eq. (2.3c). First, we see in  that the SUSY masses required to fit the (g − 2) µ anomaly are much smaller for the BHL and BHR scenarios than in the WHL case. This is a reflection of the fact that the contributions from the BHL and BHR diagrams are an order of magnitude smaller than that from the WHL diagram if all SUSY masses are taken to be equal. We also see in both BHL µ and BHR µ planes that regions with M 1 < |µ| are severely constrained by the neutralino overproduction. This is because theχ 0 1 is Bino-like in this region and the relic abundance enhances due to a small annihilation cross-section for Binos. In these two planes one can also see that regions with |µ| ∼ M 1 are disallowed by the XENON1T constraint. This can be understood since the Higgsino-Bino mixing in χ 0 1 becomes large if |µ| ∼ M 1 , which enhances the spin-independent neutralino-nucleon scattering cross-section, as discussed in the previous section. XENON1T also excludes the M 1 ∼m l L and M 1 ∼m l R regions in the BHL L (top right) and BHR L (bottom right) planes, respectively. In these regions |µ| ∼ M 1 is realised due to the condition for the µ-parameter. We see that PandaX-4t (green dashed) can explore the entire 1 σ (g − 2) µ region in both the BHL and BHR scenarios. Also, the entire 2 σ regions will be probed by the LZ and XENONnT (blue dashed) experiments.
Taking DM constraints into account, only regions with M 1 |µ| + 100 GeV are allowed. In those regions the LHC constraints are rather weak. This is because only Bino has a significant mass gap from theχ 0 1 , though the Bino cross-section at the LHC is minuscule when the squarks are decoupled. The only LHC constraint that is significant in those areas comes from the limit imposed on the compressed electroweakino simplified model provided by the ATLAS soft-analysis (see Table 2). The excluded region from this constraint only appear in the BHL µ and BHL L planes, which are shaded with purple. The regions shown in these two plots are however identical.
We now show our result for the BLR 50 (tan β = 50) and BLR 10 (tan β = 10) planes in the top and bottom panels, respectively, in Fig. 4. As mentioned in subsection 2.1, in the BLR scenario the µ-parameter is set to the maximally allowed value from the vacuum stability constraint [88,89]. Large µ is preferred to ensure a large L-R smuon mixing necessary for the BLR scenario. In this regime, however, one of the stau mass eigenstates may become very light. In order to avoid stau LSP, the Bino mass, M 1 , is placed at 20 GeV below the lighter stau mass, mτ 1 . In the left panels, the dark blue region is excluded by the LEP stau mass bound 81.9 GeV [87].
In the previous section, we have argued that the relic neutralinos tend to be overproduced in the BLR scenario since theχ 0 1 is Bino-like. As expected, a large fraction of the parameter plane is excluded by this constraint as shown in the hatched black regions. In the region allowed by the Ωχ0 1 constraint, the stau-coannihilation mechanism is operative. We see that this constraint is milder for the BLR 10 plane with tan β = 10. In this plane, |µ| is significantly larger than that in BLR 50 , and the neutralino mass is therefore slightly heavier. Since we fix the stau-Bino mass difference, this results in smaller ∆m/m and more effective coannihilation.
Since µ and M 1 are fixed implicitly at each point in the BLR planes, we show, in the right panels of Fig. 4, contours of µ (black dashed in TeV) and M 1 (blue solid in GeV). As can be seen, M 1 , which is correlated with mτ 1 , is insensitive to tan β. This is because both mτ 1 and the vacuum tunnelling rate depend on the product µ·tan β (not on µ and tan β independently) and they remain unchanged when changing tan β and µ in such a way that their product µ tan β remains constant.
This also explains that from tan β = 50 to 10, the value of µ is changed by a factor of 5, as can be checked by comparing the top right and bottom right plots in Fig. 4.
We finally comment on the LHC constraint on the BLR scenario. As can be seen, the whole (g − 2) µ region is excluded by the limit coming from the slepton simplified model bounds in CMS + − analysis. Although the Higgsino production can give a non-negligible contribution in the bottom left region of the BLR 50 plane, the main process contributing to the exclusion is pp →˜ +˜ − followed by˜ ± → ±χ0 1 , which is identical to the process targeted in CMS + − analysis. It is remarkable that the entire (g − 2) µ region is excluded in our BLR planes. We note that this strong conclusion is based on our assumption of universal slepton masses. In the BLR planes µ tan β is taken to be at the maximum value allowed by the vacuum stability constraint. For such a high value of µ tan β, the lighter stau mass eigenvalue mτ 1 becomes significantly smaller than the masses of the other sleptonsμ/ẽ when the universal mass assumption is adopted. This leads to a large mass gap betweenμ/ẽ and Bino, since the Bino mass is set 20 GeV below mτ 1 to avoid stau being the LSP. With such a large mass gap, the slepton pair production leads to the final state with two high p T leptons plus large E miss T and the preferred (g − 2) µ region with light sleptons are excluded by the CMS + − analysis. In many popular scenarios (e.g. constrained MSSM, gauge mediation, anomaly mediation, mirage mediation, etc.) in which the universal slepton masses are generated at a very high scale (e.g. the GUT scale), some mass splitting is introduced at the electroweak scale show parameter regions and appropriate experimental constraints. Region of parameter space allowed by the latest a µ experimental results [4] is depicted with green and yellow bands, corresponding to one and two sigma agreement respectively. Hatched region is excluded by dark matter abundance criterion. Red shaded region is excluded by LHC constraints. Dark blue shaded area is excluded by LEP stau mass bound [87]. Expected exclusion range by future DDMD experiments is also shown with appropriately labelled contours. Plots on the right hand side of the figure show Bino mass (blue) and µ parameter (dashed black) values, which are implicitly fixed by other parameters.
via the renormalization group evolution [117]. 10 In this case, the SUSY breaking masses for staus are smaller than for the other sleptons, and the mass gap between the slepton and Bino must be taken even larger in the BLR planes. We therefore expect the BLR contribution cannot explain the (g − 2) µ anomaly consistently with the LHC constraints also in those popular scenarios. One way to avoid this situation is to relax the universality assumption on the slepton masses and take the stau mass to be much heavier than the other sleptons. Generally, however, such a splitting leads to too large LFV and is not phenomenologically viable, unless the charged lepton Yukawa matrix and the L-and R-slepton mass matrices are tuned in such a way that they can be simultaneously diagonalised in the same basis to a very good precision. To see whether the (g − 2) µ anomaly can be explained by the BLR contribution in this case, we introduce two new planes, in which the parameters are defined similarly to BLR 50/10 . The difference from the previous scenarios is that in the new planes,m l 3 andm τ R are fixed at 10 TeV independently from the first and second generation slepton masses (mass degeneracy is still assumed among the first two generations:m l L ≡m l 1 =m l 2 andm l R ≡m e R =m µ R ). Due to the large stau masses, the vacuum stability constraint is now relaxed and it is always possible to find a value for µ tan β such that (g − 2) µ is fitted to the central value of the measurement as given in Eq. (1.1). We therefore fix µ tan β at this value, instead of the maximum value allowed by the vacuum stability. The Bino mass parameter, M 1 , is taken at 20 GeV below the smuon mass, mμ 1 to avoid the appearence of a charged LSP and to relax the LHC constraints with compressed mass spectra. Fig. 5 shows the DM and LHC constraints on the new planes with the mass splitting assumption, mτ R ,ml 3 m l L ,m l R . As mentioned above, the central value of the (g −2) µ measurement is realised everywhere across the planes. We see that the regions withm l L 220 GeV (for tan β = 50) and m l L 300 GeV (for tan β = 10) are excluded by the DM overproduction. This bound may be relaxed if one adopts even more stringent compressed mass assumption, mμ 1 − M 1 < 20 GeV, so that the coannihilation channel becomes more efficient. We also see that the constraint from the CMS + − analysis (red shaded regions) excludesm l L 300 GeV and 200 m l R /GeV 600. This bound may also be relaxed by taking more compressed mass hierarchy, mμ 1 − M 1 < 20 GeV. We see that the central value of the (g − 2) µ measurement can be fitted in a large region of the parameter space (withm l R 600 GeV,m l L 220 (300) GeV for tan β = 50 (10)) consistently with the other phenomenological constraints if the universal mass assumption is removed and the stau masses are taken to be much larger than the other slepton masses.
We also comment that for other parameter planes adopting the mass splitting assumption with heavy staus makes the LHC constraint even more stringent, since the decay branching ratios of heavier electroweakinos into the eletron and muon final states relatively enhance due to absence of staus in the low energy spectrum.

Baryonic R-parity violation
In the previous section, we have shown the SUSY (g − 2) µ regions are severely constrained by overproduction of relic neutralinos, DMDD experiments, and LHC searches in the large-E miss T channel. These constraints are a direct consequence of the assumption thatχ 0 1 is the LSP and stable. If χ 0 1 is the LSP but decays into visible SM particles, the limits fromχ 0 1 overproduction and DMDD experiments are entirely avoided and the E miss T signature at the LHC will not be available. On the other hand, the SUSY contribution to (g − 2) µ is basically unaffected unless the new particle or the operator introduced to permit theχ 0 1 decay give large contribution to (g − 2) µ . One way to makeχ 0 1 unstable is to introduce R-parity violation (RPV). Without extending the MSSM particle contents, the most general renormalisable superpotential admits four types of RPV operators (see, e.g., Ref. [119] for a review). Three of them break the lepton number (L) but not the baryon number (B) symmetry, while it is the other way around for the last one. In this section, we study the scenario with the B-breaking RPV operator. Namely, we extend the MSSM superpotential with  This scenario has three advantages. Firstly, the term in Eq. (6.1) does not generate the other three RPV operators radiatively, since they are protected by the L-symmetry. This justifies us to set the other L-violating RPV couplings to zero at all scales. Secondly, protons remain stable because, within the MSSM, proton decay requires both B-and L-violations. Lastly, these operators do not contribute to the (g − 2) µ . The one and two sigma (g − 2) µ regions found in the previous section are unchanged.
With the term in Eq. (6.1), the LSP neutralino can decay asχ 0 1 → u i d j d k andū idjdk via an off-shell squark. If i, j, k = 3, a neutralinoχ 0 1 decays into three light-flavour quarks, and the large-E miss T signature of the SUSY events at the LHC is replaced by light-flavour jets. Since jets are easily produced at the LHC, discrimination of such a signature from the SM background is challenging. Namely, this case provides a scenario that maximally opens up the allowed SUSY (g−2) µ parameter region. For simplicity, we switch off all the RPV couplings but λ 112 in our analysis. However, the result may apply for non-zero λ ijk couplings as far as i, j, k = 3. The decay rate of the neutralino is proportional to |λ 112 /m 2 q | 2 . If all RPV couplings other than λ 112 are switched off, there is no strong phenomenological constraint on λ 112 [120]. One can therefore safely assume that λ 112 is large enough so that the neutralino decay,χ 0 1 → uds (ūds), is prompt and does not yield displaced vertices.
ATLAS and CMS searches targeting RPV scenarios are rather rare and simplified model limits are not widely available. We therefore estimate the LHC constraints by a chain of Monte Carlo simulations outlined in subsection 3.3. Although all 8 and 13 TeV analyses implemented in CheckMATE (see Appendix A for the complete list) are included in our analyses, only a few provided 95 % CL exclusions. Those analyses are listed in Table 3.
The left panel of Fig. 6 shows the WHL µ plane. As in the previous plots, the green and    [4] is depicted with green and yellow bands, corresponding to one and two sigma agreement respectively. Orange shaded region corresponds to a µ bigger than experimental value by more than 2 sigma. Red, green and blue shaded regions are excluded by LHC constraints.
yellow regions correspond to the areas where the predicted (g − 2) µ agrees with the experimental value within the 1 and 2 σ accuracy, respectively. On the other hand, the red, blue and green shaded regions correspond to the 95 % CL exclusions obtained from ATLAS multijet+ [121], CMS multilepton [122] and ATLAS jets+E miss T [123] analyses, respectively. Due to the R-parity violating U DD operator, the LSP neutralino decays into three (anti-)quarks. The constraints from the neutralino overproduction and DM direct detection experiments are absent unlike the stable neutralino case.
In the WHL L plane, shown in the right of Fig. 6, our results in the region with M 2 >m l L (i.e., below the diagonal dotted line) are essentially identical to those of the M 2 > µ region of the WHL µ plane, which is becausem l L = µ + 20 GeV in the both regions and other parameters are fixed in the same way. Meanwhile, in the other half (m l L > M 2 ), an wider region is excluded by the ATLAS multijet+ (red). We note that Wino and Higgsino are substantially mixed in the mass eigenstates in this region. The exclusion is mainly driven by the pair production ofχ heavy , which ends up into soft lepton(s) plus multijet final state, where the leptons are originated from theχ heavy decay via an off-shell W (Z) boson.
The top and middle panels of Fig. 7 display the BHL and BHR planes, respectively. In the BHL µ (top left) and BHR µ (middle left) planes with |µ| > M 1 , the excluded regions from AT-LAS multijet+ (red) and CMS multilepton (blue) analyses are visible. The main SUSY process contributing to the exclusion is production of Higgsino-like states followed by the decay into l 3l3 with l 3 = (τ, ν τ ) andl 3 = (τ ,ν τ ). Thel 3 then decays into l 3 plus jets via an on-shell Bino. The region with the opposite condition, M 1 > |µ|, is not constrained by the current LHC data, since in this region, Higgsino states decay directly into jets via the RPV operator and any leptons and neutrinos (as a source of E miss T ) are unavailable. In the BHL L (top right) and BHR L (middle right), the regions on the right-hand side to the diagonal lines are identical (in the sense that they share exactly the same parameters) to the corresponding regions in the BHL µ (top left) and BHR µ (middle left), respectively. On the other hand, the regions on the left-hand side of the diagonal lines are new and unconstrained. This is again because the Higgsino states decay directly into jets in those regions. Comparing with the BHL and BHR planes in the stable neutralino scenario (see Fig. 3) more (g − 2) µ region is opened up for the RPV case due to the lack of DM direct detection constraints as well as the ATLAS soft- [109] forχ 0 2χ 0 Fig. 7 show the BLR 50 and BLR 10 planes, respectively. We see that the entire (g − 2) µ region is unconstrained for both planes. This provides a dramatic contrast with the stable neutralino case, where the whole (g − 2) µ region is excluded by the CMS + − [107] analysis as well as the neutralino overproduction. In the RPV BLR planes, the dominant SUSY process is the slepton pair production, pp →˜ +˜ − . However, small slepton cross-sections cannot lead to statistically distinguishable signature when the missing transverse energy is replaced by hadronic jets.

Gravitino LSP
The gravitino LSP scenario, often realised in gauge mediated SUSY breaking (GMSB) [124][125][126], is another well-known case in whichχ 0 1 becomes unstable even if it is the lightest among the MSSM particles. In GMSB, SUSY is spontaneously broken by the F -term of a singlet field Φ in a hidden sector. The hidden sector also contains messenger fields (Ψ, Ψ) with non-trivial representations under the SM gauge group. Assuming a direct coupling between the SUSY breaking field and messenger fields in the superpotential, W hid κΦΨΨ, the scalar and fermion masses in the messenger multiplets split and the SUSY breaking is mediated to the MSSM sector radiatively via gauge interaction. In the GMSB, the ratio between the gravitino mass, m G , and the MSSM soft SUSY breaking mass scale, m, is roughly given by where ∼ 0.1 − 0.01 is a loop factor, M P is the Planck scale and M mess is the mass scale of the messenger fields. Since we often take M mess M P , the gravitino LSP is realised.  [4] is depicted with green and yellow bands, corresponding to one and two sigma agreement respectively. Blue and red shaded regions are excluded by LHC constraints. Dark blue shaded area is excluded by LEP stau mass bound [87].
In fact, very light gravitinos have several theoretical motivations. One of them is to relax apparent fine-tuning in the Higgs sector. In GMSB, the radiative correction to the SUSY breaking Higgs mass parameter depends logarithmically on the messenger mass, δm 2 Hu ∝ ln(M mess / m), for a fixed m, and the fine-tuning may be ameliorated by taking M mess ∼ m. For example, taking M mess ∼ (10 − 100) TeV and m ∼ 1 TeV leads to m G ∼ (1 − 10) eV.
Another motivation is in cosmology. Gravitinos tend to be overproduced at the reheating era unless the reheating temperature T R is very low [127]. However, such a low T R is inconsistent with the thermal leptogenesis scenario [128,129], which typically requires T R > 10 9 GeV. This is known as the cosmological gravitino problem. Very light gravitino with m G 100 eV provides an attractive solution. Such light gravitinos are thermalised in the early Universe and the relic density becomes independent of T R and can be smaller than the observed DM density if m G 100 eV. 11 Backed by these observations, we study LHC constraints on the SUSY (g−2) µ solution assuming approximately massless gravitinos. 12 Although in the minimal GMSB model sfermion and gaugino masses are uniquely determined by the messenger scale, they are less constrained in non-minimal models [134][135][136]. Since the aim of this study is to understand how the LHC constraints change from the stable neutralino case to the gravitino LSP scenario, we take the bottom-up approach and treat the MSSM parameters as free. In the first part of this section we keep the same parameter planes defined in section 2.1 and put the approximately massless gravitino at the bottom of the spectrum. We then introduce small modifications to the definition of our planes to accommodate non-neutralino NLSP cases in the second part of the section. Since the gravitino contribution to (g−2) µ is negligible, the preferred (g−2) µ regions on the parameter planes are essentially unchanged from the previous cases. We also assume that the decay of the next-to-the-lightest SUSY particle (NLSP) is prompt and does not yield displaced vertices. As we will see below this assumption is consistent with very light gravitinos with m G 10 eV.

Neutralino NLSP
In our analysis the NLSPχ 0 1 decays into the (almost) massless gravitino together with a neutral boson (γ, Z or h). The partial decay rates are given by [137,138]: where N ij is the element of the neutralino mixing matrix, (c W , s W ) = (cos θ W , sin θ W ), (c β , s β ) = (cos β, sin β) and (7.3)  We see from this expression the neutralino decay is prompt, cτχ0 1 1 mm, for mχ0 1 100 GeV and m G 10 eV.
Branching ratios for the pure neutralino states can readily be computed by taking N 1i = δ 1i for pure-Bino, N 1i = δ 2i for pure-Wino, and N 13 = −N 14 = 1/ √ 2, N 11 = N 22 = 0 for pure-Higgsino. The branching ratios of those pure-ino states in the large tan β regime are listed in Table 4. In the heavyχ 0 In Ref. [139], CMS analysed 13 TeV, 35.9 fb −1 data in the multilepton channel and excluded NLSP Higgsinos with µ 650 GeV together with a massless gravitino for all values of Br( H → Z G). Similar conclusion has also been obtained in more recent CMS + − analysis (137 fb −1 ) [107], where the pure-Higgsino was excluded up to 650 GeV assuming Br( H → Z G) = Br( H → h G) = 0.5. In the same CMS paper the [Z G][Z G] final state is analysed and the limit σ < 3 fb is derived for mχ0 1 500 GeV. For the pure-Wino this limit translates to m W > 780 GeV.
At this point, it is already clear that the preferred (g − 2) µ region in the WHL planes is entirely excluded, since this region requires either µ 650 GeV or M 2 750 GeV (see Fig. 2 or 6). Note also that µ ∼ 650 GeV implies M 2 ∼ 750 GeV in the (g − 2) µ region, and the exclusion around this region is even stronger, since both Higgsino-and Wino-like states are produced and contribute to signal regions.
The same observation can be made for the BHL and BHR planes. In the half of those planes (|µ| < M 1 )χ 0 1 is Higgsino-like and a good (g − 2) µ fit requires mχ0 1 200 GeV, which violates the aforementioned CMS bounds, 650 GeV.
In order to constrain the Bino-like neutralino region, we use the information provided in CMS γ+E miss T analysis [140] (13 TeV, 35.9 fb −1 ). In this analysis, the final state with an energetic photon and large E miss T is selected. The signal region is defined by the following requirements: at least one photon with p T > 180 GeV, jets are separated from photons and the missing transverse energy as ∆R(γ, jet) > 0.5 and ∆φ(p miss T , p jet T ) > 0.3 for jets with p T > 100 GeV, E miss T > 300 GeV, m T ≡ 2p γ 1 T E miss T (1 − cos ∆φ(p γ 1 T , p miss T )) > 300 GeV and S γ T ≡ E miss T + i p γ i T > 600 GeV, where γ i (i = 1, 2, · · · ) are sorted in the descending order of p T . We note that the event selection is rather inclusive and neither jets nor leptons are required. The result of this analysis is interpreted for the simplified model of electroweakinos with Br(χ 0 1 → γ G) = 0.5 and Br(χ 0 1 → Z G) = Br(χ 0 1 → h G) = 0.25, and the cross-section limit is derived as a function of mχ0 1 . In particular, σ < 100 fb is derived for mχ0 1 ∼ 300 GeV. The same result is also interpreted for the squark simplified model,qq → [qχ 0 1 ][qχ 0 1 ] followed byχ 0 1 → γ G, which results in the [γ G][γ G] + jets final state. For mχ0 1 ∼ (100 − 300) GeV the cross-section limit, σ 10 fb, has been derived.
In the Bino-like neutralino region (M 1 < |µ|) of the BHL and BHR planes, a good (g − 2) µ fit requires mχ0 1 200 GeV and |µ| 500 GeV. In this regionχ 0 1 decays almost exclusively to γ G, since the Z G mode is phase-space suppressed. In this region Higgsino productions (with µ 500 GeV) have cross-section 33 fb, and they contribute to the [γ G][γ G] + jets final state via H → h B, B → γ G. This exceeds the aforementioned bound, σ 10 fb, for [γ G][γ G] + jets, and we conclude the SUSY (g − 2) µ solution in the BHL and BHR planes is incompatible with LHC constraints.
In the BLR planes the NLSP neutralino is almost pure-Bino. We see in the right panels of Fig. 4 that a good (g − 2) µ fit requires M 1 mχ0 1 250 GeV. In this region the cross-section bound, σ 10 fb, may apply to the inclusive SUSY processes giving rise to the [γ G][γ G] + jets final state.
The dominant SUSY processes in the BLR planes are found to be the left-handed slepton/sneutrino and the right-handed slepton productions, and the bound σ 10 fb translates tom l L 400 GeV andm l R 200 GeV. We see that these mass bounds exclude the whole 1 σ (g − 2) µ region in the BLR 50 plane. At the vicinity of the LEP stau mass bound, a small part of the 2 σ region in the BLR 50 and BLR 10 planes as well as of the 1 σ region in the BLR 10 plane is still allowed. More detailed analysis for the NLSP neutralino scenario is provided in Appendix B.

Slepton/Stau NLSP
In the previous subsection, we showed that the SUSY (g − 2) µ solution is incompatible with the LHC constraints in the gravitino LSP scenario when the NLSP isχ 0 1 . In this case the neutralino decays into the massless gravitino and a neutral gauge boson,χ 0 1 → G + γ/Z/h, and the final states must contain two energetic neutral bosons plus large E miss T . Such final states are much easier to be detected compared to the stable neutralino case studied in section 5. 13 In this subsection, we study the gravitino LSP scenario where the NLSP is given by the slepton, sneutrino or stau. To do this systematically, we take the WHL µ , BHL µ and BHR µ planes and simply put the slepton soft mass parameter 20 GeV below the smallest of the gaugino and Higgsino masses. Namely we set In this setup, the NLSP turns out to be sneutrinos in WHL µ and BHL µ , since the electroweak symmetry breaking introduces a small mass splitting in the left-handed slepton doublet in such a way that the sneutrino becomes lighter than the charged slepton. In BHR µ the NLSP is given by three mass-degenerate right-handed charged sleptons,ẽ R ,μ R andτ R . In BLR 50/10 , on the other hand, the NLSP becomes the lighter stau,τ 1 , in whichτ R andτ L are largely mixed. 13 We designed the parameter planes such that the two smallest mass parameters are 20 GeV apart from each other. When lighter particles are pair produced and decay into the LSP neutralinos, the momentum direction ofχ 0 1 tends not to alter significantly from the original direction of the mother particle due to the small mass difference. Consequently the two neutralinos in such an event are almost back-to-back in the transverse plane and cancel the net missing transverse energy. This cancellation is, on the other hand, lost in the gravitino LSP scenario, since the mass difference betweenχ 0 1 and the massless G cannot be small and the decayχ 0 1 → G + γ/Z/h is energetic.  Since the LHC signature of these models is unconventional and simplified model limits are not available, we estimate the LHC constraints using the MC chain outlined in subsection 3.3. Although we include in our analysis all of 8 and 13 TeV analyses listed in the tables in Appendix A, only a few give the 95 % CL exclusion. We list those relevant ATLAS and CMS analyses in Table 5.
We first look at the WHL µ plane shown in the top left panel of Fig. 8. In the plot the blue and orange shaded regions are excluded at 95 % CL by the CMS multilepton [122] and CMS soft + − [141] analyses, respectively. The CMS multilepton analysis requires the signal electrons (muons) to have p T > 25 (20) GeV. This analysis is therefore sensitive to production of heavier electroweakino states,χ heavy , since relatively hard leptons may be produced from the decay of χ heavy . On the other hand, CMS soft + − exploits soft leptons with 5 < p T /GeV < 30. Therefore, the analysis is sensitive even for production of lighter electroweakino states. In the orange shaded region in the plot the production of Wino pairs followed by W → ll ( ) (l ( ) = (e, µ, ν e , ν µ )) mainly contributes to the exclusion. We find that the charged sleptons then decay into the sneutrino and a ff pair via an off-shell W ± . The decay of NLSP sneutrino is totally invisible,ν → ν G, and the NLSP sneutrinos may effectively be treated as the invisible LSP. In the µ < M 2 region, Higgsino production cannot give a large contribution to the CMS soft + − analysis, since the Higgsinos predominantly decay into the third generation sleptons; H → l 3l ( ) 3 (l ( ) 3 = (τ, ν τ )). Due to this almost entire 1 σ (g − 2) µ region remains unconstrained in the µ M 2 region.
The top right panel of Fig. 8 displays the BHL µ plane. We see that the region withm l L M 1 < µ is excluded by the ATLAS τ + τ − (magenta). The excluded signature comes from the Higgsino pair production followed by H 0 → ττ L → τ τ G and H ± → τν L . In the bottom left region with m l L M 1 < µ, we also see the exclusion from the ATLAS soft-(purple) and the CMS multilepton (blue) analyses. The direct slepton/sneutrino production, pp →ll, followed byl → l G (with l = or ν) is responsible for these exclusions. We also see the ATLAS soft-exclusion is slightly extended to the µ M 1 region. There, the NLSP is given byχ 0 1 due to large mixing between Higgsino and Bino and sleptons decay intol → l +χ 0 1 , producing soft leptons. However, in general, the region with µ < M 1 is rather unconstrained since all SUSY particles with large production cross-sections (Higgsinos and sleptons/sneutrinos) are mass degenerate.
We turn to the BHR µ plane shown in the middle panel of Fig. 8. Two excluded regions are visible; ATLAS τ + τ − (magenta) and CMS multilepton (blue). The ATLAS τ + τ − exclusion appears for the same reason as in the BHL µ plane. It excludes the production of Higgsino pair followed by H → ττ R up to |µ| ∼ 520 GeV provided a large mass difference between H andτ R . In this plane energetic taus can also be produced from the NLSP stau decay,τ L → τ + G. We also see that CMS multilepton analysis excludes some region with 300 GeV |µ| M 1 . In this region, Bino and Higgsino are largely mixed and leptons are produced from the decay of EW bosons originated from heavier electroweakino decays,χ heavy →χ light + X (X = W ± , Z, h). As can be seen in the  , BHR µ (middle), BLR 50 (lower left) and BLR 10 (lower right) planes. Region of parameter space allowed by the latest a µ experimental results [4] is depicted with green and yellow bands, corresponding to one and two sigma agreement respectively. Blue, orange, purple and magenta shaded regions are excluded by LHC constraints. Dark blue shaded area is excluded by LEP stau mass bound [87].
Finally, we show the BLR 50 (left) and BLR 10 (right) planes in the bottom panels of Fig. 8. In the BLR scenario,τ 1 becomes significantly lighter than the other sleptons due to the large L-R mixing. Since M 1 is fixed at 20 GeV above mτ 1 , there is a mass gap between slepton and Bino. In the region excluded by the CMS multilepton (blue), the slepton production followed by˜ ± → ± B contributes, where the Bino further decays into the NLPτ 1 and a soft-τ , thenτ 1 → τ + G.
In BLR 50 , the region excluded by ATLAS τ + τ − (magenta) is also visible, where the Higgsino production followed by H → l 3l3 andτ 1 → τ + G contributes. Together with the CMS multilepton constraint, almost all 1 σ (g − 2) µ region is excluded In the BLR 10 plane, on the other hand, the exclusion from ATLAS τ + τ − is absent. This is because the Higgsinos in this plane are about five times heavier than those in BLR 50 as discussed in Section 5 (see also the right panels of Fig. 4). Because of this and the fact that the 1 σ (g − 2) µ region is more extended to largerm l L , the 1 σ (g − 2) µ region is allowed form l R <m l L in the BLR 10 plane.
We have seen that several regions (i.e. the magenta regions in BHL µ , BHR µ and BLR 50 in Fig. 8) in the slepton NLSP scenario studied in this section are excluded by the ATLAS τ + τ − analysis. We briefly comment what we expect for these regions if staus are artificially decoupled as we considered for the stable neutralino case at the end of subsection 5.2. In the BHL µ plane, the magenta region was excluded due to the Higgsino pair production followed by H 0(±) → ττ L (ν τ ). Ifτ L andν τ are decoupled, the Higgsinos predominantly decay into V +χ 0 1 (V = W, Z, h) if these modes are kinematically allowed. Otherwise, the Higgsinos decay into ll ( ) (l ( ) = e, µ, ν e , ν µ ). The latter case is severely constrained by the analyses targeting the leptonic final states [122,141,142]. In the former case, the constraint is milder and the region will be allowed as long as the Higgsino-Bino mass splitting is smaller than ∼150 GeV [107]. In the BHR µ plane, the NLSPs will beẽ R and µ R after decouplingτ R . Hard leptons are produced when they decay into the LSP gravitino and even larger region will be excluded ifτ R is decoupled. In the BLR 50 plane, if staus are decoupled, the vacuum stability constraints can easily be satisfied even for very large µ and µ can be fixed such that the SUSY (g − 2) µ fits the experimental central value. However, in this case, the NLSPs are given by light flavour sleptons. Their production and decay into the massless gravitino contributes to the final state with hard leptons plus large E miss T , which is strongly constrained by the CMS + − analysis [107]. Overall, we do not expect in general more regions to open up for SUSY (g − 2) µ by decoupling staus in the slepton NLSP scenario discussed in this subsection.

Conclusion
We have investigated in this article supersymmetric solutions to the muon (g − 2) anomaly with and without the stable neutralino. We have organised our analysis based on the four types of 1-loop diagrams generating the SUSY (g − 2) µ contribution, denoted by WHL, BHL, BHR and BLR, as shown in Fig. 1.
For significant contribution the first three types require large gaugino-Higgsino mixing in the neutralino state. However, such a neutralino is severely constrained by the DM direct detection experiments. On the other hand, when the BLR contribution is dominant, the neutralino is almost pure Bino, which suffers from the problem of neutralino overproduction. In order to fit the observed (g − 2) µ anomaly, some of the electroweakly interacting sparticles must be very light, O(500) GeV. Those light particles would be produced at the LHC and give rise to the signature with large E miss T .
In section 5 we have demonstrated these points numerically. We showed that the preferred (g − 2) µ region is in fact severely constrained by the combination of constraints from the neutralino overproduction, DM direct detection and direct BSM searches at the LHC, assuming the neutralino is stable. We nevertheless found still-allowed 1 σ (g −2) µ regions in the WHL planes with (M 2 , µ, m l L ) ∼ (1500, 300, 320) GeV and (350, 600, 370) GeV (Fig. 2) and the BHR planes with (M 1 , µ, m l R ) ∼ (300, −120, 150) GeV (Fig. 3), though these regions will be explored in the next generation DM direct detection experiments. Interestingly the whole (g − 2) µ region is excluded in the BLR planes by the LHC constraint in the dilepton plus large E miss T channel (Fig. 4). Moving to supersymmetric scenarios without stable neutralino, the DM constraints can be trivially avoided. The LHC constraints may also be relaxed though the details of the LHC signature and constraints depend on the exact scenarios.
In section 6 we studied the scenario where the neutralino decays into three (anti-)quarks via the R-parity violating U DD operator. We have shown in Figs. 6 and 7 that although the large E miss T signature is unavailable in this scenario, some regions are excluded by the LHC searches exploiting the multijet and multilepton channels. Nevertheless, the allowed (g − 2) µ region is significantly extended compared to the stable neutralino case. In particular, the BLR solution, which is excluded in the stable neutralino case, is entirely allowed in the RPV scenario.
Another well-known example where the neutralino is unstable is the scenario with gravitino LSP, which we studied in section 7. Adding the approximately massless gravitino in the spectrum, the NLSP neutralino decays into the gravitino in association with one of the neutral boson, γ/Z/h. There are several LHC analyses targeting these final states. Using the cross-section limits derived in these analyses, we showed that the SUSY (g − 2) µ solution in our 2D parameter planes is incompatible with the LHC constraints. A more detailed analysis is provided in Appendix B. This is an interesting example showing that making the neutralino unstable does not always relax the phenomenological constraints for (g − 2) µ .
In subsection 7.1 we studied a gravitino LSP scenario where the NLSP is provided by˜ /ν/τ 1 . In order to study these cases systematically we simply set the soft breaking slepton mass 20 GeV below the minimum of the gaugino and the Higgsino mass in the WHL µ , BHL µ and BHR µ planes. For the BLR 50/10 planes we set M 1 20 GeV above the lighter stau mass, mτ 1 . In this setup, the NSLP is found to be sneutrino in WHL µ and BHL µ , right-handed charged sleptons in BHR µ and the lighter stau in BLR 50/10 . As shown in Fig. 8 some regions are constrained by the LHC searches exploiting the multilepton, soft-and di-τ channels. Unlike the neutralino NLSP case, we found several 1 σ (g − 2) µ regions particularly in the WHL µ , BHR µ and BLR 10 planes.
For unstable neutralinos, the region favoured by the (g − 2) µ anomaly can only be probed by collider experiments. If the observed (g − 2) µ anomaly is indeed a sign of supersymmetry, it is very important to understand the collider signature and develop a methodology to find it for both stable and unstable neutralino scenarios. We however leave this task as an interesting future work.

Appendix A CheckMATE analyses
In the tables below, we list the ATLAS and CMS analyses used in our numerical analyses with CheckMATE. Tables 6, 7, 8, and 9 list the 13 TeV ATLAS, 13 TeV CMS, 8 TeV ATLAS and 8 TeV       allowed by the latest a µ experimental results [4] is depicted with green and yellow bands, corresponding to one and two sigma agreement respectively. Blue, green, red and magenta shaded regions are excluded by LHC constraints. Dark blue shaded area is excluded by LEP stau mass bound [87].