Likelihood analysis of the pMSSM11 in light of LHC 13-TeV data

We use MasterCode to perform a frequentist analysis of the constraints on a phenomenological MSSM model with 11 parameters, the pMSSM11, including constraints from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 36$$\end{document}∼36/fb of LHC data at 13 TeV and PICO, XENON1T and PandaX-II searches for dark matter scattering, as well as previous accelerator and astrophysical measurements, presenting fits both with and without the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(g-2)_\mu $$\end{document}(g-2)μ constraint. The pMSSM11 is specified by the following parameters: 3 gaugino masses \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{1,2,3}$$\end{document}M1,2,3, a common mass for the first-and second-generation squarks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\tilde{q}}$$\end{document}mq~ and a distinct third-generation squark mass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\tilde{q}_3}$$\end{document}mq~3, a common mass for the first-and second-generation sleptons \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\tilde{\ell }}$$\end{document}mℓ~ and a distinct third-generation slepton mass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\tilde{\tau }}$$\end{document}mτ~, a common trilinear mixing parameter A, the Higgs mixing parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}μ, the pseudoscalar Higgs mass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_A$$\end{document}MA and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tan \beta $$\end{document}tanβ. In the fit including \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(g-2)_\mu $$\end{document}(g-2)μ, a Bino-like \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^0_{1}$$\end{document}χ~10 is preferred, whereas a Higgsino-like \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^0_{1}$$\end{document}χ~10 is mildly favoured when the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(g-2)_\mu $$\end{document}(g-2)μ constraint is dropped. We identify the mechanisms that operate in different regions of the pMSSM11 parameter space to bring the relic density of the lightest neutralino, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^0_{1}$$\end{document}χ~10, into the range indicated by cosmological data. In the fit including \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(g-2)_\mu $$\end{document}(g-2)μ, coannihilations with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^0_{2}$$\end{document}χ~20 and the Wino-like \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^\pm _{1}$$\end{document}χ~1± or with nearly-degenerate first- and second-generation sleptons are active, whereas coannihilations with the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^0_{2}$$\end{document}χ~20 and the Higgsino-like \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\chi }^\pm _{1}$$\end{document}χ~1± or with first- and second-generation squarks may be important when the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(g-2)_\mu $$\end{document}(g-2)μ constraint is dropped. In the two cases, we present \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi ^2$$\end{document}χ2 functions in two-dimensional mass planes as well as their one-dimensional profile projections and best-fit spectra. Prospects remain for discovering strongly-interacting sparticles at the LHC, in both the scenarios with and without the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(g-2)_\mu $$\end{document}(g-2)μ constraint, as well as for discovering electroweakly-interacting sparticles at a future linear \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^+ e^-$$\end{document}e+e- collider such as the ILC or CLIC.


Introduction
Supersymmetric (SUSY) models of TeV-scale physics are being subjected to increasing pressure by the strengthening constraints imposed by LHC experiments [1,2] and searches for Dark Matter (DM) [3][4][5][6]. In particular, in the context of models with soft supersymmetry-breaking parameters constrained to be universal at a high unification scale, the LHC limits on sparticle masses have been in increasing tension with a supersymmetric interpretation of the anomalous magnetic moment of the muon, (g − 2) μ , which would require relatively light sleptons and electroweak gauginos [7][8][9][10][11][12][13]. This pressure has been ratcheted up by the advent of ∼ 36/fb of data from Run 2 of the LHC at a centre-of-mass energy of 13 TeV [14][15][16], 1 which probe supersymmetric models at significantly higher mass scales than was possible in Run 1 at 7 and 8 TeV in the centre of mass. In parallel, direct searches for DM scattering have also been making significant progress towards the neutrino 'floor' [17,18], in particular with the recent data releases from the LUX, PICO, XENON1T and PandaX-II experiments [3][4][5][6]. Here we analyze these constraints in the minimal supersymmetric extension of the Standard Model (MSSM), which, because of R-parity, has a stable cosmological relic particle that we assume to be the lightest neutralino,χ 0 1 [19,20]. The strengthening phenomenological, experimental and astrophysical constraints on supersymmetry (SUSY) were initially explored mainly in the contexts of models in which SUSY breaking was assumed to be universal at the GUT scale, such as the constrained MSSM (CMSSM) [7][8][9][10][11][21][22][23][24][25][26][27][28][29][30][31], non-universal Higgs models (NUHM1,2) [11,12], 2 the minimal anomaly-mediated SUSY-breaking model (mAMSB) [33], and models based on the SU(5) group [34]. These models are tractable by virtue of having a relatively limited number of parameters, though the universality assumptions they employ are not necessarily well supported in scenarios motivated by fundamental principles, such as string theory. Their limited parameter spaces are amenable to analysis, e.g., in the frequentist approach we follow, in which one constructs a global likelihood function that embodies all the information provided by the multiple constraints.
Alternatively, one may study phenomenological models in which the soft SUSY-breaking parameters are not constrained by any universality assumptions, though subject to milder constraints emanating, in particular, from upper limits on SUSY contributions to flavour-changing processes. These phenomenological MSSM (pMSSM) [35][36][37][38][39][40][41][42][43] models contain many more parameters, whose exploration is computationally demanding. There have been cut-based global analyses of variants of the pMSSM with as many as 19 parameters [44][45][46][47][48] and global fits focused on specific sectors or parameter ranges [49,50], however in the past we have restricted our frequentist attentions to a variant of the pMSSM with 10 parameters, the pMSSM10 [13,51]. These were taken to be 3 independent gaugino masses, M 1,2,3 , a common electroweak-scale mass for the first-and second-generation squarks, mq , a distinct mass for the third-generation squarks, mq 3 , a common electroweak-scale mass ml for the sleptons, a single trilinear mixing parameter A that is universal at the 1 We use here results from SUSY searches by the CMS Collaboration: the results from ATLAS [2] yield similar constraints. 2 For a recent analysis of these models in light of ∼ 13/fb of LHC data at 13 TeV, see [32]. This analysis does not include the PICO, XENON1T and most recent PandaX-II results, and has other differences that are noted later in this paper. electroweak scale, the Higgs mixing parameter μ, the pseudoscalar Higgs mass, M A and the ratio of Higgs vevs, tan β. 3 It is desirable to extend this type of analysis to more general variants of the pMSSM, for a couple of reasons. One is that the lower bounds on sparticle masses will, in general, be weaker in models with more parameters, so one should explore such models before making statements about the magnitudes of these lower bounds and prospects for discovering sparticles at the LHC or elsewhere. Another reason is that reconciling the strengthening LHC constraints with the cosmological DM density constraint requires, in general, specific relations between sparticle masses that suppress the relic density via coannihilation effects and/or rapid annihilations through direct-channel resonances. Therefore one should study models capable of accommodating these DM mechanisms [51].
Examples of DM mechanisms that have been studied extensively in the past [51] include coannihilation with the lighter stau slepton,τ 1 , the lighter chargino,χ ± 1 , or the lighter stop squark,t 1 , and rapid annihilations via the Z boson, the 125-GeV Higgs boson, h, or the heavier MSSM Higgs bosons, H/A. More recently, the possibility of coannihilation with gluinos,g, has been explored in models with nonuniversal gaugino masses [53,54], and coannihilation with the right-handed up-type squarks of the first two generations, u R /c R , emerged as a possibility in an SU(5) model with nonuniversal scalar masses m 5 , m 10 for sfermions in5 and 10 representations [34].
All of these were possibilities in the pMSSM10, but in that scenario the stau and smuon masses were fixed to be equal, putting the LHC constraints on stau coannihilation in tension with the possibility of a SUSY interpretation of (g − 2) μ , a tension that has increased with the advent of the first LHC data at 13 TeV. In this paper we study two possible resolutions of this issue. We study an extension of the parameter space of the pMSSM10 to 11 parameters by relaxing the equality between the soft SUSYbreaking contributions to the stau mass and to the (still common) masses of the smuon and selectron, the pMSSM11. In order to assess the importance of the (g − 2) μ constraint, we also consider a fit omitting the SUSY interpretation of (g − 2) μ . The principal results of this paper are comparisons between the likelihoods of different spectra in the pMSSM11 with and without (g − 2) μ , and comparisons between the likelihoods of different DM mechanisms includingτ 1 ,˜ ,q andg coannihilation, highlighting the impacts of the LHC 13 TeV and recent DM scattering data.
The layout of this paper is as follows. In Sect. 2 we specify the framework of our analysis. Section 2.1 speci- Table 1 The ranges of the pMSSM11 parameters sampled, which are divided into the indicated numbers of segments, yielding the total number of sample boxes shown in the last row. In the last column, we indicate the kind of prior used, where "soft" means a flat prior with Gaussian tails  (5) had (M Z ) [56] μ = 0.02771, σ = 0.00011 1 Gaussian Total number of boxes 384 fies the pMSSM11, establishes our notation for its parameters and describes our procedure for sampling the pMSSM11 parameter space. In Sect. 2.2 we review the MasterCode tool to construct a global χ 2 likelihood function combining constraints on model parameters, Sect. 2.3 describes our treatments of the electroweak and flavour constraints, including some updates compared with our previous analyses. In Sect. 2.4 we give details on our DM analysis, which includes constraints on both spin-independent and -dependent DM scattering [3][4][5][6]. Our implementations of the constraints from ∼ 36/fb of LHC at 13 TeV [14][15][16] are discussed in Sect. 2.5. Then, in Sect. 3.1 we present results for the global likelihood function in various parameter planes, highlighting the regions where different DM mechanisms operate and comparing results with and without the (g − 2) μ constraint being applied. Section 4 displays the one-dimensional profile likelihood functions for various masses, mass differences and other observables in these two cases, and also shows predictions for spin-independent and -dependent DM scattering. Section 5 highlights the impacts of the LHC 13-TeV data [14][15][16] and the recent direct searches for astrophysical DM [3][4][5][6]. Section 6 discusses the best-fit points, favoured and allowed spectra in these pMSSM scenarios. Finally, Sect. 7 summarizes our conclusions.
2 Analysis framework

