UvA-DARE (Digital Academic Repository) A global fit of the MSSM with GAMBIT

We study the seven-dimensional Minimal Super-symmetric Standard Model (MSSM7) with the new GAM-BIT software framework, with all parameters deﬁned at the weak scale. Our analysis signiﬁcantly extends previous weak-scale, phenomenological MSSM ﬁts, by adding more and newer experimental analyses, improving the accuracy and detail of theoretical predictions, including dominant uncertainties from the Standard Model, the Galactic dark matterhaloandthequarkcontentofthenucleon,andemploy-ing novel and highly-efﬁcient statistical sampling methods to scan the parameter space. We ﬁnd regions of the MSSM7

The most general version of the MSSM has a very large number of parameters.Assuming only C P conservation, there are in total 63 free parameters, coming mostly from soft supersymmetry-breaking terms.There are two distinct ways to approach the exploration of the MSSM.The first is to take inspiration from the apparent unification of the gauge couplings at a high scale, defining a small set of unified parameters at that scale -a Grand Unified Theory (GUT) -and then run them down to the electroweak scale in order to compare with experiment.This is done for example in the four-parameter constrained MSSM (CMSSM) with common mass parameters for the scalar and fermion soft masses [160], and various generalisations of the CMSSM, such as the non-universal Higgs mass models 1 and 2 (NUHM1 and NUHM2), which split the Higgs soft masses from the other scalar masses [161][162][163][164][165]. We have carried out extensive global fits to these GUT-motivated models in a companion paper to this one [166].
The other approach is to remain agnostic about the highscale properties of the theory, and to define all the parameters at an energy near the electroweak scale.This is the socalled 'phenomenological MSSM' (pMSSM [167]).Given the impracticality of studying the complete parameter space, it is necessary to make some physically-motivated assumptions and simplifications in order to focus upon a relevant and tractable lower-dimensional subspace of the full model.
In this paper, we perform the first global fit of the weakscale phenomenological MSSM to make use of the GAMBIT global-fitting framework [168].Our work improves upon previous pMSSM studies on the following fronts: 1.The larger number of observables in our composite likelihood function, including: sparticle searches at LEP, Run I and Run II of the LHC, observables and constraints on the Higgs sector from LEP, the Tevatron and the LHC, direct and indirect dark matter searches from a multitude of experiments, and a large number of flavour and precision observables.2. Experimental likelihoods reconstructed from event rates, where applicable, including: Monte Carlo simulation of LHC observables, kinematical acceptance-based event rates for LEP sparticle searches, indirect DM search likelihoods computed at the level of single events, and direct DM search limits based on sophisticated simulation of the relevant experiments.3. The use of the GAMBIT hierarchical model database and statistical framework, for an advanced treatment of uncertainties from dominant SM nuisance parameters across different observables, and a consistent theoretical and statistical treatment of all likelihoods.4. The use of Diver [169], a new scanner based on differential evolution, which improves sampling performance compared to techniques used previously, and allows us to more accurately locate the high-likelihood regions.5.The results in this paper are based on a public, opensource software, which allows for the reproduction and extension of our results by the interested reader. 1e begin in Sect. 2 by introducing the models, parameters and nuisances over which we scan, followed by our adopted priors and sampling methodology.In Sect.3, we then give a brief summary of the observables and likelihoods that we employ.We present our main results in Sect. 4 and their implications for future searches for the MSSM in Sect.5, and conclude in Sect.6.
All GAMBIT input files, generated likelihood samples and best-fit benchmarks for this paper are publicly available online through Zenodo [170].
2 Models and scanning framework

Model definitions and parameters
GAMBIT makes no fundamental distinction between parameters of BSM theories and nuisance parameters, scanning over each on an equal footing.Here we sample simultaneously from four different models: a 7-parameter phenomenological MSSM, and three models describing constraints on different areas of known physics relevant for calculating observables in the MSSM.These nuisance models respectively describe the SM, the Galactic DM halo, and nuclear matrix elements for different light quark flavours (relevant for direct detection of DM).

MSSM7
The most general formulation of the C P-conserving MSSM is given by the GAMBIT model MSSM63atQ.Full details of the Lagrangian can be found in Sec.5.4.3 of Ref. [168].This model has 63 free, continuous MSSM parameters: 3 gaugino masses M 1 , M 2 and M 3 , 9 parameters each from the trilinear coupling matrices A u , A d and A e , 6 real numbers associated with each of the matrices of squared soft masses L and m 2 e , and three additional parameters describing the Higgs sector.We choose to work with the explicit mass terms m 2 H u and m 2 H d for the two Higgs doublets.By swapping the Higgs bilinear couplings b and μ for the ratio of vacuum expectation values for the up-type and downtype Higgs fields tan β ≡ v u /v d , and demanding that the model successfully effect Electroweak Symmetry Breaking, we can reduce the remaining continuous freedom to a single parameter (tan β).This leaves only a free sign for μ, which constitutes an additional (64th) discrete parameter.In this definition, tan β is specified at the scale m Z , and all other parameters are defined at some other generic scale Q, usually taken to be near to the weak scale.
This parameter set is currently too large to explore in a global fit, and in any case much of the phenomenology can be captured in smaller models that incorporate simplifying assumptions.In this first paper, we explore the MSSM7atQ, a 7-parameter subspace of the MSSM63atQ.Inspired by GUT theories, we set Table 1 MSSM7 parameters, ranges and priors adopted in the scans of this paper.For a parameter x of mass dimension n, the "hybrid" prior is flat where |x| < (100 GeV) n , and logarithmic elsewhere.The "split hybrid" prior for M 2 refers to the fact that we carried out every scan twice: once with a hybrid prior over 0 ≤ M 2 ≤ 10 TeV, and again with a hybrid prior over − 10 TeV ≤ M 2 ≤ 0. In addition to the priors listed here, we also carry out additional scans of fine-tuned regions associated with specific relic density mechanisms, where we restrict models to mass spectra that satisfy various conditions.See text for details Flat, hybrid Flat, hybrid at the scale Q.We assume that all entries in A u , A d and A e are zero except for L and m 2  e to be zero, so as to suppress flavour-changing neutral currents.By setting all remaining mass matrix entries to a universal squared sfermion mass m 2 f , we reduce the final list of free parameters to H d and tan β (plus the input scale Q and the sign of μ).The MSSM7 has been studied in significant work in the previous literature, e.g.[171][172][173][174][175][176].
We assume that R-parity is conserved, making the lightest supersymmetric particle (LSP) absolutely stable, and discard all parameter combinations where the LSP is not a neutralino.This choice is discussed in more detail in Sec.2.1.1 of the companion paper [166].
In Table 1, we give the parameter ranges over which we scan the MSSM7 in this paper.We choose to define all parameters other than tan β at Q = 1 TeV, and investigate positive μ (for a definition of μ please see the superpotential given in Sec.5.4.3 of Ref. [168].).We intend to return to the μ < 0 branch of this model in future work, where we compare with less constrained subspaces of the MSSM63atQ.

