Likelihood Analysis of the Sub-GUT MSSM in Light of LHC 13-TeV Data

We describe a likelihood analysis using MasterCode of variants of the MSSM in which the soft supersymmetry-breaking parameters are assumed to have universal values at some scale $M_{in}$ below the supersymmetric grand unification scale $M_{GUT}$, as can occur in mirage mediation and other models. In addition to $M_{in}$, such `sub-GUT' models have the 4 parameters of the CMSSM, namely a common gaugino mass $m_{1/2}$, a common soft supersymmetry-breaking scalar mass $m_0$, a common trilinear mixing parameter $A$ and the ratio of MSSM Higgs vevs $\tan\beta$, assuming that the Higgs mixing parameter $\mu>0$. We take into account constraints on strongly- and electroweakly-interacting sparticles from $\sim 36$/fb of LHC data at 13 TeV and the LUX and 2017 PICO, XENON1T and PandaX-II searches for dark matter scattering, in addition to the previous LHC and dark matter constraints as well as full sets of flavour and electroweak constraints. We find a preference for $M_{in} \sim 10^5$ to $10^9$ GeV, with $M_{in} \sim M_{GUT}$ disfavoured by $\Delta \chi^2 \sim 3$ due to the ${\rm BR}(B_{s, d} \to \mu^+\mu^-)$ constraint. The lower limits on strongly-interacting sparticles are largely determined by LHC searches, and similar to those in the CMSSM. We find a preference for the LSP to be a Bino or Higgsino with $\tilde{\chi^0_1} \sim 1$ TeV, with annihilation via heavy Higgs bosons $H/A$ and stop coannihilation, or chargino coannihilation, bringing the cold dark matter density into the cosmological range. We find that spin-independent dark matter scattering is likely to be within reach of the planned LUX-Zeplin and XENONnT experiments. We probe the impact of the $(g-2)_\mu$ constraint, finding similar results whether or not it is included.

We describe a likelihood analysis using MasterCode of variants of the MSSM in which the soft supersymmetrybreaking parameters are assumed to have universal values at some scale Min below the supersymmetric grand unification scale MGUT, as can occur in mirage mediation and other models. In addition to Min, such 'sub-GUT' models have the 4 parameters of the CMSSM, namely a common gaugino mass m 1/2 , a common soft supersymmetry-breaking scalar mass m0, a common trilinear mixing parameter A and the ratio of MSSM Higgs vevs tan β, assuming that the Higgs mixing parameter µ > 0. We take into account constraints on stronglyand electroweakly-interacting sparticles from ∼ 36/fb of LHC data at 13 TeV and the LUX and 2017 PICO, XENON1T and PandaX-II searches for dark matter scattering, in addition to the previous LHC and dark matter constraints as well as full sets of flavour and electroweak constraints. We find a preference for Min ∼ 10 5 to 10 9 GeV, with Min ∼ MGUT disfavoured by ∆χ 2 ∼ 3 due to the BR(B s,d → µ + µ − ) constraint. The lower limits on strongly-interacting sparticles are largely determined by LHC searches, and similar to those in the CMSSM. We find a preference for the LSP to be a Bino or Higgsino with mχ0

Introduction
Models invoking the appearance of supersymmetry (SUSY) at the TeV scale are being sorely tested by the negative results of high-sensitivity searches for sparticles at the LHC [1,2] and for the scattering of dark matter particles [3][4][5][6]. There have been many global analyses of the implications of these experiments for specific SUSY models, mainly within the minimal supersymmetric extension of the Standard Model (MSSM), in which the lightest supersymmetric particle (LSP) is stable and a candidate for dark matter (DM). This may well be the lightest neutralino, χ 0 1 [7], as we assume here. Some of these studies have assumed universality of the soft SUSYbreaking parameters at the GUT scale, e.g., in the constrained MSSM (the CMSSM) [8][9][10][11] and in models with non-universal Higgs masses (the NUHM1,2) [9,12]. Other analyses have taken a phenomenological approach, allowing free variation in the soft SUSY-breaking parameters at the electroweak scale (the pMSSM) [13][14][15][16].
A key issue in the understanding of the implications of the LHC searches for SUSY is the exploration of regions of parameter space where compressed spectra may reduce the sensitivity of searches for missing transverse energy, / E T . These regions also have relevance to cosmology, since models with sparticles that are nearly degenerate with the LSP allow for important coannihilation processes that suppress the relic LSP number density, allowing heavier values of mχ0 1 . The accompanying heavier SUSY spectra are also more challenging for the LHC / E T searches. The CMSSM offers limited prospects for coannihilation, and examples that have been studied in some detail include coannihilation with the lighter stau slepton,τ 1 [17,18], or the lighter stop squark,t 1 [19]. Other models offer the possibilities of different coannihilation partners, such as the lighter chargino,χ ± 1 [14,20], some other slepton [16] or squark flavour [21], or the gluino [22,23]. In particular, the pMSSM allows for all these possibilities, potentially also in combination [16].
In this paper we study the implications of LHC and DM searches for an intermediate class of SUSY models, in which universality of the soft SUSY-breaking parameters is imposed at some input scale M in below the GUT scale M GUT but above the electroweak scale [24,25], which we term 'sub-GUT' models. Models in this class are well motivated theoretically, since the soft SUSYbreaking parameters in the visible sector may be induced by some dynamical mechanism such as gluino condensation that kicks in below the GUT scale. Specific examples of sub-GUT models include warped extra dimensions [26] and mirage mediation [27].
Mirage mediation can occur when two sources of supersymmetry breaking play off each other, such as moduli mediation based, e.g., on moduli stabilization as in [28] and anomaly mediation [29].
The relative contributions of each source of supersymmetry breaking can be parametrized by the strength of the moduli mediation, α, and allows one to interpolate between nearly pure moduli mediation (large α) and nearly pure anomaly mediation (α → 0). For example, gaugino masses, M i , can be written as where M s is related to the gravitino mass in anomaly mediation (m 3/2 = 16π 2 M s ), and b i , g i are the beta functions and gauge couplings. This leads to a renormalization scale, M in = M GU T e −8π 2 /α at which gaugino masses and soft scalar masses take unified values, although there is no physical threshold at M in in this model. We are not concerned here with the detailed origin of M in , simply postulating that there is a scale below the GUT scale where the supersymmetry breaking masses are unified.
Sub-GUT models are of particular phenomenological interest, since the reduction in the amount of renormalization-group (RG) running below M in , compared to that below M GUT in the CMSSM and related models, leads naturally to SUSY spectra that are more compressed [24]. These may offer extended possibilities for 'hiding' SUSY via suppressed / E T signatures, as well as offering enhanced possibilities for different coannihilation processes. Other possible effects of the reduced RG running include a stronger lower limit on mχ0 1 because of the smaller hierarchy with the gluino mass, a stronger lower limit on the DM scattering cross section because of a smaller 3 hierarchy between mχ0 1 and the squark masses, and greater tension between LHC searches and a possible SUSY explanation of the measurement of (g − 2) µ [30,31], because of the smaller hierarchies between the gluino and squark masses and the smuon andχ 0 1 masses. We use the MasterCode framework [8,9,12,14,16,21,[32][33][34][35] to study these issues in the sub-GUT generalization of the CMSSM, which has 5 free parameters, comprising M in as well as a common gaugino mass m 1/2 , a common soft SUSY-breaking scalar mass m 0 , a common trilinear mixing parameter A and the ratio of MSSM Higgs vevs tan β, assuming that the Higgs mixing parameter µ > 0, as may be suggested by (g − 2) µ 1 . Our global analysis takes into account the relevant CMS searches for strongly-and electroweakly-interacting sparticles with the full 2016 sample of ∼ 36/fb of data at 13 TeV [36][37][38], and also considers the available results of searches for long-lived charged particles [39,40] 2 . We also include a complete set of direct DM searches published in 2017, including the PICO limit on the spin-dependent scattering cross section, σ SD p [4], as well as the first XENON1T limit [5] and the most recent PandaX-II limit [6] on the spinindependent scattering cross section, σ SI p , as well as the previous LUX search [3]. We also include full sets of relevant electroweak and flavour constraints.
We find in our global sub-GUT analysis a distinct preference for M W M in M GUT , with values of M in ∼ 10 5 or ∼ 10 8 to 10 9 GeV being preferred by ∆χ 2 ∼ 3 compared to the CMSSM (where M in = M GUT ). This preference is driven principally by the ability of the sub-GUT MSSM to accommodate a value of BR(B s,d → µ + µ − ) smaller than in the Standard Model (SM), as preferred by the current data [41][42][43]. As discussed later, this effect can be traced to the different RGE evolution of A t in the sub-GUT model, which enables it have a different sign from that in the CMSSM. The lower limits on 1 We have also made an exploratory study for µ < 0 with a limited sample, finding quite similar results within the statistical uncertainties. 2 The ATLAS SUSY searches with ∼ 36/fb of data at 13 TeV [2] yield similar constraints. strongly-interacting sparticles are similar to those in the CMSSM, being largely determined by LHC searches. The favoured DM scenario is that the LSP is a Bino or Higgsino with mχ0 1 ∼ 1 TeV, with the cold DM being brought into the cosmological range by annihilation via heavy Higgs bosons H/A and stop coannihilation, or chargino coannihilation. In contrast to the CMSSM and pMSSM11, the possibility that mχ0 1 1 TeV is strongly disfavoured in the sub-GUT model, so the LHC constraints have insignificant impact. The same is true of the LHC searches for longlived charged particles.
The likelihood functions for fits with and without the (g − 2) µ constraint are quite similar, reflecting the anticipated difficulty in accounting for the (g − 2) µ anomaly in the sub-GUT MSSM. Encouragingly, we find a preference for a range of σ SI p just below the current upper limits, and within the prospective sensitivities of the LUX-Zeplin (LZ) [44] and XENONnT [45] experiments.
The outline of this paper is as follows. In Section 2 we summarize the experimental and astrophysical constraints we apply. Since we follow exactly our treatments in [16], we refer the interested reader there for details. Then, in Section 3 we summarize the MasterCode framework and how we apply it to the sub-GUT models. Our results are presented in Section 4. Finally, Section 5 summarizes our conclusions and discusses future perspectives for the sub-GUT MSSM.

Electroweak and Flavour Constraints
Our treatments of these constraints are identical to those in [16], which were based on Table 1 of [21] with the updates listed in Table 2 of [16]. Since we pay particular attention in this paper to the impact on the sub-GUT parameter space of the (g − 2) µ constraint [30], we note that we assume (1) to be the possible discrepancy with SM calculations [31] that may be explained by SUSY. As we shall see, the BR(B s,d → µ + µ − ) measure-ment [41][42][43] plays an important role in indicating a preferred region of the sub-GUT parameter space.

Higgs Constraints
In the absence of published results on the Higgs boson based on Run 2 data, we use in this global fit the published results from Run 1 [46], as incorporated in the HiggsSignals code [47].
Searches for heavy MSSM Higgs bosons are incorporated using the HiggsBounds code [48], which uses the results from Run 1 of the LHC. We also include the ATLAS limit from ∼ 36/fb of data from the LHC at 13 TeV [49].

Dark Matter Constraints and Mechanisms Cosmological density
Since R-parity is conserved in the MSSM, the LSP is a candidate to provide the cold DM (CDM). We assume that the LSP is the lightest neutralinoχ 0 1 [7], and that it dominates the total CDM density. For the latter we assume the Planck 2015 value: Ω CDM h 2 = 0.1186 ± 0.0020 EXP ± 0.0024 TH [50].

Density mechanisms
As in [16], we use the following set of measures related to particle masses to indicate when specific mechanisms are important for bringing Ω CDM h 2 into the Planck 2015 range, which have been validated by checks using Micromegas [51].

•Chargino coannihilation
This may be important if theχ 0 1 is not much lighter than the lighter chargino,χ ± 1 , and we introduce the following coannihilation measure: chargino coann. : mχ± We shade green in the 2-dimensional plots in Section 4 the parts of the 68 and 95% CL regions where (2) is satisfied.
•Rapid annihilation via direct-channel H/A poles We find that LSP annihilation is enhanced sig-nificantly if the following condition is satisfied: and shade in blue the parts of the 68 and 95% CL regions of the two-dimensional plots in Section 4 where (3) is satisfied.

•Stau coannihilation
We introduce the following measure for stau coannihilation: and shade in pink the corresponding area of the 68 and 95% CL regions of the two-dimensional sub-GUT parameter planes. We do not find regions where coannihilation with other charged slepton species, or with sneutrinos, is important.

•Stop coannihilation
We introduce the following measure for stop coannihilation: and shade in yellow the corresponding area of the 68 and 95% CL regions of the two-dimensional sub-GUT parameter planes. We do not find regions where coannihilation with other squark species, or with gluinos, is important.
•Focus-point region The sub-GUT parameter space has a focuspoint region where the DM annihilation rate is enhanced because the LSPχ 0 1 has an enhanced Higgsino component as a result of near-degeneracy in the neutralino mass matrix. We introduce the following measure to characterize this possibility: and shade in cyan the corresponding area of the 68 and 95% CL regions of the two-dimensional 5 sub-GUT parameter planes.

•Hybrid regions
In addition to regions where one of the above DM mechanisms is dominant, there are also various 'hybrid' regions where more than one mechanism is important. These are indicated in the two-dimensional planes below by shadings in mixtures of the 'primary' colours above, which are shown in the corresponding figure legends. For example, there are prominent regions where both chargino coannihilation and directchannel H/A poles are important, whose shading is darker than the blue of regions where H/A poles are dominant.

Direct DM searches
We apply the constraints from direct searches for weakly-interacting dark matter particles via both spin-independent and -dependent scattering on nuclei.
Indirect astrophysical searches for DM As discussed in [16], there are considerable uncertainties in the use of IceCube data [56] to constrain σ SD p and, as we discuss below, the global fit yields a prediction that lies well below the current PICO [4] constraint on σ SD p and the 3 We note that a recent analysis using covariant baryon chiral perturbation theory yields a very similar central value of Σ πN [54]. However, we emphasize that there are still considerable uncertainties in the estimates of σ 0 and Σ πN and hence the N |ss|N matrix element that is important for σ SI p [55].
current IceCube sensitivity, so we do not include the IceCube data in our global fit.

Stop and sbottom searches
We also implement the CMS simplified model searches with ∼ 36/fb of data at 13 TeV in the jets + 0 [36] and 1 [37] [16].

Searches for electroweak inos
We also consider the CMS searches for electroweak inos in multilepton final states with ∼ 36/fb of data at 13 TeV [38], constraining χ ± 1χ 0 2 → [Wχ 0 1 ][Zχ 0 1 ], 3 ± + 2χ 0 1 via˜ ± /ν intermediate states, and 3τ ± + 2χ 0 1 viaτ ± intermediate states using Fastlim [57] as described in [16]. These analyses can also be used to constrain the production of electroweak inos in the decays of coloured sparticles, since these searches do not impose conditions on the number of jets. However, as we discuss below, in the sub-GUT model the above-mentioned searches for stronglyinteracting sparticles impose such strong limits on the mχ0 1 and mχ± 1 that the searches for electroweak inos do not have significant impact on the preferred parameter regions.

Searches for long-lived or stable charged particles
We also consider a posteriori the search for long-lived charged particles published in [39], which are sensitive to lifetimes ns, and the search for massive charged particles that escape from the detector without decaying [40]. How-ever, these also do not have significant impact on the preferred parameter regions, as we discuss in detail below, and are not included in our global fit.

Model Parameters
As mentioned above, the five-dimensional sub-GUT MSSM parameter space we consider in this paper comprises a gaugino mass parameter m 1/2 , a soft SUSY-breaking scalar mass parameter m 0 and a trilinear soft SUSY-breaking parameter A 0 that are assumed to be universal at some input mass scale M in , and the ratio of MSSM Higgs vevs, tan β. Table 1 displays the ranges of these parameters sampled in our analysis, as well as their divisions into segments, which define boxes in the five-dimensional parameter space.

Parameter
Range # of segments M in (10 3 , 10 16 ) GeV 6 m 1/2 (0, 6) TeV 2 m 0 (0, 6) TeV 2 A 0 (−15, 10) TeV 2 tan β ( 1 , 60) 2 Total # of boxes 96 Table 1 The ranges of the sub-GUT MSSM parameters sampled, together with the numbers of segments into which they are divided, together with the total number of sample boxes shown in the last row. This sample is for positive values of the Higgs mixing parameter, µ. As already noted, a smaller sample for µ < 0 gives similar results. Note that our sign convention for A is opposite to that used in SoftSusy [58].

Sampling Procedure
We sample the boxes in the five-dimensional sub-GUT MSSM parameter space using the MultiNest package [59], choosing for each box a prior such that 80% of the sample has a flat distribution within the nominal box, and 20% of the sample is in normally-distributed tails extending outside the box. This eliminates features associated with the boundaries of the 96 boxes, by providing a smooth overlap between them. In total, our sample includes ∼ 112 million points with ∆χ 2 < 100.

The MasterCode
The MasterCode framework [8,9,12,14,16,21,[32][33][34][35], interfaces and combines consistently various private and public codes using the SUSY Les Houches Accord (SLHA) [60]. This analysis uses the following codes: SoftSusy 3.7.2 [58] for the MSSM spectrum, FeynWZ [61] for the electroweak precision observables, SuFla [62] and SuperIso [63]  The top left panel of Fig. 1 displays the onedimensional profile χ 2 likelihood function for M in , as obtained under various assumptions 4 . In this and subsequent one-dimensional plots, the solid lines represent the results of a fit including results from ∼ 36/fb of data from the LHC at 13 TeV (LHC13), whereas the dashed lines omit these results, and the blue lines include (g − 2) µ , whereas the green lines are obtained when this constraint is dropped.
We observe in the top left panel of Fig sub-GUT sub-GUT sub-GUT Here and in subsequent one-dimensional plots, the solid lines include the constraints from ∼ 36/fb of LHC data at 13 TeV and the dashed lines drop them, and the blue lines include (g − 2) µ , whereas the green lines drop these constraints. Here and in subsequent two-dimensional plots, the red (blue) (green) contours are boundaries of the 1-, 2-and 3-σ regions, and the shadings correspond to the DM mechanisms indicated in the legend. 13-TeV data and (g −2) µ are both included (solid blue line), falling to 5.9×10 5 GeV when the 13-TeV data are dropped (dashed blue line). There is little difference between the global χ 2 values at these two minima, but values of M in < 10 5 GeV are strongly disfavoured. The rise in ∆χ 2 when M in increases to ∼ 10 6 GeV and the LHC 13-TeV data are included (solid lines) is largely due to the contribution of BR(B s,d → µ + µ − ). At lower M in , the H → τ + τ − constraint allows a larger value of tan β, which leads (together with an increase in the magnitude of A) to greater negative interference in the supersymmetric contribution to BR(B s,d → µ + µ − ), as preferred by the data.
For both fits including the LHC 13-TeV data (solid lines), the ∆χ 2 function ∼ 1 for most of the range M in ∈ (10 5 , 10 11 ) GeV, apart from localized dips, whereas ∆χ 2 rises to 2 for M in 10 12 GeV. As already mentioned and discussed in more detail later, the reduction in the global χ 2 function for M in 10 12 GeV arises because for these values of M in the sub-GUT model can accommodate better the measurement of BR(B s,d → µ + µ − ), whose central experimental value is somewhat lower than in the SM.
When the (g − 2) µ constraint is dropped, as shown by the green lines in top left panel of Fig. 1, there is a minimum of χ 2 around M in 1.6 × 10 5 GeV, whether the LHC 13-TeV constraint is included, or not. The values of the other input parameters at the best-fit points with and without these data are also very similar, as are the values of ∆χ 2 . On the other hand, the values of ∆χ 2 for M in ∈ (10 5 , 10 8 ) GeV are generally smaller when the LHC 13-TeV constraints are dropped, the principal effect being due to the H/A → τ + τ − constraint.
In contrast, when M in 10 9 GeV the ∆χ 2 function in the top left panel of Fig. 1 is quite similar whether the LHC 13-TeV and (g − 2) µ constraints are included or not, though ∆χ 2 0.5 lower when the (g − 2) µ constraint is dropped, as seen by comparing the green and blue lines. This is because the tension between (g − 2) µ and LHC data is increased when M 3 /M 1 is reduced, as occurs because of the smaller RGE running when M in < M GUT . Conversely, lower M in is relatively more favoured when (g − 2) µ is dropped, leading to this increase in ∆χ 2 at high M in though the total χ 2 is reduced.
We list in Table 2 the parameters of the best-fit points when we drop one or both of the (g − 2) µ and LHC13 constraints, as well as the values of the global χ 2 function at the best-fit points. We see that the best-fit points without (g − 2) µ are very similar with and without the LHC 13-TeV constraint. On the other hand, the best-fit points with (g − 2) µ have quite different values of the other input parameters, as well as larger values of M in , particularly when the LHC 13-TeV data are included.
The top right panel of Fig. 1 displays the (m 0 , m 1/2 ) plane when the (g − 2) µ and LHC13 constraints are applied. Here and in subsequent planes, the green star indicates the best-fit point, whose input parameters are listed in Table 2: it lies in a hybrid stop coannihilation and rapid H/A annihilation region.
This parameter plane and others in Fig. 1 and subsequent figures also display the 68% CL (1σ), 95% CL (2-σ) and 99.7% (3-σ) contours in the fit including both (g − 2) µ and the LHC13 data as red, blue and green lines, respectively. We note, here and subsequently, that the green 3-σ contours are generally close to the blue 2-σ contours, indicating a relatively rapid increase in χ 2 , and that the χ 2 function is relatively flat for m 0 , m 1/2 1 TeV. The regions inside the 95% CL contours are colour-coded according to the dominant DM mechanisms, as shown in the legend beneath Fig. 1 5 . Similar results for this and other planes are obtained when either or both of the (g − 2) µ and LHC13 constraints are dropped.
We see that chargino coannihilation is important in the upper part of the (m 0 , m 1/2 ) plane shown in the top right panel of Fig. 1, but rapid annihilation via the H/A bosons becomes important for lower m 1/2 , often hybridized with other mechanisms including stop and stau coannihilation. We also note smaller regions with m 1/2 ∼ 1.5 to 3 TeV where stop coannihilation and focus-point mechanisms are dominant.
The middle left panel of Fig. 1 shows the cor-  Table 2 Values of the sub-GUT input parameters at the best-fit points with and without (g − 2) µ and the LHC 13-TeV data.
responding (M in , m 0 ) plane, where we see a significant positive correlation between the variables that is particularly noticeable in the 68% CL region. In most of this and the 95% CL region with M in 10 13 GeV the relic LSP density is controlled by chargino coannihilation, though with patches where rapid annihilation via the A/H bosons is important, partly in hybrid combinations. In contrast, the (M in , m 1/2 ) plane shown in the middle right panel of Fig. 1 does not exhibit a strong correlation between the variables. We see again the importance of chargino coannihilation, with the A/H mechanism becoming more important for lower m 1/2 and larger M in , and for all values of m 1/2 for M in 10 14 GeV.
Also visible in the middle row of planes are small regions with M in ∼ 10 13 to 10 14 GeV where stau coannihilation is dominant, partly hybridized with stop coannihilation. The reduction in the global χ 2 function for M in 10 12 GeV visible in the top left panel of Fig. 1 is associated with the 68% CL regions in this range of M in visible in the two middle planes of Fig. 1.
The one-dimensional profile likelihood functions for m 0 and m 1/2 are shown in the bottom panels of Fig. 1. We note once again the similarities between the results with/without (g − 2) µ (blue/green lines) and the LHC13 constraints (solid/dashed lines). The flattening of the χ 2 function for m 0 at small values reflects the extension to m 0 = 0 of the 95% CL region in the top right panel of Fig. 1. On the other hand, the χ 2 function for m 1/2 rises rapidly at small values, reflecting the close spacing of the 95 and 99.7% CL contours for m 1/2 ∼ 1 TeV seen in the same plane. The impact of the LHC13 constraints is visible in the differences between the solid and dashed curves at small m 0 , in particular. The (g − 2) µ constraint has less impact, as shown by the smaller differences between the green and blue curves. We see that the χ 2 function for m 0 rises by 1 at large mass values, whereas that for m 1/2 falls monotonically at large values. The χ 2 function for m 1/2 exhibits a local maximum at m 1/2 ∼ 3 TeV, which corresponds to the separation between the two 68% CL regions in the top right plane of Fig. 1. These are dominated by chargino coannihilation (larger m 1/2 , green shading) and by rapid annihilation via A/H bosons (smaller m 1/2 , blue shading) and other mechanisms, respectively.

Squarks and gluinos
The various panels of Fig. 2 show the limited impact of the LHC 13-TeV constraints on the possible masses of strongly-interacting sparticles in the sub-GUT model, comparing the solid and dashed curves. The upper left panel shows that the 95% CL lower limit on mg ∼ 1.5 TeV, whether the LHC 13-TeV data and the (g − 2) µ constraint are included or not. However, the bestfit value of mg increases from ∼ 2 TeV to a very large value when (g − 2) µ is dropped, although the ∆χ 2 price for mg ∼ 2 TeV is ∼ 1. The upper right panel shows similar features in the profile likelihood function for mq R (that for mq L is similar), with a 95% CL lower limit of ∼ 2 TeV, which is again quite independent of the inclusion of (g − 2) µ and the 13-TeV data. The lower panels of Fig. 2 show the corresponding profile likelihood functions for mt 1 (left panel) and mb 1 (right panel). We see that these could both be consid-  sub-GUT sub-GUT sub-GUT erably lighter than the gluino and the first-and second-generation squarks, with 95% CL lower limits mt 1 ∼ 900 GeV and mb 1 ∼ 1.5 TeV, respectively.

The lightest neutralino and lighter chargino
The top left panel of Fig. 3 shows the profile likelihood function for mχ0 1 , and the top right panel shows that for mχ± 1 . We see that in all the cases considered (with and without the (g − 2) µ and LHC13 constraints), the value of ∆χ 2 calculated using the LHC constraints on stronglyinteracting sparticles is larger than 4 for mχ0 1 750 GeV and mχ± 1 800 GeV. Therefore, the LHC electroweakino searches [38] have no impact on the 95% CL regions in our 2-dimensional projections of the sub-GUT parameter space, and we do not include the results of [38] in our global fit.
We now examine the profile likelihood functions for the fractions of Bino, Wino and Higgsino in theχ 0 1 composition: which are shown in Fig. 4. As usual, results from an analysis including the 13-TeV data are shown as solid lines and without them as dashed lines, with (g − 2) µ as blue lines and without it as green   lines. The top left panel shows that in the LHC 13-TeV case with (g − 2) µ an almost pureB composition of theχ 0 1 is preferred, N 11 → 1, though the possibility that this component is almost absent is only very slightly disfavoured. Conversely, before the LHC 13-TeV data there was a very mild preference for N 11 → 0, and this is still the case if (g − 2) µ is dropped. The upper right panel shows that a smallW 3 component in theχ 0 1 is strongly preferred in all cases. Finally, the lower panel confirms that smallH u,d components are preferred when the LHC 13-TeV and (g −2) µ constraints are applied, but largeH u,d components are preferred otherwise.
Theχ 0 1 compositions favoured at the 1-, 2-and 3-σ levels (blue, yellow and red) are displayed in Fig. 5 for fits including LHC 13-TeV data with (without) the (g−2) µ constraint in the left (right) panel. We see that these regions are quite similar in the two panels, and correspond to small Wino admixtures. On the other hand, the Bino fraction N 2 11 and the Higgsino fraction N 2 13 + N 2 14 are relatively unconstrained at the 95% CL. The best-fit points are indicated by green stars, and the left panel shows again that in the fit with (g − 2) µ the LSP is an almost pure Bino, whereas an almost pure Higsino composition is favoured in the fit without (g −2) µ , as also seen in Table 3. These two extremes have very similar χ 2 values in each of the fits displayed.BW  Table 3 Composition of theχ 0 1 LSP at the best-fit points with and without (g − 2) µ and the LHC 13-TeV data.
The global χ 2 function is minimized for mχ0 1 1.0 TeV, which is typical of scenarios with a Higgsino-like LSP whose density is brought into the Planck 2015 range by coannihilation with a nearly-degenerate Higgsino-like charginoχ ± 1 . Indeed, we see in the top right panel of Fig. 3 that χ 2 is minimized when also mχ± 1 mχ0 1 1.0 TeV. Table 3 displays the LSP composition of the sub-GUT model at the best-fit points with and without (g − 2) µ and the LHC 13-TeV data. We see again that theχ 0 1 LSP is mainly a Higgsino with almost equalH u andH d components, except in the fit with both LHC 13-TeV data and (g − 2) µ included, in which case it is an almost pure Bino.
Looking at the middle left panel of Fig. 3, we see that the best-fit point has a chargino-LSP mass difference that may be O(1) GeV or ∼ 200 to 300 GeV, with similar χ 2 in all the cases considered, namely with and without the (g−2) µ and LHC13 constraints. As seen in the middle right panel of Fig. 3, in the more degenerate case the preferred chargino lifetime τχ± 1 ∼ 10 −12 s. The current LHC searches for long-lived charged particles [39] therefore do not impact this chargino coannihilation region, and are also not included in our global fit.
The top right panel of Fig. 3 displays an almostdegenerate local minimum of χ 2 with mχ± 1 ∼ 1.3 TeV, corresponding to a second, local minimum of χ 2 where mχ± 1 − mχ0 1 ∼ 200 to 300 GeV, as seen in the middle left panel. In this region the relic density is brought into the Planck 2015 range by rapid annihilation through A/H bosons, as can be inferred from the bottom left panel of Fig. 3, where we see that at this secondary minimum M A 2 TeV 2mχ0 1 . Theχ ± 1 lifetime in this region is too short to appear in the middle and bottom right panels of Fig. 3, and too short to have a separated vertex signature at the LHC.

Sleptons
The upper left panel of Fig. 6 shows the profile likelihood function for mμ R (that for mẽ R is indistinguishable, theμ L andẽ L are slightly heavier). We see that in the sub-GUT model small values of mμ R were already disfavoured by earlier LHC data (dashed lines), and that this tendency has been reinforced by the LHC 13-TeV data (compare the solid lines). The same is true whether the (g − 2) µ constraint is included or dropped (compare the blue and green curves).
The upper right panel Fig. 6 shows the corresponding profile likelihood function for mτ 1 , which shares many similar features. However, we note that the χ 2 function for mτ 1 is generally lower than that for mμ R ∈ (1, 2) TeV, though the 95% lower limits on mτ 1 and mμ R are quite similar, and both are 1 TeV when the LHC 13-TeV constraints are included in the fit.
The lower left panel of Fig. 6 shows that very small values of mτ 1 − mχ0 1 in the stau coannihilation region are allowed at the ∆χ 2 ∼ 1 level in all the fits with the (g − 2) µ constraint, rising to ∆χ 2 2 for mτ 1 − mχ0 1 20 GeV when the LHC 13-TeV data are included.
The lower right panel of Fig. 6 shows the (mτ 1 , ττ 1 ) plane, where we see that ττ 1 ∈ (10 −7 , 10 3 ) s is allowed at the 68% CL, for 1600 GeV mτ 1 2000 GeV and at the 95% CL also for mτ 1 ∼ 1100 GeV. This region of parameter space is close to the tip of the stau coannihilation strip. Lowerτ 1 masses are strongly disfavoured by the LHC constraints, particularly at 13 TeV, as seen in the upper right panel of Fig. 6. The heavierτ 1 masses with lower ∆χ 2 seen there do not lie in the stau coannihilation strip, and have larger mτ 1 −mχ0 1 and hence smaller lifetimes that are not shown in the lower right panel of Fig. 6. Because of the lower limit on mτ 1 seen in this panel, neither the LHC search for longlived charged particles [39] nor the LHC search for (meta-)stable massive charged particles that exit the detector [40] are relevant for our global fit.
In view of this, and the fact that the search for long-lived particles [39] is also insensitive in the chargino coannihilation region, as discussed above, the results of [39,40] are not included in the calculation of the global likelihood function.

(g − 2) µ
We see in the left panel of Fig. 7 that only a small contribution to (g − 2) µ is possible in sub-GUT models, the profile likelihood functions with and without the LHC 13-TeV data and (g − 2) µ being all quite similar. This is because in the sub-GUT model with low M in the LHC searches for strongly-interacting sparticles constrain thẽ µ mass more strongly than in the GUT-scale CMSSM. The dotted line shows the ∆χ 2 contribution due to our implementation of the (g − 2) µ constraint alone. We see that in all cases it contributes ∆χ 2 9 to the global fit.

The (M A , tan β) Plane
The right panel of Fig. 7 shows the (M A , tan β) plane when the LHC 13-TeV data and the (g−2) µ constraint are included in the fit. We see that M A 1.3 TeV at the 95% CL and that, whereas tan β ∼ 5 is allowed at the 95% CL. Larger values tan β 30 are favoured at the 68% CL, and the best-fit point has tan β 36. (This increases to tan β ∼ 45 if either the LHC 13-TeV and/or (g − 2) µ constraint is dropped.) As in the previous two-dimensional projections of the sub-GUT parameter space, the 99.7% (3-σ) CL contour lies close to that for the 95% CL.

B Decay Observables
We see in the left panel of Fig. 8  sub-GUT 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5   that we have studied can accommodate comfortably the preference seen in the data (dotted line) for such a small value of BR(B s,d → µ + µ − ) 6 , which is not the case in models such as the CMSSM that impose universal boundary conditions on the soft supersymmetry-breaking parameters at the GUT scale, if µ > 0. The right panel of Fig. 8 shows how the contributions of the flavour (blue shading) and other observables to the global likelihood function depend on M in for values between 10 4 and 10 16 GeV. This variation in the flavour contribution (which is dominated by BR(B s,d → µ + µ − )) is largely responsible for the sub-GUT preference for M in < M GUT seen in the top left panel of Fig. 1. Values of M in ∈ (10 5 , 10 12 ) GeV can accommodate very well the experimental value of BR(B s,d → µ + µ − ). This preference is made possible by the different RGE running in the sub-GUT model, which can change the sign of the product A t µ that controls the relative signs of the SM and SUSY contributions to the B s,d → µ + µ − decay amplitudes, permitting negative interference that reduces BR(B s,d → µ + µ − ). As already discussed, the reduction in BR(B s,d → µ + µ − ) and the global χ 2 function for 10 8 GeV M in 10 12 GeV is associated with the blue 68% CL regions with M in 10 12 GeV seen in the middle panels of Fig. 1. On the other hand, we see in Fig. 9 that sub-GUT models favour values of BR(b → sγ) that are close to the SM value.
The contributions to the global χ 2 function of other classes of observables as functions of M in are also exhibited in the right panel of Fig. 8. In addition to the aforementioned reduction in the flavour contribution when M in 10 12 GeV (blue shading), there is a coincident (but smaller) increase in the contribution of the electroweak precision observables (orange shading) related to tension in the electroweak symmetry-breaking conditions. The other contributions to the global χ 2 function, namely the nuisance parameters (red shading), Higgs mass (light green), (g − 2) µ (teal) and DM (red), vary smoothly for M in ∼ 10 12 GeV.

Higgs Mass
We see in Fig. 10 that the profile likelihood function for M h lies within the contribution of the direct experimental constraint convoluted with the uncertainty in the FeynHiggs calculation of M h (dotted line). We infer that there is no tension between the direct experimental measurement of M h and the other observables included in our global fit. We have also calculated (not shown) the branching ratios for Higgs decays into γγ, ZZ * and gg (used as a proxy for gg → h production), finding that they are expected to be very similar to their values in the SM, with 2-σ ranges that lie well within the current experimental uncertainties.

Searches for Dark Matter Scattering
The left panel of Fig. 11 shows the nominal predictions for the spin-independent DM scattering cross-section σ SI p obtained using the SSARD code [53]. We caution that there are considerable uncertainties in the calculation of σ SI p , which are taken into account in our global fit. Thus points with nominal values of σ SI p above the experimental limit may nevertheless lie within the 95% CL range for the global fit. We see that sub-GUT models favour a range of σ SI p close to the present limit from the LUX, XENON1T and PandaX-II experiments 7 . Moreover, at the 95% CL, the nominal sub-GUT predictions for σ SI p are within the projected reaches of the LZ and XENON1/nT experiments. However, they are subject to the considerable uncertainty in the σ SI p matrix element, and might even fall below the neutrino 'floor' shown as a dashed orange line in [69].
We see in the right panel of Fig. 11 that the sub-GUT predictions for the spin-dependent DM scattering cross-section σ SD p lie somewhat below the present upper limit from the PICO direct DM search experiment. Spin-dependent DM scattering is also probed by indirect searches for neutrinos produced by the annihilations of neutralinos trapped inside the Sun after scattering on protons in its interior. If the neutralinos annihilate into τ + τ − , the IceCube experiment sets the strongest    Figure 11. Left panel: Two-dimensional profile likelihood function for the nominal value of σ SI p calculated using the SSARD code [53] in the (mχ0 1 , σ SI p ) plane, displaying also the upper limits established by the LUX [3], XENON1T [5] and PandaX-II Collaborations [6] shown as solid black, blue and green contours, respectively. The projected future 90% CL sensitivities of the LUX-Zeplin (LZ) [67] and XENON1T/nT [68] experiments are shown as dashed magenta and blue lines, respectively, and the neutrino background 'floor' [69] is shown as a dashed orange line with yellow shading below. Right panel: Two-dimensional profile likelihood function for the nominal value of σ SD p calculated using the SSARD code [53] in the (mχ0 1 , σ SD p ) plane, showing also the upper limit established by the PICO Collaboration [4].We also show the indirect limits from the Icecube [56] and Super-Kamiokande [70] experiments, assumingχ 0 1χ 0 1 → τ + τ − dominates, as well as the 'floor' for σ SD p calculated in [71].
such indirect limit [56], and we also show the constraint from Super-Kamiokande [70]. These constraints are currently not sensitive enough to cut into the range of the (mχ0 1 , σ SD p ) plane allowed in our global fit. We also show the neutrino 'floor' for σ SD p , taken from [71]: wee that values of σ SD p below this floor are quite possible in the sub-GUT model.

Impacts of the LHC 13-TeV and New Direct Detection Constraints
We show in Fig. 12 some two-dimensional projections of the regions of sub-GUT MSSM parameters favoured at the 68% (red lines), 95% (blue lines) and 99.7% CL (green lines), comparing the results of fits including the LHC 13-TeV data and recent direct searches for spin-independent dark matter scattering (solid lines) and discarding them (dashed lines). The upper left panel shows the (mq R , mg) plane, the upper right plane shows the (mq R , mχ0 1 ) plane, the lower left plane shows the (mg, mχ0 1 ) plane, and the lower right panel shows the (mχ0 1 , σ SI p ) plane. We see that in the upper panels that the new data restrict the favoured parameter space for mq R ∼ 2 TeV, the two left panels show a restriction for mg ∼ 1.3 TeV, and the right and lower panels show that the new data also restrict the range of mχ0 1 to 800 GeV. However, the lower right panel does not show any new restriction on the range of possible values of σ SI p . sub-GUT w/ LHC13: best fit, 1σ, 2σ, 3σ sub-GUT w/o LHC13: best fit, 1σ, 2σ, 3σ Figure 12. Two-dimensional projections of the global likelihood function for the sub-GUT MSSM in the (mq R , mg) plane (upper left panel), the (mq R , mχ0 1 ) plane (upper right panel), the (mg, mχ0 1 ) plane (lower left panel), and the (mχ0 1 , σ SI p ) plane (lower right panel). In each panel we compare the projections of the sub-GUT parameter regions favoured at the 68% (red lines), 95% (blue lines) and 99.7% CL (green lines) in global fits with the LHC 13-TeV data and results from LUX, XENON1T, and PandaX-II [3,5,6] (solid lines), and without them (dashed lines).

Best-Fit Points, Spectra and Decay Modes
The values of the input parameters at the bestfit points with and without the (g − 2) µ and LHC 13-TeV constraints have been shown in Table 2. The best fits have M in between 1.6 × 10 5 and 4.1 × 10 8 GeV, and we note that the input parameters are rather insensitive to the inclusion of the 13-TeV data when (g − 2) µ is dropped. Table 4 displays the mass spectra obtained as outputs at the best-fit point including the 13-TeV data (quoted to 3 significant figures) and including (left column) or dropping (right column) the (g − 2) µ constraint. As could be expected, the sparticle masses are generally heavier when (g − 2) µ is dropped. However, the differences are small in the cases of theχ 0 1 ,χ 0 2 andχ ± 1 , being generally < 10 GeV. We also give in the nextto-last line of Table 4 the values of the global χ 2 function at these best-fit points, dropping the HiggsSignals contributions, as was done previously [21,33] to avoid biasing the analysis.
The contributions of different observables to the global likelihood function at the best-fit points including LHC13 data are shown in Fig. 13. We compare the contributions when (g−2) µ is included (pink histograms) and without (g − 2) µ (blue histograms). We note, in particular, that the contribution of BR(B s,d → µ + µ − ) is very small in both cases, which is a distinctive feature of sub-GUT models.
The last line of Table 4 shows the p-values for the best fits with and without (g − 2) µ , which were calculated as follows. In the case with (without) (g − 2) µ , setting aside HiggsSignals so as to avoid biasing the analysis [21,33], the number of constraints making non-zero contributions to the global χ 2 function (not including nuisance parameters) is 29 (28), and the number of free parameters is 5 in each case. Hence the numbers of degrees of freedom are 24 (23) in the two cases. The values of the total χ 2 function at the best-fit points, dropping the HiggsSignals contribution, are 28.9 (18.0) and the corresponding p-values are 23% (76%). The qualities of the global fits with and without (g−2) µ are therefore both good. and the fit including (g − 2) µ is not poor enough to  Table 4 The spectra at the best-fit points including the LHC 13-TeV data and including (left column) or dropping (right column) the (g − 2) µ constraint. The masses are quoted in GeV. The three bottom lines give the values of the χ 2 function dropping HiggsSignals, the numbers of degrees of freedom (d.o.f.) and the corresponding p-values.
reject this fit hypothesis. The spectra for the best fits are displayed graphically in Fig. 14, including the (g − 2) µ constraint (upper panel) and dropping it (lower panel). Also shown are the decay modes with branching ratios > 5%, as dashed lines whose intensities increase with the branching ratios. The heavy Higgs bosons decay predominantly to SM final states, hence no dashed lines are shown. We see that in both cases the squarks and gluino are probably too heavy to be discovered at the LHC, and the sleptons are too heavy to be discovered at any planned e + e − collider. The best prospects for sparticle discovery may be forχ ± 1 andχ 0 2 production at CLIC running at E CM 2 TeV [75].   Figure 13. Contributions to the global χ 2 function at the best-fit points found in our sub-GUT analysis including LHC 13-TeV data, in the cases with and without the (g − 2) µ constraint (pink and blue histograms, respectively).  Figure 14. The spectra of Higgs bosons and sparticles at the best-fit points in the sub-GUT model including LHC 13-TeV data, including the (g − 2) µ constraint (upper panel) and dropping it (lower panel), with dashed lines indicating the decay modes with branching ratios > 5%. These plots were made using PySLHA [76].
The global likelihood function is quite flat at large sparticle masses, and very different spectra are consistent with the data, within the current uncertainties. The 68 and 95% CL ranges of Higgs and sparticle masses are displayed in Fig. 15 as orange and yellow bands, respectively, with the best-fit values indicated by blue lines. The upper panel is for a fit including the (g − 2) µ constraint, which is dropped in the lower panel. At the 68% CL there are possibilities for squark and gluino discovery at the LHC and theτ 1 ,μ R andẽ R become potentially discoverable at CLIC if it operates at E CM = 3 TeV [75].

Summary and Perspectives
We have performed in this paper a frequentist analysis of sub-GUT models in which soft supersymmetry-breaking parameters are assumed to be universal at some input scale M in < M GUT . The best-fit input parameters with and without (g − 2) µ and the LHC 13-TeV data are shown in Table 2. The physical sparticle masses including the LHC data, with and without (g − 2) µ , are shown in Table 4 and in Fig. 14, where decay patterns are also indicated. As seen in the bottom line of Table 4, the p-values for the fits with and without (g−2) µ are 23% and 76%, respectively.
Compared to the best fits with M in = M GUT , we have found that the minimum value of the global χ 2 function may be reduced by ∆χ 2 ∼ 2 in the sub-GUT model, with the exact amount depending whether the (g − 2) µ constraint and/or LHC13 data are included in the fit. Whether these observables are included, or not, the global χ 2 minimum occurs for M in ∼ 10 7 GeV, and is due to the sub-GUT model's ability to provide a better fit to the measured value of BR(B s,d → µ + µ − ) than in the CMSSM. Although intriguing, this improvement in the fit quality is not very significant, but it will be interesting to monitor how the experimental measurement of BR(B s,d → µ + µ − ) evolves.
In all the scenarios studied (with/without (g − 2) µ and/or LHC13), the profile likelihood function for mg (mq) varies by 1 for mg 1.9 TeV (mq 2.2 TeV). The corresponding slowly-varying ranges of χ 2 for mt 1 (mb 1 ) start at ∼ 1 TeV (∼ 1.6 TeV), respectively. On the other hand, we find a more marked preference for mχ0 1 ∼ 1 TeV, with theχ ± 1 andχ 0 2 being slightly heavier and large mass values being disfavoured at the ∆χ 2 ∼ 3 level. The best-fit point is in a region where rapid annihilation via H/A poles is hybridized with stop coannihilation, with chargino coannihilation and stau coannihilation also playing roles in both the 68 and 95% CL regions. Within the 95% CL region, the chargino lifetime may exceed 10 −12 s, and the stau lifetime may be as long as one second, motivating continued searches for long-lived sparticles at the LHC.
Taking the LHC13 constraints into account, we find that the spin-independent DM cross-section, σ SI p , may be just below the present upper limits from the LUX, XENON1T and PandaX-II experiments, and within the reaches of the planned XENONnT and LZ experiments. On the other hand, the spin-dependent DM cross-section, σ SD p , may be between some 2 and 5 orders of magnitude below the current upper limit from the PICO experiment.
Within the sub-GUT framework, therefore, we find interesting perspectives for LHC searches for strongly-interacting sparticles via the conventional missing-energy signature. Future / E T searches for electroweakly-interacting sparticles and for long-lived massive charged particles may also have interesting prospects. The best-fit region of parameter space accommodates the observed deviation of BR(B s,d → µ + µ − ) from its value in the SM, and it will be interesting to see further improvement in the precision of this measurement. A future e + e − collider with centre-of-mass energy above 2 TeV, such as CLIC [75], would have interesting perspectives for discovering and measuring the properties of electroweakly-interacting sparticles. There are also interesting perspectives for direct DM searches via spin-independent scattering.