Model parameters
As mentioned above, in this paper we consider a pMSSM scenario with eleven parameters, namely where q 1,2 ≡ u, d, s, c, we assume soft SUSY-breaking parameters for left-and right-handed sfermions, and the sneutrinos have the same soft SUSY-breaking parameter as the corresponding charged sfermions. All of these parameters are specified at a renormalisation scale M SUSY given by the geometric mean of the masses of the scalar top eigenstates, M SUSY ≡ √ mt 1 mt 2 , which is also the scale at which electroweak symmetry breaking conditions are imposed. We allow the sign of the mixing parameter μ to be either positive or negative. The important difference from the pMSSM10 scenario we studied previously [13] is that the first-and second-generation slepton mass m˜ and the stau mass mτ are decoupled in the pMSSM11. 4 The ranges of these parameters sampled in our analysis are displayed in Table 1. In each case, we indicate in the third column of Table 1 how the ranges of most of these parameters are divided into segments, much as we did previously for our analysis of the pMSSM10 [13].
These segments define boxes in the eleven-dimensional parameter space, which we sample using the MultiNest package [57][58][59]. In order to ensure a smooth overlap between boxes and eliminate features associated with their boundaries, we choose 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. An initial scan over all mass parameters with absolute values ≤ 4 TeV showed that non-trivial behaviour of the global likelihood function was restricted to |M 1 | 1 TeV and m˜ 1 TeV. In order to achieve high resolution efficiently, we restricted the range of m˜ to < 2 TeV in the full scan. 5 To study properly the impact of the (g −2) μ , we performed separate sampling campaigns with and without it. On the other hand, during the sampling phase the constraints coming from LHC13 results have not been included. Since their impact consists in providing lower bounds to the sparticle masses, this choice allows for a proper assessment of their impact on the full parameter space. Moreover, we also performed dedicated scans for various DM annihilation mechanisms, in such a way to improve the quality of the sample in the description of the fine-tuned spectrum configurations that characterize them. The data sets from the various campaigns have been merged into a single set on which the likelihood is computed dynamically including or excluding the (g − 2) μ and/or the LHC13 constraints according to our interest. The total number of points in our pMSSM11 parameter scan is ∼ 2 × 10 9 .

MasterCode
We perform a global likelihood analysis of the pMSSM11 including constraints from direct searches for SUSY particles at the LHC, measurements of the Higgs boson mass and signal strengths, LHC searches for SUSY Higgs bosons, precision electroweak observables, flavour constraints from Band K -physics observables, the cosmological constraint on the overall cold dark matter (CDM) density, and upper limits on spin-independent and -dependent LSP-nuclear scattering. We treat (g −2) μ as an optional constraint, presenting results from global fits with and without it, and we treat m t , α s and M Z as nuisance parameters.