Nuisance parameters
We make use of three different nuisance models in our scans: the SM as defined in SLHA2 [168,177], a model of the Galactic DM halo that follows a truncated Maxwell- Boltzmann velocity distribution [168,178], and a model containing nuclear matrix elements required for calculating DMnucleon couplings [168,178].We vary a total of five parameters across these models: the strong coupling at scale m Z in the M S scheme, the top pole mass, the local density of DM, and the strange and up/down quark contents of the nucleon.The treatment of these models and parameters here is identical to the treatment in the companion paper [166], where we scan the nuisance parameters using flat priors, and impose constraints on them within the likelihood function.We refer the reader to Sects.2.1.2-2.1.4and Sect.3.1 of Ref. [166] for details.For ease of reference however, here we reproduce in Table 2 the parameter ranges and values used in our scans.

Scanning methodology
Our scanning methodology is similar to that detailed in Sect.2.2 of the companion paper [166]; here we give only a short summary, focussing in particular on any points of difference.
We use a number of different priors, two different samplers,2 and a range of different sampler configurations to scan the parameter space of the MSSM7.We then combine the results of all scans into a single set of likelihood samples, and use it to generate profile likelihoods for the MSSM7 parameters and observables.Using multiple scanners, priors and sampling settings allows for accurate determination of both the highest-likelihood point and profile likelihood contours.We do not consider Bayesian posteriors in the current paper, as our preliminary investigations indicate that Bayesian results in weak-scale MSSM models are dominated by the choice of priors.This suggests that a careful choice of prior (based on e.g.fine-tuning considerations) is needed for later interpretation, which is beyond the scope of the current paper.
We primarily rely on the differential evolution scanner Diver 1.0.0 [169], but also perform some supplementary scans with the nested sampler MultiNest 3.10 [179].For a performance comparison of Diver, MultiNest and other modern scanning algorithms, please see Sec. 11 of Ref. [169].Unlike in the companion paper, we repeat all Diver scans with both self-adaptive jDE [180] and its λjDE variant [169].The λjDE algorithm is more aggressive in finding the maximum likelihood, whereas jDE ensures a more thorough exploration of large regions of degenerate likelihood.Carrying out jDEonly scans (as opposed to retaining λjDE for all scans, as in Ref. [166]) is more beneficial in the MSSM7 than in the NUHM2 or its subspaces, as the additional freedom of the MSSM7 means that individual regions satisfying the bound set by the DM relic density, whilst more numerous, occupy a more fragmented and isolated volume of the parameter space than in 'smaller' models like the CMSSM.
We carry out scans with effectively logarithmic priors using both Diver and MultiNest, where all parameters except tan β follow so-called "hybrid" priors.The hybrid prior we use for parameter x of mass dimension n is flat in the region |x| < (100 GeV) n , and logarithmic for larger |x|.We supplement these scans with additional Diver runs using flat priors on the parameters most sparsely sampled at large values in log-prior scans ( A u 3 , A d 3 and M 2 H u ), and additional Multi-Nest scans using flat priors on all dimensionful parameters.To further improve the completeness of the sampling across the parameter space, we also subdivide every run into separate scans for M 2 > 0 and M 2 < 0.
As we show in Sect.4, chargino co-annihilation is the dominant mechanism for producing an acceptable DM abundance in large parts of the allowed parameter space.To ensure that we also explore the narrow parameter regions where sfermion co-annihilation or resonant annihilation is responsible for depleting the relic density to or below the observed value, we perform four 'directed' scans with different additional conditions on the particle spectrum: For all four scans, we also require that , to avoid parameter combinations where chargino co-annihilation dominates.We carry out these scans exclusively with Diver, using pure jDE.In these directed scans, we allow M 2 to vary freely over its full positive and negative range (Table 1).
The mass conditions effectively act as priors, allowing us to obtain a starting population of points roughly in the right area, before refining the fit to the best-fit likelihood region corresponding to each mechanism for depleting the relic den-sity.Note that the bounds that we use for this process are quite generous compared to the actual definitions of these regions that we use later in this paper, as they are only designed to direct the scanner to the correct vicinity of parameter space in which to look for good fits, not to act as cut for physical interpretation.To generate an initial population of parameter points satisfying the cut in question, we simply draw randomly from the overall parameter space, and assign zero likelihood to all parameter combinations that do not satisfy the cut.We then take the set of resulting points satisfying the cut, and use them with Diver to discover new points with higher likelihood values, continuing until the algorithm indicates convergence (suggesting that the best fit has been found, to within some tolerance).The looser mass cut on sleptons compared to squarks ensures that we are able to generate an initial population within the required cut in a reasonable time.
Taking into account all sampling configurations, priors and parameter space partitionings, this leads to a total of 2 sgn(M 2 ) × 2 priors × (1 MultiNest + 2 Diver configurations) + 4 directed scans = 16 scans.Each of the 16 scans took on the order of a few days to run on 2400 modern (Intel Core i7) cores.In total, these scans resulted in 1.8 × 10 8 valid samples, of which 1.76 × 10 7 (1.37 × 10 7 ) were within the 2D 95% (68%) CL region of the global best-fit point.After the scans we ran all samples through a postprocessing step, using the postprocessor scanner of ScannerBit [169], in order to correct for a bug in the flavour likelihood and apply Run II LHC searches which had just recently become available in ColliderBit [181].Following this postprocessing step 1.67×10 7 (2×10 5 ) of the original samples ended within the 95% (68%) CL region of the global best-fit point.
For the main Diver scans, we set the population size NP to 19 200, and the convergence threshold convthresh to 10 −5 ; for the directed scans, we instead set NP = 6000 and convthresh = 10 −4 .For MultiNest, we use nlive = 5000 live points and a stopping tolerance of tol = 0.1.These are relatively loose choices, but this results from the fact that we are only using MultiNest to bulk out the dominant likelihood mode of each scan rather than to locate the best fit point with high accuracy (the later is performed by Diver).

Observables and likelihoods
We compute a wide range of collider, flavour, DM and precision observables, and combine them with results from the latest experimental searches to produce an extensive set of likelihood components, which we then use to construct a composite likelihood function for driving the samplers in this paper.The composite likelihood includes -the DM relic density measurement c h 2 = 0.1188 ± 0.0010 from Planck [182] (implemented as an upper limit, so as to permit an additional non-neutralino DM component) -Fermi-LAT Pass 8 dwarf limits on DM annihilation [183], -IceCube 79-string limits on DM annihilation in the Sun [47,184], -direct DM searches by LUX [185][186][187], Panda-X [188], PICO [189,190], XENON100 [191], SuperCDMS [192] and SIMPLE [193], -the anomalous magnetic moment of the muon [194,195], -MSSM contributions to the mass of the W boson, -59 different flavour observables measured by LHCb, Babar and Belle, -13 different ATLAS and CMS analyses from Run I and Run II (as in the companion paper [166], we apply the Run II searches as a postprocessing step), -Higgs limits, mass measurements and signal strengths from LEP and the LHC, and -limits from LEP on sparticle production and decay in 18 different channels.-nuisance likelihoods associated the local density of DM [178], the quark contents of the nucleon [178], the top mass [196] and the strong coupling [196].
The theoretical treatments, experimental data and final likelihood functions for all these observables match those described in Sec. 3 of the companion paper [166], so we refer the reader to that paper, and the GAMBIT module papers [178,181,196,197], for details.The only exception is the DM relic density calculation, which we perform for the MSSM7 with micrOMEGAs 3.6.9.2 [198] (with settings fast = 1, Beps = 10 −5 ), rather than DarkSUSY 5.1.3[173], because the former is faster for some highly degenerate sfermion coannihilation models.