Electroweak and flavour constraints
Our treatments of many of these constraints follow those we have used previously, which were summarized most recently in Table 1 in [34]. Table 2 summarizes the updates we make in this paper. As noted there, the only change in the electroweak sector is in M W . 9 Here we follow [96] in combining naively the recent ATLAS measurement M W = 80.370±0.019 GeV with the previous world average value M W = 80.385 ± 0.015 GeV, obtaining M W = 80.379 ± 0.012 GeV. 10 Since one of our objectives in this paper is to emphasize the impact on the pMSSM11 parameter space of the (g − 2) μ constraint, for reference we also include in Table 2 the implementation of this constraint that we use as an option. 11 As can be seen in Table 2, we have also updated a number of flavour constraints. In particular, we have updated the global analysis of BR(B s,d → μ + μ − ) to include the latest 6 We use here an updated version of FeynWZ (not yet publicly available) in which the M W evaluation is based on [65] and is identical to that implemented in FeynHiggs, which gives more reliable results in parameter regions with larger SUSY masses or small SUSY mass splittings. The other EWPO are treated in the same way as in [63,64]. 7 We note that FeynHiggs incorporates resummation effects in Higgs mass calculations that are not included in the MSSM FlexibleSUSY generator [72] used in [32,52], although these are available through other FlexibleSUSY generators, HSSUSY/SplitSUSY [73][74][75] and FlexibleEFT [76]. It should also be noted that FeynHiggs has recently been improved for higher SUSY mass scales [77,78]. 8 The SSARD computation of the scattering cross-section follows the computations detailed in [86][87][88]. The uncertainties in the crosssections are derived from a straightforward propagation of errors in the input quantities which determine the cross-section. The dominant uncertainties are discussed below in more detail. 9 We emphasize that, although they are not displayed in Table 2 because they have not changed since [34], we use a complete set of electroweak constraints, not restricted to M W as used in [32,52]. We also note that the FeynWZ code we use to calculate M W incorporates 2-loop corrections that are not included in the FlexibleSUSY code [72] used in [32,52]. 10 In so doing, we neglect correlations in the uncertainties due to PDFs, QED and boson p T modelling, but our results are relatively insensitive to the details of this combination. 11 The (g − 2) μ evaluation in FeynHiggs contains less sophisticated two-loop corrections than GM2CALC [97]. However, the difference is small compared with other uncertainties in our analysis.  Table 1 in [34]. We indicate separately the experimental and applicable theoretical errors in the SM and SUSY (sometimes in combination, labelled "MSSM"). The contribution of the τ (B s → μ + μ − ) constraint to the global χ 2 likelihood function is essentially constant across the relevant region of the pMSSM11 parameter space, and it is not included in the fit.  [3,4,6] Combined likelihood in the (mχ0 [5] Likelihood in the (mχ0 [14,15] Combined likelihood in the (mg, mχ0 Run 2 result from LHCb [112] as well as the Run 1 results of CMS, LHCb [110] and ATLAS [111]. We assume minimal flavour violation (MFV) when combining the BR(B d → μ + μ − ) constraint with that from BR(B s → μ + μ − ) into the quantity R μμ [11], and take into account the correlation between the theoretical calculations of f B s and f B d .
The LHCb Collaboration has also published [112] a first determination of the effective B s lifetime as measured in B s → μ + μ − decays, providing a constraint on the quantity A via the relation where [114] y s = τ B s s 2 = 0.0675 ± 0.004, where τ B s is the inclusive B s decay lifetime, the complex numbers p, q specify the relation between the mass eigenstates of the B 0 s −B 0 s system and the flavour eigenstates [114], and A(B 0 s → μ + μ − ) and A(B 0 s → μ + μ − ) are the B 0 s andB 0 s decay amplitudes. In the Standard Model (SM), A = 1 so that τ (B s → μ + μ − )| SM = τ B s /(1 − y s ) = 1.619 ± 0.009 ps. On general grounds, A ∈ [−1, 1]. The LHCb measurement τ (B s → μ + μ − ) = 2.04±0.44(stat.)± 0.05(syst.) ps corresponds formally to A = 7.7 ± 10.0, implying that the current LHCb result does not constrain significantly the pMSSM11 parameter space, and we do not include it in our fit. However, in the later discussion of our fit results we present for information the χ 2 profile likelihood functions we find for A and τ (B s → μ + μ − ).
We have also updated our implementations of b → sγ , B → τ ν, B → X s , M B s and M B d to take account of updated theoretical calculations within the SM. For the same reason, in the kaon sector we have also updated our implementations of K → μν and K → πνν. 12 Since there are, in general, supersymmetric contributions to the observables commonly used in global fits to CKM parameters, we remove these contributions and make a global fit to the CKM parameters without them.
In general, we treat the electroweak precision observables, (g − 2) μ and all B-and K -physics observables (except for BR(B s,d → μ + μ − )) as Gaussian constraints, combining in quadrature the experimental and applicable SM and SUSY theory errors.

Dark matter constraints and mechanisms
Cosmological density Since we work in the framework of the MSSM, R-parity is conserved, so that the lightest SUSY particle (LSP) is a candidate to provide the CDM. We assume that the LSP is the lightest neutralinoχ 0 1 [19,20], and that it is the dominant component of the CDM. As in our recent papers [33,34], we use the Planck 2015 constraint on the total CDM density: CDM h 2 = 0.1186 ± 0.0020 EXP ± 0.0024 TH [128].

Density mechanisms
As one of the primary objectives in our analysis is to investigate the relevances of various mechanisms for bringing the relicχ 0 1 density into the range allowed by astrophysics and cosmology, we introduce a set of measures related to particle masses that were found in our previous analyses [51] to indicate when specific mechanisms were dominant. 13 These may be grouped as follows.

• Coannihilation with an Ino
This may be important if theχ 0 1 is not much lighter than the lighter chargino,χ ± 1 , and the second neutralino,χ 0 2 , or the gluino,g. For these cases we introduce the coannihilation measures Ino coann. : We find that chargino andχ 0 2 coannihilation is important in our analysis, and in our 2-dimensional plots we shade green the regions where (4) is satisfied when the Ino is the lighter chargino,χ ± 1 (which is almost degenerate with thẽ χ 0 2 ). On the other hand, we find that gluino coannihilation is not important in the pMSSM11 when the (g − 2) μ constraint is imposed. This is due to the fact that (g − 2) μ forces the neutralino mass to values for which a gluino of equivalent mass would be excluded by current LHC results.
• Coannihilation with sleptons 13 We have checked specifically the validity of these measures using Micromegas, finding good consistency in most cases. However, in certain hybrid regions where more than one mechanism satisfied the criteria we found that just one mechanism dominates. Moreover, it might also happen that some regions of the parameter space are not classified by a given measure even if the corresponding mechanism is active.
In the version of the pMSSM that we study here, the two stau mass eigenvalues are similar, since the soft SUSYbreaking parameters are specified at the TeV scale and the left-right mixing ∝ m τ is relatively small, but the stau masses are not degenerate with the selectron and smuon masses, in general. We find that smuon and selectron coannihilation are in general more important than stau coannihilation, thanks to the greater multiplicity of near-degenerate states. We introduce the following coannihilation measure: and shade in yellow (pink) the regions of our twodimensional plots where (5) is satisfied for = μ, e (τ ), respectively.

• Coannihilation with squarks
Similarly, this may be important for squarksq that are not much heavier than theχ 0 1 . The case considered most often has beenq =t 1 , but here we consider all possibilities, including coannihilations with first-and secondgeneration squarks, which we find to be important when the LHC 13-TeV constraint or (g − 2) μ is dropped. We introduce the coannihilation measurẽ q coann. : mq mχ0 and we use the following colours in our plots for the regions where (6) is satisfied:q =d/s/ũ/c L ,R cyan,t 1 grey,b 1 purple.

• Annihilation via a direct-channel boson pole
When there is a massive boson B with mass M B ∼ 2mχ0 1 , χ 0 1χ 0 1 annihilation is enhanced along a 'funnel' in parameter space. We have found that such a mechanism is likely to dominate if the following condition is satisfied: We have considered the cases B = h, Z and H/A, and use blue shading for the regions of our subsequent plots where (7) is satisfied when B = H/A. We comment later on a small region where rapid annihilation via the h and Z poles is important.

• Enhanced Higgsino component
We have also considered a somewhat different possibility, namely that theχ 0 CMSSM: Regions where the condition (8) is satisfied generally satisfy the chargino coannihilation condition with a Higgsino-like LSP, and are also shaded green.

• Hybrid regions
In addition to the 'primary' regions where only one of the conditions (4)-(8) is satisfied, there are also 'hybrid' regions where more than one condition is satisfied. These are indicated in the following by mixtures of the corresponding primary colours.

Direct DM searches
We implement experimental constraints from direct searches for supersymmetric DM via both spin-independent and -dependent scattering on nuclei. We use the LUX [4], XENON1T [6] and PandaX-II [3] constraints on the spinindependent DM scattering cross section σ SI p , which we implement via a combined two-dimensional likelihood function in the (mχ0 1 , σ SI p ) plane. Our treatment of the spin-independent nuclear scattering matrix element follows that in our previous work [12] and is based on SSARD [85]. As reviewed, for example, in [88] the largest uncertainties in the matrix element are those associated with the pion-nucleon σ -term, π N , and the SU(3) octet symmetry-breaking contribution to the nucleon mass, σ 0 . These may be expressed as follows in terms ofqq matrix elements in the nucleon: from which we see that thess matrix element It is well known that σ SI p is sensitive to the value of y, and hence to the values of σ 0 and π N . We follow [129] in interpreting the measured octet baryon mass differences as yielding σ 0 = 36 ± 7 MeV, 14 and we follow our previous work in assuming here that π N = 50 ± 7 MeV, 15 corresponding to a central value of y = 0.28. For comparison, two recent determinations of π N give somewhat larger values that are, however, compatible with the value we assume, within the quoted uncertainties: π N = 59.1 ± 3.5 MeV (from pionic atoms) [132] and 58 ± 5 MeV (from π -nucleon scattering) [133] (see also [134], which found the value π N = 59 ± 7 MeV). On the other hand, lattice calculations [135][136][137][138] yield systematically smaller values of π N that are in tension with these data-driven estimates, as discussed in [133]. Our value of π N is intermediate and relatively conservative in that it implies a smaller value of y than the data-driven estimates of π N . 16 We also implement in this paper the PICO [5] constraint on the spin-dependent DM scattering cross section σ SD p , also using the SSARD code [85]. As discussed in [139], the spindependentχ 0 1 p scattering matrix element is determined by the light quark contributions to the proton spin, which we take to be [88] u = +0.84 ± 0.03, where the uncertainties are dominated by those in measurements of polarized deep-inelastic scattering, and hence are correlated: the uncertainty in the combination u− d (from g A ) is very small, and that in u + d −2 s (from semileptonic octet baryon decays) is also somewhat smaller. 17 Indirect astrophysical searches for DM These include searches for γ -rays from DM annihilations near the Galactic centre and in dwarf galaxies, and for energetic neutrinos produced by the annihilations of DM particles trapped inside the Sun. There are large astrophysical uncertainties in estimates of the possible γ -ray flux from the Galactic centre, and other studies have indicated that the available limits on the fluxes from dwarf galaxies do not yet impose competitive constraints on supersymmetric modelssee, for example, [140] and [32]. The strongest constraints on energetic solar neutrinos are those provided by the IceCube Collaboration [141]. Their impact depends on the annihilation final states, being strongest for annihilations into τ + τ − , somewhat weaker for W + W − , and much weaker forbb final states. The capture of dark matter particles in the Sun is often assumed to be dominated by energy loss due to spindependent scattering on protons, in which case an upper limit on the neutrino flux may be used to constrain the spindependent cross-section σ SD p , as done by the IceCube Col- 16 For comparison, a similar value of π N = 59 ± 9 MeV is assumed in [32], but with σ s ≡ m s N |ss|N = 43 ± 8 MeV inferred from lattice calculations. This corresponds to π N − σ 0 = (m u + m d )σ s /m s ∼ 3.5 MeV, implying a value of σ 0 different from the value we use, which is based on octet baryon masses. 17 The values (11) of the q that we use are similar to those used in [32].
laboration [141]. However, the interpretation of this constraint [141] depends on the importance of spin-independent scattering on 4 He and heavier nuclei inside the Sun, and whether the DM density inside the Sun is in equilibrium between capture and annihilation [142]. As discussed in Sect. 4.10, we have found in an exploratory study that the IceCube constraint has little impact once the more recent PICO constraint [5] on σ SD p is taken into account. In view of the fact that it has fewer uncertainties, we use the PICO result in our global fit, setting aside the IceCube result [141]. 18

13 TeV LHC constraints
The LHC constraints we consider are those from searches for coloured sparticles in events with missing transverse energy, / E T , accompanied by jets and possibly leptons, searches for electroweak inos in events with multiple leptons, searches for long-lived charged particles, measurements of the 125 GeV Higgs boson h, and searches for the heavier SUSY Higgs bosons H, A, H ± . Our principal focus in this paper is on the implications of Run-2 LHC searches with ∼ 36/fb of data at 13 TeV, though we also make comparisons with the situation before these constraints were released. Our implementations of the constraints from LHC Run 1 at energies of 7 and 8 TeV used in our previous analysis of the pMSSM10 model were described in [13], and our implementations of / E T searches with ∼ 13/fb of data at 13 TeV in the gluino and squark production channels were described in [34], as were our implementations of searches for long-lived charged particles and for H, A, H ± with similar data sets. We refer the reader to these publications for details of those implementations, focusing here on our implementations of the Run 2 searches with ∼ 36/fb of data.

Searches for gluinos and squarks
We consider the constraints from CMS simplified model searches using events with / E T and jets but no leptons released in [14] and events with / E T and jets and a single lepton released in [15].
In the approach taken, e.g., by CheckMATE [143], ColliderBit [144] and MadAnalysis 5 [145], Monte Carlo simulations are used to estimate the signal yield from a model point after the event selection and to test it by comparing it with the upper bound given by an experimental collaboration. However, such a method is time-consuming and computationally prohibitive for our purpose. To circumvent this issue, we take the Fastlim [146] approach 19 and consider the implications of [14] for the following supersymmetric topologies:gg → [qqχ 0 1 ] 2 and [bbχ 0 1 ] 2 , and qq → [qχ 0 1 ][qχ 0 1 ], and the implications of [15] for the 18 In contrast, [32] uses the IceCube result, but not the PICO result. 19 The SmodelS code [147,148] takes a similar approach, as described in [146].
topologygg → [ttχ 0 1 ] 2 . The kinematics of each of these topologies depends on a reduced subset of sparticle masses, e.g., (mg, mχ0 1 ) in the case of thegg → [qqχ 0 1 ] 2 topology, and the CMS publications [14,15] provide in Root files 95% CL upper limits σ UL on the cross sections in the corresponding parameter planes. For each point in the main pMSSM11 sample, we calculate for thegg initial state and various final states contributions to the global χ 2 likelihood function of the form where SM denotes the Standard Model particles considered in each topology, SM ≡ qq, bb and tt, and analogously for theqq → [qχ 0 1 ][qχ 0 1 ] topology, where SM ≡ q andq. We use NLL-fast [149,150] to compute the cross sections for coloured sparticle pair-production up to NLO+NLL level.
If gluino and squarks have comparable masses, associated gluino-squark production may be sizeable. In the mg mq region, a fraction of the gq →gq process where the gluino decays intoq +q may be regarded as the production of a squark-antisquark pair with a soft quark jet. Ignoring this soft jet, we can constrain this process by considering the qq →qq simplified model limit. In the analyses we consider, jets are treated inclusively and this extra quark jet tends to slightly increase the acceptance. Ignoring the soft jet therefore results in underestimation of the signal acceptance, leading to a conservative limit. In order to constrain the gq →gq →qqq process in the same way as qq →qq, we rescale the squark cross-section as σqq → σqq + σgq · BRg →qq before applying squark simplified model limit.
Similarly, in the mq mg region we rescale the gluino cross-section as σgg → σgg + σgq · BRq →qg to constrain the gq →gq →ggq process using the gluino simplified model limit.

Stop and sbottom searches
Our treatment of LHC 13 TeV limits on stops and sbottoms is similar in principle to our implementation of the gluino and squark constraints described above. It is based on CMS simplified model searches in the jets + 0 [14] and 1 [15] lepton final states, where the results are interpreted as limits on the following topologies: We also use Fastlim to implement the CMS constraints in all these channels, following the same procedure as described above for gluinos and squarks, and estimating the corresponding contributions to the global χ 2 likelihood function as where SM = t, c and bW + forq 3 =t 1 and SM = b for q 3 =b 1 , respectively. In a significant part of the pMSSM11 parameter space, the neutralino relic abundance is brought into the observed range by Wino or Higgsino coannihilation mechanisms. In these regions,χ ± 1 andχ 0 1 are highly mass degenerate, with a mass difference that is typically smaller than 5 GeV. Since the decay products of theχ ± 1 →χ 0 1 transition are too soft to affect the signal acceptance, we can replaceχ ± 1 byχ 0 1 in the simplified topology. This approximation allows us to constrain thet 1 Thus, in the Wino and Higgsino coannihilation regions, we replace, e.g., the numerator in (13) by σt , enhancing the sensitivity.

Searches for electroweak inos
The CMS Collaboration has also released results from searches for electroweak ino production at the LHC in multilepton final states with ∼ 36/fb of data at 13 TeV [16]. The As in the cases of searches for strongly-interacting sparticles described above, we use Fastlim to compare the cross-section times branching ratio with the 95% CL upper limit released by CMS [16]. We obtain the corresponding contributions to the global χ 2 likelihood function as where SM ≡ W or Z , one or two ± and one or two τ ± , respectively. One complication compared to the previous coloured sparticle cases is that σχ± 1χ 0 2 depends on many MSSM parameters: and it is not feasible to tabulate the cross section directly in a multi-dimensional look-up table. We have therefore used the code EWK-fast [151], which is based on the observation that σ ( pp →χ ± 1χ 0 2 ) factorizes mathematically (whereχ i andχ j represent any chargino and/or neutralino): where T a (U) is a function of the mixing matrices U = {U, V, N } that can be calculated analytically. The factor F a (mχ i , mχ j , m a ) captures the kinematics and the effect of the parton distribution function and is tabulated in 3-dimensional look-up tables as a function of mχ i , mχ j and m a , where m a = mq L , mũ R or md R . The electroweak ino analyses described above can be extended to constrain models in which electroweak inos can be produced in the decays of coloured sparticles. This is because these searches do not impose conditions on the number of jets and the final states in such events resemble those arising from the direct production of electroweak inos associated with initial-state QCD radiation. In order to constrain this class of events we include an extra contribution to the electroweak ino cross-section, much as we discussed above in the case of theqg constraint. For example, in order to constrainqq →χ iχ j + jets, we rescale the cross-section: σχ iχ j → σχ iχ j + σq˜q BRq → jχ i BR˜q → jχ j before applying the electroweak ino simplified limit. 20 2.6 Combination of contributions to global χ 2 function from LHC sparticle searches The total contribution of LHC Run-2 sparticle searches is obtained by adding the contributions from the coloured sparticle (12) and (13) and electroweak ino searches (14): where the sum is over all the distinct SM final states mentioned above. The simple sum is justified because event samples with different final states are statistically independent, so that their correlations are not important for our analysis.
We summarise the simplified model limits we use in our scan in Table 3.

Measurements of the h(125) boson
These are incorporated via the HiggsSignals code [90,91], which implements the information from ATLAS and CMS measurements from LHC Run 1, as summarized in the joint ATLAS and CMS publication [152].  [16] 2.8 Searches for heavy MSSM Higgs bosons These are incorporated via the HiggsBounds code [92][93][94][95], which implements the information from ATLAS and CMS measurements from LHC Run 1, supplemented by the constraint from ∼ 36/fb of data from the LHC at 13 TeV provided by ATLAS [127].

Searches for long-lived or stable charged particles
The CMS Collaboration has published a search for charged particles with lifetimes 3 ns [123], and a search for massive charged particles that leave the detector without decaying [122]. We do not include the results of these searches in our global likelihood analysis, but comment later on their potential impacts. The only constraint that we impose on long-lived charged sparticles a priori is to require the lifetime to be smaller than 10 3 s so as to avoid modifying the successful predictions of cosmological nucleosynthesis calculations [153][154][155][156][157][158][159].

Global fit results
The input parameter values for our best-fit points with and without (g −2) μ are shown in the second and fourth columns of Table 4, and the spectra and dominant decays shown in Fig. 1. The third and fifth columns show input values for other points of interest that we discuss below. Lower rows of Table 4 show the total χ 2 per degree of freedom (d.o.f.) for each point, dropping the contributions from HiggsSignals that are shown in the last line. We also show the corresponding p-values, as calculated using the prescription described in [34] to estimate the number of degrees of freedom. 21 We ignored the contribution to the likelihood 21 In previous studies (see, e.g., the first paper in [7][8][9][10]) we have validated our naive p-value approximation with toy experiments, and found that it provides a reasonably accurate and conservative estimate of the coming from the nuisance parameters, and we removed the contribution to the likelihood from HiggsSignals, so as to avoid biasing our results by giving too much importance to the Higgs signal rates. Since all the other constraints contribute significantly to χ 2 function somewhere in the pMSSM11, we include them all in the d.o.f. count. However, we merged into a single constraint the LHC direct searches for sparticle production at 8 and 13 TeV, and also combined the 8-and 13-TeV limits on heavy Higgs bosons from A/H → τ + τ − searches. This results in totals of 31 and 30 constraints for the cases with and without (g − 2) μ , respectively. Since the number of free parameters is 11, this yields 20 and 19 for the numbers of d.o.f. in the two cases, as stated in Table 4. We note that the p-values are all comfortably high, whether (g − 2) μ is included, or not.

Parameter planes
We now display results from our global fits with and without (g − 2) μ in pairs of 2-dimensional pMSSM11 parameter planes. We indicate the locations of the best-fit points in these two-dimensional projections by green stars, We also show in these planes the χ 2 = 2.30, 5.99 and 11.3 contours, corresponding approximately to the boundaries of the regions preferred/allowed/possible at the 1-/2-/3-σ levels (68%, 95% and 99.7% CL), as red, blue and green solid lines, respectively. Within the 2-σ contours, we use colour coding to indicate the dominant DM mechanisms, as discussed in Sect. 2.4, for the parameter sets that minimize χ 2 at each point in the plane.

Squarks and gluinos
The top row of plots in Fig. 2 show (mq , mg) planes, where mq is an average over the masses of the left-and right-handed first-and second-generation squarks, which are very similar underlying p-value of the likelihood distribution. This was confirmed by a study in the last paper in [21][22][23][24][25][26][27][28][29][30][31], which compared for different scenarios the naive p-value calculation with that obtained from toys. in the pMSSM11. 22 In the top left panel, where (g − 2) μ is included, we see 95% CL lower bounds mq 2000 GeV and mg 1400 GeV, with regions favoured at the 68% CL appearing at slightly larger masses. We note that the best-fit point, denoted by the green star, is at large mq > 4000 GeV and mg ∼ 3900 GeV. The full set of pMSSM parameter values at this point, as well as the value of the global χ 2 function, are listed in the second column of Table 4. Important sparticle production cross-sections and decay modes at this best-fit point are shown in the top panel of Table 5.
Within the 2-σ contour, the dominant DM mechanism is slepton coannihilation, with stau coannihilation also playing a role for mq ∼ 2.5 TeV, andχ ± 1 coannihilation playing a role at mg ∼ 1500 GeV and when mg 2500 GeV and mq 2800 GeV. Finally, we observe that at the 3-σ level much smaller values of mq are allowed, and that there is also a peninsula at small mg and larger mq that appears at the same level. These regions avoid the LHC exclusion searches in virtue of the same mechanisms which allow lower masses when the (g − 2) μ constraint is not applied and which will be described more in detail below. However, they are not able to satisfy the (g − 2) μ and this is why they take a χ 2 11 penalty which makes them allowed only at 3-σ .
We also note a 'nose' feature corresponding to a reduction in the lower bounds when mq ∼ 2.2 TeV and 0 < mq − mg 200 GeV. We have verified that this is due to a loss of search sensitivity whenq R →g + q, the q jet is soft, and g → qq +χ * , whereχ * denotes any electroweak ino other than the LSP, compared to a high sensitivity forq R → qχ 0 third column of Table 4. The upper panel of Fig. 3 displays relevant sparticle masses and the most important sparticle decay chains at this point, and numerical values are given in the second panel of Table 5. We see that the right-handed squarks decay into a variety of final states involving heavier neutralinos and charginos via intermediate gluinos due to mg < mq , reducing the effectiveness of / E T -based searches in this 'nose' region, compared to simpleq → q +χ 0 1 decays. We see significant differences in the top right panel where (g − 2) μ is dropped. The best-fit in this case is close to the 68% CL boundary at (mq , mg) ∼ (1000, 1600) GeV, with the parameters and χ 2 value shown in the fourth column of Table 4. As we discuss later, BR(B s,d → μ + μ − ) and the DM density constraint play important roles in preferring a relatively low value of mq . The dominant particle production and decay modes for this best-fit point are shown in the third panel of Table 5. It is notable that the 95% CL lower limits on mq and mg are reduced to ∼ 1000 GeV, and a lesspronounced 'nose' feature now appears when mq ∼ 1 TeV and 0 < mg − mq 200 GeV. Again, we have verified that this reflects a loss of search sensitivity wheng →q +q, theq jet is soft, andq → q +χ * (χ 0 1 ), whereχ 0 1 is much heavier than in the fit with (g − 2) μ (for which a large SUSY contribution requires mχ0 1 to be small), since the direct decaỹ g → qqχ * (χ 0 1 ) in the mq > mg case is more sensitive than the above cascade decay in the compressed spectrum. The lower panel of Fig. 3 shows the most important sparticle decay chains at the representative point in this region whose parameters are listed in the fourth column of Table 4, and the numerical values of branching ratios are given in the bottom panel of Table 5.
The differences between the fits with and without the (g − 2) μ constraint are driven primarily by the fact that Fig. 1 Higgs and sparticle spectra for the best-fit points for the pMSSM11 with (top) and without the (g − 2) μ constraint (bottom), showing also decay paths with branching ratios > 5%, the widths of the lines being proportional to the branching ratios. These plots were prepared using the code presented in [160] the fit with (g − 2) μ prefers small mχ0 1 , in which case the LHC 13-TeV searches require large mq and mg, whereas the fit without (g − 2) μ favours a region with larger mχ0 1 . In this case, the loss of search efficiency due to a compressed spectrum allows mq and mg to be smaller than in the fit with (g − 2) μ . As we see later, in this compressed region the LSP is mainly a neutral Higgsino, and coannihilations with a nearby charged Higgsino and theχ 0 2 are important in determining the relic neutralino density. Coannihilations with first-and second-generation squarks are also relevant here and in a band with mq ∼ 1 TeV mg (coloured cyan), whereas coannihilations with gluinos are important along a band with (1 TeV, 2 TeV) mg mq (coloured magenta). In this plane the 1-, 2-and 3-σ contours lie relatively close to each other.
In the middle row of Fig. 2 we display the corresponding (mq , mχ0 1 ) planes. We see a preference for mχ0 1 550 GeV in the left panel, where the (g − 2) μ constraint is included, whereas much larger values of mχ0 1 are allowed at the 3-σ level. These larger values of mχ0 1 appear within the 1-and 2-σ contours in the middle right panel where the (g − 2) μ constraint is dropped. We also see again that larger values of mq are favoured when (g − 2) μ is included, whereas a small mq − mχ0 1 mass difference is preferred when the (g − 2) μ constraint is dropped. In both the middle panels the dominant DM mechanisms are slepton andχ ± 1 coannihilation, with the rapid annihilation via the heavy H/A Higgs bosons becoming important at large masses when (g − 2) μ is dropped. Similar features are seen in the (mg, mχ0 1 ) planes displayed in the bottom row of Fig. 2. ) planes (bottom panels), including the (g − 2) μ constraint (left panels) and dropping it (right panels)   ) planes in the lower panels, again including the (g −2) μ constraint in the left panels and dropping it in the right panels. We see that both the third-generation squark masses may be considerably smaller than those in the first two generations. Specifically, an isolated, low stop-mass region where (mt 1 , mχ0 1 ) ∼ (500, 300) GeV is allowed at the 95% CL 23 in both the cases with and without (g − 2) μ , which is connected in the latter case to the rest of the 95% CL region at the 3-σ level. The low stop-mass island is allowed and defined by different physics mechanisms. First, the third-generation- 23 For relatively low stop masses, large values of X t /M S A t /M S √ 6 are required to avoid tension with the Higgs mass measured at the LHC. Constraints from Charged and Color Breaking (CCB) minima can be relevant [162] in such a case, but we have not taken these in account [163] in our analysis. This is because our best-fit point region is characterized by relatively small values of A t /M S , as it can be seen from Table 4, for which this issue is not relevant. squark spectra are sufficiently compressed to allow the points to bypass the LHC13 constraints. Moreover, it is characterized by compressed-slepton spectra as well, which explains the fact that the region is shaded in yellow in the plots. We also note that it can not be extended to lower stop masses because otherwise it would be disallowed by sbottom searches, since in our scenario the masses of the stop and sbottom squarks are defined by a single soft SUSY-breaking mass term and the sbottoms would not be sufficiently compressed to be allowed by LHC searches. LHC constraints also limit its extensions in the direction of lower neutralino (too light third-generation squarks) or higher stop masses (due to the loss of compression). Finally, at heavier neutralino masses slepton coannihilation is insufficient to reduce the relic density into the allowed range.4 When (g−2) μ is dropped, extended 95% CL regions with mχ0 1 500 GeV appear when mt 1 1100 GeV and mb 1 1250 GeV. When (g − 2) μ is included, there are extended regions with mχ0 1 500 GeV that appear at the 3-σ level. Within the 1-and 2-σ contours, the dominant DM mechanisms are slepton andχ ± 1 coannihilation, with rapid annihilation via the heavy H/A Higgs bosons again becoming important at large mχ0 1 when (g − 2) μ is dropped. The same mechanism is also active inside the white regions between 800 GeV (1 TeV) mt 1 (mb 1 ) 1.1 (1.2) TeV and 400 GeV mχ0 1 600 GeV, the blue shading being absent due to the proxy-measure being not sufficiently descriptive in this parameter space region. Stop and sbottom coannihilation are also important for small mt 1 − mχ0 1 and mb 1 − mχ0 1 .

Sleptons
As was to be expected, there are large differences between the (mμ R , mχ0 1 ) planes with and without the (g − 2) μ constraint, shown in the upper panels in Fig. 5. We see in the upper left plane a preference for mμ R 550(750) GeV and mχ0 1 500(550) GeV at the 68 (95)% CL, enforced by the (g − 2) μ constraint, with larger masses allowed at the 3-σ level. There is also a 68% CL region with similar ranges of mμ R and mχ0 1 in the case without (g − 2) μ (upper right panel), but the 95% CL region extends to much larger values of mμ R and mχ0 1 , and there is also a second, extended 68% CL region that is separated by a band of points with only slightly higher χ 2 . In both these plots, we see a very narrow strip where slepton-χ 0 1 coannihilation is important, whereasχ ± 1 coannihilation dominates in most of the regions allowed at the 95% CL, supplemented by annihilation via the H/A bosons at large mχ0 1 when (g −2) μ is dropped. We do not display the corresponding (mμ L , mχ0 1 ) and (mẽ L ,R , mχ0 1 ) planes, which are very similar because we impose universality on the soft SUSY-breaking masses of the first two slepton generations.
However, in the pMSSM11 the soft SUSY-breaking stau masses are allowed to be different, with the result seen in the lower panels of Fig. 5 that large values of mτ 1 are allowed at the 68 and 95% CL even when (g − 2) μ is imposed. The 1 TeV at larger tan β. We also note that the ) planes (upper panels) and the (mτ 1 , mχ0 1 ) planes (lower panels), including the (g − 2) μ constraint (left panels) and dropping it (right panels) dominant DM mechanisms display significant differences. Chargino coannihilation is important in both planes, but slepton coannihilation appears only in the case where (g − 2) μ is included. In this case annihilation via the H/A poles appears only when M A 1 TeV, but it appears also at larger M A when (g − 2) μ is dropped. We see in both cases a limited region with M A ∼ 2 TeV and tan β 10 where stau coannihilation dominates. In our previous pMSSM10 analysis [13] the interplay of the LHC electroweak searches, (g − 2) μ and the DM constraints, heavily relying on the fact that only one independent slepton mass parameter was allowed, led to a region with 25 tan β 45 being preferred at the 68% CL. However, in the pMSSM11, dropping the restriction mτ = m˜ now allows values of tan β < 5 for a wide range of M A values. Also, despite the updated (stronger) con-straints on H/A → τ τ , values down to M A ∼ 500 GeV are still allowed at the 95% CL.

One-dimensional likelihood functions
In this section we present the profile χ 2 likelihood functions corresponding to various one-dimensional projections of the results from our global fits, again comparing those with and without the (g − 2) μ constraint. In the following series of plots, results including the LHC 13-TeV constraints are shown as solid lines, and those using only 8-TeV results are shown as dashed lines. Results obtained including (g − 2) μ are shown in blue and those obtained without (g − 2) μ are shown in green. ) planes (upper panels) and the (M A , tan β) planes (lower panels), including the (g − 2) μ constraint (left panels) and dropping it (right panels) 4.1 (g − 2) μ As a preliminary, Fig. 7 shows the one-dimensional profile likelihood functions for (g − 2) μ with (blue) and (green) without applying the (g − 2) μ constraint a priori. Comparing the solid and dashed lines, we see very little difference between the results using and discarding the LHC 13-TeV data. The results including (g − 2) μ (blue lines) largely reflect our implementation of the (g − 2) μ constraint shown in Table 2. Interestingly, when this constraint is not applied a priori (green lines), whilst a very small SUSY contribution to (g −2) μ is preferred, a wide range of values of (g −2) μ are found to be allowed at the χ 2 ∼ 2 level and the experimental value can be accommodated at the 1.5-σ level. Although the other data certainly do not favour a large SUSY contribution to (g − 2) μ , neither do they exclude it.