Best fits
In much of the parameter space of the MSSM7 (and indeed the MSSM more generally), the annihilation cross-section of heavy neutralino DM is so small that the thermal relic 3 FlexibleSUSY gets model-dependent information from SARAH [199,200] and uses some numerical routines from SOFTSUSY [201,202].
density greatly exceeds the value measured by Planck.Such models are robustly ruled out.The only way for a model to respect this upper limit is to exhibit one or more specific mechanisms for depleting the thermal abundance, typically associated with co-annihilation with another supersymmetric species, or resonant annihilation via a neutral boson 'funnel'.
Five such mechanisms play a role within the final 95% confidence level (CL) regions of our scans.In the following discussion, we classify samples and colour regions according which mechanism(s) they display: -chargino co-annihilation: where 'heavy' may be H 0 or A 0 , and 'light' may be h 0 or Z 0 , and a parameter combination qualifies as a member of a region if either condition is satisfied.Indeed, this is the strategy we adopt in general: if a model fulfils one of these conditions, we include it in the region, even if it ends up becoming a member of multiple regions, and even if some dominate over others.For clarity, we do not make any attempt to identify hybrid regions, or determine which of the mechanisms dominates (as to do so would require assumptions about relative temperature dependences and interferences between different partial annihilation rates).The union of these regions contain the full set of models allowed at 95% CL.
In Table 3, we show the details of the best-fit point in each of these five regions, breaking down the final log-likelihood into contributions from the different observables included in the fit.The overall best fit occurs in the chargino coannihilation region, where the lightest two neutralinos and the lightest chargino are all dominantly Higgsino, and thus highly degenerate in mass.All pairwise annihilations and co-annihilations between any of these three species can thus contribute significantly to the final relic density in this region.In Fig. 1 we give a visual representation of the mass spectrum of this point, where one can see clearly that we have some very light neutralinos and charginos in this model.The masses are around 260 GeV, making them potentially interesting targets for future LHC searches (Sect.5.1).
We also define a so-called 'ideal' reference likelihood in Table 3.This is the best likelihood that a model could realistically achieve were it to predict all observed quantities precisely, and predict no additional contribution beyond the expected background in searches that have produced only limits.We compute this for most likelihood components by assuming that the model prediction is either equal to the measured value or the background-only prediction.For some highly composite observables, where many different chan- nels enter and the SM or background-only prediction can in principle be improved upon by introducing a BSM contribution, we take the ideal case to be the best fit achievable in a more general, effective phenomenological framework.The two likelihoods that we apply this treatment to are those associated with LHC measurements of Higgs properties, and the angular observables of the B 0 → K * 0 μ + μ − decay observed by LHCb.In the former case, we take the ideal likelihood to be the best fit obtainable by independently varying the mass, width and branching fractions of a single Higgs in order to fit the LHC data contained in HiggsSignals.For the latter, we take the ideal likelihood from the best-fit point that we found in a flavour EFT global fit, discussed in Sec.6.2 of the FlavBit paper [197].
The log-likelihood difference between the global best fit and the ideal case is ln L BF = 36.345.Compared to the CMSSM ( ln L BF = 36.820;4 BSM parameters + 1 sign), NUHM1 ( ln L BF = 36.702;5 BSM parameters + 1 sign) and NUHM2 ( ln L BF = 36.362;6 BSM parameters + 1 sign) [166], we see a fairly mild improvement from moving to the MSSM7 (7 BSM parameters).It is possible to use ln L BF to estimate the overall goodness of fit, but this requires knowledge of the effective number of degrees of freedom (see Sec. 4.1 of Ref. [166] for details).Guessing this to be between 30 and 50 (on the basis of the number of active observables in the fit) would lead to a p value of between 2 × 10 −5 and 0.02.Comparing the specific case of e.g.37 degrees of freedom to the equivalent calculation for the NUHM2, NUHM1 and CMSSM with 38, 39 and 40 degrees of freedom, respectively, the p value for the MSSM7 would be 5.9 × 10 −4 , compared to 5.9 × 10 −4 (NUHM2), 7.1 × 10 −4 (NUHM1) and 9.4 × 10 −4 (CMSSM).This comparison is not entirely fair, given that we have allowed sgn(μ) to vary in the other three theories but not in the MSSM7.Nonetheless, it does suggest that the likelihood improve-ment in the MSSM7 is not sufficient to overcome the p value penalty arising from the larger number of free parameters compared to the three GUT-scale models.
As with the CMSSM, NUHM1 and NUHM2, the individual contributions to ln L BF = 36.345come almost entirely from known anomalies unexplainable in either the SM or MSSM7.In particular, these include the D and D * meson decay ratios R D and R D * (contained in the tree-level B and D decay likelihood; see Sec. 5.1 of Ref. [197]), the magnetic moment of the muon (a μ ; see Sec. 4.2.2 of Ref. [196]) and the angular observables of the electroweak penguin decay B 0 → K * 0 μ + μ − (see Sec. 5.2 of Ref. [197]).We refer the reader to Sects.4.1-4.3 of the companion paper [166] for further discussion.
The best fits possible by relic density mechanisms other than chargino co-annihilation are not drastically worse than the global best fit.The best models in the four other regions all lie within ln L < 1.2 of the global optimum.Compared to the best fit with chargino co-annihilation, the best stop co-annihilation model has the light stop (m t1 = 760 GeV) needed to co-annihilate with the neutralino, and therefore a light sfermion spectrum more generally, due to the universality of m f at the weak scale.This point thus ends up being penalised by both the LHC Higgs likelihood and B → X s γ , but advantaged by B 0 → K * 0 μ + μ − .In contrast, the best-fit sbottom co-annihilation point has a heavier spectrum, with all sfermions masses above 1 TeV, and hence only suffers on B 0 → K * 0 μ + μ − relative to the global best fit.Both the light and heavy funnel best fits are hybrids with chargino co-annihilation, showing light charginos and neutralinos.The best-fit A/H funnel point is only marginally worse than the global best, improving slightly on it in terms of B 0 → K * 0 μ + μ − but losing out due to slightly worse fits to a μ , B → X s γ and the LHC Higgs likelihood.The spectrum of the best fit in the Z / h funnel region is split, with heavy sfermions and gluinos, but light charginos and neutralinos.The latter lead to significant SUSY loop corrections to the W self energy.This model is also slightly worse than the best fit in terms of B 0 → K * 0 μ + μ − and a μ , but recovers somewhat by making a smaller contribution to B → X s γ .