Squarks and gluinos
The profile likelihood functions for squarks and gluinos are shown in Fig. 8. The left panel is for mq , where we see that when the 13-TeV LHC data and (g − 2) μ constraint are included (solid blue line), there is a monotonic decrease in χ 2 as mq increases, with mq 1.9 TeV at the 95% CL (horizontal dotted line). This constraint is much stronger than that obtained with 8-TeV data alone (dashed blue and green lines): mq 1.0 TeV at the 95% CL. In particular, the 13- Fig. 7 One-dimensional profile likelihood functions for (g − 2) μ in the pMSSM11, with (blue) and without (green) applying the (g − 2) μ constraint a priori and with (solid) and without (dashed) applying the constraints coming from the LHC run at 13 TeV. Also shown as a dotted line is the experimental constraint [108,109], taking into account the theoretical uncertainty [100][101][102][103][104][105][106][107] within the Standard Model TeV data exclude a squark coannihilation strip that had been allowed by the 8-TeV data. When (g − 2) μ is dropped but the 13-TeV data retained (solid green line), the χ 2 function exhibits a global minimum at mq ∼ 1 TeV, with a plateau at χ 2 1.5 at larger mq . Important roles in the location of this global minimum are played by the BR(B s,d → μ + μ − ) constraint as discussed in Sect. 4.4, whose contribution to the global χ 2 function at this point is ∼ 1.1 lower than at large mq , and by the relic DM density constraint, which is satisfied thanks to multiple coannihilation processes as discussed in Sect. 4.6.
In the right panel of Fig. 8 for mg, we see that with both the LHC 13-TeV data and (g−2) μ included mg 1.8 TeV (solid blue line), whereas without (g − 2) μ we find mg 1.0 TeV (solid green line). On the other hand, in the absence of the LHC 13-TeV data (dashed lines), mg 500 GeV would have been allowed at the 95% CL, whether (g − 2) μ is included, or not. The LHC 13-TeV run has excluded a region of gluino coannihilation that was allowed by the 8-TeV data. Third-generation squarks An analogous pair of plots showing the profile likelihood functions for the masses of thet 1 andb 1 are shown in the left and right panels of Fig. 9. When the LHC 13-TeV data are included we see in the left panel a well-defined local minimum of the χ 2 function in a compressed-stop region with χ 2 ∼ 2.3 for mt 1 ∼ 400 GeV. This is followed by a local maximum that exceeds χ 2 > 9 for mt 1 ∼ 800 GeV when (g − 2) μ is included (solid blue line) but is lower when (g − 2) μ is dropped (solid green line). This is followed in both cases by a monotonic decrease for larger mt 1 and a global minimum of χ 2 for mt 1 ∼ 1800 GeV.
In the case of mb 1 (right panel of Fig. 9). when the 13-TeV LHC data and (g −2) μ are included (solid blue line) there are some irregularities in the χ 2 function for mb 1 ∼ 1000 GeV, but no hint of a compressed-sbottom region when (g − 2) μ is dropped (dashed blue line). Comparing with the situation when only LHC 8-TeV used, we see that the 13-TeV data have increased significantly the pressure on scenarios with mb 1 1.5 TeV. At larger masses the χ 2 functions mb 1 are very similar to those for mt 1 , whether (g − 2) μ is included or not. Sleptons Figure 10 displays analogous plots of the profile likelihood functions for mμ R (left panel, those for mμ L and mẽ L ,R are very similar) and mτ 1 (right panel, that for mτ 2 is quite similar). When the (g − 2) μ constraint is implemented (blue lines), the χ 2 function for mμ R exhibits the expected welldefined minimum at mμ R ∼ 200 to 500 GeV when the LHC 13-TeV data are included. In the absence of the (g − 2) μ constraint (green lines), this is replaced by a plateau with χ 2 ∼ 2 that extends to mμ R ∼ 900 GeV, where the profile likelihood function drops to very small values for larger mμ R . The drop occurs because this fit prefers mχ0 1 ∼ 900 to 1000 GeV, and any heavierμ R can decay into aχ 0 1 in this mass range.
We see in the right panel of Fig. 10 that when (g − 2) μ is included (blue lines) the profile likelihood function for mτ 1 is quite different from that for mμ R , thanks to the decoupling between their soft SUSY-breaking masses in the pMSSM11. The χ 2 function falls monotonically to a local minimum when mτ 1 ∼ 300 GeV and remains small for larger mτ 1 , whether the LHC 13-TeV data are included (solid line), or not (dashed line). However, when (g − 2) μ is dropped (green lines), the profile likelihood function for mτ 1 is quite similar to that for mμ R , also exhibiting a plateau with χ 2 ∼ 2 and falling to small values for mτ 1 900 GeV when the LHC 13-TeV data are included. This feature appears because, in order to avoid a charged LSP, a smaller value of mτ 1 would require a smaller value of mχ0 1 , which is disfavoured as seen in the left panel of Fig. 11 and discussed below. Figure 11 shows the profile likelihood functions for the lightest neutralinoχ 0 1 (left panel) and the lighter charginõ χ ± 1 (right panel). When the (g − 2) μ constraint is applied (blue lines), the χ 2 function for mχ0 1 including 13-TeV data exhibits a well-defined but broad minimum at mχ0 1 ∼ 100 to 400 GeV. This preference for small mχ0 1 was already seen in the upper boundaries of the 68% and 95% CL regions in the planes involving mχ0   On the other hand, when the (g−2) μ constraint is dropped (green lines) we see a preference for mχ0 1 ∼ 950 GeV. Despite the fact that the LSP is a nearly-pure Higgsino at this best-fit point, this mass of ∼ 950 GeV is below the ∼ 1.1 TeV mass expected for a Higgsino dark matter candidate. This arises because, at the best-fit point, several of the squark masses lie close to the LSP mass, making multiple coannihilation important. Due to the relatively large number of states with masses close to the Higgsino, their density actually increases the final LSP relic density, 24 thereby pushing the mass of the Higgsino below its nominal ∼ 1.1 TeV value.

Electroweak inos
Turning now to the profile likelihood functions for the lighter charginoχ ± 1 (right panel of Fig. 11), we see that when (g − 2) μ is taken into account (blue lines) the χ 2 function also features a well-defined minimum for mχ± 1 ∼ 200 to 500 GeV (that forχ 0 2 is very similar), reflecting the importance ofχ ± 1 −χ 0 1 coannihilation. This minimum is followed by a rise to a local maximum at mχ± 1 ∼ 600 GeV, which is more pronounced when the 13-TeV data are included (solid blue), followed by a slow decrease as mχ± 1 increases further. When the (g − 2) μ constraint is dropped and the LHC 13-TeV data are included (solid green line), the χ 2 functions for mχ±

Neutralino composition
It is interesting also to examine the profile likelihood functions for the amplitudes N 1i characterizing theχ 0 1 composition: which are shown in Fig. 12, again for the analysis with the 13-TeV data as solid lines and without them as dashed lines, and with (g − 2) μ as blue lines and without it as green lines.
The top left panel shows that, when (g − 2) μ is included, an almost pureB composition of theχ 0 1 is preferred, N 11 → 1, though the possibility that this component is almost absent is also allowed at the level χ 2 ∼ 4. On the other hand, when the constraint from (g − 2) μ is removed, there is a mild ( χ 2 ∼ 1) preference for N 11 → 0. The reason for this is again the preference for a largeH u,d components in the latter case, where the neutralino mass is allowed to be larger, due to flavor constraints slightly favoring a 1 TeV neutralino as a solution to the observed DM relic density. The upper right panel shows that a smallW 3 component in theχ 0 1 is preferred in all cases. 25 Finally, the lower panel confirms that small H u,d components are preferred by χ 2 4 when (g − 2) μ is included, whereas there would have been a preference for these components to dominate in the absence of the (g − 2) μ constraint. Figure 13 displays information about the preferred and disfavouredχ 0 1 compositions in two triangular panels. Both are for fits including LHC 13-TeV data (those dropping these data are quite similar), the left panel includes the (g − 2) μ constraint and the right panel drops it. The χ 2 for the bestfit points at each location in the triangles are colour-coded as indicated. We see in the left panel that in the case with (g − 2) μ a small Wino fraction N 2 12 < 0.1 is strongly favoured, while the relative proportions of the Bino fraction N 2 11 and the Higgsino fraction N 2 13 + N 2 14 are relatively unconstrained at the 95% CL. On the other hand, the right panel shows that almost all binary combinations of Bino, Wino and Higgsino (along the edges of the triangle) are allowed at the 95% CL, but three-way mixtures (in the interior of the triangle) are strongly disfavoured. Table 6 compares the composition of the LSPχ 0 1 found at the best-fit points in our present pMSSM11 analysis based on LHC 13-TeV data (with and without the (g − 2) μ constraint) with the composition at the best-fit point from our previous pMSSM10 analysis that also applied the (g − 2) μ constraint [13]. We see that both the pMSSM11 and pMSSM10 analyses with (g − 2) μ prefer an almost pureB composition. On the other hand, when the (g − 2) μ constraint is dropped the pMSSM11 analysis prefers an almost equal mixture of H u andH d components with a small admixture ofB and again a very small admixture ofW 3 because we only scan m˜ and mτ , hence mχ0 1 < 2 TeV. Table 6 also displays the composition of the second-lightest neutralino,χ 0 2 , and we see that its content is mainlyW 3 in the fit to the pMSSM11 with (g − 2) μ and in the pMSSM10 fit, but is mainly Higgsino in the fit to the pMSSM11 without (g − 2) μ . Figure 14 displays the one-dimensional profile likelihood functions for BR(B s,d → μ + μ − ) in the pMSSM11 (left panel) and the BR(B s → X s γ ) branching ratio (right panel), with and without the LHC 13-TeV data and the (g −2) μ constraint. We see in the left panel that a value of BR(B s,d → μ + μ − ) close to the SM value is preferred if both these 25 This is because we only scan m˜ and mτ < 2 TeV, hence mχ0  constraints are applied, though deviations at the level of ± ∼ 10% are allowed at the level of χ 2 = 4 (2 σ ), corresponding to the 95% CL. On the other hand, if (g − 2) μ is dropped, a larger range of BR(B s,d → μ + μ − ) is allowed, with a larger deviation at the level of ± ∼ 30% becoming allowed at the level of χ 2 = 4. In particular, when the LHC13 data are included but (g − 2) μ is dropped, the global χ 2 function is minimized at a value of BR(B s,d → μ + μ − ) below the SM value, as hinted by the present experimental data, with the SM value being mildly disfavoured by χ 2 1. It will be interesting to see how measurements of BR(B s,d → μ + μ − ) evolve.

B-physics observables
The analogous curves for BR(B s → X s γ ) in the right panel of Fig. 14 show preferences for values close the SM predictions, with 2 σ ranges that are ±20%. Discriminating between the SM and the pMSSM11 would require significant reductions in both the theoretical and experimental uncertainties in BR(B s → X s γ ).
As already mentioned in Sect. 2.3, the LHCb Collaboration has recently announced the first experimental measurement of τ (B s → μ + μ − ), which is related to the quantity A that takes the value +1 in the SM, but may be different in a SUSY model such as the pMSSM11. Figure 15 displays the profile likelihood functions for A (left panel) and τ (B s → μ + μ − )/τ B s (right panel), in our pMSSM11 fits with and without the LHC 13-TeV data and (g − 2) μ . We restrict our attention to positive values of A , corresponding to τ (B s → μ + μ − )/τ B s > 0. 94. We see that all the fits favour values of A close to unity, with that dropping both the LHC 13-TeV data and (g − 2) μ allowing the widest range. Values of τ (B s → μ + μ − )/τ B s close to unity are also favoured, with χ 2 9 for τ (B s → μ + μ − )/τ B s = 0.94. Table 6 The amplitudes characterizing the decomposition of the LSP χ 0 1 and of theχ 0 2 into interaction eigenstates at the best-fit points in our present pMSSM11 analysis including LHC 13-TeV data, with and without the (g−2) μ constraint, compared with the composition at the best-fit point found in our previous pMSSM10 analysis that also included the (g − 2) μ constraint, but only LHC 8-TeV data [13] Model StateBW The new LHCb measurement [112] does not challenge any of these model predictions. Figure 16 shows similar plots of M h (upper left panel), and of the ratios of the branching ratios for h → γ γ, Z Z * and h → gg (treated as a proxy for σ (gg → h)) to their values in the SM in the upper right, lower left and lower right panels, respectively. Taking into account the theoretical uncertainties in the calculation of M h in a supersymmetric model [66][67][68][69][70][71], which we take to be ±3 GeV, 26 there is no tension with the global fits. These also favour values of the decay branching ratios that are similar to those in the SM whether (g − 2) μ is included in the fit, or not, though with uncertainties that are typically ± ∼ 20%. As discussed in [34], the global combination of ATLAS and CMS measurements using LHC Run-1 data has significantly larger uncertainties.

Dark matter measures
In Sect. 2.4 we introduced various possible mechanisms for bringing the relicχ 0 1 density into the range allowed by Planck and other data, proposing measures of their prospective importance that we portrayed using different colours in the two-dimensional parameter planes shown in Sect. 3. We emphasized there and in the subsequent discussions of onedimensional profile likelihood functions earlier in Sect. 4 the roles played by certain of these DM mechanisms. In this subsection we display profile likelihood functions for the most interesting of these DM measures, discussing the χ 2 levels at which they become relevant. As in the previous sections, we compare results for the analysis in which the (g − 2) μ constraint is applied with those when (g − 2) μ is discarded. Figure 17 displays the profile likelihood functions for the selected DM measures. The top left panel shows the first-and second-generation slepton measure, and we see that χ 2 is generally small throughout this region. Theτ 1 measure is shown in the top right panel, and we see that with (g − 2) μ included, whether or not the LHC 13-TeV results are included, the χ 2 function has a shallow minimum within the region where this mechanism may dominate (shown as the vertical pink band), but very small values of theτ 1 coan- nihilation measure are disfavoured, and larger values of this measure also appear with a negligible likelihood price. On the other hand, when (g − 2) μ is dropped we find that χ 2 ∼ 2 is almost independent of mτ 1 /mχ0 The χ 2 function rises as mτ 1 /mχ0 1 → 1 when (g − 2) μ is included, because this constraint prefers small values of mχ0 1 , for which the relic density constraint cannot be satisfied when mτ 1 /mχ0 1 → 1. However, since the first-and second-generation slepton masses are independent of mτ 1 in the pMSSM11 there is no such obstacle disfavouring mμ R /mχ0 1 → 1. Therefore the profile χ 2 function for the first-and second-generation DM measure does not rise in this limit, as seen in the top left panel of Fig. 17.
In the case of theχ ± 1 coannihilation measure shown in the middle left panel of Fig. 17, we see that the best-fit pMSSM11 points lie within this shaded band, whether the LHC 13-TeV data and/or (g − 2) μ are included or not. In the case with (g − 2) μ , the best-fit point has mχ± 1 /mχ0 1 ∼ 1.1 whether the LHC 13-TeV data are included or not, whereas when (g−2) μ is dropped there is a strong preference for mχ± 1 /mχ0 1 close to unity, which is possible in the case because the LSP is Higgsino-like. As in the case of theτ 1 DM measure, the relic density constraint disfavours mχ± 1 /mχ0 1 → 1 when (g − 2) μ is included. We find some parameter sets with mχ± 1 − mχ0 1 10 MeV that have χ 2 4, which occur when M 1 is negative, near the border of a region where the LSP would be theχ ± 1 . In the case of the A/H measure shown in the middle right panel, we see that χ 2 > 3 in this region when the (g − 2) μ and LHC 13-TeV constraints are both used. However, the χ 2 price of rapid annihilation through the A/H poles is reduced if either of these constraints is dropped. Indeed, including the (g − 2) μ constraint forces the neutralino mass to be at most 500 GeV, in which case the funnel condition implies an upper bound on M A 1 TeV, well within the reach of LHC 13-TeV searches for tan β 15.
The bottom left panel of Fig. 17 displays the profile likelihood function for the squark coannihilation measure mq L /mχ0 1 − 1. We see that before the LHC-13 data the bestfit point with (g − 2) μ included was in the squark coannihilation region with mq L /mχ0 1 < 1.1, though this feature was absent when (g − 2) μ was dropped. Including the LHC 13-TeV data, the best-fit points with and without (g − 2) μ have mq L mχ0 1 , but there is still a vestige of the squark coannihilation region with χ 2 < 4 when (g −2) μ is dropped. The reason for this is that lifting the (g − 2) μ constraint allows for a heavier neutralino, which in turn implies heavier squark masses still allowed by LHC-13 TeV data Finally, the bottom right panel of Fig. 17 shows the gluino coannihilation measure, and we see that this may also play a role when χ 2 < 4, unless both the LHC 13-TeV data and (g − 2) μ are included.