Preferred regions
We begin by giving the 1D profile likelihoods for each of our input parameters in Fig. 2. For simplicity, we refer to m f ≡ (m 2 f ) 1/2 rather than the input parameter m 2 f .We also give 1D profile likelihoods for the derived parameters μ and M 1 .The GUT-inspired relation (Eq. 1) means that M 1 ≈ 0.48M 2 ≈ 0.18M 3 , while |μ| is determined from the EWSB conditions.With M 1 < M 2 it follows that M 1 and μ are the main mass parameters controlling the composition of the lightest neutralino.In our results we show M 1 and μ at the scale where the spectrum is calculated, Q = √ m t1 m t2 ≡ M SUSY .Due to the central role played by the μ parameter, it is more instructive to discuss the results connected to the Higgs sector in terms of μ and m A 0 than m 2 H u and m 2 H d .In Fig. 2 and throughout this paper, we show profile likelihood regions coloured according to the different coannihilation and funnel mechanisms contributing to keeping the neutralino relic density low enough to evade the Planck bound.These are: chargino co-annihilation (yellow), stop coannihilation (red), sbottom co-annihilation (blue), the A/H funnel (orange) and the h/Z funnel (purple); the definitions of these classifications can be found in the previous subsection.
Figure 3 shows the 2D joint profile likelihood for M 1 and μ (top) and M 2 and m f (bottom).The edges of the coloured regions here correspond to 95% CL relative to the best fit of the entire sample, not relative to the best fits of each coloured region.Here we see that the parameter space allowed at 95% CL encompasses three distinct regions, each expressing a different composition for the lightest neutralino and chargino: after which the wino component dominates.
Due to Eq. 1, a purely wino-dominated χ 0 1 is not possible in the MSSM7.
For Regions 1 and 2, the masses of the lightest chargino and the two lightest neutralinos are nearly degenerate, and all very close to μ.The neutralino relic density is therefore depleted by all pairwise annihilations and co-annihilations between the three species, which we collectively refer to simply as 'chargino co-annihilation'.In Region 1, where the lightest neutralino is essentially a pure Higgsino, the relic density constraint implies μ 1.2 TeV.The A/H -funnel also contributes across most of Regions 1 and 2, except in the case of very low μ or μ |M 1 |, where the dependence of m A 0 on |μ| makes it difficult to satisfy the funnel relation In Region 3, a small mass difference between the lightest neutralino and chargino is no longer automatic.The dominant relic density mechanisms in this parameter region are stop and sbottom co-annihilation, supported by annihilation through the A/H funnel.The tuning required in the former to get the lightest neutralino and lightest squark nearly degenerate in mass shows up as strongly-correlated bands in the M 2 -m f plane (lower panels of Fig. 3).Because the MSSM7 employs a common sfermion soft-mass parameter m 2 f at the input scale (Q = 1 TeV in our case), mass splittings among 123 Fig. 2 1D profile likelihood ratio for the input parameters A d3 , A u3 , M 2 , tan β, m H d , m Hu and m f , as well as the derived parameters M 1 and μ different sfermions are mostly generated by varying amounts of mixing.In comparison, the contribution from RGE running from Q = 1 TeV to Q = M SUSY , which splits m 2 f into individual soft masses, is generally subdominant.
In the tree-level stop mass matrix the off-diagonal element is vy t (A u 3 sin β − μ cos β), while it is vy b,τ (A d 3 cos β − μ sin β) in the sbottom and stau mass matrices, where y t,b,τ are the fermion Yukawa couplings and v ≈ 246 GeV.Because increased left-right mixing reduces the mass of the lighter of the two mass eigenstates, the large top Yukawa ensures that t1 is the lightest sfermion across most of the allowed parameter space (including for models that exhibit sbottom co-annihilation).With 3 ≤ tan β ≤ 70 the terms A u 3 sin β (stop) and μ sin β (sbottom and stau) dominate the sfermion mixing in large regions of parameter space.The dependence on large μ to obtain a sbottom mass significantly lower than the mass set by the common m f parameter explains why the sbottom co-annihilation region does not extend as far to small μ as the stop co-annihilation region in Fig. 3. Also, since y b ≈ 2.5y τ , the lightest stau remains heavier than the lightest sbottom in the regions of parameter space with large mixing for the down-type sfermions, which explains the absence of any region dominated by stau co-annihilation in our results.The requirement that all sfermions are heavier than the lightest neutralino excludes large regions of parameter space at m f 1.3|M 2 | ≈ 2.6|M 1 | in the bottom panels of Fig. 3.The steep slope of the exclusion boundary can broadly be understood as a consequence of the μ-dependent mixing in the sfermion sector.The region close to the boundary at m f ≈ 1.3|M 2 | is part of Regions 2 and 3 (μ |M 1 |) in the μ-M 1 plane, so that increasing M 2 in this region pushes up both M 1 and μ.To keep the lightest sfermion heavier than the neutralino, m f must therefore increase enough to compensate both the increase in neutralino mass from M 1 and the potential decrease in the light sfermion mass due to the μdependent left-right mixing.We come back to this interplay between the parameters of the neutralino and sfermion sector when discussing the μ-tan β plane in Fig. 5.
The region of small |M 1 | in the upper row of Fig. 3 (and therefore also small |M 2 | in the lower row) is strongly constrained by LHC limits.Direct LHC searches are also strongly constraining at low m f (lower panels).Gluino searches are particularly effective, as Eq. ( 1) implies that the gluino mass parameter is M 3 ≈ 7M 1 at a scale of 1 TeV.Given that simplified gluino mass limits reach up to 1.9 TeV [216], this disfavours bino masses in the MSSM7 of up to ∼ 300 GeV.Indeed, this is the main reason that we do not find the same preference for very light binos observed in Ref. [217], where each of the gaugino masses was allowed to vary independently.The common sfermion mass parameter means that, for light 3rd generation squarks, the 1st and 2nd generation squarks are not necessarily decoupled.Thus, LHC searches for 1st and 2nd generation squarks also constrain how far down towards low neutralino masses (small μ or |M 1 |; upper panels of Fig. 3) and low sfermion masses (lower panels of Fig. 3) the stop and sbottom co-annihilation regions extend.Measurements of the 125 GeV Higgs, limits from DM direct detection experiments, flavour physics and precision measurements of the W mass also contribute to disfavouring low gaugino masses in our fits.
At low sfermion masses, we also see a weak preference for positive M 2 , stemming from the (g−2) μ likelihood.Because we assume μ > 0 for our model, having M 2 (and thus M 1 ) positive ensures a positive SUSY contribution to (g − 2) μ .
In Fig. 4 we explore the impacts of the relic density constraint on the MSSM7 in more detail, investigating the profile likelihood of h 2 and m A as a function of the mass of the neutralino.The behaviour of Higgsino DM follows a relatively well-known pattern, seen also in the CMSSM and NUHM [166]: Higgsino DM co-annihilates steadily less efficiently as the neutralino mass increases, passing through the observed value of the relic density at m χ 0 1.2 TeV, the efficiency of Higgsino coannihilation makes for sub-dominant Higgsino DM, as seen in the diagonal chargino co-annihilation region in the upperright panel of Fig. 4.This can be tempered by fine-tuning the Higgsino-bino mixture, bringing up the relic density to the observed value, but such combinations are now very strongly constrained by direct detection, where mixed gaugino-Higgsino DM maximises both the spin-dependent and spin-independent neutralino-nucleon scattering crosssections.
At very low masses, the chargino co-annihilation region reaches down far enough that resonant annihilation via the SM Higgs further boosts the annihilation cross-section, leading to a region of hybrid chargino co-annihilation-h funnel models with neutralino masses as low as 61 GeV. 4 The best 4 For dedicated analyses of scenarios with a very light neutralino in MSSM parameterisations without a GUT relation on the gaugino masses, see for instance Refs.[218][219][220][221][222]. χ − 1 at LEP.Indeed, such masses would naively seem to be in contradiction with published limits, e.g.m χ ± 1 > 94 GeV [223,224].However, this particular limit assumes m χ ± 1 − m χ 0 1 > 3 GeV, and does not strictly apply to our best fit.The GAMBIT implementation of LEP limits in ColliderBit, detailed in Sec.2.2 of Ref. [181], takes into account the mass-dependent signal efficiency for the chargino and neutralino searches.These are quite important for cases where the spectrum has some degenerate masses, as in our best fit.In this case, the relevant search is the one for leptonic decays of the chargino at L3, with results shown in Fig. 2b of Ref. [225].Our treatment is a significant improvement on the hard lower limits that have often been used in the past.
Figure 4 shows that the heavy Higgs funnel can work for a wide range of neutralino masses in the MSSM7, from ∼ 200 GeV up to many TeV.The lower limit here comes from the lower limit on the mass of the CP-odd Higgs boson, seen in the bottom-left corner of the m A 0 -m χ 0 1 plane (Fig. 4).This arises due to penalties from the flavour physics likelihoods and the LHC Higgs likelihood.Because A 0 is close H + = m 2 A 0 + m 2 W at tree level), having a light A 0 causes tension with the B R(B → X s γ ) likelihood, which in isolation requires m H + 570 GeV at 95% CL for type-II two Higgs doublet models such as the MSSM [226].For large tan β, the likelihoods for tree-level leptonic and semi-leptonic B and D decays also penalise low A 0 masses.The tension with these likelihoods at low masses is to some extent compensated for by an improvement in the fit to the electroweak penguin decay B 0 → K * 0 μ + μ − , but for m 0 A 400 GeV, the combined restrictions imposed by flavour physics and measurements of the 125 GeV Higgs push the likelihood below the 95% CL, as evident in Fig. 4.
In this paper we have allowed neutralinos to be a subdominant component of DM.Were we to instead require that they constitute all of DM, our fits would be concentrated in the area around the horizontal line in the upper panels of Fig. 4.This would restrict the Higgsino-dominated DM models of the chargino co-annihilation region to m χ 0 1 1 TeV, moving the best-fit point to the A/H funnel and a mass of m χ 0 1 = 416 GeV.In terms of the neutralino mass itself, this would rule out m χ 0 1 < 250 GeV at 95% CL (1D).As we discuss later in this section, the absence of light charginos would also degrade the (already poor) fit to a μ .
In Fig. 5, we show the preferred regions and relic density mechanisms active in the μ-tan β and A d 3 -A u 3 planes.
The shape of the allowed region in the μ-tan β plane can be understood as follows.For the scenario in Region 1 of the upper panels of Fig. 3, μ M 1 and the lightest neutralino is dominantly Higgsino.This leads to the relic density bound μ 1.2 TeV.In Region 2, where the lightest neutralino is a mixture of bino and Higgsino, this upper bound on μ increases to ∼ 2.5 TeV.This limit is where we see the edge of the chargino co-annihilation and A/H -funnel regions at intermediate and large tan β in the upper panels of Fig. 5.
In Region 3, μ > |M 1 | and the lightest neutralino is dominantly bino, so there is no upper bound on μ from the relic density.In this case, the viable relic density mechanisms are stop/sbottom co-annihilation and the heavy Higgs funnel.Stop/sbottom co-annihilation can only work if the bino mass (M 1 ) is similar to the mass of the lightest squark.At large tan β, the left-right mixing in the sbottom mass matrix is proportional to μ, meaning that to keep the sbottom from becoming tachyonic, the diagonal entry (m 2 f ) must be increased as μ is increased.Pulling up m 2 f therefore pulls up the mass of the lightest squark, which in turn requires pulling |M 1 | up in order to stay in the stop/sbottom co-annihilation region.This is a delicate game, as |M 1 | needs to be kept below μ in order to remain in the bino LSP region (Region 3) at all.Whether or not this is possible depends on small corrections 123 from other parameters.At smaller values of tan β, the leftright mixing picks up an additional contribution proportional to A d 3 , and the adjustment can be pulled off with the help of some additional tuning in A d 3 .The net result is that |M 1 | remains less than μ, but not by more than a factor of a few.Because the heavy Higgs bosons receive their dominant mass contribution from |μ|, this sets their masses to be a factor of a few times that of the lightest neutralino, making stop/sbottom co-annihilation at higher μ in Region 3 appear mostly as a hybrid with the A/H funnel.
At large μ and large tan β, models in Region 3 are also impacted significantly by the Higgs likelihood.As discussed in Refs.[227,228], the bottom Yukawa coupling receives important SUSY corrections proportional to μ tan β, coming from gluino-sbottom and charged Higgsino-stop loops.For large μ and tan β, this increases the decay rate (h 0 → bb), which reduces the signal strengths for all other Higgs channels.The gluino-sbottom contribution is generally dominant, and for μ > 0 it is always positive.On the other hand, the Higgsino-stop contribution is proportional to A u 3 , so that for large and negative A u 3 it can compensate the gluino-sbottom correction.Thus, the good-fit region extending out towards large μ is dominantly associated with large, negative A u 3 .
Large |A u 3 | may cause the scalar potential of the MSSM to develop a minimum that breaks gauge invariance.We checked this in the same way as described in Sec.4.1 of the companion paper [166], finding even less impact in the MSSM7 than in the CMSSM or NUHM: whilst a small number of individual points are potentially affected by colour-or charge-breaking vacua, the overall preferred regions of the model remain unaffected.We naively carried out the same tests for |A d 3 | as well, swapping all up-type parameters for their down-type equivalents.We found that a few more models were affected than in the up-type tests, in particular those at large μ and small tan β discussed in the context of Fig. 5 above, where A d 3 helps to prevent the sbottoms becoming tachyonic.However, the impact was still quite isolated and had no impact on the overall inference.
In Fig. 6, we show the profile likelihood for the SUSY contribution a μ to the magnetic moment of the muon, compared with the experimental likelihood function for the observed discrepancy a μ,obs −a μ,SM = (28.7±8.0)×10−10 .Chargino co-annihilation models give the largest SUSY contributions, as they exhibit lighter charginos than other models.However, due to the relatively large values preferred for m f , which governs the masses of μ and νμ , it is essentially impossible to fit a μ simultaneously with all other observables even in the chargino co-annihilation region.
Compared to the MSSM10 results discussed in Ref. [217], we see broadly similar and consistent phenomenology, up to differences expected from the slightly different models being scanned.Both studies find the light Higgs funnel, chargino co-annihilation and stop/sbottom co-annihilation in essen-Fig.6 1D profile likelihood ratio for the SUSY contribution a μ to the anomalous magnetic moment of the muon.In green we show a Gaussian likelihood for the observed value a μ,obs −a μ,SM = (28.7±8.0)×10−10 , where we have combined the experimental and Standard Model (SM) theoretical uncertainties in quadrature tially the same areas.As already discussed, we find that the MSSM7 does not permit stau co-annihilation, and we see a preference for larger neutralino and sfermion masses than Ref. [217], a consequence of the unified gaugino and sfermion mass parameters in the MSSM7 and our inclusion of constraints from Run II of the LHC.We also see stop/sbottom co-annihilation extend to higher masses than in Ref. [217], reflecting either a lower likelihood for such models relative to the best fit in the MSSM10 than in the MSSM7, or improved sampling in the current paper.Unlike in the MSSM10, we find that it is not possible to consistently explain a μ in the MSSM7.