NLSP lifetimes
We display in Fig. 18 the one-dimensional profile likelihood for the NLSP lifetime, τ NLSP , including all possible NLSP species. There is little difference between the χ 2 functions with (g − 2) μ , whether or not the LHC 13-TeV data are included (blue curves). In both cases, we find that χ 2 4 for τ NLSP 10 −10 s. On the other hand, when the (g − 2) μ constraint is dropped (green curves), we see that values of τ NLSP 10 3 s are allowed at the χ 2 4 level, again whether or not the LHC 13-TeV data are included (green curves). As already mentioned, we exclude from our scan parameter sets with NLSP lifetimes exceeding 10 3 s, as they could alter the successful predictions of standard Big Bang nucleosynthesis [153][154][155][156][157][158][159].
The upper panels of Fig. 19 display the χ 2 distributions for chargino (left) and stau lifetimes (right) between 10 −7 s and 10 3 s, for the fits omitting (g − 2) μ (fits including (g − 2) μ give χ 2 outside the displayed range). We see that, whereas shorter lifetimes are favoured, lifetimes as long as 10 3 s are allowed at the 95% CL for both sparticle species when (g − 2) μ is dropped, whether or not the LHC 13-TeV data are included. The lower panels of Fig. 19 display the corresponding mass-lifetime planes for the chargino and stau. We see that a long-lived chargino would have a mass mχ± 1 ∼ 1.1 TeV, and a long-lived stau would have a mass mτ 1 ∼ 1.5 TeV, both beyond the reaches of current LHC searches for long-lived charged particles. We have also checked the possible lifetimes of other NLSP candidates, finding that squarks and gluinos generally have lifetimes 10 −17 (10 −10 ) s at the 95% CL in fits including LHC 13-TeV with (without) the (g − 2) μ constraint, with just a few points having longer lifetimes. Hence they also do not offer good prospects for LHC searches for long-lived particles.

Spin-independent scattering cross section
We now discuss the prospects for direct detection ofχ 0 1 DM via spin-independent elastic scattering. Fig. 20 shows (mχ0 1 , σ SI p ) planes, including the LHC 13-TeV data, with (left panel) and without (right panel) the (g − 2) μ constraint. The values of σ SI p displayed are the nominal values calculated using the central values of the matrix elements in the SSARD code. The pale green shaded region is that excluded by the combined LUX [4], XENON1T [6] and PandaX-II [3] limit, which is shown as a solid black line. 27 The yellow shaded region lies below the neutrino 'floor', which is shown as an orange dashed line. We see that mχ0 1 100 GeV in both the cases with and without the (g − 2) μ constraint, with upper limit mχ0 1 550 at the 95% CL when (g − 2) μ is included. When this constraint is dropped, the 95% CL range extends up to 2 TeV, the upper limit for which our analysis is applicable, because we have limited our scan to slepton masses ≤ 2 TeV.
We see that the nominal prediction for σ SI p at the bestfit point is at the level of the sensitivities projected for the planned LUX-Zeplin (LZ) and XENON1T/nT experiments (solid purple line) when the (g − 2) μ constraint is Fig. 19 Upper panels: One-dimensional profile likelihood plots for the lifetime of theχ ± 1 (left) and theτ 1 (right). Lower panels: The corresponding mass-lifetime planes for theχ ± 1 andτ 1 , with the 95% CL regions shaded according to the dominant DM mechanisms dropped, and somewhat higher if (g − 2) μ is included. However, we emphasize that there are considerable uncertainties in the estimate of σ SI p , which are reflected in the fact that the range of nominal SSARD predictions extends above the current combined limit from the LUX [4], XENON1T [6] and PandaX-II [3] experiments. There is no incompatibility when the uncertainties in the σ SI p estimate are taken into account. The 68 and 95% CL ranges of the nominal values of σ SI p extend slightly below the neutrino 'floor' in the case with (g − 2) μ included, and much lower in the case where (g − 2) μ is dropped. In both cases, large values of σ SI p occur in the chargino coannihilation region (green shaded area), with other DM mechanisms including squark coannihilation yielding large values of σ SI p for mχ0 1 1 TeV. However, this and the other DM mechanisms indicated also allow much smaller values of σ SI p . As in the case of the pMSSM10 studied in [13], we expect that points with very small values of σ SI p would, in general, have similarly small values for the spin-independent scattering cross section on neutrons.
4.9 Spin-dependent scattering cross section Figure 21 displays the corresponding planes of (mχ0 1 , σ SD p ) with (left panel) and without (right panel) the (g − 2) μ constraint applied. Here the neutrino 'floor' is taken from [170]. As in the σ SI p case, we see that the allowed ranges of mχ0 1 extend from ∼ 100 GeV to ∼ 550 GeV when (g − 2) μ is included and up to the sampling limit of 2 TeV when (g −2) μ is dropped. The uncertainties in the calculation of σ SD p are significantly smaller than those for σ SI p , and we see that the ranges of the 68 and 95% regions in the nominal σ SD p calculations lie below the upper limit from the PICO experiment [5] (solid purple line). In both the left and right panels, the nominal predictions for the best-fit points lie some ∼ 3 orders of magnitude below the current PICO limit. For completeness, we also show the upper limits from SuperKamiokande [171] and IceCube [141] searches for energetic solar neutrinos, assuming that the LSPs annihilate predominantly into τ + τ − (which is not always the case in the pMSSM11) and neglect-  [4], XENON1T [6] and PandaX-II [3] Collaborations are shown as green, magenta and blue contours, respectively, and the combined limit is indicated by a black line with green shading above. The projected future 90% CL exclusion sensitivities of the LUX-Zeplin (LZ) [168] and XENON1T/nT [169] experiments are shown as solid purple and dashed blue lines, respectively, and the neutrino background 'floor' is shown as a dashed light-blue line with a shading of the same colour below  [85]. The upper limit established by the PICO Collaboration [5] is shown as a purple contour, with green shading above. The neutrino 'floor' for σ SD p is taken from [170]. We also show the indicative upper limits from SuperKamiokande [171] and IceCube [141] searches for energetic solar neutrinos obtained assuming that the LSPs annihilate predominantly into τ + τ − , which are subject to the caveats discussed in the text ing the uncertainties in interpretation mentioned earlier: see the discussion in the following section.
We see in the left panel of Fig. 21 (when (g − 2) μ is included) that points with chargino coannihilation as the dominant DM mechanism yield nominal predictions for σ SD that extend over many orders of magnitude below the current PICO limit and well below the τ + τ − floor. Points for which slepton coannihilation is the dominant DM mechanism do not reach so close to the PICO limit, but may also lie many orders of magnitude below it. We see in the right panel (when (g − 2) μ is dropped) similar ranges of nominal σ SD p values. We also see that when mχ0 1 1 TeV many competing DM mechanisms come into play, and may give small values of σ SD p . However, in the case of squark coannihilation σ SD p may lie within ∼ 3 orders of magnitude of the PICO upper limit.

Indirect astrophysical searches for dark matter
We have explored the possible impact of indirect searches for DM via annihilations into neutrinos inside the Sun. If the DM inside the Sun is in equilibrium between capture and annihilation, the annihilation is quadratically sensitive to the local Galactic DM density. However, as discussed earlier, equilibrium is not always a good approximation. We note also that the capture rate is not determined solely by spindependent scattering on protons in the Sun, but also depends on the amount of spin-independent scattering on Helium and heavy nuclei. As we have seen, the σ SI p matrix element is more uncertain than that for σ SD p , and this uncertainty should be propagated into the constraint on σ SD p . Finally, we note that the greatest sensitivity of the IceCube search for energetic neutrinos from the Sun [141] is for annihilations into τ + τ − and W + W − , which are not always the dominant final states in the pMSSM11 models of interest.

Fig. 23
Higgs and sparticle spectrum for the pMSSM11 with and without the (g − 2) μ constraint applied (upper and lower panels, respectively). The values at the best-fit points are indicated by blue lines, the 68% CL ranges by orange bands, and the 95% CL ranges by yellow bands Using the nominal values of the matrix elements from SSARD and neglecting the astrophysical uncertainties, we have calculated the signals in the IceCube detector for a subset of our pMSSM11 points that are consistent with the PICO constraint [5]. We find that the IceCube W + W − constraint [141] has negligible impact on these parameter sets, and that only a fraction are affected by the IceCube τ + τ − constraint. In view of this and the uncertainties in the interpretation of the IceCube searches, we have not included them in our fits.

Impacts of the LHC 13-TeV and new direct detection constraints
In this section we illustrate the impact of the LHC 13-TeV data and the recent updates from the Xenon-based direct detection experiments LUX, XENON1T, and PandaX-II [3,4,6] on relevant pMSSM11 parameter planes. In the left panel of Fig. 22 we display the impact of the new results on the (mq , mg) plane: the solid red, blue and green lines are the current 68%, 95% and 99.7% CL contours, and the dashed lines are those for the corresponding 68, 95% and 99.7% CL contours in a global fit omitting the LHC 13-TeV constraints and those from the Xenon-based direct detection experiments. The right panel of Fig. 22 makes a similar comparison Fig. 24 The χ 2 pulls at the best-fit points in the pMSSM11 including (left) and without the (g − 2) μ constraint (right). In the rightmost plot, the χ 2 pull from (g − 2) μ is shown (hatched orange bar), but its penalty is not included in the fit of the 68, 95 and 99.7% CL regions in the (mχ0 1 , σ SI p ) plane found in global fits including LHC 13-TeV and Xenon-based detector data (solid lines) and omitting these data (dashed lines).
We see in the upper left panel of Fig. 22 that the LHC 13-TeV constraints exclude bands of parameter space at low mq and mg, disallowing in particular a squark coannihilation region at mq ∼ 500 GeV and large mg and a gluino coannihilation strip at mg ∼ 500 GeV that were allowed by the LHC 8-TeV data. The impact on the gluino and squark coannihilation strips can also be appreciated from the upper right and lower left panels, where they appear as dashed-blue islands along the diagonal where the mass is degenerate with the neutralino that disappear completely after the inclusion of the LHC 13-TeV constraints. The bottom right panel of Fig. 22 shows that low values of σ SI p that would have been allowed in a fit without the LHC 13-TeV data are now disallowed. This effect is in addition to the downwards pressure on σ SI p exerted by the new generation of Xenon-based direct detection experiments.

Best-fit points, spectra and decays
Following our previous discussions of some two-dimensional projections of the pMSSM11 parameter space and various one-dimension profile likelihood functions, we now discuss in more detail the best-fit points in the pMSSM11 fits incorporating the LHC 13-TeV data, both with and without the (g − 2) μ constraint, whose input pMSSM11 parameter values were given in the first and third columns of Table 4. We note, however, that the likelihood functions are very flat for larger masses, so these best-fit points should not be taken as definite predictions. Figure 1 displays the spectra of Higgs bosons and sparticles at the best-fit points for the pMSSM11 including (upper panel) and excluding (lower panel) the (g − 2) μ constraint. 28 In each case we also show the decay paths with branching ratios > 5%, the widths of the lines being proportional to the branching ratios. The heavier Higgs bosons H, A, H ± , are lighter in the case without (g − 2) μ , whereas the sleptons and the electroweak inos are heavier. The branching ratio patterns differ in the two cases, with the Higgs bosons mainly decaying to SM particles when (g − 2) μ is not imposed. We note that the first-and second-generation sleptons are much lighter than the third-generation sleptons in the case with (g − 2) μ . The third-generation squarks are also heavier when (g − 2) μ is dropped, whereas the gluino and the firstand second-generation squarks are lighter in this case. In both cases, the third-generation squarks may lie within reach of future LHC runs, whereas the first-and second generation squarks would be accessible only if (g − 2) μ is dropped. The gluino would also be accessible in this case, and possibly also if (g − 2) μ is included.
We re-emphasize that the remarks in the previous paragraph apply to the best-fit points, and that the spectra might differ significantly, as the likelihood functions are quite flat for large masses. The 68 and 95% CL ranges are displayed in Fig. 23 as orange and yellow bands, respectively, with the best-fit values indicated by blue lines. We see that for most sparticles the 95 and even 68% CL ranges extend into the ranges accessible to future LHC runs. As was to be expected, the best prospects for measuring sparticles at a linear e + e − collider such as ILC [172,173] or CLIC [174] are offered by first-and second-generation sleptons and the lighter electroweak inosχ 0 1 ,χ 0 2 andχ ± 1 in the case with the (g − 2) μ constraint applied. Figure 24 displays the breakdowns of the global χ 2 functions in the cases with (left panel) and without (right panel) the (g − 2) μ constraint. 29 The different classes of observables are grouped together and colour-coded. We see that M W makes only a small contribution, and that the total contribution to the global χ 2 function of the precision electroweak observables are quite similar in the two cases. The total contribution of the flavour sector is slightly reduced when (g − 2) μ is dropped: χ 2 ∼ −1.2, largely because of a better fit to BR(B s → μ + μ − ), but this improvement is not very significant. The contributions of the Higgs, LEP, LHC and DM sectors are again very similar in the fits with and without (g − 2) μ .

Conclusions
In this paper we have used the MasterCode tool to analyze the constraints on the parameter space of the pMSSM11 model, in which the soft SUSY-breaking contributions to the masses of the first-and second-generation sleptons are allowed to vary independently from the third-generation slepton mass. We have taken into account the available constraints on strongly-and electroweakly-interacting sparticles from ∼ 36/fb of LHC data at 13 TeV [14][15][16] and the most recent limits from the LUX, PICO, XENON1T and PandaX-II experiments [3][4][5][6] searching directly for DM scattering. In addition, we have updated the constraint from the measurement of M W and some constraints from flavour observables, as described in Table 2. We have presented the results from two global fits, one including the (g − 2) μ constraint and without it. We have also made various comparisons with fits without the LHC 13-TeV data. Comparing with our earlier fit to the pMSSM10 [13], we note that the freedom for m˜ = mτ plays an important role in best fits. Furthermore, there is a big difference between M 1 and M 2 at the best-fit point without (g − 2) μ .
The most visible impact of the LHC 13-TeV constraints has been on the masses of the strongly-interacting sparticles: see the left panels of Figs. 8 and 9 and compare the solid and dashed curves. On the other hand, the impact of the LHC constraints on electroweak inos has been less marked: see Fig. 11. As was to be expected, the importance of the (g − 2) μ constraint is seen in the likelihood functions for charged slepton masses and electroweak inos: compare the blue and green curves in Figs. 10 and 11. The composition of the LSPχ 0 1 is also different in the cases with and without (g − 2) μ : as seen in Fig. 12 and Table 6, aB LSP is preferred when (g − 2) μ is included, whereas aH LSP is preferred when (g − 2) μ is dropped. Moreover, the inclusion of the (g − 2) μ constraint also has significant indirect implications for the squark masses, as also seen in Figs. 8 and 9. This anal-ysis reinforces the importance of clarifying the interpretation of the difference between the experimental measurement and the SM calculation of (g − 2) μ . We therefore welcome the advent of the Fermilab (g − 2) μ experiment [175] and continued efforts to refine the SM calculation.
We have also analyzed in this paper the importances of different mechanisms for bringing the relic LSP density into the range favoured by Planck 2015 and other data: see the shadings in Figs. 2, 4, 5, 6, 19, 20 and 21, and the profile χ 2 functions for the DM measures in Fig. 17. As we see there, important roles are played by chargino coannihilation, slepton coannihilation and rapid annihilation via direct-channel H/A boson exchange, though other mechanisms such as stau and squark coannihilation may be important in limited regions of parameter space. 30 In the case where the (g − 2) μ constraint is dropped, there is a preference for a region where mχ0 1 ∼ mχ± 1 ∼ mq ∼ mg where multiple coannihilation processes play a role, and the compressed spectrum reduces the sensitivity of the LHC sparticle searches.
In general, our analysis favours quite small deviations from the SM predictions for electroweak, flavour and Higgs observables: see Figs. 12 and 14, in particular. We have also analyzed the pMSSM11 predictions for the A and τ (B s → μ + μ − ) observables recently measured for the first time by the LHCb Collaboration [112]. As seen in Fig. 13, the pMSSM11 predictions for these observables are very similar to those in the SM, deviating by much less than the current experimental uncertainties. Accordingly, we do not include A and τ (B s → μ + μ − ) in our global fits.
We find that current LHC searches for long-lived particles do not impact our scan of the pMSSM11 parameter space. However, the pMSSM11 still offers significant prospects for the discovery of long-lived particles. When the (g − 2) μ constraint is imposed, we find that χ 2 4 for τ NLSP 10 −10 s. However, when the (g − 2) μ constraint is dropped, values of τ NLSP as long as 10 3 s (the limit we impose in order to maintain successful Big Bang nucleosynthesis) are allowed at the χ 2 4 level, As seen in Figs. 20 and 21, the pMSSM11 offers interesting prospects for the detection of supersymmetric DM. In both the spin-independent and -dependent cases, cross sections close to the present experimental upper limits are favoured at the 68% CL, whether or not (g − 2) μ is included in the set of constraints. Interestingly, in the case of σ SI p with (g − 2) μ included, there is a lower limit that is not far below 30 Compared to the pMSSM7 analysis in [52], we find that stop coannihilation is less prominent, and that rapid annihilation through the Z and the light Higgs boson is of very limited importance. In these respects the more general realization of the MSSM with four additional free parameters yields substantially different results. the neutrino 'floor', 31 whereas σ SI p may be much lower when (g − 2) μ is dropped, and low values of σ SD p are allowed in both cases.
We turn finally to the prospects for discovering sparticles in future runs of the LHC, or with a future linear e + e − collider. As seen in Fig. 21, whether or not (g − 2) μ is included in the global fit, the third-generation squarks may well be within reach of future LHC runs, and the first-and secondgeneration squarks and the gluino may also be accessible if the (g − 2) μ constraint is dropped. If it is included, on the other hand, there are also good prospects for discovering electroweakly-interacting sparticles at an e + e − collider, in particular theẽ,μ,χ 0 1 ,χ 0 2 andχ ± 1 . It is often said that the night is darkest just before dawn, and the same may be true for supersymmetry. the TU Munich for hospitality during the final stages of this work and has been partially supported by the DFG cluster of excellence is supported by XuntaGal. The work of K.A.O. is supported in part by DOE grant de-sc0011842 at the University of Minnesota. The work of G.W. is also supported in part by the European Commission through the "HiggsTools" Initial Training Network PITN-GA-2012-316704. During part of this work we used the middleware suite udocker [176] to deploy MasterCode on clusters, developed by the EC H2020 project INDIGO-Datacloud (RIA 653549). We are particularly grateful to Jorge Gomes for his kind support. We also thank DESY and especially the DESY IT department for making us available the computational resources of the BIRD/NAF2 cluster, which have been used intensively to carry out this work.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative 31 However, we repeat that the uncertainties in the calculation of σ SI p are large, and these remarks apply within the framework of a calculation of σ SI p using SSARD.
Commons license, and indicate if changes were made. Funded by SCOAP 3 .