LHC
In Fig. 7 we show the 1D profile likelihoods for the masses of χ 0 1 , χ ± 1 , g, t1 , b1 and τ1 .The 2σ preferred region for the gluino mass extends upwards from ∼ 2 TeV, which is on the border of exclusion by current LHC searches for 0-lepton final states, to ∼ 20 TeV, well beyond the reach of the LHC.Similarly for m τ1 , where the small, weak production crosssection ensures that the predicted mass range is currently unobservable at the LHC.
More interesting are the χ 0 1 and χ ± 1 profile likelihoods, which are both peaked at low values.Given that these are naively within range of both LEP and the LHC Run I analyses, it is worth examining the properties of these low mass points in detail.Figure 8 shows our profile likelihood function in the χ ± 1 -χ 0 1 mass plane, zoomed into the low-mass region, , m g , m t1 , m b1 and m τ1 .We show separate distributions for each mechanism that allows the models to obey the relic density constraint Fig. 8 Left: profile likelihood in the χ± 1 − χ0 1 mass plane.Centre: subregions within the 95% CL profile likelihood region, coloured according to mechanisms by which the relic density constraint is satisfied.The regions shown correspond to neutralino co-annihilation with charginos, stops or sbottoms, and resonant annihilation through the light or heavy Higgs funnels.Superimposed in red is the latest CMS Run II simplified model limit for χ± 1 χ0 1 production and decay with decoupled slep-tons [229].This limit should be interpreted with caution (see main text for details).Right: the same information as the central plot, but zoomed into the low-mass region.Note that, although the CMS limit appears to have excluded part of the chargino co-annihilation region, this is a binning effect.One should instead refer to the plot of the χ ± 1 − χ0 1 mass difference in Fig. 7, which provides finer resolution on the mass difference in this region along with colour-coding indicating which mechanisms help to satisfy the relic density constraint.For the part of our 2σ region with m χ ± 1 275 GeV, an acceptable relic density is mostly generated via chargino co-annihilation, leading to very degenerate χ ± 1 and χ 0 1 masses.This explains the lack of exclusion by the LEP and LHC analyses included in our scan Fig. 9 1D profile likelihood ratios for the χ± 1 − χ0 1 mass difference (top) and the t1 − χ0 1 mass difference (bottom).Left: separate distributions for each mechanism allowing models to obey the relic density constraint.The regions correspond to neutralino co-annihilation with charginos, stops or sbottoms, and resonant annihilation through the light or heavy Higgs funnels.Right: as per the left, but zoomed in to small mass differences (which lose sensitivity for compressed spectra).Notably, our more careful treatment of the LEP limits than in previous studies has allowed models within the naive LEP reach to emerge unscathed.
One might wonder if other LHC analyses will soon (or have probed this low-mass The most EW gaugino limits are from CMS [231][232][233][234], using 36 fb −1 of 13 TeV data.A detailed study of the impact of these results would the addition of the relevant analyses to the liderBit module, and the calculation a complete likelihood similar to the equivalent Run I analyses already included in ColliderBit.In the present case, however, we can already obtain insight from a basic analysis of the simplified model limits presented by CMS Collaboration.CMS their results for each state in a range of simplified models of chargino and neutralino production, in which they set the branching fractions for specific decays to 100%, fixed the gaugino content, and set a 95% CL exclusion limit in the χ ± 1 -χ 0 1 mass plane.Figure 7 demonstrates that the sleptons are heavy across our entire preferred 2σ region, which is a natural consequence of having a unified sfermion mass in our parameterisation of the MSSM7.At least one stop mass must be high to induce large radiative corrections to the Higgs mass, which has the effect of dragging up the sfermion mass scale.In addition, and as mentioned previously, the τ1 will typically be heavier than the t1 and b1 due to the smaller Yukawa coupling.Thus, the relevant CMS simplified models are those featuring decoupled sleptons [229].We caution that these limits do not apply in general, and do not directly translate to limits on our model points without a detailed check that the neutralino and chargino mixing matrices and decay branching fractions match the CMS assumptions.One can, however, treat the CMS limits as the most optimistic possible exclusion in the χ ± 1 -χ 0 1 plane, to obtain a rough guide to the CMS sensitivity. Proceeding in this spirit, we see that the current CMS limits just barely touch our 2σ contour in regions where the spectrum is not compressed (Fig. 8).Indeed, the highest likelihood region looks to be out of reach in the near future.Note that if the GUT-inspired constraint on M 2 is relaxed, more solutions would fall within the CMS exclusion limit, so these searches will be important for global fits with more parameters.For compressed spectra, the details are less clear, as the ability of the CMS soft dilepton search to exclude the lightest models depends crucially on the precise χ ± 1 -χ 0 1 mass splitting.This is shown in the top of Fig. 9, where it is apparent that the chargino co-anihilation points appear as a peak in the likelihood at χ ± 1 -χ 0 1 mass differences of less than 10 GeV.This is too small to be probed by the recent CMS results.The chargino co-annihilation region remains free from LHC exclusion, assuming prompt decays of the chargino.We note, however, that for very small mass differences (approaching the pion mass), long-lived particle searches might provide additional constraints.We defer a detailed analysis of these to future work.
We now look at whether it is possible to probe the squark sector of the MSSM7 at the LHC in the near future.The lightest squarks are the t1 and b1 .Figure 7 shows that the peak of the sbottom profile likelihood lies out of reach of the LHC in the near future, and that masses below ∼ 800 GeV are disfavoured at the 2σ level.Figure 10 shows the b1 − χ 0 1 mass plane, revealing that the lower sbottom masses are associated with a small b1 − χ 0 1 mass difference.This arises from the fact that stop and/or sbottom co-annihilation often account for the generation of an acceptable relic density in this low-mass region.However, there are also low-mass regions in which resonant A/H annihilation or chargino co-annihilation contribute to DM annihilation, giving a wider range of mass differences.As above, comparison with recent CMS simplified model limits provides some insights into the ability of the LHC to probe these models in the near future.A variety of CMS searches for sbottom production have been interpreted in the context of a simplified model of sbottom pair production and decay to a bottom quark and the lightest neutralino [235][236][237].We again treat these limits as a rough guide to the most favourable possible LHC exclusion potential, and compare our results to the CMS summary plot given in Reference [238].The current analyses have potentially probed a small region of Fig. 10 (with χ 0 1 masses below 600 GeV and b1 masses below ≈ 1 TeV).However, almost our entire 2σ preferred region remains unconstrained.Directly ruling out sbottom co-annihilation as a viable contributor to an acceptable relic density would require probing compressed spectra in sbottom decays up to a mass of ∼ 4 TeV, an impossible task at the LHC.Nonetheless, the fact that current limits are nearing the tip of the stop co-annihilation strip means that discovery prospects even in the next run of the LHC are quite promising (although more so for models that exhibit only stop co-annihilation than those that display both stop and sbottom co-annihilation).
The stop mass has a marginally higher likelihood at lower masses (Fig. 7). Figure 11 shows the profile likelihood ratio in the t1 -χ 0 1 mass plane, along with colour-coded regions illustrating the relevant relic density mechanisms.As for the sbottom mass, points with a t1 mass below 1 TeV show a strong mass correlation with the lightest neutralino, as they lie in the stop co-annihilation region.Comparison with the most recent CMS Run II simplified model results [236,237,[242][243][244] Fig. 11 Left: the profile likelihood ratio in the t1 − χ0 1 mass plane.Centre: colour-coding shows mechanism(s) that allow models within the 95% CL region to avoid exceeding the observed relic density of DM.The regions shown correspond to neutralino co-annihilation with charginos, stops or sbottoms, and resonant annihilation through the light or heavy Higgs funnels.Superimposed in red is the latest CMS Run II simplified model limit for stop pair production [230].Right: the same information as the central plot, zoomed into the low-mass region Fig. 12 Spin-independent neutralino-proton scattering cross-sections in the MSSM7, rescaled by the fraction f of the observed relic density predicted by each model.Left: profile likelihood, showing 68 and 95% CL contours.Right: colour-coding shows mechanism(s) that allow models within the 95% CL region of the profile likelihood to avoid exceeding the observed relic density of DM, corresponding to neutralino co-annihilation with charginos, stops or sbottoms, and resonant annihilation through the light or heavy Higgs funnels.Overplotted are 90% CL constraints from LUX, [187], and projections for the reach of XENON1T after two years of exposure, XENONnT/LZ, assuming 1-3 years of data and an exposure of 20 tonne-years [239], and DARWIN, assuming 3-4 years of data and 200 tonne-years of exposure [240].The dashed grey line indicates the "neutrino floor" where background events from coherent neutrino scattering start to limit the experimental sensitivity [241].The exact placement of this limit is subject to several caveats; see [241] for further details reveals that the lowest-mass points in the stop co-annihilation region remain unprobed, as do the chargino co-annihilation and A/H -funnel points.The t1 − χ 0 1 mass difference is shown in the bottom panels of Fig. 9.Although this is of course small for the stop and sbottom co-annihilation region points, it is not, contrary to the chargino case, sharply peaked at sufficiently low values that decay products can be assumed to be hard to reconstruct at the LHC.This offers hope that LHC searches for compressed spectra (sensitive to smallish mass differences) can eventually tackle these models.

Direct detection of dark matter
In this section, we examine the preferred spin-independent (SI; Fig. 12) and spin-dependent (SD; Fig. 13) neutralinoproton scattering cross-sections in the MSSM7.Here we rescale the scattering cross-sections by the ratio f of the predicted to the observed relic density, so as to ease comparison with various experimental limits and projections.Figure 12 shows that SI limits from direct detection are already highly constraining, with many models with high likelihoods lying  [47,184], assuming that dark matter annihilates exclusively via the bb or τ + τ − channel, PICO-60 [245], and projections for the reach of PICO-250 [246] just below the current sensitivity of LUX [187], and very soon to be probed by XENON1T [239] and its successors.Eventually, DARWIN [240] looks set to probe the entirety of the light Higgs funnel and the chargino co-annihilation region, as well as large parts of the heavy funnel and stop/sbottom co-annihilation regions.
In the SD sector, IceCube already constrains mixed gaugino-Higgsino models in the MSSM, as noted in Refs.[47,175,247].PICO [245] is not competitive for MSSM models, but its future upgrades appear set to make significant inroads into both Higgs funnels and the chargino co-annihilation region.However, it remains to be seen if XENON1T will probe such models on a shorter timescale.Future neutrino telescopes such as KM3NeT [248] and proposed upgrades to IceCube [249,250] may also offer significantly improved sensitivity to models in the MSSM7, but to date the expected sensitivity to DM masses above 100 GeV is not known.Whilst not particularly constraining in terms of SD proton scattering, LUX [251] already provides constraints on the SD neutralino-neutron cross-section, which are just beginning to touch on the allowed parameter space of the MSSM7 (not shown, but included in our scans via DDCalc [178]).
Although models exist down to SI and SD cross-sections of 10 −55 cm 2 in the stop/sbottom co-annihilation and A/H funnel regions of our fits, the large cancellations required to produce such cross-sections may be spoilt by loop corrections [254,255].This raises hope that future direct detection experiments will discover neutralino DM in the MSSM7 or a similar model.However, specific investigations in the MSSM7 suggest that this is not necessarily expected for all parameter combinations, so some parts of the parameter space should still be expected to lie well below any future sensitivity, even after applying higher-order corrections [172].

Indirect detection of dark matter
Let us finally address discovery prospects of the MSSM7 with indirect DM searches.To this end, we show in Fig. 14 the velocity-weighted annihilation cross section, σ v, in the limit of vanishing relative velocity of the annihilating DM particles, as a function of the lightest neutralino mass.We rescale this quantity by the square 5 of the fraction f of the calculated neutralino relic density to the observed DM abundance, thereby accounting for the possibility of the lightest neutralino making up only a fraction of DM.In the left panel, we show the profile likelihood, while in the right panel we indicate the mechanism(s) responsible for increasing the (co-)annihilation rate in the early Universe, and hence decreasing the present neutralino relic density to or below the observed DM abundance.For comparison, we also indicate the same current and projected future limits from selected indirect detection experiments as in Fig. 21 of the companion paper [166], namely present Fermi-LAT [183] limits for bb and τ + τ − final states from observations of 15 dwarf galaxies, projected Fermi-LAT limits for bb, and the projected sensitivity of the Chrerenkov Telescope Array (CTA) for bb, assuming 500 h of Galactic halo observations [253].
Across almost the entire neutralino mass range, we find models within the 95% CL region of the profile likelihood that exhibit present-day annihilation rates above the canon-Fig.14 Zero-velocity neutralino self-annihilation cross-sections in the MSSM7, rescaled by the square of the fraction f of the observed relic density predicted by each model.Left: profile likelihood, showing 68 and 95% CL contours.Right: mechanism(s) that allow models within the 95% CL region of the profile likelihood to avoid exceeding the observed relic density of DM.Overplotted are 95% CL constraints from the search for dark matter annihilation in 15 dwarf spheroidal galaxies by the Fermi-LAT Collaboration [183].These limits are based on 6 years of Pass 8 data, and are given for two different assumed annihilation final states ( bb and τ + τ − ).We also show the projected improvement in the bb channel after 15 years, if the number of known dwarfs were to quadruple in that time [252].The final curve is the best-case projected sensitivity of the Cherenkov Telescope Array to annihilation in the Galactic halo, computed assuming bb final states, neglecting systematic errors, and assuming 500 h of observation [253] ical thermal value of 3 × 10 −26 cm 3 s −1 .Those models are a subset of the A/H funnel region, where the pseudoscalar Higgs is almost exactly twice as heavy as the lightest neutralino, m A 2 χ 0 1 .This leads to resonant enhancement of the annihilation rate as v → 0, as is the case today -but not in the early Universe, where thermal effects mean that v = 0 in general.For some models in this part of the parameter space, current Fermi limits are already quite constraining.Projected Fermi limits, assuming 15 years of data on 60 dwarf galaxies (vs.6 years and 15 dwarfs for the current limits), will start to cut into the (current) 68% CL region.For neutralino masses above around 300 GeV, CTA will be even more sensitive than this.Large parts of the MSSM7 parameter space, however, will be impossible to probe with any planned indirect detection experiment; this includes, unfortunately, both the global best fit point of the MSSM7 and the best-fit points of all of the individual parameter regions corresponding to different mechanisms of lowering the relic density.
We emphasise that even though the CTA limits shown here are rather optimistic, in that they neglect the effect of systematic uncertainties [256], the above discussion somewhat underestimates the prospects of indirect DM searches.One reason is that we have neglected in our discussion other detection channels than gamma-rays, in particular charged cosmic rays.As discussed in some more detail in Sec.4.4.3 of the companion paper [166], radiative corrections in particular, e.g.[257][258][259], as well as Sommerfeld enhancement [73,[260][261][262], are further effects that we have not taken into account here.For parts of the parameter space this leads to increased annihilation rates and/or distinct spectral features, which are much easier to constrain or identify with experi-ments than the featureless gamma-ray spectra from the final states that we have considered here.A full discussion of these effects, and their impact on indirect DM searches within the MSSM7, is beyond the scope of this study, although we plan to return to this in future work.

Conclusions
We have carried out an extensive global fit of the 7-parameter, weak-scale phenomenological MSSM, using the newlyreleased GAMBIT global fitting framework.Our fit takes into account updated experimental data, improved theoretical calculations and more advanced statistical sampling methods than previous studies of similar models.We have also considered leading uncertainties from the Standard Model, the dark matter halo of the Milky Way, and the quark content of the nucleon, fully scanning over the relevant parameters and profiling them out in the final fit.Finally, we have explored the full range of experimentally-allowed parameters, by allowing neutralinos to constitute any fraction of the observed cosmological dark matter.
The MSSM7 shows quite a rich selection of phenomenology across its parameter space, ranging from Higgsinodominated dark matter annihilating through co-annihilations with other Higgsinos in the early Universe, to resonant annihilation via the light and heavy Higgs funnels, to coannihilation of neutralinos with both light stops and sbottoms.We find a preference for light, Higgsino-dominated neutralinos, with m χ 0 1 750 GeV at 68% CL and m χ 0 1 2.5 TeV at 95% CL.We have shown that stop/sbottom co-123 annihilation models lie just out of reach of current LHC searches, with stops and sbottoms as light at 500 GeV.This makes the prospects for probing at least some such models at the LHC in the near future quite promising.Both direct and indirect searches for dark matter place significant constraints on the allowed parameter ranges in the MSSM7, and the next generation of these experiments will probe large parts of the highest-likelihood areas of its parameter space.
The current study is essentially a starting point for detailed, modular scans of supersymmetric models defined at the weak scale with GAMBIT.GAMBIT's hierarchical model database already contains many generalisations of the MSSM7, which would themselves make very interesting targets for global analyses similar to this one.
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 Commons license, and indicate if changes were made.Funded by SCOAP 3 .

Fig. 1
Fig. 1 Sparticle mass spectrum of the best-fit point

Fig. 3
Fig. 3 Left: joint profile likelihoods in the μ-M 1 (top) and M 2 -m f planes (bottom).Stars indicate the point of highest likelihood in each plain, and white contours correspond to the 1σ and 2σ CL regions with respect to the best-fit point.Right: coloured regions indicating in which parts of the 2σ best-fit region different co-annihilation and funnel mechanisms contribute to keeping the relic density low.The best-fit point in each region is indicated by a star with the corresponding colour

Fig. 4 1 ∼ 1 . 2 1 .
Fig. 4 Left: joint profile likelihoods in the mass of the lightest neutralino and its relic density h 2 (top), and in the masses of the lightest neutralino and the CP-odd Higgs A 0 (bottom).Right: coloured regions indicating in which parts of the 2σ best-fit region different co-annihilation and funnel mechanisms contribute to the relic density.The best-fit point in each region is indicated by a star with the corresponding colour

1 = 1 = 71 . 6 2 =
fit in this region (Table3) has m χ 0 69.2 GeV, m χ ± GeV and m χ 0 73.7 GeV, while the other sparticles are fairly heavy.This leads to considerable cross sections for direct pair production of χ

Fig. 5
Fig. 5 Left: joint profile likelihoods in μ-tan β (top) and A u3 -A d3 planes.Right: coloured regions indicating in which parts of the 2σ best-fit region different co-annihilation and funnel mechanisms contribute to the relic density.The best-fit point in each region is indicated by a star with the corresponding colour.In the bottom right plot the yellow chargino co-annihilation region covers the entire plane and the orange A/H funnel region spans the entire plane below A u3 ∼ 5 TeV

Fig. 10
Fig. 10 Left: The profile likelihood ratio in the b1 − χ0 1 mass plane.Centre: Colour-coding shows mechanism(s) that allow models within the 95% CL region to avoid exceeding the observed relic density of DM.The regions shown correspond to neutralino co-annihilation with

Fig. 13
Fig. 13 Spin-dependent neutralino-proton scattering cross-sections in the MSSM7, rescaled by the fraction f of the observed relic density predicted by each model.Left: profile likelihood, showing 68 and 95% CL contours.Right: mechanism(s) that allow models within the 95% CL region of the profile likelihood to avoid exceeding the observed

Table 2
Standard Model, dark matter halo and nuclear nuisance parameters and ranges.

Table 3
[170]fit points in the A/H -funnel, b co-annihilation, t coannihilation, χ± 1 co-annihilation and Z / h funnel regions.For each point, we show the individual likelihood contributions, parameter values (including nuisance parameters) and derived quantities crucial for interpreting the mass spectrum.Other SM and astrophysical parameters are set to the fixed values given in Table2.SLHA1 and SLHA2 files for the best-fit points can be found in the online data associated with this paper[170]