# Global fits of GUT-scale SUSY models with GAMBIT

Peter Athron, Csaba Balázs, Torsten Bringmann, Andy Buckley, Marcin Chrząszcz, Jan Conrad, Jonathan M. Cornell, Lars A. Dal, Joakim Edsjö, Ben Farmer, et al.

► To cite this version:

Peter Athron, Csaba Balázs, Torsten Bringmann, Andy Buckley, Marcin Chrząszcz, et al. Global fits of GUT-scale SUSY models with GAMBIT. European Physical Journal C: Particles and Fields, 2017, 77 (12), pp.824. 10.1140/epjc/s10052-017-5167-0. hal-01703823

HAL Id: hal-01703823

https://hal.science/hal-01703823v1

Submitted on 22 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Global fits of GUT-scale SUSY models with GAMBIT

The GAMBIT Collaboration: Peter Athron<sup>1,2,a</sup>, Csaba Balázs<sup>1,2</sup>, Torsten Bringmann<sup>3</sup>, Andy Buckley<sup>4</sup>, Marcin Chrzaszcz<sup>5,6</sup>, Jan Conrad<sup>7,8</sup>, Jonathan M. Cornell<sup>9</sup>, Lars A. Dal<sup>3</sup>, Joakim Edsjö<sup>7,8</sup>, Ben Farmer<sup>7,8,b</sup>, Paul Jackson<sup>2,10</sup>, Abram Krislock<sup>3</sup>, Anders Kvellestad<sup>11,c</sup>, Farvah Mahmoudi<sup>12,13,f</sup>, Gregory D. Martinez<sup>14</sup>, Antje Putze<sup>15</sup>, Are Raklev<sup>3</sup>, Christopher Rogan<sup>16</sup>, Roberto Ruiz de Austri<sup>17</sup>, Aldo Saavedra<sup>2,18</sup>, Christopher Savage<sup>11</sup>, Pat Scott<sup>19,d</sup>, Nicola Serra<sup>5</sup>, Christoph Weniger<sup>20</sup>, Martin White<sup>2,10,e</sup>

Received: 22 May 2017 / Accepted: 24 August 2017 / Published online: 4 December 2017
© The Author(s) 2017. This article is an open access publication

**Abstract** We present the most comprehensive global fits to date of three supersymmetric models motivated by grand unification: the Constrained Minimal Supersymmetric Standard Model (CMSSM), and its Non-Universal Higgs Mass generalisations NUHM1 and NUHM2. We include likelihoods from a number of direct and indirect dark matter searches, a large collection of electroweak precision and flavour observables, direct searches for supersymmetry at LEP and Runs I and II of the LHC, and constraints from Higgs observables. Our analysis improves on existing results not only in terms of the number of included observables, but also in the

level of detail with which we treat them, our sampling techniques for scanning the parameter space, and our treatment of nuisance parameters. We show that stau co-annihilation is now ruled out in the CMSSM at more than 95% confidence. Stop co-annihilation turns out to be one of the most promising mechanisms for achieving an appropriate relic density of dark matter in all three models, whilst avoiding all other constraints. We find high-likelihood regions of parameter space featuring light stops and charginos, making them potentially detectable in the near future at the LHC. We also show that tonne-scale direct detection will play a largely complementary role, probing large parts of the remaining viable parameter space, including essentially all models with multi-TeV neutralinos.

<footnote n="1">¹ School of Physics and Astronomy, Monash University, Melbourne, VIC 3800, Australia</footnote>

<footnote n="2">² Australian Research Council Centre of Excellence for Particle Physics at the Tera-scale, Australia, http://www.coepp.org.au/</footnote>

<footnote n="3">³ Department of Physics, University of Oslo, 0316 Oslo, Norway</footnote>

<footnote n="4">⁴ SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, UK</footnote>

<footnote n="5">⁵ Physik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zurich, Switzerland</footnote>

<footnote n="6">⁶ Polish Academy of Sciences, H. Niewodniczański Institute of Nuclear Physics, 31-342 Kraków, Poland</footnote>

<footnote n="7">⁷ Oskar Klein Centre for Cosmoparticle Physics, AlbaNova University Centre, 10691 Stockholm, Sweden</footnote>

<footnote n="8">⁸ Department of Physics, Stockholm University, 10691 Stockholm, Sweden</footnote>

<footnote n="9">⁹ Department of Physics, McGill University, 3600 rue University, Montréal, QC H3A 2T8, Canada</footnote>

<footnote n="10">¹⁰ Department of Physics, University of Adelaide, Adelaide, SA 5005, Australia</footnote>

<footnote n="11">¹¹ NORDITA, Roslagstullsbacken 23, 10691 Stockholm, Sweden</footnote>

<footnote n="12">¹² ENS de Lyon, Centre de Recherche Astrophysique de Lyon UMR5574, Univ Lyon, Univ Lyon 1, CNRS, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="13">¹³ Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland</footnote>

<footnote n="14">¹⁴ Physics and Astronomy Department, University of California, Los Angeles, CA 90095, USA</footnote>

<footnote n="15">¹⁵ LAPTh, Université de Savoie, CNRS, 9 chemin de Bellevue, B.P.110, 74941 Annecy-le-Vieux, France</footnote>

<footnote n="16">¹⁶ Department of Physics, Harvard University, Cambridge, MA 02138, USA</footnote>

<footnote n="17">¹⁷ Instituto de Física Corpuscular, IFIC-UV/CSIC, Valencia, Spain</footnote>

<footnote n="18">¹⁸ Faculty of Engineering and Information Technologies, Centre for Translational Data Science, School of Physics, The University of Sydney, Camperdown, NSW 2006, Australia</footnote>

<footnote n="19">¹⁹ Department of Physics, Blackett Laboratory, Imperial College London, Prince Consort Road, London SW7 2AZ, UK</footnote>

<footnote n="20">²⁰ GRAPPA, Institute of Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands</footnote>

<footnote>a e-mail: peter.athron@coepp.org.au</footnote>

<footnote>b e-mail: benjamin.farmer@fysik.su.se</footnote>

<footnote>c e-mail: anders.kvellestad@nordita.org</footnote>

<footnote>d e-mail: p.scott@imperial.ac.uk</footnote>

<footnote>e e-mail: martin.white@adelaide.edu.au</footnote>

<footnote>f Also Institut Universitaire de France, 103 boulevard Saint-Michel, 75005 Paris, France.</footnote>
## Contents

1 Introduction .................................................. 2

2 Models and scanning framework .......................... 3

2.1 Model definitions and parameters .......................... 3

2.1.1 SUSY models ........................................ 3

2.1.2 Standard model .................................... 4

2.1.3 Dark matter halo model ............................ 4

2.1.4 Nuclear model ...................................... 5

2.2 Scanning methodology .................................... 5

3 Observables and likelihoods ................................ 7

3.1 Nuisance likelihoods .................................... 7

3.1.1 Standard model .................................... 7

3.1.2 Local halo model .................................... 7

3.1.3 Nuclear matrix elements ............................ 7

3.2 Spectrum calculation .................................... 7

3.3 Relic density of dark matter ............................ 7

3.4 Gamma rays from dark matter annihilation ............ 8

3.5 High-energy neutrinos from dark matter annihilation in the Sun .................................................. 8

3.6 Direct detection of dark matter .......................... 9

3.7 Electroweak precision observables ........................ 9

3.8 Flavour physics likelihoods ............................ 9

3.9 Searches for superpartners at LEP ........................ 10

3.10 Searches for supersymmetry at the LHC .............. 10

3.11 Higgs physics ............................................ 12

4 Results ........................................................ 12

4.1 CMSSM ........................................................ 12

4.2 NUHM1 ........................................................ 19

4.3 NUHM2 ........................................................ 22

4.4 Discovery prospects ........................................ 24

4.4.1 LHC ........................................................ 24

4.4.2 Direct detection ...................................... 32

4.4.3 Indirect detection .................................... 33

5 Conclusions .................................................... 35

References .......................................................... 36

## 1 Introduction

Although the Standard Model (SM) of particle physics has long provided a spectacularly successful description of physics at and below the electroweak scale, it remains incomplete. Explaining dark matter (DM), the asymmetry between matter and antimatter, the hierarchy between the Planck and electroweak scales, the origin of the fundamental forces and charges, or anomalies in low-energy precision and flavour measurements, requires extending the SM by adding one or more new particles.

The Minimal Supersymmetric extension of the SM (MSSM) offers solutions to many of these shortcomings, with substantial implications for dark matter [1–40], the cosmic matter-

antimatter asymmetry [41–43], Higgs physics [44–60], the unification of gauge forces [61–76], the stability of the electroweak vacuum [77–80], cosmological inflation [81–90], precision measurements [91–99] and flavor physics [100–102].

Even though the MSSM framework is predictive, its Lagrangian terms responsible for softly breaking SUSY contain over a hundred new parameters. This impairs the practical predictivity of the model. Mediation of supersymmetry breaking by Planck-scale physics is a popular and viable motivation for reducing the free parameters to a small number at the Grand Unified Theory (GUT) scale [103–108]. Here we analyse three scenarios motivated by gravity mediation: the Constrained MSSM (CMSSM) [109] and two of its Non-Universal Higgs Mass (NUHM1, NUHM2) extensions [110–114].

Opinions differ regarding the phenomenological feasibility of these models, especially the CMSSM. Global fits of the CMSSM after Run I of the Large Hadron Collider (LHC) indicate that its experimentally-viable parameter space has been pushed to regions with superpartners heavier than 1 TeV. The relic density of DM in these scenarios is set by neutralino-stau co-annihilation, resonant annihilation through a heavy Higgs, or a large Higgsino component [115–143]. Relaxing the assumption of scalar soft-mass universality at the GUT scale provides much more flexible phenomenology [120, 134, 144–154]. Global fits of the NUHM1 show that the tension that exists in the CMSSM between the measured Higgs mass and precision/flavour observables is reduced due to the decoupling of the Higgs sector from the squark and slepton sectors, and a new region of chargino co-annihilation opens up for dark matter [117, 119, 120, 126, 134, 145, 148]. Extending the parameter space to the NUHM2 [117, 145, 148] relaxes the constraints on the scalar masses even further.

In this work we use the GAMBIT framework [155–160] to scan and assess the viability of the parameter spaces of each of these three GUT-scale scenarios in detail. We also carry out a detailed comparison of our results with previous ones, to understand the impact of the improved theoretical calculations and updated experimental data that we include, and as a verification of our new computational framework. We have also carried out similar analyses of scalar singlet DM [161] and ‘phenomenological’ (weak-scale) SUSY models [162] with the GAMBIT framework.

There are several important features of our study that make it the most definitive exploration of the CMSSM, NUHM1 or NUHM2 to date:

1. We apply the DM relic density constraint as an upper bound only. This requires that the cosmological density of the lightest neutralino does not exceed the observed
density of DM. This is a conservative option from the point of view of excluding a light mass spectrum, as it introduces more possibilities for light Higgsino and light Higgsino-bino DM than in studies where a lower bound is also applied.

2. We include a significantly higher number of observables in our combined likelihood than has been done before. These include rates in multiple direct and indirect searches for DM, a wide range of LHC sparticle searches and Higgs observables, and a up-to-date set of flavour physics observables and electroweak precision measurements.

3. In addition to improving the quantity of data included in the fit, we have also improved the quality of the typical simulation treatments, including direct Monte Carlo simulation of LHC observables during the global fit, event-level indirect search likelihoods, and direct DM search limits based on rigorous simulation of the relevant experiments.

4. Using **GAMBIT** allows us to pursue a thorough theoretical and statistical approach, where theoretical assumptions are consistently treated across different observables and experimental searches. This includes the accurate treatment, via nuisance parameters, of uncertainties associated with the local DM distribution, nuclear matrix elements relevant for direct detection, and SM parameters.

5. **GAMBIT** includes an interface to *Diver* [160], a new scanner based on differential evolution, which provides significantly improved sampling performance compared to conventional techniques. This allows us to more accurately locate and more comprehensively map small regions of high likelihood.

6. The public, open-source nature of **GAMBIT**<fnref n="1" /> makes our study transparent, reproducible and extendible by the reader.

In Sect. 2, we introduce the CMSSM, NUHM1 and NUHM2, along with their parameters, the ranges and priors over which we vary those parameters, and the algorithms and settings that we use for sampling them. Section 3 contains a summary of the experimental data, observables and likelihood calculations that go into each fit. We then present our results in Sect. 4, before looking at the implications of our scans for future searches for the models in question, and concluding in Sect. 5.

All input files, samples and best-fit benchmarks produced for this paper are publicly accessible from Zenodo [163].

# 2 Models and scanning framework

## 2.1 Model definitions and parameters

From a statistical standpoint, there is no fundamental difference between models that describe SM physics (and astrophysics) and physics beyond the SM (BSM). **GAMBIT** therefore treats BSM models on exactly the same footing as models that describe nuisance parameters, which are designed to quantify uncertainties on better-constrained quantities. The only difference is that nuisance models are generally more strongly constrained by the likelihood than BSM models.

In this paper, we simultaneously sample from four models in each scan: one GUT-scale SUSY model (CMSSM, NUHM1 or NUHM2; Sect. 2.1.1), and three specific nuisance models. The first nuisance model includes the parameters of the SM (Sect. 2.1.2), the second parameterises the density and velocity distribution of the DM halo (Sect. 2.1.3), and the third encapsulates the nuclear uncertainties relevant for DM direct detection (Sect. 2.1.4).

### 2.1.1 SUSY models

The definitions of the MSSM superpotential and soft-breaking Lagrangian that we use are specified in Sect. 5.4.3 of Ref. [155], and we follow the conventions established there. All the BSM models that we investigate in this paper are subsets of the **GAMBIT** model MSSM63atMGUT [155], which is the most general formulation of the CP-conserving MSSM, with the soft masses defined at the scale where the gauge couplings $g_1$ and $g_2$ unify (the GUT scale).

The complexity of the MSSM can be reduced considerably if one makes simplifying assumptions about the values of the soft masses at the GUT scale:

**CMSSM** The soft mass parameters at the GUT scale are fixed to a universal scalar mass $m_0$, a universal gaugino mass $m_{1/2}$ and a universal trilinear coupling $A_0$. The diagonal elements in the sfermion mass-squared matrices $\mathbf{m}_Q^2, \mathbf{m}_u^2, \mathbf{m}_d^2, \mathbf{m}_L^2$ and $\mathbf{m}_e^2$ are set to $m_0^2$, all off-diagonal elements are set to zero, and the scalar Higgs mass-squared parameters $\mathbf{m}_{H_u}^2$ and $\mathbf{m}_{H_d}^2$ are set to $m_0^2$. The gaugino masses $M_1, M_2$ and $M_3$ are set to $m_{1/2}$ and all trilinear couplings, $(A_u)_{ij}, (A_d)_{ij}$ and $(A_e)_{ij}$ are set to $A_0$. Electroweak symmetry breaking (EWSB) conditions fix the soft-breaking bilinear $b$ and the magnitude of the superpotential bilinear $\mu$ at the SUSY scale. The remaining free parameters in the Higgs sector are the sign of $\mu$ and the ratio of the vacuum expectation values of the two Higgs doublets $\tan \beta \equiv v_u/v_d$, which is defined at the

<footnote n="1">¹ http://gambit.hepforge.org.</footnote>
scale $m_Z$. The CMSSM is defined by $m_0, m_{1/2}, A_0, \tan \beta(m_Z)$, and sgn($\mu$).

**NUHM1** The GUT-scale constraint on the soft scalar Higgs masses is relaxed, introducing the additional free parameter $m_H$. The soft Higgs masses $m_{H_u}$ and $m_{H_d}$ are not set equal to $m_0$, but instead obey the relation $m_{H_u} = m_{H_d} \equiv m_H$ at the GUT scale. Here $m_H$ is treated as a real dimension-one parameter, ignoring scenarios where $m_H^2 < 0$. This means that at the GUT scale we require $m_{H_u}^2 = m_{H_d}^2 > 0$, which acts as the boundary condition under which the correct shape of the Higgs potential must be radiatively generated at the electroweak scale. The parameters of the NUHM1 are $m_0, m_{1/2}, A_0, \tan \beta(m_Z)$, sgn($\mu$) and $m_H$.

**NUHM2** The constraint on the soft Higgs masses is further relaxed so that $m_{H_u}$ and $m_{H_d}$ become independent, real, dimension-one parameters at the GUT scale. As in the NUHM1, $m_{H_u}^2$ and $m_{H_d}^2$ are always positive at the GUT scale, and the correct shape of the Higgs potential at the electroweak scale must be radiatively generated. The parameters are thus $m_0, m_{1/2}, A_0, \tan \beta(m_Z)$, sgn($\mu$), $m_{H_u}$ and $m_{H_d}$.

We assume throughout that R-parity is conserved, making the lightest supersymmetric particle (LSP) stable. In this paper we consider only the possibility of neutralino LSPs, assigning zero likelihood to all parameter combinations where this is not the case. Sneutrino DM in the MSSM [164] is now essentially ruled out by direct detection, though it remains viable in MSSM extensions (see Ref. [165] for a review). Gravitino LSP scenarios (e.g. [27,166]) are still viable even in the CMSSM, so adding such models to the results that we present here would be an interesting future extension.

The parameter ranges that we scan over for the CMSSM, NUHM1 and NUHM2 can be found in Table 1. We allow the magnitudes of all dimensionful parameters to vary between 50 GeV and 10 TeV. The lower cutoff is motivated by the constraints on sparticle masses from existing searches. The upper cutoff is somewhat arbitrary, but designed to encompass the mass range interesting for solving the hierarchy problem, and for leading to potentially-observable phenomenology. We consider both positive and negative $\mu$, and the full range of $\tan \beta$ over which particle spectra can be consistently calculated and EWSB achieved in such models.

### 2.1.2 Standard model

Here we define the SM as per SLHA2 [167], sampling from the GAMBIT model StandardModel_SLAVE2 [155]. We identify the strength of the strong coupling at the scale of the Z mass, $\alpha_s(m_Z)$, and the top quark pole mass, $m_t$, as the

<table><caption>**Table 1** CMSSM, NUHM1 and NUHM2 parameters, ranges and priors adopted in the different scans contributing to the final results of this paper. The “hybrid” prior for $A_0$ is flat where $|A_0| < 100$ GeV, and logarithmic elsewhere. The “binary” prior for sgn($\mu$) indicates that we repeated every scan for each sign. In addition to the listed priors, we also performed supplementary scans restricted to models with either $m_{\tilde{\tau}_1} < 1.5 m_{\tilde{\chi}_1^0}$ or $m_{\tilde{\tau}_1} < 1.5 m_{\tilde{\chi}_1^0}$. Details can be found in Sec. 2.2.</caption><thead><tr><th>Parameter</th><th>Minimum</th><th>Maximum</th><th>Priors</th></tr></thead><tbody><tr><td><strong>CMSSM</strong></td><td></td><td></td><td></td></tr><tr><td>$m_0$</td><td>50 GeV</td><td>10 TeV</td><td>flat, log</td></tr><tr><td>$m_{1/2}$</td><td>50 GeV</td><td>10 TeV</td><td>flat, log</td></tr><tr><td>$A_0$</td><td>-10 TeV</td><td>10 TeV</td><td>flat, hybrid</td></tr><tr><td>$\tan \beta$</td><td>3</td><td>70</td><td>flat</td></tr><tr><td>$\text{sgn}(\mu)$</td><td>-</td><td>+</td><td>binary</td></tr><tr><td><strong>NUHM1</strong> – as per CMSSM plus</td><td></td><td></td><td></td></tr><tr><td>$m_H$</td><td>50 GeV</td><td>10 TeV</td><td>flat, log</td></tr><tr><td><strong>NUHM2</strong> – as per CMSSM plus</td><td></td><td></td><td></td></tr><tr><td>$m_{H_u}$</td><td>50 GeV</td><td>10 TeV</td><td>flat, log</td></tr><tr><td>$m_{H_d}$</td><td>50 GeV</td><td>10 TeV</td><td>flat, log</td></tr></tbody></table>

most relevant nuisance parameters within this model. Both affect the running of soft-breaking masses from the GUT scale. The mass of the SM-like Higgs boson is also very sensitive to the top quark mass, and has a strong influence on the scan through the Higgs likelihood (see Sect. 3.11).

In all our fits, we allow both these parameters to vary within $\pm 3\sigma$ of their observed central values [168,169]. The resulting parameter ranges are shown in Table 2. We adopt flat priors on both $\alpha_s$ and $m_t$; their values are sufficiently well-determined that the prior has no impact on results. The values of other SM parameters that we keep fixed in our scans are also shown in Table 2.

### 2.1.3 Dark matter halo model

The density and velocity distributions that characterise the DM halo of the Milky Way constitute an important source of uncertainty for astrophysical observations, particularly direct and indirect searches for DM. In this paper, we employ the GAMBIT model Halo_gNFW_rho0 [155] to describe the halo. This consists of a generalised NFW [170] spatial profile, tied to a locally Maxwell-Boltzmann velocity distribution by a specific input local density $\rho_0$.

Because we do not employ any observables in our fits that depend on the Milky Way density profile, the spatial part of this model plays no role. The local distribution of DM velocities $v$ is given by

$$ \tilde{f}(\mathbf{v}) = \frac{1}{N_{\text{esc}}} (\pi v_0^2)^{-3/2} e^{-v^2/v_0^2}, \quad (1) $$

where $v_{\text{esc}}$ is the local Galactic escape velocity, $v_0$ is the most probable particle speed and
<table><caption>**Table 2** Standard Model, dark matter halo and nuclear nuisance parameters and ranges. We vary each of the parameters in the first section of the table simultaneously with CMSSM, NUHM1 or NUHM2 parameters in all of our fits, employing flat priors on each. The parameters listed in the second section of the table are constant in all scans</caption><thead><tr><th>Parameter</th><th></th><th>Value(±Range)</th></tr></thead><tbody><tr><td colspan="3"><strong>Varied</strong></td></tr><tr><td>Strong coupling</td><td>α<sub>s</sub><sup>MS</sup>(m<sub>Z</sub>)</td><td>0.1185(18)</td></tr><tr><td>Top quark pole mass</td><td>m<sub>t</sub></td><td>173.34(2.28) GeV</td></tr><tr><td>Local DM density</td><td>ρ<sub>0</sub></td><td>0.2–0.8 GeV cm<sup>-3</sup></td></tr><tr><td>Nuclear matrix el. (strange)</td><td>σ<sub>s</sub></td><td>43(24) MeV</td></tr><tr><td>Nuclear matrix el. (up + down)</td><td>σ<sub>l</sub></td><td>58(27) MeV</td></tr><tr><td colspan="3"><strong>Fixed</strong></td></tr><tr><td>Electromagnetic coupling</td><td>1/α<sub>s</sub><sup>MS</sup>(m<sub>Z</sub>)</td><td>127.940</td></tr><tr><td>Fermi coupling × 10<sup>5</sup></td><td>G<sub>F,5</sub></td><td>1.1663787</td></tr><tr><td>Z pole mass</td><td>m<sub>Z</sub></td><td>91.1876 GeV</td></tr><tr><td>Bottom quark mass</td><td>m<sub>b</sub><sup>MS</sup>(m<sub>b</sub>)</td><td>4.18 GeV</td></tr><tr><td>Charm quark mass</td><td>m<sub>c</sub><sup>MS</sup>(m<sub>c</sub>)</td><td>1.275 GeV</td></tr><tr><td>Strange quark mass</td><td>m<sub>s</sub><sup>MS</sup>(2 GeV)</td><td>95 MeV</td></tr><tr><td>Down quark mass</td><td>m<sub>d</sub><sup>MS</sup>(2 GeV)</td><td>4.80 MeV</td></tr><tr><td>Up quark mass</td><td>m<sub>u</sub><sup>MS</sup>(2 GeV)</td><td>2.30 MeV</td></tr><tr><td>τ pole mass</td><td>m<sub>τ</sub></td><td>1.77682 GeV</td></tr><tr><td>CKM Wolfenstein parameters</td><td>λ</td><td>0.22537</td></tr><tr><td></td><td>A</td><td>0.814</td></tr><tr><td></td><td>ρ̄</td><td>0.117</td></tr><tr><td></td><td>η̄</td><td>0.353</td></tr><tr><td>Most probable halo speed</td><td>v<sub>0</sub></td><td>235 km s<sup>-1</sup></td></tr><tr><td>Local disk circular velocity</td><td>v<sub>rot</sub></td><td>235 km s<sup>-1</sup></td></tr><tr><td>Local escape velocity</td><td>v<sub>esc</sub></td><td>550 km s<sup>-1</sup></td></tr><tr><td>Up contribution to proton spin</td><td>Δ<sub>u</sub><sup>(p)</sup></td><td>0.842</td></tr><tr><td>Down contrib. to proton spin</td><td>Δ<sub>d</sub><sup>(p)</sup></td><td>-0.427</td></tr><tr><td>Strange contrib. to proton spin</td><td>Δ<sub>s</sub><sup>(p)</sup></td><td>-0.085</td></tr></tbody></table>

$$N_{\text{esc}} \equiv \text{erf}\left(\frac{v_{\text{esc}}}{v_0}\right) - \frac{2v_{\text{esc}}}{\sqrt{\pi}v_0} \exp\left(-\frac{v_{\text{esc}}^2}{v_0^2}\right), \quad (2)$$

is the normalisation factor induced by truncating the distribution at $v_{\text{esc}}$.

In the Earth's rest frame, DM particles have a velocity distribution given by:

$$f(\mathbf{u}, t) = \tilde{f}(\mathbf{v}_{\text{obs}}(t) + \mathbf{u}), \quad (3)$$

where $\mathbf{v}_{\text{obs}}$ is the velocity of the Earth relative to the Milky Way DM halo. This is given by:

$$\mathbf{v}_{\text{obs}} = \mathbf{v}_{\text{LSR}} + \mathbf{v}_{\odot, \text{pec}} + \mathbf{V}_{\oplus}(t), \quad (4)$$

where $\mathbf{v}_{\odot, \text{pec}} = (11, 12, 7) \text{ km s}^{-1}$ is the peculiar velocity of the Sun, which is known with very high precision [171]. The Local Standard of Rest (LSR) in Galactic coordinates moves with a velocity $\mathbf{v}_{\text{LSR}} = (0, v_{\text{rot}}, 0)$, while $\mathbf{V}_{\oplus}(t) =$

29.78 km s⁻¹ [172] denotes the speed of the Earth in the solar rest frame.

For an NFW profile $v_0$ is within 10% of $v_{\text{rot}}$. As shown in Table 2, we set both these parameters to 235 km s⁻¹ [156, 173, 174] in all our scans. Similarly, we adopt a fixed value of 550 km s⁻¹ for the local escape speed [175].

Because it has a substantial impact on direct detection and high-energy solar neutrino signals from DM, we vary the local density of DM as a nuisance parameter in all scans (Table 2). Here we adopt an asymmetric range of +0.4–0.2 GeV cm⁻³ around the canonical value of $\rho_0 = 0.4 \text{ GeV cm}^{-3}$, reflecting the log-normal form of the likelihood that we apply to this parameter (see Sect. 3.1.2). The prior on $\rho_0$ has no impact because it is sufficiently well-constrained by the associated nuisance likelihood; we choose to make it flat.

See Refs. [156,176] for further discussion and details of the DM halo model, parameters and uncertainties.

### 2.1.4 Nuclear model

A final class of uncertainty relevant for direct detection and neutralino capture by the Sun is due to the effective nuclear couplings in WIMP-nucleon cross-sections. For spin-independent interactions, these depend on the light-quark composition of the proton and the neutron. We scan the GAMBIT model nuclear_params_sigmas_signal, parameterising the 6 individual hadronic matrix elements in terms of just two nuclear matrix elements

$$\sigma_l \equiv \frac{1}{2}(m_u + m_d)\langle N | \bar{u}u + \bar{d}d | N \rangle, \quad (5)$$

$$\sigma_s \equiv m_s \langle N | \bar{s}s | N \rangle, \quad (6)$$

which we take to be identical for $N = p$ and $N = n$ [177]. These two parameters respectively describe the light-quark and strange-quark contents of the nucleus. We vary $\sigma_l$ and $\sigma_0$ over their $\pm 3\sigma$ ranges in all fits. Discussion of the values and uncertainties of these parameters can be found in Sect. 3.1.3 and the DarkBit paper [156]. Like all other nuisance parameters listed in Table 2, the nuclear matrix elements are sufficiently well constrained that the prior is irrelevant, so we choose it to be flat.

The spin-dependent couplings are described by the spin content of the proton and neutron $\Delta_q^{(N)}$ for each light quark $q \in \{u, d, s\}$. As the values for the proton and neutron are related, only three of these parameters are independent. As listed in Table 2, we specify the values for the proton, and set them to the central values discussed in Ref. [156].

## 2.2 Scanning methodology

In this paper we carry out a number of different scans of each of the three GUT-scale models, employing multiple pri-
<table><caption>**Table 3** Samplers and their settings for the different scans of this paper</caption><thead><tr><th>Scanner</th><th>Parameter</th><th>Setting</th></tr></thead><tbody><tr><td rowspan="2">Diver</td><td>NP</td><td>19200</td></tr><tr><td>convthresh</td><td>10<sup>-5</sup></td></tr><tr><td rowspan="2">Diver<br>(co-annihilation)</td><td>NP</td><td>6000</td></tr><tr><td>convthresh</td><td>10<sup>-4</sup></td></tr><tr><td rowspan="2">MultiNest</td><td>nlive</td><td>5000</td></tr><tr><td>tol</td><td>0.1</td></tr></tbody></table>

ors, sampling algorithms and settings. We then merge the results of all scans for each model, in order to obtain the most complete sampling of the profile likelihood possible. We leave discussion and presentation of Bayesian posteriors for a future paper, as they remain strongly dominated by the choice of prior even in such low-dimensional versions of the MSSM, and a detailed analysis of their implications for fine-tuning and naturalness (e.g. [178,179]) is beyond the scope of the current paper.

The parameter ranges and priors that we employ in scans of each model are listed in Table 1. We repeat every scan for positive and negative values of $\mu$. We carry out scans with both flat and logarithmic priors on all dimensionful parameters. In the case of the trilinear coupling $A_0$, which may be positive or negative, in our log-prior scans we employ a hybrid prior (log_flat_join in the language of Ref. [160]), consisting of a symmetric logarithmic prior at large $|A_0|$, truncated to a flat prior at $|A_0| < 100$ GeV.

We use two different samplers for our scans: Diver 1.0.0 [160] and MultiNest 3.10 [180]. The settings that we use for each can be found in Table 3.

Diver is a self-adaptive sampler based on differential evolution [181]. It samples the profile likelihood far more efficiently than traditional algorithms [160], allowing high-quality profile likelihoods to be computed in a fraction of the time of previous SUSY global fits, using significantly fewer likelihood samples. As a result, the majority of our results are driven by the Diver scans. Following the extensive tests discussed in Ref. [160], for most scans<fnref n="2" /> we choose a population size of NP = 19,200 and a convergence threshold of convthresh = $10^{-5}$. The latter is defined in terms of the smoothed fractional improvement in the mean likelihood across the entire population. Other than these two parameters, we employ Diver with the default settings defined in Scan-nerBit [160]. In particular, this includes the $\lambda$jDE version (introduced in Ref. [160]) of the self-adaptive jDE algorithm [182].

MultiNest is an implementation of nested sampling [183], a method optimised for the calculation of the Bayesian evidence. As a by-product, it also produces posterior samples, which it obtains via likelihood evaluation. It can therefore be very useful for sampling profile likelihoods as well, especially for smoothly mapping isolikelihood contours. However, it typically requires a rather long runtime to properly find the global best fit and any highly-localised likelihood modes [160, 184, 185]. We employ it here mainly to bulk out our sampling of the main likelihood mode of each scan a little, in regions where the profile likelihood is comparatively flat. For this purpose, we run MultiNest with relatively loose settings, choosing nlive = 5000 live points and a stopping tolerance of tol = 0.1. The tolerance is given in terms of the estimated fractional remaining unsampled evidence.

For more details on the performance of the two scanning algorithms, and comparisons to others, please see the Scan-nerBit paper [160].

To more densely sample the narrow strips in parameter space where neutralino-sfermion co-annihilations play an important role in determining the relic density of DM, we also carry out two specially-targeted versions of each log-prior scan. In these scans, we restrict the mass of either the lightest slepton or the lightest squark to within 50% of the mass of the lightest neutralino, i.e. $m_{\tilde{l}_1} \le 1.5 m_{\tilde{\chi}_1^0}$ or $m_{\tilde{q}_1} \le 1.5 m_{\tilde{\chi}_1^0}$. Although these additional scans are not necessary for finding the sfermion co-annihilation regions (Diver typically uncovers these regions anyway in untargeted scans), they are useful for ensuring that the boundaries of these regions are mapped thoroughly.

We carry out the additional scans using Diver only, building the initial population exclusively from models that satisfy the mass-ratio cut, before evolving it as usual. As it takes many random draws from the prior to successfully build such an initial population, we run these scans with a reduced population of NP = 6000, and a looser convergence criterion (convthresh = $10^4$) than the untargeted equivalents.

This results in a total of 3 models × 2 sgn($\mu$) × (2 priors × 2 scanners + 2 targeted co-annihilation scans) = 36 separate scans. For the CMSSM, NUHM1 and NUHM2, this results in a total of 71, 94 and 117 million viable samples, respectively. Each of these 36 scans typically took 1–3 days to run on 2400 modern (Intel Core i7) supercomputer cores.

In all profile likelihood plots that we show in the paper, we sort the samples of our scans into 60 bins across the range of data values that they cover in each direction. We then interpolate with a bilinear scheme to a finer resolution of 500 when plotting [186]. This expressly avoids any smoothing of the resulting profile likelihoods, which would amount to manipulation of the likelihoods of our samples. The binning

<footnote n="2">² For the special case of flat-prior CMSSM scans, where less stringent parameters already provide quite sufficient sampling, we use NP = 14,400, convthresh = $10^{-4}$.</footnote>
and interpolation process can produce some cosmetic artefacts, in particular a sawtooth pattern in regions where the likelihood drops off sharply. Using a fixed number of bins across the data range (rather than the plot range) can also sometimes produce surprising effects when plotting multiple regions on the same axes, as the same region will typically appear smaller and ‘smoother’ for subsets of the data where all samples lie within a small range of parameter values than for subsets with samples spread across the entire plot plane (as the size of an individual bin is much larger in the latter case). This should be kept in mind especially when viewing plots of mass correlations for multiple co-annihilation mechanisms and comparing preferred regions to LHC sensitivity curves (e.g. Fig. 15).

# 3 Observables and likelihoods

## 3.1 Nuisance likelihoods

### 3.1.1 Standard model

We include independent Gaussian likelihoods for each of the two SM nuisance parameters in our scans. We evaluate the strong coupling $\alpha_s$ at the scale $\mu = m_Z$ in the $\overline{MS}$ scheme, and compare with $\alpha_s(m_Z) = 0.1185 \pm 0.0006$ from lattice QCD [169]. We interpret the quoted uncertainty as a $1\sigma$ confidence interval, and do not incorporate any additional theoretical uncertainty.

For the top quark pole mass $m_t$, we compare with the combined measurements of experiments at the Tevatron and LHC: $m_t = 173.34 \pm 0.27(\text{stat}) \pm 0.71(\text{syst})$ GeV, with a total uncertainty of 0.76 GeV [168]. We do not assign any separate systematic error to our interpretation of the experimental result as the top pole mass.

### 3.1.2 Local halo model

The canonical local density of DM extracted from fits to stellar kinematic data is $\bar{\rho}_0 = 0.4 \text{ GeV cm}^{-3}$ (see e.g. [187,188]). Because arbitrarily small or negative densities are unphysical, we adopt a log-normal distribution for the likelihood of $\rho_0$,

$$ \mathcal{L}_{\rho_0} = \frac{1}{\sqrt{2\pi}\sigma'_{\rho_0}\rho_0} \exp\left(-\frac{\ln(\rho_0/\bar{\rho}_0)^2}{2\sigma'^2_{\rho_0}}\right), \quad (7) $$

where $\sigma'_{\rho_0} = \ln(1 + \sigma_{\rho_0}/\rho_0)$ and $\sigma_{\rho_0}$ is taken to be 0.15 GeV/cm$^3$. We refer the reader to Ref. [155] for additional implementation details of the GAMBIT log-normal likelihood, and Refs. [156, 176, 189] for a more extended discussion of the central value and uncertainty on this parameter.

### 3.1.3 Nuclear matrix elements

We constrain the nuclear matrix elements $\sigma_s$ and $\sigma_l$ using Gaussian likelihood functions, with central $\pm 1\sigma$ values of 43 $\pm$ 8 and 58 $\pm$ 9 MeV, respectively. The former is based on lattice calculations [190], whereas the latter is a weighted average of a number of different results in the literature [156].

## 3.2 Spectrum calculation

We use FlexibleSUSY 1.5.1 [191] to compute the mass spectrum of the MSSM. This code obtains model-dependent information from SARAH [192,193], and borrows some numerical routines from SOFTSUSY [194,195]. FlexibleSUSY employs full three-family, two-loop renormalisation group equations (RGEs) and full one-loop self-energies and tadpoles. In addition, it computes the Higgs mass using two-loop corrections at $O(\alpha_t\alpha_s)$, $O(\alpha_b\alpha_s)$, $O(\alpha_t^2)$, $O(\alpha_b^2)$, $O(\alpha_t^2)$ and $O(\alpha_t\alpha_b)$ [196–199].

Large logarithms appear when the supersymmetric spectrum is very heavy. To improve precision, these can be resummed using techniques from effective field theory (EFT) [44, 46, 200–203]. This method has been implemented in several public codes [44, 46, 202]. However, the hierarchical spectrum assumed in the EFT calculation only appears in small subspaces of the models over which we scan in this paper. The public codes that implement this calculation, and that are suitable for cluster-scale parameter scans, are rather new.$^3$ It was also pointed out in Ref. [44] that the accuracy of the fixed-order Higgs mass prediction in FlexibleSUSY is much better than one would naively expect at large sparticle masses, due to accidental cancellations. We have therefore retained the fixed-order calculations in FlexibleSUSY for calculating the Higgs mass in this paper.

## 3.3 Relic density of dark matter

The thermal relic abundance of the lightest neutralino is a strong constraint on the MSSM. Many parameter combinations lead to more DM than the cosmological abundance observed by Planck, which is $\Omega_{\text{c}}h^2 = 0.1188 \pm 0.0010$ [204]. For the relic density not to exceed this value, if the lightest neutralino is heavier than $\sim 100$ GeV, typically one or more specific depletion mechanisms must be active. These include co-annihilation of light sfermions or charginos with the lightest neutralino, and resonance or ‘funnel’ effects, where the lightest neutralino has a mass very close to half that of another neutral species.

<footnote>$^3$ For example, during this work we found a bug in the resummation that affected results at low masses when testing with FeynHiggs 2.12.0 – though this should be corrected in later versions.</footnote>
We compute the relic density of each model taking into account DM annihilation to all two-body final states, including full co-annihilation [205,206], thermal and resonance effects, using the native DarkBit relic density calculator [156], connected to various subroutines of DarkSUSY 5.1.3 [207]. We obtain the effective annihilation rate $W_{\text{eff}}$ from DarkSUSY, passing all spectrum, decay and SM information from GAMBIT, and considering co-annihilations with particles up to 60% heavier than the lightest neutralino. We also employ the DarkSUSY Boltzmann solver, setting the option fast = 1. This ensures that the relic density calculation for most models takes less than a second. This setting controls the convergence criteria of the Boltzmann solver, and is the recommended option unless accuracy of better than 1% is required.

The likelihood that we employ penalises only models that predict more than the observed relic density. The likelihood function is a half Gaussian (see Ref. [156] and Sect. 8.3 of Ref. [155]), centered on the Planck observation but treating it as an upper limit. Consistent with our choice of the fast parameter for the Boltzmann solver, we retain the DarkBit default theoretical uncertainty of 5% on the relic density, adding it in quadrature to the observational error. Further discussion of this number in the context of higher-order corrections can be found in Refs. [3, 18, 156, 208–212].

## 3.4 Gamma rays from dark matter annihilation

Neutralino annihilation in astrophysical objects would produce a variety of final states, leading to both prompt gamma rays and those produced as final decay products. Dwarf spheroidal galaxies are particularly important targets, as they are strongly dominated by dark rather than visible matter, and exhibit little or no astrophysical gamma-ray emission. Limits from gamma-ray observations of dwarf galaxies have therefore played an increasingly important role in global fits, e.g. [132,213–215].

For most neutralino masses, the most stringent gamma-ray limits on DM annihilation come from joint analyses of multiple Milky Way satellite galaxies [216–220] using data from the Fermi Large Area Telescope (Fermi-LAT). We employ likelihoods from the analysis of six years of Pass 8 data in the direction of 15 dwarf spheroidal galaxies [220], as implemented in gamLike [156].

gamLike constructs a composite likelihood

$$ \ln \mathcal{L}_{\text{exp}} = \sum_{k=1}^{N_{\text{dSph}}} \sum_{i=1}^{N_{\text{ebin}}} \ln \mathcal{L}_{ki}(\Phi_i \cdot J_k), \quad (8) $$

from gamma-ray data sorted into $N_{\text{dSph}}$ fields of view (one for each dwarf) and $N_{\text{ebin}}$ energy bins. The partial likelihoods $\mathcal{L}_{ki}$ describe the likelihood of obtaining the observed number

of photons in the $i$th energy bin from the $k$th dwarf. The energy-dependent factor

$$ \Phi_i = \frac{\langle \sigma v \rangle_0}{8\pi m_{\chi}^2} \int_{\Delta E_i} dE \frac{dN_{\gamma}}{dE} \quad (9) $$

depends on the MSSM model, whereas the astrophysical factor

$$ J_k = \int_{\Delta\Omega_k} d\Omega \int_{\text{l.o.s.}} ds \, \rho_{\chi}(s)^2 \quad (10) $$

is a model-independent property of each dwarf galaxy. Here the differential gamma-ray multiplicity per annihilation is $dN_{\gamma}/dE$, the zero-velocity annihilation cross-section is $\langle \sigma v \rangle_0 \equiv \sigma v|_{v \to 0}$, $m_{\chi}$ is the DM mass, and $\rho_{\chi}(s)$ is the DM density along the line of sight parameter $s$ in a given dwarf. $\Delta E_i$ is the width of the $i$th energy bin, and $\Delta\Omega_k$ is the solid angle around the position of the $k$th dwarf over which gamma-ray data are being considered.

As in [220], gamLike profiles over the $J_k$-factors as nuisance parameters, giving a final likelihood of

$$ \ln \mathcal{L}_{\text{dwarfs}}^{\text{prof.}}(\Phi_i) = \max_{J_1 \dots J_k} (\ln \mathcal{L}_{\text{exp}} + \ln \mathcal{L}_J), \quad (11) $$

where

$$ \ln \mathcal{L}_J = \sum_{k=1}^{N_{\text{dSph}}} \ln \mathcal{N}(\log_{10} J_k | \log_{10} \hat{J}_k, \sigma_k). \quad (12) $$

Here the probability distribution for each $J_k$ is assumed to follow an independent log-normal distribution with mean $\hat{J}_k$ and width $\sigma_k$.

We compute the predicted spectrum $dN_{\gamma}/dE$ for each model by combining tabulated two-body annihilation spectra from DarkSUSY [207] with yields computed on the fly with the DarkBit Fast Cascade Monte Carlo [156]. To this, we add the dominant contribution from photon internal bremsstrahlung [221]. For each parameter combination, we rescale the expected gamma-ray flux by the squared ratio of the predicted relic density to the observed value, allowing for the fact that neutralinos may only be a fraction of DM. We limit this scaling factor to 1, not rescaling signals when the predicted relic density is greater than the observed value.

## 3.5 High-energy neutrinos from dark matter annihilation in the Sun

The Sun is expected to capture neutralinos from the local halo by nuclear scattering. Subsequent neutralino annihilation and interaction of the annihilation products in the solar core would produce GeV-energy neutrinos, which may
be detectable at the Earth. The rate-limiting step for all MSSM models is the capture process, which depends sensitively on both the spin-dependent and spin-independent nuclear scattering cross-sections. The dominant constraints on intermediate and high-mass neutralino annihilation in the solar interior currently come from the IceCube experiment [222,223]. We use the DarkBit interface [156] to nulike 1.0.4 [7,224], which computes an unbinned likelihood from the event-level energy and angular information contained in the three independent event selections of the 79-string IceCube dataset [222]. We predict the neutrino spectra at Earth using DarkSUSY 5.1.3, which contains tabulated results previously obtained from WimpSim [225].<fnref n="4" /> Slightly stronger limits are also available from the 86-string dataset [223], but not in a format that allows them to be accurately applied to MSSM models.

### 3.6 Direct detection of dark matter

The dominant direct DM constraints on the models in this paper come from the LUX [226–228], Panda-X [229] and PICO [230,231] experiments. We also include likelihoods from XENON100 [232], SuperCDMS [233] and SIMPLE [234]. A new analysis from PICO-60 [235] appeared after much of this paper was already finalised, but the majority of MSSM models susceptible to that limit are already probed in our scans by the IceCube 79-string likelihood (Sect. 3.5). We do not include the recent XENON1T result [236], but given that its sensitivity improvement relative to LUX is smaller than the error in our likelihood approximation [156], this will not impact our results.

For each experimental search and combination of MSSM, halo and nuclear parameters, we use the likelihood functions contained in DDCalc [156] to compute a Poisson likelihood,

$$ \mathcal{L}_i(N_{p,i} | N_{o,i}) = \frac{(b_i + N_{p,i})^{N_{o,i}} e^{-(b_i + N_{p,i})}}{N_{o,i}!} \quad (13) $$

Here $N_{o,i}$ is the number of observed events in the analysis region of the $i$th experiment, $b_i$ is the expected number of background events, and $N_{p,i}$ is the expected number of signal events. DDCalc computes the latter by interpolating in pre-computed efficiency tables, which include both detector and acceptance effects. The signal prediction takes into account both the spin-dependent and spin-independent interactions expected from each MSSM model. We compute the DM-nucleon couplings for each MSSM model using DarkSUSY 5.1.3 [207].

We scale the direct detection yields for each parameter combination by the ratio of the predicted relic density to the

value observed by Planck [204], allowing for the fact that neutralinos may not constitute all of DM. We do not rescale direct detection rates when the predicted relic density is larger than the observed value.

### 3.7 Electroweak precision observables

We include likelihoods from PrecisionBit [159] for the W mass and the anomalous magnetic moment of the muon $a_μ$. These functions construct a basic Gaussian likelihood based on the difference between the calculated and measured value, and combine theoretical and experimental uncertainties in quadrature.

The W mass must be recalculated using the details of the SUSY spectrum. In the present scans, the value of $m_W$ comes from FlexibleSUSY. SpecBit assigns a theoretical uncertainty of 10 MeV to this quantity, based on the size of two-loop corrections [159]. PrecisionBit compares these to $m_W = 80.385 \pm 0.015$ GeV [169], based on mass measurements and uncertainties from the Tevatron and LEP experiments.

For $a_μ$, we assume an SM contribution of $a_{μ,SM} = (11659180.2 \pm 4.9) \times 10^{-10}$, which comes from theoretical calculations based on $e^+e^-$ data [237]. We evaluate the supersymmetric contribution using GM2Calc 1.3.0 [92], which determines an uncertainty on its result by estimating the magnitude of neglected higher-order corrections using the two-loop Barr-Zee corrections [238]. The total predicted value is the sum of the SM and MSSM contributions, and the total uncertainty the sum in quadrature of their individual uncertainties. We compare this with the experimental measurement of $a_μ = (11659208.9 \pm 6.3) \times 10^{-10}$ [239,240], where the experimental error is the sum in quadrature of the systematic $(3.3 \times 10^{-10})$ and statistical $(5.4 \times 10^{-10})$ contributions.

### 3.8 Flavour physics likelihoods

Scans in this paper include 59 flavour observables from FlavBit [158]. These are sorted into four different categories for likelihood calculation:

1. Tree-level leptonic and semi-leptonic B and D meson decays (8 observables). Branching fractions for $B \to D\mu\nu_\mu$, $B \to D^*\mu\nu_\mu$, $B \to \tau\nu_\tau$, $D \to \mu\nu_\mu$, $D_s \to \mu\nu_\mu$ and $D_s \to \tau\nu_\tau$, as well as ratios $R_D \equiv B(B \to D\tau\nu_\tau)/B(B \to D\ell\nu_\ell)$ and $R_{D^*} \equiv B(B \to D^*\tau\nu_\tau)/B(B \to D^*\ell\nu_\ell)$. Here either $\mu$ or $e$ may be substituted for $\ell$, as both are effectively massless in the B-meson system.

2. Electroweak penguin decays (48 observables). Eight observables ($AFB, FL, S_3, S_4, S_5, S_7, S_8$ and $S_9$) for the

<footnote n="4">⁴ WimpSim is available at www.fysik.su.se/~edsjo/wimpsim/.</footnote>
decay $B^0 \to K^{*0} \mu^+ \mu^-$, each in six different angular $(q^2)$ bins.

3. Rare leptonic B decays (2 observables). Branching fractions for $B^0 \to \mu^+ \mu^-$ and $B_s^0 \to \mu^+ \mu^-$.

4. The branching fraction for $B \to X_s \gamma$, for photon energies $E_\gamma > 1.6$ GeV (1 observable).

All observable predictions draw on **SuperIso 3.6** [241,242]. We have not included the $B_s-\bar{B}_s$ meson mass difference $\Delta M_s$, owing to the fact that it is only calculable within **FlavBit** via **FeynHiggs**, which we otherwise avoided for the scans of this paper in the interests of speed, and due to worries about its most recent versions' accuracy in parts of the parameter space (some details of which have been mentioned earlier in this section).

Recent LHCb results in the exclusive modes have already provided substantial additional constraints as compared to the available inclusive results from the B factories. In particular, several angular observables in the $B^0 \to K^{*0} \mu^+ \mu^-$ decay have been measured for the first time.

We construct a separate likelihood function for observables in each of the four categories above, including correlated uncertainties on observables within each category wherever warranted. The likelihood functions consider correlations between experimental measurement errors separately from correlations between theoretical errors (arising from e.g. common scale or form factor uncertainties), and then sum them to obtain the final covariance matrix. **FlavBit** then computes the likelihood within each category using a $\chi^2$ approximation,

$$ \ln \mathcal{L} = -\frac{1}{2} \chi^2 = -\frac{1}{2} \sum_{i,j=1}^{N} (y_i - x_i)(V^{-1})_{ij}(y_j - x_j), \quad (14) $$

where $x_i$ and $y_i$ are the experimental measurements and theoretical predictions, respectively, and $V$ is the covariance matrix.

In the first likelihood category, **FlavBit** includes experimental measurements, correlations and combinations from Refs. [243–245] and theoretical uncertainties from Refs. [246,247], supplemented by our own additional calculations with a beta version of **SuperIso 3.7**. The experimental measurements and correlated uncertainties of $B^0 \to K^{*0} \mu^+ \mu^-$ angular observables come from LHCb [248], and the theoretical errors and correlations from Refs. [249,250]. Data for rare leptonic decays are the latest from LHCb and CMS [251,252], and theoretical uncertainties come from Ref. [253]. For $B \to X_s \gamma$, we use the latest average [254] of measurements by Belle [255,256] and Babar [257–259], and a theoretical uncertainty of 7% [260,261]. More details can be found in the **FlavBit** paper [158].

<table><caption>**Table 4** Results from LEP on sparticle pair production used in the scans of this paper. Here $\bar{l} = \tilde{e}_L, \tilde{e}_R, \tilde{\mu}_L, \tilde{\mu}_R, \tilde{\tau}_1$ or $\tilde{\tau}_2$</caption><thead><tr><th>Production</th><th>Decay</th><th>Experiment</th></tr></thead><tbody><tr><td>$\bar{l}\bar{l}^*$</td><td>$\bar{l} \to l\tilde{\chi}_1^0 + \text{c.c.}$</td><td>ALEPH [262], L3 [263]</td></tr><tr><td>$\tilde{\chi}_i^0 \tilde{\chi}_1^0$<br>(i = 2, 3, 4)</td><td>$\tilde{\chi}_i^0 \to q\bar{q} \tilde{\chi}_1^0$</td><td>OPAL [264]</td></tr><tr><td>$\tilde{\chi}_i^+ \tilde{\chi}_i^-$<br>(i = 1, 2)</td><td>$\tilde{\chi}_i^+ \tilde{\chi}_i^- \to \bar{l}\bar{l} \tilde{\chi}_1^0$<br>$\tilde{\chi}_i^+ \tilde{\chi}_i^- \to q\bar{q}'q\bar{q}' \tilde{\chi}_1^0 \tilde{\chi}_1^0$<br>$\tilde{\chi}_i^+ \tilde{\chi}_i^- \to q\bar{q}'l\nu \tilde{\chi}_1^0 \tilde{\chi}_1^0$<br>$\tilde{\chi}_i^+ \tilde{\chi}_i^- \to l\nu l\nu \tilde{\chi}_1^0 \tilde{\chi}_1^0$</td><td>L3 [265]<br>OPAL [264]<br>OPAL [264]<br>OPAL [264], L3 [265]</td></tr></tbody></table>

## 3.9 Searches for superpartners at LEP

Even though they are typically overshadowed by constraints from the LHC, LEP searches can have a significant impact in some parts of the parameter spaces that we consider in this paper. This is especially true for light, highly-degenerate spectra. Direct limits on sparticle production at LEP have typically taken the form of hard lower limits on sparticle masses, at e.g. 95% CL, computed with model-dependent assumptions [207,266]. **ColliderBit** instead uses the individual cross-section limits on pair production of neutralinos, charginos and sleptons from the ALEPH, L3 and OPAL experiments, as a function of the sparticle masses.

For each MSSM parameter combination, we compute the pair-production cross-sections at LEP for the processes given in Table 4, using the cross-section calculations included in **ColliderBit** and based on the results of Refs. [267–270]. We take the relevant sparticle decay branching fractions from **DecayBit** (choosing to obtain widths from a suitably-patched SUSY-HIT 1.5 [159,271]), and calculate the product of the cross-section and branching fraction for each process. This number can then be compared to digitised LEP cross-section limits in the plane of $m_{\tilde{\chi}_1^0}$ and the mass of the directly-produced sparticle, interpolating when the masses do not fall exactly on a grid point. This takes care of the mass-dependent experimental acceptance for each parameter point. We then calculate the likelihood of the experimental result assuming a Gaussian form, accounting for the dominant theoretical uncertainty on the signal prediction by varying the mass of the pair-produced sparticles within the uncertainties provided by **SpecBit**. Finally, we multiply the likelihoods from the various experiments and channels, taking them as independent measurements.

Further details of the cross-section and likelihood calculations can be found in the **ColliderBit** paper [157].

## 3.10 Searches for supersymmetry at the LHC

Many searches for supersymmetric particles have been performed at the LHC by ATLAS and CMS, in a variety of final
states arising from proton–proton collisions at $\sqrt{s} = 7, 8$ and 13 TeV [272,273]. The results of all searches to date are consistent with the predictions of the SM, placing strong, model-dependent constraints on the mass spectrum of the MSSM. Taking into account the complete list of LHC searches is impractical; here we implement the most constraining analyses:

1. 0-lepton supersymmetry searches (ATLAS & CMS, Run I & Run II). These provide the best constraints on models with a light gluino or one or more light squarks. The analyses look for an excess of events in final states with jets and missing energy, using a variety of kinematic variables [274–276].

2. Third generation squark searches (ATLAS & CMS, Run I). These searches target stop pair production, with subsequent decay to either a top quark and the lightest neutralino, or to a b quark and a chargino. We include the results of ATLAS searches in 0-, 1- and 2-lepton final states [277–279], and CMS searches for 1- and 2-lepton final states [280,281]. We also include the ATLAS search for direct sbottom production in final states with b-jets and missing energy [282].

3. Multilepton supersymmetry searches (ATLAS and CMS, Run I). We include 2- and 3-lepton searches by ATLAS [283,284] and the 3-lepton search by CMS [285]. These are typically the most constraining searches for direct production of charginos and neutralinos, and the 2 lepton search is also sensitive to slepton pair production and decay.

4. Dark matter searches (CMS, Run I). We include the CMS monojet search [286], which constrains supersymmetric particle production in the case of compressed mass spectra.

We use ColliderBit to calculate the expected signal yield for each combination of model parameters, in each analysis region, using the external Monte Carlo (MC) event generator Pythia 8 [287,288], the native ColliderBit detector parameterisation BuckFast [157], and the ColliderBit implementation of the analysis cuts applied in each LHC paper. ColliderBit contains a number of code optimisations of the Pythia 8 routines, including parallelisation of the main event loop via OpenMP. These modifications make it feasible to run 20,000 MC events per parameter combination during the global fit itself, as we do here. Due to the computational cost of calculating next-to-leading order (NLO) cross sections, we normalise the signal yields using leading-order (LO) plus leading-log (LL) cross-sections only, as provided by Pythia 8. For a more exhaustive discussion of this choice see the ColliderBit paper [157].

In a specific signal region with a predicted number of signal events s and an expected number of background

events b, the likelihood of observing n events is described in ColliderBit by a marginalised form of the Poisson likelihood [213,224,289],

$$ \mathcal{L}(n|s, b) = \int_{0}^{\infty} \frac{[\xi(s+b)]^n e^{-\xi(b+s)}}{n!} P(\xi) d\xi, \quad (15) $$

where $\xi$ is a scaling variable with a probability distribution centred on 1, designed to describe the effective rescaling of the signal + background prediction due to systematic uncertainties. Marginalising over $\xi$ this way, it is possible to include the effects of fractional systematic uncertainties on both the signal prediction ($\sigma_s$) and the background estimate ($\sigma_b$).<fnref n="5" />

We assume a log-normal distribution for $\xi$,

$$ P(\xi|\sigma_{\xi}) = \frac{1}{\sqrt{2\pi}\sigma_{\xi}} \frac{1}{\xi} \exp\left[-\frac{1}{2}\left(\frac{\ln\xi}{\sigma_{\xi}}\right)^2\right], \quad (16) $$

where $\sigma_{\xi}^2 = \sigma_s^2 + \sigma_b^2$. We compute this integral using the highly-optimised implementation in nulike 1.0.4 [7,224].

The analyses listed above are statistically independent, either because they use a completely independent dataset (based on collisions at ATLAS versus CMS, or during Run I versus Run II), or because they utilise signal regions that have no overlap with the signal regions of any of the other searches. This allows us to simply multiply the likelihoods of all analyses in order to arrive at a combined likelihood.

However, within each analysis, signal regions are not always orthogonal, i.e. some contain events or significant systems in common. Given that there is no public information describing the correlations across these signal regions,<fnref n="6" /> we calculate the likelihood for an analysis based on the signal region that is expected to give the strongest limit. We determine the expected limit from each signal region by computing the expected ratio between the signal plus background and background-only likelihoods, in the hypothetical scenario where the observed number of events is exactly equal to the background expectation,

$$ \Delta \ln \mathcal{L}_{\text{pred}} = \ln \mathcal{L}(n = b|s, b) - \ln \mathcal{L}(n = b|s = 0, b). \quad (17) $$

<footnote n="5">⁵ Due to our use of LO cross-sections, including a signal systematic associated with finite MC statistics is in practice rather pointless, as with 20,000 simulated events this is basically always dwarfed by the systematic error associated with neglecting NLO corrections. Considering that the LO cross-sections in the MSSM are known to almost always lie significantly below the NLO cross-section, our approach is in any case very conservative. In the present scan we have thus set $\sigma_s = 0$. For details, see the ColliderBit paper [157].</footnote>

<footnote n="6">⁶ This is at least true for most analyses; Refs. [290,291] are notable recent exceptions.</footnote>
Taking the difference with respect to the background log-likelihood prevents erroneous model-to-model jumps in the likelihood function (see Ref. [157] for more details).

Given the absence of published correlations between the yields (and uncertainties) in the various signal regions, this is arguably the best possible treatment, and it has the added merit of giving conservative results. Because no significant excess has been observed in any of the LHC searches that we include, we restrict the combined LHC Run I and combined Run II log-likelihood each to a maximum of 0, i.e. forbidding mildly better fits than the SM (which are achievable via statistical fluctuations in the data or Monte Carlo simulation, at a little less than the $1\sigma$ level).

We included all Run I searches listed above directly in our main scans of the CMSSM, NUHM1 and NUHM2. We then applied the likelihoods associated with the 13 TeV, 13 fb⁻¹ Run II ATLAS and CMS 0-lepton searches in a postprocessing step, using the ScannerBit postprocessor scanner (see Sect. 6 of Ref. [160]). These searches uncovered no excesses, and therefore do not change the regions preferred by our scans except to disfavour a strip of additional models (compared to the Run I searches) at sparticle masses of a few hundred GeV. The accuracy of our sampling is therefore unaffected by their inclusion via postprocessing rather than in the original scans.<fnref n="7" />

### 3.11 Higgs physics

We use likelihoods from HiggsBounds 4.3.1 [292–294] and HiggsSignals 1.4.0 [295], as interfaced via ColliderBit [157]. These provide two likelihood terms: one based on limits from LEP, and the other on measurements of Higgs masses and signal strengths at the LHC (plus some subdominant contributions from the Tevatron).

The combined LEP Higgs likelihood is an approximate Gaussian likelihood, valid in the asymptotic limit. HiggsBounds constructs this from the full CL$_{s+b}$ distribution, accounting for the effect of varying production cross-sections and Higgs masses by interpolating in a grid of pre-calculated values.

The LHC Higgs likelihood is based on mass and signal-strength measurements reported by ATLAS and CMS. The mass and signal-strength data contribute separate $\chi^2$ terms to the overall LHC Higgs log-likelihood. For each channel where a mass measurement is available, a $\chi^2$ contribution is calculated for the hypothesis that each neutral Higgs particle is responsible for the observed 125 GeV boson [296,297]. Only the minimum value enters the final likelihood. This min-

imisation allows for the possibility that multiple resonances exist at 125 GeV with near-degenerate masses. The signal-strength contribution to the $\chi^2$ uses a covariance matrix that contains all published experimental uncertainties on all measurements of signal strengths, including their correlations.

As discussed in Sect. 3.2, we obtain theoretical predictions of Higgs masses from FlexibleSUSY, adopting an uncertainty of 2 GeV on the mass of the lightest neutral Higgs, and 3% on all other Higgses [159]. We compute Higgs decay rates and branching fractions using SUSY-HIT 1.5 [271] via DecayBit [159]. To obtain the neutral Higgs boson production cross sections, we employ an effective coupling approximation, assuming that the BSM-to-SM ratios of Higgs production cross sections are equal to the ratios of the relevant squared couplings. We determine the coupling ratios using the partial width approximation, in which the ratios of squared BSM-to-SM couplings are taken to be equal to the ratios of the equivalent partial decay widths. To obtain branching fractions for SM-like Higgs bosons of equivalent mass to those in our MSSM models, we use lookup tables computed with HDECAY 6.51 [298,299]. More details can be found in the DecayBit paper [159].

## 4 Results

### 4.1 CMSSM

In the left panel of Fig. 1, we show the joint profile-likelihood ratio for the mass of lightest neutralino and the relic density in the CMSSM. In the right panel, we show the same 95% CL regions colour-coded according to the possible mechanisms by which different models may avoid exceeding the observed relic density of DM. We classify these regions as follows:

- stau co-annihilation: $m_{\tilde{\tau}_1} \le 1.2 m_{\tilde{\chi}_1^0}$,

- stop co-annihilation: $m_{\tilde{\tau}_1} \le 1.2 m_{\tilde{\chi}_1^0}$,

- chargino co-annihilation: $\tilde{\chi}_1^0 \ge 50\%$ Higgsino,

- A/H-funnel: $1.6 m_{\tilde{\chi}_1^0} \le m_{\text{heavy}} \le 2.4 m_{\tilde{\chi}_1^0}$,

where ‘heavy’ may be $A^0$ or $H^0$, i.e. a model qualifies if either Higgs is in range.

We emphasise that this classification is not exclusive. The labels that we give to these regions are merely a convenient shorthand for the precise mass/composition relations that we give above. In particular, they should not be interpreted as definitive indications that a specific mechanism is solely (nor even predominantly) responsible for setting the relic density of the neutralino. These relations indicate necessary but not sufficient conditions for a given mechanism to play a significant role in setting the relic density. The colour-coding in Fig. 1 (right) is done on the basis of the subset of the points

<footnote n="7">⁷ We applied the Run II searches this way not for reasons of computational speed, but just as a matter of practicality, given when supercomputing time, Run II results and different components of GAMBIT respectively became available.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Profile likelihoods and confidence regions for the CMSSM, in terms of the mass and thermal relic abundance ($\Omega_{\chi}h^2$) of the lightest neutralino. Left: The profile likelihood ratio, plotted with $1\sigma$ and $2\sigma$ contour lines drawn in white, and the best fit point indicated by a star. Right: Mechanisms for ensuring that the relic density of DM does not exceed the measured value, through either chargino co-annihilation, res-</figcaption></figure>

onant annihilation via the A/H-funnel, or stop co-annihilation. Other potential mechanisms (e.g. stau co-annihilation) are not shown, as they do not lie within $2\sigma$ of the best-fit point of the entire sample. $2\sigma$ contours for each mechanism are plotted using darker lines, and best-fit points are indicated by a correspondingly coloured star

in the $2\sigma$ region of the full scan that fulfil each of the mass/-composition relations, and the resulting shading of regions is overlaid. In many cases, as we will show, single parameter combinations can satisfy two or more of the mass/composition conditions, and can thus be classified as members of multiple regions. In these cases, one of the mechanisms sometimes dominates over the others. Hybrid sub-regions also exist where the relic density is controlled by two or more mechanisms. For clarity, we make no attempt to show any of these cases as separate regions, nor to colour according to which (if any) mechanism dominates in overlapping regions. For specific cases of interest, we do however attempt to clarify these finer issues in our discussion of the results that we show.

Even within individual regions, readers should be wary of the need for nuance in interpreting the “relic density mechanism” labels. Points labelled “chargino co-annihilation” will typically exhibit co-annihilation of the lightest neutralino with both the lightest chargino and the next-to-lightest neutralino, as small $\tilde{\chi}_1^0-\tilde{\chi}_2^0$ and $\tilde{\chi}_1^0-\tilde{\chi}_1^\pm$ mass splittings are an automatic consequence of a predominantly Higgsino LSP. Nevertheless, both these co-annihilation processes are outweighed in many models simply by boosted $\tilde{\chi}_1^0-\tilde{\chi}_1^0$ annihilation, brought about by the dominance of the Higgsino component in the lightest neutralino. Similarly, A/H-funnel points will exhibit resonant annihilation through both the CP-odd Higgs, $A^0$, and the heavy CP-even Higgs, $H^0$, which are close to degenerate in mass in the CMSSM (and NUHM models). The CP-odd Higgs resonance dominates at the present day however, as s-channel annihilation via the CP-even state is velocity suppressed.

In contrast to previous studies of the CMSSM, we apply the relic density measurement as an upper limit only, allowing for the possibility that thermal neutralinos do not constitute all of DM. This has important consequences for the resulting phenomenology.

Higgsino LSPs are automatically nearly degenerate with the lightest chargino and next-to-lightest neutralino, leading to efficient co-annihilation and an under-abundant relic density for $m_\chi \lesssim 1$ TeV. In isolation, this effect naturally gives the observed relic density at neutralino masses of about a TeV, and lower and higher values at smaller and larger neutralino masses, respectively.<fnref n="8" /> This effect can be seen in the low-mass yellow strip in Fig. 1. If the LSP is instead a “well-tempered” [300] admixture of Higgsino and bino,<fnref n="9" /> then the efficiency of the co-annihilation effect can be tuned to give the exact observed relic density, even at very low neutralino masses. Such scenarios are however heavily constrained by recent LUX [227,228] and Panda-X [229] limits on the spin-independent scattering cross-section [308–310]. As we see in the low-mass section of Fig. 1 however, relaxing the demand that the neutralino must explain all of DM allows models to be more Higgsino-dominated, leading to subdominant neutralino DM. The reduced relic density also helps Higgsino models avoid limits from spin-dependent nuclear scattering, which would otherwise prove rather constraining.

<footnote n="8">⁸ Note that the Sommerfeld effect can be important in the context of pure Higgsino DM; see Sect. 4.4.3 for details.</footnote>

<footnote n="9">⁹ In the CMSSM, this well-tempered mixture is realised within the “focus point” region [301–307].</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Left: The profile likelihood ratio in the CMSSM, for $m_0$ and $m_{1/2}$ (top) and $\tan \beta$ and $A_0$ (bottom), with explicit 68% and 95% CL contour lines drawn in white, and the best fit point indicated by a star. Right: Colour-coding shows the mechanisms active in models within the 95% CL contour for avoiding thermal overproduction of neutralino</figcaption></figure>

dark matter, through either chargino co-annihilation, resonant annihilation via the A/H funnel, or stop co-annihilation. Other potential mechanisms (e.g. stau co-annihilation) are not present, as they do not lie within the 95% CL contour

Similarly, at masses above 1 TeV, the not-quite-efficient-enough Higgsino co-annihilation can be supplemented by additional resonant annihilation through the heavy Higgs funnel, bringing the relic density down to the observed value, or lower. These models can be seen as overlapping yellow and orange regions at $m_\chi \gtrsim 1$ TeV in the right panel of Fig. 1.

We now see that relaxing the relic density constraint to an upper limit opens up a much richer set of phenomenologically-viable scenarios, with lighter Higgsino or mixed Higgino-bino LSPs. From the perspective of global fits, treating the relic density as an upper bound is a conservative approach, and allows us to test whether the preference for heavy spectra found in recent studies [115,145,311] persists even when a greater variety of light LSPs is permitted.

The right panel of Fig. 1 shows that at 95% CL, all of the identified annihilation mechanisms (stop co-annihilation, A/H-funnel and chargino co-annihilation) permit solutions where the measured relic density is fully accounted for, as well as scenarios where only a very small fraction of the DM relic abundance is explained in the CMSSM. The fit does not

demonstrate any clear preference for the relic density to be under-abundant or very close to the measured value. Looking at the top of this plot, we indeed see the established picture for chargino co-annihilation discussed above, where a pure Higgsino DM candidate should have a mass of around 1 TeV to fit the observed relic density.

In Fig. 2, we show 2D CMSSM joint profile likelihoods for $m_0$ and $m_{1/2}$, as well as for $\tan \beta$ and $A_0$. Here the plots include both positive and negative $\mu$, and are again coloured by relic density mechanism. We see a large region of high likelihood at large $m_0$ and $m_{1/2}$, consisting of overlapping chargino co-annihilation and A/H-funnel points. The A/H-funnel region is concentrated at high $\tan \beta$, as is well known from previous studies of the CMSSM (e.g. Ref. [312]). The chargino co-annihilation region disfavours large negative $A_0$, in agreement with existing results in the literature.<fnref n="10" />

<footnote n="10">¹⁰ See for example Fig. 2d of Ref. [311], and the middle panels of Fig. 2 of Ref. [115].</footnote>
At lower $m_0$ and $m_{1/2}$, a stop co-annihilation region appears, with a light stop very close in mass to the lightest neutralino. Due to constraints from direct searches, as well as Higgs-mass measurements at the LHC, which push up the sfermion masses, these scenarios can only be obtained through very large stop mixing. This restricts the stop co-annihilation region to very large and negative $A_0$ values, and low-to-moderate tan $\beta$, as can be seen in the bottom panels of Fig. 2. This region has not been seen in most of the recent global fit literature, as revealing it requires not only consideration of large, negative $A_0$ values, but also very careful scanning of the parameter space.<fnref n="11" />

The preference for large and negative $A_0$ in stop co-annihilation could lead to colour- or charge-breaking minima in the scalar potential. We have investigated the presence of such problems for points in the stop co-annihilation region, using several conditions that have been proposed in the literature:

1. $A_t^2 < 3.0(\mathbf{m}_{Q_{3,3}}^2 + \mathbf{m}_{\mu_{3,3}}^2 + \mu^2 + m_{H_u}^2)$ [313],

2. $A_t^2 < 7.5(\mathbf{m}_{Q_{3,3}}^2 + \mathbf{m}_{\mu_{3,3}}^2) - 3\mu^2$ [314], and

3. $A_t^2 < 3.4(\mathbf{m}_{Q_{3,3}}^2 + \mathbf{m}_{Q_{3,3}}^2) + 60(m_{H_u}^2 + \mu^2)$, based on the results in Ref. [79].

We found that whilst some points in this region do violate one or more of these conditions, removing all points that do so neither modifies the shapes of the likelihood contours in our plots, nor the fact that the best-fit occurs in the stop co-annihilation region. This question could in principle be investigated further by calculating the tunnelling probability for each point, e.g. using Vevacious [315]. However, it is not possible to do this in a reasonable amount of time with the large number of points in our scans. Even though the conditions above are not definitive, being neither necessary nor sufficient to establish that the vacuum of the theory breaks gauge invariance, neither is studying stability with tools such as Vevacious, due to the large number of scalar fields in the MSSM and the resulting difficulty of finding all relevant minima of the potential. We therefore leave detailed investigation of such issues for a future paper.

Looking at the lower-right panel of Fig. 2, the stop co-annihilation region undoubtedly extends to even lower values of $A_0$ than we have considered here. Combined with possible impacts of Sommerfeld enhancement on the relic density [316], this would have the effect of allowing stop co-annihilation to extend to very large values of $m_0$ (Ref. [316] found stop co-annihilation models with $m_0$ as large as 13 TeV). However, as $A_0$ becomes more negative, colour- and charge-breaking vacua become an ever-increasing concern.

In contrast with previous results, we do not find a stau co-annihilation region inside the 95% CL region surrounding our overall best fit. We do find stau co-annihilation solutions in the same region of parameter space as seen in the literature<fnref n="12" /> when we look at $4\sigma$ confidence regions. In addition to this, we also see small islands of stau co-annihilation appear inside the $2\sigma$ contours if we remove the LHC Run II likelihood (leaving only Run I analyses), although these are much smaller than seen in the previous literature. Therefore, the likelihood of the stau co-annihilation strip is suppressed in our results relative to those in the literature. Beyond the LHC Run II likelihood, which has also been shown to impact this region in Ref. [115], the suppression compared to other regions of good fit comes mostly from the LHC Higgs likelihood. This is influenced by the following differences in our analysis compared to existing analyses (for recent examples see Refs. [115,145,311]): (1) relaxation of the relic density constraint to an upper bound, allowing light Higgsino DM scenarios and consequentially relaxing the constraint on $\mu$; (2) differences in the Higgs mass calculation (FlexibleSUSY rather than FeynHiggs) and branching ratio calculations (HDECAY rather than FeynHiggs); (3) a wider prior mass range than some previous scans; and (4) an improved scanning technique, which finds a modestly better fit in the other regions, relative to the stau co-annihilation region.

The plots in Fig. 2 combine scans for $\mu < 0$ and $\mu > 0$. As the sign of $\mu$ is a discrete parameter, it is useful to also investigate each sign independently. These plots are shown in Fig. 3. A preference for positive $\mu$, from the anomalous magnetic moment of the muon $a_\mu$, has been reported previously (e.g. [145]). In our results Higgs observables and LHC sparticle search likelihoods (and the large allowed range for dimensionful parameters) push up the mass scale of the preferred sparticle spectrum, minimising the impact of the $a_\mu$ likelihood and removing the preference for positive $\mu$. We see a mild preference for $\mu < 0$, which has a best fit log-likelihood of -263.75, as compared to -265.00 for $\mu > 0$. The negative $\mu$ results also exhibit an enhanced stop co-annihilation region at low mass, and a reduced $A/H$-funnel region at higher mass, relative to the positive $\mu$ results.

The larger, better-fitting stop co-annihilation region at $\mu < 0$ is driven entirely by the Higgs signal likelihood, in particular the fit to the gauge boson signal strengths. Positive $\mu$ in this region suppresses the $h \to WW, ZZ, \gamma\gamma$ branching fractions to below the observed values, leading to a best-fit likelihood worse than the $\mu < 0$ equivalent by $\Delta \ln \mathcal{L} = 1.3$. Indeed, the $\mu < 0$ fit is actually slightly better than the fit of the SM to the Higgs data, by $\Delta \ln \mathcal{L} = 0.9$ units. Although the implied preference for $\mu < 0$ over $\mu > 0$ is weak, at just $1.1\sigma$ (in 2D), this demonstrates that precision Higgs physics

<footnote n="11">¹¹ As this manuscript was undergoing final editing, an updated version of Ref. [115] was released, showing a stop co-annihilation region in good agreement with ours.</footnote>

<footnote n="12">¹² See for example Fig. 3 of Ref. [311], Fig. 2 of Ref. [115] and Fig. 1 of [145].</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Left: Profile likelihood ratio in the $m_0 - m_{1/2}$ plane of the CMSSM, for $\mu \ge 0$ (top) and $\mu < 0$ (bottom). Right: Colour-coding showing the mechanisms for avoiding a relic density of DM that exceeds the observed value</figcaption></figure>

has now reached the stage where it can directly constrain the parameters of supersymmetric scenarios.

The vastly different size of the heavy Higgs funnel for $\mu < 0$ and $\mu > 0$ is due to differences in LSP composition. For $\mu > 0$, the A/H funnel contains many Higgsino LSP points, which combine with the chargino co-annihilation mechanism to form an extensive hybrid region. In constrast, the bino solutions that do exist in this region are somewhat more concentrated in the $m_0-m_{1/2}$ plane. For $\mu < 0$, the heavy Higgs funnel region is almost exclusively bino, leaving the chargino co-annihilation to exist mostly as a pure mechanism, and resulting in an upper limit on the mass of the LSP of ~1.2 TeV. Although there are some relatively isolated points in the $\mu < 0$ scan exhibiting hybrid funnel-chargino co-annihilation behaviour, it seems difficult to obtain valid solutions to the RGEs with such spectra when $\mu$ is negative.

The best-fit points for each relic density mechanism (with positive and negative $\mu$ results combined) are given in Table 5. These are also shown in the figures of this section, as stars coloured by their corresponding region. We also give the mass spectrum for the global best-fit CMSSM point (column 5 of Table 5) in Fig. 4, demonstrating that the only

light superpartners are the lightest stop, the lightest two neutralinos and the lightest chargino, which is almost exactly degenerate in mass with the $\tilde{\chi}_2^0$. The $\tilde{\chi}_1^0$ is a pure bino for this point, whereas the $\tilde{\chi}_2^0$ and $\tilde{\chi}_1^\pm$ are pure wino. The point generates a relic density within the allowed range through stop co-annihilation, but with a $\tilde{t}_1 - \tilde{\chi}_1^0$ mass difference of $\simeq 40$ GeV. This mass difference should ensure prompt stop decay and potential visibility in future compressed spectrum searches at the LHC.

In Table 5, we also give a detailed breakdown of the likelihood contributions from the different searches discussed in Sect. 3, and compare to an ‘ideal’ reference likelihood. The ideal likelihood is defined as the best likelihood that a model could be expected to achieve, were it to perfectly predict all detections, and make no additional contribution beyond that predicted from background for all other searches. Computing this is straightforward for most likelihood components, as it follows directly from setting the model prediction to either the observed value (e.g. $m_W$, $\Omega_c h^2$, $a_\mu$, any nuisance parameters) or the background-only prediction (e.g. direct DM, LHC and neutrino searches). In some cases however, where multiple sub-observables are involved and the background-
<table><caption>**Table 5** Best-fit points in the CMSSM, for each of the regions characterised by a specific mechanism for suppressing the relic density of dark matter. Here we show the likelihood contributions, parameter values at each point, and some quantities relevant for the interpretation of mass spectra at the different best fits. We also give likelihood components</caption><thead><tr><th>Likelihood term</th><th>Ideal</th><th>A/H-funnel</th><th>τ̃ co-ann.</th><th>t̃ co-ann.</th><th>χ̃<sub>1</sub><sup>±</sup> co-ann.</th><th>Δ ln L<sub>BF</sub></th></tr></thead><tbody><tr><td>LHC sparticle searches</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>LHC Higgs</td><td>-37.734</td><td>-37.960</td><td>-41.296</td><td>-38.042</td><td>-38.069</td><td>0.308</td></tr><tr><td>LEP Higgs</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH selectron</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH smuon</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH stau</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 selectron</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 smuon</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 stau</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 neutralino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 chargino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino hadronic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino semi-leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL neutralino hadronic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>B<sub>(s)</sub> → μ<sup>+</sup>μ<sup>-</sup></td><td>0.000</td><td>-1.939</td><td>-2.739</td><td>-2.029</td><td>-1.939</td><td>2.029</td></tr><tr><td>Tree-level B and D decays</td><td>0.000</td><td>-15.515</td><td>-15.491</td><td>-15.283</td><td>-15.610</td><td>15.283</td></tr><tr><td>B<sup>0</sup> → K*0μ<sup>+</sup>μ<sup>-</sup></td><td>-184.260</td><td>-196.506</td><td>-197.469</td><td>-196.088</td><td>-196.309</td><td>11.828</td></tr><tr><td>B → X<sub>s</sub>γ</td><td>9.799</td><td>9.258</td><td>9.525</td><td>9.106</td><td>9.184</td><td>0.693</td></tr><tr><td>a<sub>μ</sub></td><td>20.266</td><td>13.915</td><td>14.556</td><td>13.977</td><td>13.903</td><td>6.289</td></tr><tr><td>W mass</td><td>3.281</td><td>3.084</td><td>3.093</td><td>3.050</td><td>3.095</td><td>0.231</td></tr><tr><td>Relic density</td><td>5.989</td><td>5.989</td><td>5.984</td><td>5.989</td><td>5.989</td><td>0.000</td></tr><tr><td>PICO-2L</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>0.000</td></tr><tr><td>PICO-60 F</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>SIMPLE 2014</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>0.000</td></tr><tr><td>LUX 2015</td><td>-0.640</td><td>-0.676</td><td>-0.642</td><td>-0.640</td><td>-0.727</td><td>0.000</td></tr><tr><td>LUX 2016</td><td>-1.467</td><td>-1.539</td><td>-1.472</td><td>-1.467</td><td>-1.646</td><td>0.000</td></tr><tr><td>PandaX 2016</td><td>-1.886</td><td>-1.936</td><td>-1.889</td><td>-1.886</td><td>-2.009</td><td>0.000</td></tr><tr><td>SuperCDMS 2014</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>0.000</td></tr><tr><td>XENON100 2012</td><td>-1.693</td><td>-1.675</td><td>-1.692</td><td>-1.693</td><td>-1.651</td><td>0.000</td></tr><tr><td>IceCube 79-string</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>γ rays (Fermi-LAT dwarfs)</td><td>-33.244</td><td>-33.421</td><td>-33.393</td><td>-33.381</td><td>-33.394</td><td>0.137</td></tr><tr><td>ρ<sub>0</sub></td><td>1.142</td><td>1.141</td><td>1.142</td><td>1.141</td><td>1.141</td><td>0.001</td></tr><tr><td>σ<sub>s</sub> and σ<sub>t</sub></td><td>-6.115</td><td>-6.115</td><td>-6.116</td><td>-6.115</td><td>-6.117</td><td>0.000</td></tr><tr><td>α<sub>s</sub>(m<sub>Z</sub>)(M̄S)</td><td>6.500</td><td>6.487</td><td>6.479</td><td>6.481</td><td>6.479</td><td>0.019</td></tr><tr><td>Top quark mass</td><td>-0.645</td><td>-0.645</td><td>-0.645</td><td>-0.649</td><td>-0.645</td><td>0.004</td></tr><tr><td>Total</td><td>-226.927</td><td>-264.273</td><td>-268.287</td><td>-263.747</td><td>-264.546</td><td>36.820</td></tr></tbody></table><table><thead><tr><th>Quantity</th><th>A/H-funnel</th><th>τ̃ co-ann.</th><th>t̃ co-ann.</th><th>χ̃<sub>1</sub><sup>±</sup> co-ann.</th></tr></thead><tbody><tr><td>A<sub>0</sub></td><td>9924.435</td><td>-1227.154</td><td>-9965.036</td><td>9206.079</td></tr><tr><td>m<sub>0</sub></td><td>9136.379</td><td>1476.893</td><td>4269.402</td><td>9000.628</td></tr><tr><td>m<sub>1/2</sub></td><td>2532.163</td><td>2422.340</td><td>1266.043</td><td>2256.472</td></tr><tr><td>tan β</td><td>49.048</td><td>48.594</td><td>14.857</td><td>49.879</td></tr><tr><td>sgn(μ)</td><td>-</td><td>+</td><td>-</td><td>-</td></tr></tbody></table>

for a canonical 'ideal' likelihood (see text), along with its offset from the global best-fit point. SLHA1 and SLHA2 files corresponding to the best-fit point in each region can be found in the online data associated with this paper [163]
<table><caption>**Table 5** continued</caption><thead><tr><th>Quantity</th><th>A/H-funnel</th><th>τ̃ co-ann.</th><th>τ̃ co-ann.</th><th>χ̃<sub>1</sub><sup>±</sup> co-ann.</th></tr></thead><tbody><tr><td>m<sub>t</sub></td><td>173.366</td><td>173.358</td><td>173.267</td><td>173.329</td></tr><tr><td>α<sub>s</sub>(m<sub>Z</sub>)(M&#773;S)</td><td>0.119</td><td>0.119</td><td>0.119</td><td>0.119</td></tr><tr><td>ρ<sub>0</sub></td><td>0.394</td><td>0.401</td><td>0.403</td><td>0.394</td></tr><tr><td>σ<sub>s</sub></td><td>42.950</td><td>43.031</td><td>42.975</td><td>43.503</td></tr><tr><td>σ<sub>l</sub></td><td>57.976</td><td>58.544</td><td>57.887</td><td>58.155</td></tr><tr><td>M<sub>1</sub></td><td>1140.417</td><td>1089.994</td><td>556.554</td><td>1011.999</td></tr><tr><td>μ</td><td>-1409.433</td><td>2621.118</td><td>-4073.398</td><td>-983.112</td></tr><tr><td>m<sub>τ̃<sub>1</sub></sub></td><td>6554.967</td><td>3594.650</td><td>592.052</td><td>6279.661</td></tr><tr><td>m<sub>τ̃<sub>2</sub></sub></td><td>6590.901</td><td>1076.748</td><td>4071.458</td><td>6407.136</td></tr><tr><td>m<sub>A</sub></td><td>2292.366</td><td>2182.200</td><td>5612.268</td><td>1953.735</td></tr><tr><td>m<sub>h</sub></td><td>124.896</td><td>124.054</td><td>125.007</td><td>124.797</td></tr><tr><td>m<sub>χ̃<sub>1</sub><sup>0</sup></sub></td><td>1133.191</td><td>1076.738</td><td>565.069</td><td>973.418</td></tr><tr><td>(%bino, %Higgsino)</td><td>(99, 1)</td><td>(100, 0)</td><td>(100, 0)</td><td>(44, 56)</td></tr><tr><td>m<sub>χ̃<sub>2</sub><sup>0</sup></sub></td><td>1432.774</td><td>1999.921</td><td>1083.062</td><td>-1005.489</td></tr><tr><td>(%bino, %Higgsino)</td><td>(1, 98)</td><td>(0, 1)</td><td>(0, 0)</td><td>(0, 100)</td></tr><tr><td>m<sub>χ̃<sub>1</sub><sup>±</sup></sub></td><td>1430.811</td><td>2000.084</td><td>1083.224</td><td>1002.018</td></tr><tr><td>(%wino, %Higgsino)</td><td>(1, 99)</td><td>(99, 1)</td><td>(100, 0)</td><td>(1, 99)</td></tr><tr><td>m<sub>g̃</sub></td><td>5545.587</td><td>5017.077</td><td>2926.857</td><td>5002.109</td></tr><tr><td>Ωh<sup>2</sup></td><td>6.88 × 10<sup>-2</sup></td><td>1.06 × 10<sup>-1</sup></td><td>4.62 × 10<sup>-2</sup></td><td>4.00 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Particle mass spectrum of the CMSSM best-fit point</figcaption></figure>

only or SM prediction can be improved on by including a BSM contribution, a more nuanced calculation is required. This is the case for the LHC Higgs and electroweak penguin ($B^0 \to K^*0\mu^+\mu^-$) likelihoods. For these components, we define the ideal likelihood to be the highest value possible in a more general phenomenological scenario. In the flavour sector, we use the $B^0 \to K^*0\mu^+\mu^-$ likelihood at the best fit point of the scan of the flavour EFT shown in Ref. [158]. In the Higgs sector, we take the best-fit likelihood obtainable by allowing the mass, width and decay branching fractions of a single scalar to vary freely in order to fit the full set of data contained in HiggsSignals.

The log-likelihood difference of the best fit in the CMSSM to the ideal likelihood is $\Delta\mathcal{L}_{BF} = 36.820$. This difference is largely driven by known anomalies that cannot be explained by either the SM or the MSSM, including the magnetic moment of the muon ($\Delta\mathcal{L} = 6.289$; see Ref. [159]) and the $B^0 \to K^*0\mu^+\mu^-$ angular observables ($\Delta\mathcal{L} = 11.828$; see Ref. [158]). The largest contribution to $\Delta\mathcal{L}_{BF}$ comes from anomalies in tree-level B and D decays, in particular the partially-correlated branching fraction ratios $R_D$ and $R_{D*}$. These have values of 0.308 and 0.248 respectively at the best-fit point of the scan. For comparison, the SM predictions are 0.300 and 0.252, whereas the observed values are $0.403 \pm 0.047$ and $0.310 \pm 0.017$. Further discussion and details can be found in the FlavBit paper [158].

The log-likelihood difference of a point relative to the ideal log-likelihood can be used to give some indication of the goodness of fit, as its definition is very similar to half of the “likelihood $\chi^2$” of Baker and Cousins [317]. The likelihood $\chi^2$ is known to follow a $\chi^2$ distribution in the asymptotic limit. The main difficulty in using this fact is estimating the effective degrees of freedom of the fit, as carrying out simulations to find the true distribution of our test statistic is computationally intractable.<fnref n="13" /> Given the number of observ-

<footnote n="13">¹³ We note that a similar thing has been done in the CMSSM [311], but using a likelihood function far quicker to compute than ours, based on interpolation in a 2D grid of LHC signal yields rather than explicit simulation for each parameter combination.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Left: Profile likelihood ratio in the planes of the NUHM1 parameters $m_0$ and $m_{1/2}$ (top), and tan $\beta$ and $A_0$ (bottom). Explicit contour lines for 68% and 95% CL are drawn in white and the best fit point is</figcaption></figure>

indicated with a star. Right: Colour-coding shows the mechanisms to avoid exceeding the observed relic density of DM

ables that actively constrain the fit, a reasonable guess for the effective degrees of freedom is probably something in the range of 30–50, leading to a $p$ value of between $2 \times 10^{-5}$ and 0.02 for the CMSSM; neither a particularly good fit nor catastrophically bad, given the uncertainties involved in the estimate of the $p$ value. Taking 40 as a canonical estimate of the degrees of freedom, for the sake of later comparison with the NUHM1 and NUHM2, the $p$ value would be $9.4 \times 10^{-4}$.

## 4.2 NUHM1

The main results from the NUHM1 scan are shown in Figs. 5 and 6. Figure 5 shows results in the $m_0–m_{1/2}$ and tan $\beta–A_0$ planes, with plots of the profile likelihood ratio on the left and the DM annihilation mechanisms, defined as in the previous subsection, on the right. In comparison to the CMSSM equivalent, Fig. 2, one can see that the additional freedom in the NUHM1 substantially extends the likelihood contours, so that much of the parameter space is now allowed.

In particular, we now find a stau co-annihilation region, which was absent in the CMSSM results. The extra freedom present in the Higgs sector in the NUHM1 avoids the penalty

from the LHC Higgs likelihood seen in stau co-annihilation models in the CMSSM. Furthermore, we now see chargino co-annihilation solutions within the $2\sigma$ contours that extend to arbitrarily low $m_0$. These are both a consequence of the fact that once $m_0$ is decoupled from $m_H$, the former can be pushed low without impacting EWSB. This allows light staus to exist, making stau co-annihilation viable, and also means that $|\mu|$ can be low at arbitrarily small $m_0$, leading to Higgsino LSPs.

Such low values of $m_0$ in chargino-coannihilation scenarios suggests that the first- and second-generation squarks may be light enough to be constrained directly by collider searches. However, a detailed examination reveals that their masses remain above 2 TeV, and out of reach of LHC limits, for all models within our $2\sigma$ contours.

A similar expansion of the chargino co-annihilation region<fnref n="14" /> has been seen in the previous literature comparing the CMSSM and NUHM1 models (see e.g. Fig. 6 of Ref.

<footnote n="14">¹⁴ The stau co-annihilation region also extends to arbitrarily low $m_0$, but this is because our definition of stau co-annihilation admits the possibility of an under-abundant Higgsino DM candidate with hybrid stau and chargino co-annihilation.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 6 As per Fig. 5, but for the $m_{1/2}-m_H$ (top) and $m_0-m_H$ (bottom) planes</figcaption></figure>

[129] and Fig. 1 of Ref. [145], although the contours do not reach arbitrarily low $m_0$ for all $m_{1/2}$ in those studies. This difference can be explained by the additional freedom associated with only applying the relic density measurement as a one-sided limit. We checked that demanding neutralinos make up all of DM removes some low-$m_0$ scenarios from the $2\sigma$ contours, such that a better agreement with the results in the literature is obtained.

Another interesting feature of this is that these low $m_0$ values mean that the NUHM1 admits significantly lighter squarks within the chargino co-annihilation region than in the CMSSM.

We give results in the $m_H-m_{1/2}$ and $m_H-m_0$ planes for the NUHM1 in Fig. 6. These plots show a sharp cut-off in the likelihood near the diagonals $m_0 \approx m_H$ and $m_{1/2} \approx m_H$, such that $m_H$ must be greater than both $m_0$ and $m_{1/2}$. This structure emerges from the combination of several effects. Reducing $m_H$ with respect to $m_0$ and $m_{1/2}$ leads to $m_{H_u}^2$ running to a more negative value, and this in turn leads to a larger $\mu$ value through the EWSB conditions. Past this boundary in Fig. 6, $\mu$ is then always significantly larger than the bino mass, $M_1$. As a result the neutralino LSP is always bino in these scenarios and requires either an A/H-funnel

or sfermion co-annihilation mechanism to reduce the relic density to the measured value or below. The A/H-funnel mechanism also requires $\mu$ to be small, as $\mu^2$ gives a contribution to the pseudoscalar mass. As a result, if $\mu$ is much larger than $2M_1$, then the relation for the A/H-funnel mechanism, $m_A \approx 2m_{\tilde{\chi}_1^0} \approx 2M_1$, cannot be achieved. Although we do find majority-bino LSPs annihilating through an A/H-funnel, these have smaller values of $\mu$ than can be achieved when $m_H$ is less than either $m_0$ or $m_{1/2}$. Finally, sfermion co-annihilation can be effective in this region, but only for lower values of $m_0$ and $m_{1/2}$. In those scenarios, if $m_H \lesssim m_0, m_{1/2}$ the likelihood is suppressed by the LHC Higgs likelihood, because it is difficult to fit the 125 GeV Higgs there.

We investigated charge- and colour-breaking minima in the NUHM1 in the same way as in the CMSSM (Sect. 4.1). As in the CMSSM, a number of models within our 95% CL regions are affected by one or more of the three proposed conditions, but removing all such parameter combinations does not move the best fit to a different region, nor substantially change the regions of parameter space preferred by our fits.

As we did for the CMSSM, in Table 6 we show the best-fit points for each mechanism for depleting the relic density of
<table><caption>**Table 6** Best-fit points in the NUHM1, for each of the regions characterised by a specific mechanism for suppressing the relic density of dark matter. Here we show the likelihood contributions, parameter values at each point, and some quantities relevant for the interpretation of mass spectra at the different best fits. We also give likelihood components</caption><thead><tr><th>Likelihood term</th><th>Ideal</th><th>A/H-funnel</th><th>τ̃ co-ann.</th><th>τ̃ co-ann.</th><th>χ̃<sub>1</sub><sup>±</sup> co-ann.</th><th>Δ ln L<sub>BF</sub></th></tr></thead><tbody><tr><td>LHC sparticle searches</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>LHC Higgs</td><td>-37.734</td><td>-38.646</td><td>-38.182</td><td>-38.271</td><td>-38.531</td><td>0.537</td></tr><tr><td>LEP Higgs</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH selectron</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH smuon</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH stau</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 selectron</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 smuon</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 stau</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 neutralino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 chargino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino hadronic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino semi-leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL neutralino hadronic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>B<sub>(s)</sub> → μ<sup>+</sup>μ<sup>-</sup></td><td>0.000</td><td>-1.985</td><td>-2.033</td><td>-2.032</td><td>-2.043</td><td>2.032</td></tr><tr><td>Tree-level B and D decays</td><td>0.000</td><td>-15.703</td><td>-15.286</td><td>-15.286</td><td>-15.282</td><td>15.286</td></tr><tr><td>B<sup>0</sup> → K*<sup>0</sup>μ<sup>+</sup>μ<sup>-</sup></td><td>-184.260</td><td>-196.553</td><td>-195.323</td><td>-194.855</td><td>-194.825</td><td>10.595</td></tr><tr><td>B → X<sub>S</sub>γ</td><td>9.799</td><td>9.272</td><td>8.696</td><td>8.430</td><td>8.351</td><td>1.369</td></tr><tr><td>a<sub>μ</sub></td><td>20.266</td><td>14.158</td><td>13.837</td><td>13.819</td><td>13.836</td><td>6.447</td></tr><tr><td>W mass</td><td>3.281</td><td>3.095</td><td>3.062</td><td>3.075</td><td>3.096</td><td>0.206</td></tr><tr><td>Relic density</td><td>5.989</td><td>5.989</td><td>5.989</td><td>5.989</td><td>5.989</td><td>0.000</td></tr><tr><td>PICO-2L</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>0.000</td></tr><tr><td>PICO-60 F</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>-0.001</td><td>0.000</td></tr><tr><td>SIMPLE 2014</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>0.000</td></tr><tr><td>LUX 2015</td><td>-0.640</td><td>-0.666</td><td>-0.646</td><td>-0.659</td><td>-0.676</td><td>0.019</td></tr><tr><td>LUX 2016</td><td>-1.467</td><td>-1.519</td><td>-1.479</td><td>-1.504</td><td>-1.539</td><td>0.037</td></tr><tr><td>PandaX 2016</td><td>-1.886</td><td>-1.921</td><td>-1.894</td><td>-1.912</td><td>-1.936</td><td>0.026</td></tr><tr><td>SuperCDMS 2014</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>0.000</td></tr><tr><td>XENON100 2012</td><td>-1.693</td><td>-1.680</td><td>-1.690</td><td>-1.684</td><td>-1.675</td><td>0.009</td></tr><tr><td>IceCube 79-string</td><td>0.000</td><td>-0.014</td><td>0.000</td><td>0.000</td><td>-0.135</td><td>0.000</td></tr><tr><td>γ rays (Fermi-LAT dwarfs)</td><td>-33.244</td><td>-33.384</td><td>-33.364</td><td>-33.373</td><td>-33.398</td><td>0.129</td></tr><tr><td>ρ<sub>0</sub></td><td>1.142</td><td>1.141</td><td>1.141</td><td>1.140</td><td>1.141</td><td>0.002</td></tr><tr><td>σ<sub>s</sub> and σ<sub>l</sub></td><td>-6.115</td><td>-6.115</td><td>-6.135</td><td>-6.124</td><td>-6.117</td><td>0.009</td></tr><tr><td>α<sub>s</sub>(m<sub>Z</sub>)(M̄S)</td><td>6.500</td><td>6.491</td><td>6.488</td><td>6.493</td><td>6.494</td><td>0.007</td></tr><tr><td>Top quark mass</td><td>-0.645</td><td>-0.647</td><td>-0.673</td><td>-0.655</td><td>-0.645</td><td>0.010</td></tr><tr><td>Total</td><td>-226.927</td><td>-264.907</td><td>-263.712</td><td>-263.629</td><td>-264.115</td><td>36.702</td></tr></tbody></table><table><thead><tr><th>Quantity</th><th>A/H-funnel</th><th>τ̃ co-ann.</th><th>τ̃ co-ann.</th><th>χ̃<sub>1</sub><sup>±</sup> co-ann.</th></tr></thead><tbody><tr><td>A<sub>0</sub></td><td>9084.348</td><td>-7798.283</td><td>-7016.861</td><td>-6439.114</td></tr><tr><td>m<sub>0</sub></td><td>5139.563</td><td>1659.858</td><td>2042.775</td><td>1472.445</td></tr><tr><td>m<sub>1/2</sub></td><td>5266.693</td><td>2656.510</td><td>2245.476</td><td>2319.968</td></tr><tr><td>m<sub>H</sub></td><td>6954.864</td><td>5407.626</td><td>4990.078</td><td>5034.071</td></tr><tr><td>tan β</td><td>53.263</td><td>19.430</td><td>18.128</td><td>11.840</td></tr></tbody></table>

for a canonical 'ideal' likelihood (see text), along with its offset from the global best-fit point. SLHA1 and SLHA2 files corresponding to the best-fit point in each region can be found in the online data associated with this paper [163]
<table><caption>**Table 6** continued</caption><thead><tr><th>Quantity</th><th>A/H-funnel</th><th>τ&#773; co-ann.</th><th>t&#773; co-ann.</th><th>χ&#773;<sub>1</sub><sup>±</sup> co-ann.</th></tr></thead><tbody><tr><td>sgn(μ)</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>m<sub>t</sub></td><td>173.393</td><td>173.522</td><td>173.451</td><td>173.362</td></tr><tr><td>α<sub>s</sub>(m<sub>Z</sub>)(M&#773;S)</td><td>0.119</td><td>0.118</td><td>0.119</td><td>0.119</td></tr><tr><td>ρ<sub>0</sub></td><td>0.403</td><td>0.398</td><td>0.408</td><td>0.396</td></tr><tr><td>σ<sub>s</sub></td><td>42.776</td><td>43.646</td><td>43.747</td><td>42.478</td></tr><tr><td>σ<sub>t</sub></td><td>57.737</td><td>56.355</td><td>57.132</td><td>58.024</td></tr><tr><td>M<sub>1</sub></td><td>2419.401</td><td>1184.390</td><td>994.971</td><td>1023.177</td></tr><tr><td>μ</td><td>836.283</td><td>-1753.895</td><td>-1462.491</td><td>-351.100</td></tr><tr><td>m<sub>t&#773;<sub>1</sub></sub></td><td>7902.945</td><td>1198.127</td><td>1032.608</td><td>1012.967</td></tr><tr><td>m<sub>t&#773;<sub>1</sub></sub></td><td>2231.113</td><td>1295.803</td><td>1819.486</td><td>1513.479</td></tr><tr><td>m<sub>A</sub></td><td>1805.767</td><td>5428.634</td><td>5002.455</td><td>5122.233</td></tr><tr><td>m<sub>h</sub></td><td>125.026</td><td>124.544</td><td>124.531</td><td>124.903</td></tr><tr><td>m<sub>χ&#773;<sub>1</sub><sup>0</sup></sub></td><td>856.207</td><td>1179.991</td><td>993.716</td><td>358.905</td></tr><tr><td>(%bino, %Higgsino)</td><td>(0, 100)</td><td>(100, 0)</td><td>(100, 0)</td><td>(0, 100)</td></tr><tr><td>m<sub>χ&#773;<sub>2</sub><sup>0</sup></sub></td><td>-858.645</td><td>1760.580</td><td>1467.989</td><td>-364.815</td></tr><tr><td>(%bino, %Higgsino)</td><td>(0, 100)</td><td>(0, 98)</td><td>(0, 98)</td><td>(0, 100)</td></tr><tr><td>m<sub>χ&#773;<sub>1</sub><sup>±</sup></sub></td><td>857.791</td><td>1760.608</td><td>1467.887</td><td>362.366</td></tr><tr><td>(%wino, %Higgsino)</td><td>(0, 100)</td><td>(2, 98)</td><td>(2, 98)</td><td>(0, 100)</td></tr><tr><td>m<sub>g&#773;</sub></td><td>10470.041</td><td>5462.593</td><td>4705.842</td><td>4823.285</td></tr><tr><td>Ωh<sup>2</sup></td><td>7.03 × 10<sup>-2</sup></td><td>5.24 × 10<sup>-2</sup></td><td>9.29 × 10<sup>-2</sup></td><td>1.59 × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Particle mass spectrum of the NUHM1 best-fit point</figcaption></figure>

DM. The best-fit point in the chargino co-annihilation region has small $\tilde{\chi}_1^0$, $\tilde{\chi}_2^0$ and $\tilde{\chi}_1^\pm$ masses, but escapes LHC exclusion due to the highly compressed mass spectrum for these sparticles. As in the CMSSM, the overall best-fit point lies in the stop co-annihilation region. Its mass spectrum is shown in Fig. 7. There are important differences to the CMSSM case, however. Firstly, the stop is heavier, now sitting just above 1 TeV in mass. The $\tilde{t}_1 - \tilde{\chi}_1^0$ mass difference is once again roughly 40 GeV, ensuring prompt decay of the stop. The heavier stop mass is accompanied by a heavier mass spectrum in general, with no sparticles lighter than 800 GeV in mass. The $\tilde{\chi}_1^0$ is pure bino, but the $\tilde{\chi}_2^0$, $\tilde{\chi}_3^0$ and $\tilde{\chi}_1^\pm$ are now

predominantly Higgsino in character, leaving the $\tilde{\chi}_4^0$ and $\tilde{\chi}_2^\pm$ to be mostly wino. Discovery of this point would be very challenging at the LHC in the near future, due to the heavy weakly-coupled states, and the lack of light coloured states that have a large mass splitting with the $\tilde{\chi}_1^0$.

For the NUHM1, $\Delta \ln \mathcal{L}_{\text{BF}} = 36.702$, slightly better than what we found in the CMSSM. For the sake of comparison with the CMSSM ($p = 9.4 \times 10^{-4}$ if computed with 40 degrees of freedom), we can compute a p-value assuming one less degree of freedom, i.e. 39. This gives $7.1 \times 10^{-4}$, slightly worse than the CMSSM. We see that despite the improvement in the fit, the fact that it has not delivered a sufficiently large improvement in $\Delta \ln \mathcal{L}_{\text{BF}}$ means that this is not enough to outweigh the penalty associated with the introduction of the additional parameter.

## 4.3 NUHM2

The NUHM2 results are shown in Figs. 8 and 9. Figure 8 shows results in the $m_0-m_{1/2}$ and $\tan \beta-A_0$ planes, with plots of the profile likelihood ratio on the left and the annihilation mechanism, defined at the start of Sect. 4.1, on the right.

In comparison to the NUHM1 results, the stau co-annihilation region is significantly extended, covering higher values of $m_{1/2}$, and lower $A_0$ and $\tan \beta$. This is due to modification of the RG flow for the soft scalar stau masses that
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Left: Profile likelihoods in the NUHM2, in terms of the $m_0 - m_{1/2}$ and $A_0 - \tan \beta$ planes. Right: corresponding mechanisms to avoid exceeding the observed relic density of DM</figcaption></figure>

occurs when the soft Higgs masses are split at the GUT scale. A similar effect has been observed and discussed for this model in Ref. [148]. As in the NUHM1 despite the low values of $m_0$ our models within the $2\sigma$ contours do not generate first and second generation squark masses below ~2 TeV.

In Fig. 9 we show the structure of the $m_{H_u}-m_{1/2}$ and $m_{H_d}-m_{1/2}$ planes in the same format as the $m_0-m_{1/2}$ and $\tan \beta-A_0$ planes. The $m_{H_u}-m_{1/2}$ plot is quite similar to the $m_H-m_{1/2}$ plot (Fig. 6) for the NUHM1 model discussed in Sect. 4.2, while in contrast there is not much structure in the $m_{H_d}-m_{1/2}$ plane.<fnref n="15" />

This could be anticipated from the NUHM1 results (again Fig. 6), as the structure was caused by the fact that smaller $m_{H_u}$ at the GUT scale leads to $\mu \gg M_1$ at the SUSY scale, making bino DM the only possibility. As in the NUHM1 case, the A/H-funnel mechanism for the bino again does not work because $\mu$ is too large to allow $m_A \approx 2m_{\tilde{\chi}_1^0} \approx 2M_1$. However, the extra freedom in the Higgs sector from splitting $m_{H_u}$ and $m_{H_d}$ at the GUT scale does allow a better LHC Higgs

likelihood at smaller $m_0$ and $m_{1/2}$, so that the stau and stop co-annihilation regions can be found when $m_{H_u} \lesssim m_0, m_{1/2}$.

As with the NUHM1, we checked the three charge- and colour-breaking conditions mentioned in Sect. 4.1. The results were as in the NUHM1: some individual parameter combinations are affected, but the overall inference is not.

We give a table of best-fit NUHM2 points in Table 7, with the mass spectrum for the overall best fit shown in Fig. 10. Once again, the overall best fit is obtained for a point that satisfies the relic density bound through stop co-annihilation. The $\tilde{\tau}_1$ mass is 950 GeV, but the mass difference with the $\tilde{\chi}_1^0$ is now less than 20 GeV, making this very difficult to resolve at the LHC. The neutralino-chargino sector features a pure bino $\tilde{\chi}_1^0$, wino-dominated $\tilde{\chi}_2^0$ and $\tilde{\chi}_1^{\pm}$, and Higgsino-dominated $\tilde{\chi}_3^0$, $\tilde{\chi}_4^0$ and $\tilde{\chi}_2^{\pm}$. The large mass of the $\tilde{\chi}_2^0$ adds to the problem of the $\tilde{\tau}_1 - \tilde{\chi}_1^0$ mass difference, making this a particularly challenging scenario for collider searches.

For the NUHM2, $\Delta \ln \mathcal{L} = 36.362$, indicating a better fit than either the CMSSM or NUHM1. This is expected to some extent, as the NUHM2 has one more free parameter than the NUHM1. Indeed, accounting for the extra freedom in the fit (i.e. adopting a canonical degree of freedom of 38 instead

<footnote n="15">¹⁵ For brevity we omit plots showing the $m_{H_u}-m_0$ and $m_{H_d}-m_0$ planes, which exhibit the same behaviour.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 9** As per Fig. 8, but for the $m_{H_u}-m_{1/2}$ (top) and $m_{H_d}-m_{1/2}$ (bottom) planes</figcaption></figure>

of 39 or 40 – see previous subsections), and computing the implied *p*-value, the result is just $5.9 \times 10^{-4}$. This actually disfavours the NUHM2 compared to the NUHM1 ($p = 7.1 \times 10^{-4}$ for 39 dof) and CMSSM ($p = 9.4 \times 10^{-4}$ for 40 dof), because its additional parameter does not provide a sufficiently large improvement to the overall fit.

We have not commented so far on the ability of any of the models to explain the large discrepancy between the measured value of $a_\mu \equiv (g-2)_\mu/2$ and that predicted by the SM [237,318]. This is because, with the heavy spectra found, a sizable supersymmetric contribution to $\Delta a_\mu$ is not expected. However, in Ref. [148] (see right panel of Fig. 12 in that paper), it was found that although the best fit for the NUHM2 predicts a very small $a_\mu$, there are points within the $2\sigma$ contours that predict significantly larger values of around $2 \times 10^{-9}$, which may give some grounds for optimism. In contrast, the MSSM contribution to $a_\mu$ within the $2\sigma$ confidence regions of our CMSSM, NUHM1 and NUHM2 fits is below $5 \times 10^{-10}$. We show this visually in Fig. 11. Therefore, with the latest data and using GM2Calc to obtain the most precise calculation available of the supersymmetric contributions to the anomalous magnetic moment of the muon, we

find that none of the GUT-scale models that we consider can make a significant contribution to resolving this discrepancy.

A similar point can be made about the flavour anomalies associated with the angular observables in $B^0 \to K^{*0}\mu^+\mu^-$ decays and the ratios $R_D$ and $R_{D*}$; none of the best-fits or 95% CL regions of our scans indicate any ability for these data to be explained within the CMSSM, NUHM1 or NUHM2.

## 4.4 Discovery prospects

In the following we discuss the discovery prospects of the CMSSM, NUHM1 and NUHM2, given current constraints. We first address prospects at the LHC, followed by direct and indirect detection of DM.

### 4.4.1 LHC

In Fig. 12, we show the 1D profile likelihood ratio for the masses of the gluino, lightest (third generation) squarks, lightest stau, lightest chargino and lightest neutralino in the CMSSM, NUHM1 and NUHM2. The likelihood is generally low for coloured sparticles light enough to be in reach of LHC
<table><caption>**Table 7** Best-fit points in the NUHM2, for each of the regions characterised by a specific mechanism for suppressing the relic density of dark matter. Here we show the likelihood contributions, parameter values at each point, and some quantities relevant for the interpretation of mass spectra at the different best fits. We also give likelihood components</caption><thead><tr><th>Likelihood term</th><th>Ideal</th><th>A/H-funnel</th><th>τ̃ co-ann.</th><th>τ̃ co-ann.</th><th>χ̃<sub>1</sub><sup>±</sup> co-ann.</th><th>Δ ln L<sub>BF</sub></th></tr></thead><tbody><tr><td>LHC sparticle searches</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>LHC Higgs</td><td>-37.734</td><td>-38.563</td><td>-37.928</td><td>-37.980</td><td>-38.484</td><td>0.246</td></tr><tr><td>LEP Higgs</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH selectron</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH smuon</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>ALEPH stau</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 selectron</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 smuon</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 stau</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 neutralino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>L3 chargino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino hadronic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino semi-leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL chargino leptonic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>OPAL neutralino hadronic</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>B<sub>(s)</sub> → μ<sup>+</sup>μ<sup>-</sup></td><td>0.000</td><td>-1.972</td><td>-2.037</td><td>-2.033</td><td>-2.030</td><td>2.033</td></tr><tr><td>Tree-level B and D decays</td><td>0.000</td><td>-15.553</td><td>-15.283</td><td>-15.283</td><td>-15.290</td><td>15.283</td></tr><tr><td>B<sup>0</sup> → K*<sup>0</sup>μ<sup>+</sup>μ<sup>-</sup></td><td>-184.260</td><td>-195.596</td><td>-195.475</td><td>-195.043</td><td>-194.415</td><td>10.783</td></tr><tr><td>B → X<sub>S</sub>γ</td><td>9.799</td><td>8.865</td><td>8.797</td><td>8.550</td><td>8.077</td><td>1.249</td></tr><tr><td>a<sub>μ</sub></td><td>20.266</td><td>14.086</td><td>13.756</td><td>13.842</td><td>13.876</td><td>6.424</td></tr><tr><td>W mass</td><td>3.281</td><td>3.060</td><td>3.078</td><td>3.074</td><td>3.097</td><td>0.207</td></tr><tr><td>Relic density</td><td>5.989</td><td>5.989</td><td>5.989</td><td>5.989</td><td>5.989</td><td>0.000</td></tr><tr><td>PICO-2L</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>-1.000</td><td>0.000</td></tr><tr><td>PICO-60 F</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>SIMPLE 2014</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>-2.972</td><td>0.000</td></tr><tr><td>LUX 2015</td><td>-0.640</td><td>-0.657</td><td>-0.641</td><td>-0.641</td><td>-0.671</td><td>0.001</td></tr><tr><td>LUX 2016</td><td>-1.467</td><td>-1.501</td><td>-1.468</td><td>-1.470</td><td>-1.529</td><td>0.003</td></tr><tr><td>PandaX 2016</td><td>-1.886</td><td>-1.909</td><td>-1.887</td><td>-1.888</td><td>-1.929</td><td>0.002</td></tr><tr><td>SuperCDMS 2014</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>-2.248</td><td>0.000</td></tr><tr><td>XENON100 2012</td><td>-1.693</td><td>-1.685</td><td>-1.693</td><td>-1.692</td><td>-1.678</td><td>0.001</td></tr><tr><td>IceCube 79-string</td><td>0.000</td><td>-0.021</td><td>0.000</td><td>0.000</td><td>-0.108</td><td>0.000</td></tr><tr><td>γ rays (Fermi-LAT dwarfs)</td><td>-33.244</td><td>-33.398</td><td>-33.371</td><td>-33.369</td><td>-33.398</td><td>0.125</td></tr><tr><td>ρ<sub>0</sub></td><td>1.142</td><td>1.141</td><td>1.137</td><td>1.141</td><td>1.131</td><td>0.001</td></tr><tr><td>σ<sub>s</sub> and σ<sub>t</sub></td><td>-6.115</td><td>-6.115</td><td>-6.116</td><td>-6.115</td><td>-6.116</td><td>0.000</td></tr><tr><td>α<sub>s</sub>(m<sub>Z</sub>)(M̄S)</td><td>6.500</td><td>6.447</td><td>6.499</td><td>6.496</td><td>6.496</td><td>0.004</td></tr><tr><td>Top quark mass</td><td>-0.645</td><td>-0.652</td><td>-0.661</td><td>-0.646</td><td>-0.645</td><td>0.001</td></tr><tr><td>Total</td><td>-226.927</td><td>-264.255</td><td>-263.524</td><td>-263.289</td><td>-263.855</td><td>36.362</td></tr></tbody></table><table><thead><tr><th>Quantity</th><th>A/H-funnel</th><th>τ̃ co-ann.</th><th>τ̃ co-ann.</th><th>χ̃<sub>1</sub><sup>±</sup> co-ann.</th></tr></thead><tbody><tr><td>A<sub>0</sub></td><td>7337.758</td><td>-6666.073</td><td>-7706.626</td><td>-8213.109</td></tr><tr><td>m<sub>0</sub></td><td>4945.237</td><td>1582.304</td><td>3294.531</td><td>2697.314</td></tr><tr><td>m<sub>1/2</sub></td><td>4981.246</td><td>2265.444</td><td>2085.463</td><td>2607.561</td></tr><tr><td>m<sub>H<sub>u</sub></sub></td><td>6845.748</td><td>3714.036</td><td>5196.468</td><td>6282.001</td></tr><tr><td>m<sub>H<sub>d</sub></sub></td><td>93.459</td><td>9285.571</td><td>8990.311</td><td>4005.580</td></tr></tbody></table>

for a canonical 'ideal' likelihood (see text), along with its offset from the global best-fit point. SLHA1 and SLHA2 files corresponding to the best-fit point in each region can be found in the online data associated with this paper [163]
<table><caption>**Table 7 continued**</caption><thead><tr><th>Quantity</th><th>A/H-funnel</th><th>τ̃ co-ann.</th><th>τ̃ co-ann.</th><th>χ̃<sup>±</sup><sub>1</sub> co-ann.</th></tr></thead><tbody><tr><td>tan β</td><td>28.221</td><td>22.567</td><td>23.345</td><td>18.075</td></tr><tr><td>sgn(μ)</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>m<sub>t</sub></td><td>173.246</td><td>173.479</td><td>173.388</td><td>173.328</td></tr><tr><td>α<sub>s</sub>(m<sub>Z</sub>)(M&#773;S)</td><td>0.119</td><td>0.119</td><td>0.119</td><td>0.119</td></tr><tr><td>ρ<sub>0</sub></td><td>0.396</td><td>0.388</td><td>0.405</td><td>0.381</td></tr><tr><td>σ<sub>s</sub></td><td>43.162</td><td>42.562</td><td>43.121</td><td>43.323</td></tr><tr><td>σ<sub>f</sub></td><td>57.980</td><td>58.022</td><td>57.890</td><td>57.764</td></tr><tr><td>M<sub>1</sub></td><td>2277.442</td><td>1004.143</td><td>925.176</td><td>1157.614</td></tr><tr><td>μ</td><td>537.021</td><td>-2480.773</td><td>-1928.496</td><td>-382.757</td></tr><tr><td>m<sub>τ̃<sub>1</sub></sub></td><td>7589.989</td><td>1030.595</td><td>948.763</td><td>1217.299</td></tr><tr><td>m<sub>τ̃<sub>1</sub></sub></td><td>4633.573</td><td>1083.376</td><td>3001.595</td><td>2261.195</td></tr><tr><td>m<sub>A</sub></td><td>1176.568</td><td>9151.605</td><td>8624.785</td><td>3808.674</td></tr><tr><td>m<sub>h</sub></td><td>125.377</td><td>124.398</td><td>125.173</td><td>125.414</td></tr><tr><td>m<sub>χ̃<sup>0</sup><sub>1</sub></sub></td><td>553.377</td><td>1004.076</td><td>930.008</td><td>391.009</td></tr><tr><td>(%bino, %Higgsino)</td><td>(0, 100)</td><td>(100, 0)</td><td>(100, 0)</td><td>(0, 100)</td></tr><tr><td>m<sub>χ̃<sup>0</sup><sub>2</sub></sub></td><td>-555.848</td><td>1868.405</td><td>1734.260</td><td>-396.274</td></tr><tr><td>(%bino, %Higgsino)</td><td>(0, 100)</td><td>(0, 1)</td><td>(0, 6)</td><td>(0, 100)</td></tr><tr><td>m<sub>χ̃<sup>±</sup><sub>1</sub></sub></td><td>554.943</td><td>1868.573</td><td>1734.450</td><td>394.095</td></tr><tr><td>(%wino, %Higgsino)</td><td>(0, 100)</td><td>(99, 1)</td><td>(94, 6)</td><td>(0, 100)</td></tr><tr><td>m<sub>g̃</sub></td><td>9979.887</td><td>4715.895</td><td>4471.116</td><td>5436.877</td></tr><tr><td>Ωh<sup>2</sup></td><td>3.06 × 10<sup>-2</sup></td><td>6.76 × 10<sup>-2</sup></td><td>4.49 × 10<sup>-2</sup></td><td>1.81 × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 10** Sparticle mass spectrum of the NUHM2 best-fit point</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11** 1D profile likelihood ratio for $\Delta a_{\mu}$ in the CMSSM (red), NUHM1 (blue) and NUHM2 (purple). For comparison we show a Gaussian likelihood for the observed discrepancy $a_{\mu,\text{obs}} - a_{\mu,\text{SM}} = (28.7 \pm 8.0) \times 10^{-10}$ (green), adding the experimental and theoretical uncertainties in quadrature</figcaption></figure>

Run II, but there is an interesting peak of high likelihood at low stop masses for all three models, centred on the best-fit masses of 592, 1030 and 950 GeV for the CMSSM, NUHM1 and NUHM2 respectively. At least naively, this appears worthy of further investigation for each model, in terms of the potential for discovery at the LHC.

Concentrating first on the profile likelihood for $m_{\tilde{\tau}_1}$ in the CMSSM, the first consideration is the mass difference $m_{\tilde{\tau}_1} - m_{\tilde{\chi}_1^0}$ for models with a low stop mass, as experimental prospects generally deteriorate rapidly for more compressed spectra. The CMSSM 1D profile likelihood ratio for the mass

difference $m_{\tilde{\tau}_1} - m_{\tilde{\chi}_1^0}$ is shown in the top panels of Fig. 13 in red, while Fig. 14 shows the 2D profile likelihood in the $\tilde{\tau}_1 - \tilde{\chi}_1^0$ mass plane. The low-mass stop solutions all satisfy the relic density constraint through stop co-annihilation, giving stop-neutralino mass differences below $\sim 50$ GeV. For very small mass differences, below the mass of the b quark,
<figure><img src="image_14.png" /><figcaption>**Fig. 12** 1D profile likelihoods for the masses of $\tilde{g}$, $\tilde{t}_1$, $\tilde{b}_1$, $\tilde{\tau}_1$, $\tilde{\chi}_1^0$ and $\tilde{\chi}_1^{\pm}$ in the CMSSM (red), NUHM1 (blue) and NUHM2 (purple)</figcaption></figure>

these points could be probed by long-lived particle searches at the LHC. We defer a detailed study of this to future work.

If the stop decays promptly, however, this region can in principle be probed by LHC compressed spectra searches, particularly in the recent Run II updates that were not included in our initial scan. Although we plan a detailed analysis of the full range of recent LHC results in a forthcoming paper, some insight can be gained by examining the

recent 36 fb⁻¹ simplified model limits presented by the CMS experiment [291,321–324] at 13 TeV. They carried out stop searches in a variety of final states, and interpreted them in terms a model in which stop pair production is immediately followed by decay to a (possibly off-shell) top quark and the lightest neutralino. Although this is not necessarily the case for our models, the simplified model limit acts as a guide to the strongest possible exclusion potential of these
<figure><img src="image_15.png" /><figcaption>**Fig. 13** Left: 1D profile likelihoods for the mass differences $m_{\tilde{t}_1} - m_{\tilde{\chi}_1^0}$ (top) and $m_{\tilde{\chi}_1^+} - m_{\tilde{\chi}_1^0}$ (bottom), in the CMSSM (red), NUHM1 (blue) and NUHM2 (purple). Right: as per left-hand plots, but zoomed to focus on the smallest mass differences</figcaption></figure>

Run II searches. We show this limit in Fig. 14 as a red line. The low-mass part of our 2$\sigma$ best-fit region remains out of reach of the latest CMS search. We have also checked that the models in this region emerge almost unscathed when compared to recent ATLAS limits on compressed stop scenarios [325–327],<fnref n="16" /> but there is some hope that at least the lower parts of this region will be probed in the near future. Completely excluding the stop co-annihilation region in the CMSSM would require probing compressed spectra in lightest stop decays up to a stop mass of approximately 900 GeV. Although finding such models is challenging at the LHC, stop pair-production is within the kinematic reach of a multi-TeV linear collider for the whole region, and dedicated analysis, similar to searches for Higgsino-dominated neutralinos, should be effective in constraining such models.

This picture changes in the NUHM1, which is most easily seen by examining which mechanism for obeying

the relic density constraint is active in each region of the $\tilde{t}_1 - \tilde{\chi}_1^0$ mass plane. Figure 15 shows that, whereas the entire CMSSM 95% CL region at low stop masses arises from stop co-annihilation, the extra freedom in the NUHM1 model allows the existence of points with low stop mass that generate the required relic density through either the stau co-annihilation or chargino co-annihilation mechanisms (or indeed some combination thereof). There is hence a region with stop masses below 1 TeV that would exhibit larger $\tilde{t}_1 - \tilde{\chi}_1^0$ mass differences, making future discovery at the LHC an easier prospect. Indeed, comparison with the most recent CMS simplified model limits demonstrates that part of this chargino co-annihilation region may already have been probed [319]. Still, the region of highest likelihood is the stop co-annihilation region, with $m_{\tilde{t}_1} - m_{\tilde{\chi}_1^0} \lesssim 50$ GeV. This can be seen in the top panels of Fig. 13 in blue. Excluding the stop co-annihilation mechanism entirely in the NUHM1 is more difficult than in the CMSSM, requiring the ability to probe compressed spectra for $\tilde{t}_1$ masses up to approximately 1700 GeV, as seen in Fig. 15. The situation in the NUHM2 model (not shown) is qualitatively similar.

<footnote n="16">¹⁶ We note that the ATLAS limit assumes a 100% branching fraction for the process $\tilde{t}_1 \to c\tilde{\chi}_1^0$. We have checked that this agrees closely with the branching fractions returned by DecayBit and SUSY-HIT for our best-fit stop co-annihilation point.</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 14** 2D profile likelihoods for the CMSSM, plotted in the $\tilde{\tau}_1 - \tilde{\chi}_1^0$ mass plane. Left: the full range of neutralino masses present in the combined sample. Right: as per the lefthand panle, but zoomed in to focus on the low-mass region. Superimposed in red is the latest CMS Run II simplified model limit for $\tilde{\tau}_1$ pair production, followed by decay to $t$ quarks and the lightest neutralino [319]. This limit should be interpreted with caution (for details see main text)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15** 95% CL 2D profile likelihoods in the $\tilde{\tau}_1 - \tilde{\chi}_1^0$ mass plane, coloured according to the mechanism(s) active in depleting the relic density. Left: the CMSSM. Right: the NUHM1. Superimposed in red is the latest CMS Run II simplified model limit for $\tilde{\tau}_1$ pair production and subsequent decay to $t$ quarks and the lightest neutralino [319]. This limit should be interpreted with caution (for details see main text)</figcaption></figure>

Figure 12 also shows the presence of a region with relatively small $\tilde{\tau}_1$ masses, particularly in the NUHM1 and NUHM2. These masses are, however, already too large to lead to substantial stau production, given the small direct production cross-section [329] at 13 TeV.

Finally, we consider the prospects for discovery of charginos and neutralinos in the CMSSM, NUHM1 and NUHM2, as there are high-likelihood model points with relatively low $\tilde{\chi}_1^0$ and $\tilde{\chi}_1^{\pm}$ masses in all three models. In Figure 16, we show the profile likelihood ratio in the $\tilde{\chi}_1^{\pm} - \tilde{\chi}_1^0$ mass plane for the CMSSM, as well as a version colour-coded by the mechanism for depleting the relic density. For low masses, there is always a strict correlation between the $\tilde{\chi}_1^{\pm}$ and $\tilde{\chi}_1^0$ masses in the CMSSM. The $m_{\tilde{\chi}_1^{\pm}} \sim m_{\tilde{\chi}_1^0}$ correlation is a consequence of a Higgsino-dominated lightest neutralino, which is always accompanied by a Higgsino-dominated chargino with a similar mass and leads to chargino co-annihilation.

In the stop co-annihilation region, the neutralino is dominantly bino, and the chargino is mostly wino, with a mass about twice that of the neutralino. This comes from the approximate 2:1 ratio between the low-scale wino and bino mass parameters, produced by RGE running from the common GUT-scale input value $m_{1/2}$. We also show the envelope of the latest CMS simplified model interpretations for $\tilde{\chi}_1^{\pm} \tilde{\chi}_1^0$ production and decay with decoupled sleptons [330–333]. This should only be used as an indicator of the optimum CMS exclusion power, as we have not performed a detailed examination of our model points to check that the EW gaugino mixing matrices and decay branching ratios match the CMS assumptions. Only the low-mass tip of the stop co-annihilation part of our $2\sigma$ region has been probed by the most recent CMS analyses, and the best fit point is far beyond the current LHC reach.
<figure><img src="image_18.png" /><figcaption>**Fig. 16** Left: Profile likelihood for the CMSSM in the $\tilde{\chi}_1^{\pm} - \tilde{\chi}_1^0$ mass plane. Right: Colour-coding shows the mechanism(s) that allow models within the 95% CL region to avoid exceeding the observed relic density of DM. Superimposed in red is the latest CMS Run II simplified</figcaption></figure>

model limit for $\tilde{\chi}_1^{\pm} \tilde{\chi}_1^0$ pair production and decay with decoupled sleptons [320]. This limit should be interpreted with caution (for details see main text)

<figure><img src="image_19.png" /><figcaption>**Fig. 17** Left: Profile likelihood for the NUHM1 in the $\tilde{\chi}_1^{\pm} - \tilde{\chi}_1^0$ mass plane. Right: Colour-coding shows the mechanism(s) that allow models within the 95% CL region to avoid exceeding the observed relic density</figcaption></figure>

of DM. Superimposed in red is the latest CMS Run II simplified model limit for $\tilde{\chi}_1^{\pm} \tilde{\chi}_1^0$ pair production and decay via sleptons [328]. This limit should be interpreted with caution (for details see main text)

Figure 17 shows the profile likelihood ratio in the $\tilde{\chi}_1^{\pm} - \tilde{\chi}_1^0$ mass plane for the NUHM1. The low-mass region now sees a contribution from stau co-annihilation in addition to chargino co-annihilation. This is interesting for LHC searches, as the assumption of decoupled sleptons clearly no longer applies. The recent CMS simplified model interpretations include a model where the sleptons are not decoupled, but the interpretation is even more fraught than that of the previous simplified models we have considered. The slepton masses are fixed in these scenarios, and one can generically expect the strength of the exclusions to decrease as one departs both from the mass assumptions, and from the branching ratio assumptions. Nevertheless, we show this limit in Fig. 17 in order to demonstrate the most optimistic possible exclusion,

compared to our 2D profile likelihood. Almost the entire region with compressed spectra remains unprobed. As the bottom right plot of Fig. 13 shows, the highest likelihood in the degenerate region is obtained for chargino-neutralino mass differences less than 15 GeV, which is small enough to escape the CMS searches. However, the part of our low-mass 95% CL region without degenerate $\tilde{\chi}_1^{\pm} - \tilde{\chi}_1^0$ masses may be within current LHC reach. Furthermore, there is hope that the LHC would prove capable of exploring part of the 68% CL region in the near future. The situation in the NUHM2 is shown in Fig. 18. The main difference with the NUHM1 is that the low-mass region potentially explored by the recent CMS analysis updates falls within our 68% CL preferred region.
<figure><img src="image_20.png" /><figcaption>**Fig. 18** Left: Profile likelihood for the NUHM2 in the $\tilde{\chi}_1^\pm - \tilde{\chi}_1^0$ mass plane. Right: Colour-coding shows mechanism(s) that allow models within the 95% CL region to avoid exceeding the observed relic density</figcaption></figure>

of DM. Superimposed in red is the latest CMS Run II simplified model limit for $\tilde{\chi}_1^\pm \tilde{\chi}_1^0$ pair production and decay via sleptons [328]. This limit should be interpreted with caution (for details see main text)

<figure><img src="image_21.png" /><figcaption>**Fig. 19** The spin-independent neutralino-proton cross-section. Upper left: Profile likelihood in the CMSSM. Lower left: Colour-coding shows the active mechanism(s) by which CMSSM models avoid exceeding the observed relic density of DM, through either chargino co-annihilation, the A/H funnel, or stop co-annihilation. Top right: Colour-coded regions in the NUHM1, now also featuring stau co-annihilation (blue). Bottom right: Colour-coded regions of the NUHM2. 90% CL exclusion limits are overlaid from the complete LUX exposure [228],</figcaption></figure>

the projected reach of XENON1T with two years of exposure, the projected reach of XENONnT/LZ with 20 tonne-years of exposure [334] (around 1–3 years of data), and the projected reach of DARWIN with 200 tonne-years of exposure [335] (around 3–4 years of data). The “neutrino floor”, where the coherent neutrino background starts to limit the experimental sensitivity, is indicated by the dashed grey line [336]. The exact position of this limit is subject to several caveats; see [336] for further details
<figure><img src="image_22.png" /><figcaption>**Fig. 20** The spin-dependent neutralino-proton cross-section. Upper left: Profile likelihood in the CMSSM. Bottom left: Colour-coding shows the active mechanism(s) by which CMSSM models avoid exceeding the observed relic density of DM, through either chargino co-annihilation, the A/H funnel, or stop co-annihilation. Top right: Colour-coded regions in the NUHM1, now also featuring stau co-annihilation</figcaption></figure>

(blue). Bottom right: Colour-coded regions of the NUHM2. 90% CL exclusion limits are overlaid from the 79-string IceCube search for DM [7,222], assuming dark matter annihilation in the Sun to $\bar{b}b$ (yellow solid) and $\tau^+\tau^-$ (red solid) final states, from PICO-60 [235] (green solid), and projected limits from PICO-250 [337] (green dashes)

### 4.4.2 Direct detection

Now we turn to a discussion of the discovery prospects at future direct DM detection experiments. In Figs. 19 and 20, we show the spin-independent (SI) and spin-dependent (SD) nuclear scattering cross-sections of the lightest neutralino as a function of its mass, scaled for the fraction of the local density of DM in neutralinos. We give the full profile likelihood for the CMSSM only. The other panels show $2\sigma$ confidence regions for each model, with color-coded mechanisms for reducing the relic density to or below the observed value.

In the CMSSM, the chargino co-annihilation and A/H-funnel regions largely overlap, predicting SI cross-sections in the range $10^{-46}-10^{-44}$ cm$^2$ and neutralino masses from 200 to 2500 GeV. This region will be almost fully probed by XENON1T after two years of data-taking (long dashed curve), and could be completely excluded by XENONnT or LZ with 1–3 years of data (short dashed curve). The stop

region has significantly lower SI cross-sections, with an LSP that can be almost pure bino in nature. Most of this region is outside of the projected reach of future multi-tonne detectors, although the proposed ~50-tonne DARWIN experiment may probe it slightly (dotted curve). As discussed in Sect. 4.4.1, this region is also difficult to see at the LHC, but may be within the reach of a future linear collider.

The NUHM1 and NUHM2 display similar properties, with large parts of the chargino co-annihilation and A/H-funnel regions able to be tested in the near future, including models with very heavy LSPs (in the NUHM2 case, up to 4.5 TeV). Some of the chargino co-annihilation region at relatively low LSP masses (< 1250 GeV) will remain untested by future direct detection experiments until DARWIN. Parts of both the stop and stau co-annihilation regions will escape all direct detection, even DARWIN, although other parts of this region at higher masses will be easily detected or excluded by XENON1T.
In Fig. 20 we show the rescaled spin-dependent neutralino-proton scattering cross-sections for the CMSSM, NUHM1 and NUHM2. Here we overplot current PICO limits not included in our scan [235], and sensitivity estimates for PICO-250, a scaled-up version of PICO-60 [337]. We also show the IceCube 79-string limits from Ref. [7], for two different annihilation final states. We can see that the preferred regions are relatively far from the current limits, so future direct detection experiments are unlikely to probe them further. However, the proposed neutrino telescopes IceCube-PINGU [338] and KM3NeT-ORCA [339] may have sufficient additional sensitivity to test the models with largest SD cross-sections. Although estimates of the expected sensitivity of these experiments to $\sigma_{\text{SD}}$ exist, those estimates do not (yet) extend above DM masses of 100 GeV, so at present they can tell us little about prospects for discovery of the CMSSM, NUHM1 or NUHM2.

Figure 19 shows the $2\sigma$ allowed region for the SI cross-section extending to substantially lower values in the CMSSM than the NUHM1 or NUHM2. This seems surprising, as the CMSSM is a subspace of the NUHM1 and NUHM2, so all viable CMSSM models are indeed also viable NUHM1 and NUHM2 models. The improvement in the best-fit likelihood in the NUHM1 compared to the CMSSM is not sufficient to explain this effect. The smallest scattering cross-sections are caused by cancellations in the tree-level matrix elements, which can be tuned to essentially arbitrary accuracy. A consequence of this is that models become steadily more fine-tuned as the cross-section asymptotically approaches zero, and therefore steadily more difficult to find for sampling algorithms. What we see here is evidence of the additional numerical difficulty of finding such points in the NUHM1 and NUHM2, due to the additional challenge of dealing with more dimensions, and a more diverse set of viable regions of parameter space. However, in models where the mass parameters unify at a high scale, loop corrections [340,341] are expected to spoil such carefully-tuned cancellations anyway, holding cross-sections well above the lowest values that we see in the CMSSM [342]. The fact that we have found scattering cross-sections as low as $10^{-60} \text{ cm}^2$ in the CMSSM, but not at quite such low values in the NUHM1 and NUHM2, is therefore ultimately of little physical significance. Even if this isn’t physically significant however, getting as low as $10^{-60} \text{ cm}^2$ in the CMSSM is nonetheless quite a remarkable numerical feat, made possible only by our use of *Diver*. This increases our confidence in the completeness of our sampling in the rest of the parameter space, and in fits of weak-scale MSSM models [162].

### 4.4.3 Indirect detection

To assess the discovery prospects for future indirect searches for DM, in Fig. 21 we show the rescaled zero-velocity annihila-

hilation cross-section $f^2 \cdot \langle \sigma v \rangle_0$, as a function of the mass of the lightest neutralino. Here $f$ is again the ratio of the neutralino relic density in the model to the observed relic density of DM. Note that we implicitly assume here that if neutralinos are not all of the DM, the other component(s) of DM cluster in the same way as neutralinos, leading to the same $f$ cosmologically, in dwarf galaxies and in the local halo. Although this needn’t be true in general, the general requirements that DM be cold and (almost) non-interacting mean that this should be a reasonably good approximation.

The upper left panel of Fig. 21 shows the profile likelihood for the CMSSM, and the remaining panels show the mechanisms by which models in the CMSSM (bottom left), NUHM1 (top right) and NUHM2 (bottom right) avoid producing too much thermal DM. In the same figure we also indicate, for comparison, current limits from dwarf galaxy observations by the Fermi-LAT [220], assuming photon spectra for DM annihilation to $\bar{b}b$ and $\tau^+\tau^-$ final states. We also show projected Fermi limits for $\bar{b}b$ final states [343], assuming 15 years of data on 60 dwarf galaxies (vs 6 years and 15 dwarfs in the current limits). Lastly, we show the projected sensitivity of CTA after 500 h of observation of the Galactic halo, also assuming $\bar{b}b$ final states [344]. We note that the actual (projected) limits depend on the final state, and hence the lines shown are only indicative for general points in the CMSSM parameter space. However, as long as the final states are hadronic, the expected variations remain within a factor of about three [345].

In general, the largest annihilation cross-sections are expected for the $A/H$ funnel region, where resonant annihilation boosts $\sigma v$. All models with annihilation cross-sections above the canonical thermal value ($3 \times 10^{-26} \text{ cm}^3 \text{ s}^{-1}$) exhibit resonant annihilation through the $A$ funnel (note that in the zero-velocity limit, due to the $CP$ properties of the initial state, only the pseudoscalar resonances can contribute). As one would expect, all regions above this value in Fig. 21 are indeed identified as being part of the $A/H$ funnel region, indicated by the fact that they are shaded orange. Some parts of these regions are also shaded yellow and/or blue, as some of the parameter points identified as belonging to the $A/H$ funnel region also satisfy the necessary mass/composition conditions to be counted as part of the stau and/or chargino co-annihilation regions. However, in all regions of overlap above $f^2 \cdot \langle \sigma v \rangle_0 \sim 3 \times 10^{-26} \text{ cm}^3 \text{ s}^{-1}$, resonant annihilation via the heavy Higgs bosons is the dominant mechanism in setting the relic density. Most such models exhibit a relic density below the observed value.<fnref n="17" />

<footnote n="17">¹⁷ The fact that $\sigma v$ is set by a resonance in the funnel region means that the present-day annihilation cross-section can be somewhat higher or lower than during freeze-out; models where $\langle \sigma v \rangle_0 > 3 \times 10^{-26} \text{ cm}^3 \text{ s}^{-1}$ but $\Omega_{ch}^2$ matches the observed value exhibit this effect. This allows the observed relic density to be achieved up to $\langle \sigma v \rangle_0 \sim 7 \times 10^{-26} \text{ cm}^3 \text{ s}^{-1}$</footnote>
<figure><img src="image_23.png" /><figcaption>**Fig. 21** The present-day neutralino self-annihilation cross-section. Upper left: Profile likelihood in the CMSSM. Bottom left: Colour-coding shows the active mechanism(s) by which CMSSM models avoid exceeding the observed relic density of DM, through either chargino co-annihilation, the A/H funnel, or stop co-annihilation. Top right: Colour-coded regions in the NUHM1, now also featuring stau co-annihilation (blue). Bottom right: Colour-coded regions of the NUHM2. 95% CL exclusion limits are overlaid from the 6-year *Fermi*-LAT search for DM</figcaption></figure>

annihilation in 15 satellite dwarf galaxies [220], assuming dark matter annihilation to $\bar{b}b$ (yellow solid) and $\tau^+\tau^-$ (red solid) final states. We also show the projected improvement for $b\bar{b}$ final states with 15 years of LAT data and four times as many dwarfs [343] (dashed yellow), and an optimistic projection of the sensitivity to $b\bar{b}$ final states of a Galactic halo search for DM annihilation by the upcoming Cherenkov Telescope Array, assuming 500 h of observations and no systematic uncertainties [344] (green dashes)

Indeed, most of the high-mass models identified in our scans as having stau and/or chargino co-annihilation also exhibit resonant annihilation through the heavy Higgs funnel.<fnref n="18" /> Indeed, as discussed at the beginning of Sec. 4 in the context of Fig. 1, above DM masses of around a TeV, chargino co-annihilation in the CMSSM is unable to deplete the relic density to the thermal value or below without additional assistance from the A/H resonance – so in fact, all chargino co-annihilation models above about a TeV are hybrid models of some kind. This point is confirmed by careful study of

all other CMSSM, NUHM1 and NUHM2 plots in this section: for $m_\chi \gtrsim 1$ TeV, chargino co-annihilation regions only appear in the presence of the heavy Higgs funnel, indicating that all high-mass chargino co-annihilation models are in fact either hybrids with, or completely dominated by, the heavy Higgs funnel. The NUHM1 and NUHM2 plots also show that the same situation holds for stau co-annihilation at $m_\chi \gtrsim 1.5$ TeV in the NUHM, which at such masses appears only in combination with the heavy Higgs funnel.

In the CMSSM, *Fermi* will generally only probe the low-likelihood tails of the A/H funnel region and its co-annihilation hybrid, with the exception of ~1 TeV Higgsinos (see below). Taking the optimistic predictions of Ref. [344] at face value, CTA will significantly cut into this region. However, this ignores the impact of detector and background systematics. Adding a systematic uncertainty of 1%, the sensi-

<footnote>Footnote 17 continued
in the CMSSM, $\langle\sigma v\rangle_0 \sim 10^{-25} \text{ cm}^3 \text{ s}^{-1}$ in the NUHM1, and $\langle\sigma v\rangle_0 \sim 3 \times 10^{-25} \text{ cm}^3 \text{ s}^{-1}$ in the NUHM2.</footnote>

<footnote n="18">¹⁸ We remind the reader that shaded models are all those that exhibit a given relic density mechanism – not those that *only* exhibit that mechanism.</footnote>
tivity would degrade by a factor of ~6 [346], and hence probe a significantly smaller part of the funnel and/or chargino co-annihilation region. The stop co-annihilation region, on the other hand, will remain largely unconstrained in the CMSSM, even with future indirect detection missions.

Let us stress, however, that indirect detection prospects will generally be better than indicated by this general discussion. One aspect not taken into account here is the impact of radiative corrections to the annihilation rate, which are particularly relevant for large neutralino masses [221,347,348]. In the stop co-annihilation region, for example, the inclusion of Higgs-Strahlung off fermion final states can increase the total annihilation rate by a factor of a few compared to what we have implemented [348]. Also antiprotons can be an efficient complementary probe of compressed mass spectra [349]. For heavy neutralinos with mass-degenerate charginos, another example is a distinct feature from $W^{+}W^{-}\gamma$ final states [350], which adds to the already large monochromatic line signal from such models [351–353]; such a signal is much more easily distinguished from astrophysical backgrounds than the spectrum from $\bar{b}b$ final states assumed for CTA in Fig. 21. This signal would also appear in observations of the Galactic Centre well before any dwarf observations shown in Fig. 21 or taken into account in our scans. Last but not least, let us mention the Sommerfeld effect [354–356], which leads to a large enhancement in particular for ~1 TeV Higgsino DM [35, 356, 357], and which we have not (yet) included in GAMBIT.

In the NUHM1 and NUHM2, Fermi appears to be just beginning to constrain A-funnel models at masses of around 1.7 TeV. The various mechanisms to suppress the relic density are not as well-separated in the $\langle\sigma v\rangle_0 - m_{\chi}$ plane in NUHM models as in the CMSSM, however. As a consequence, all these mechanisms can be (partially) tested with CTA. We note that this includes the stau co-annihilation region. The sensitivity curves shown here are overly conservative for such models, because $\tau^{+}\tau^{-}\gamma$ final states will be much more constraining than $\bar{b}b$ spectra [358]. Still, even with the most optimistic assumptions, large parts of the viable parameter spaces of all of the GUT-scale models that we consider here will remain impossible to probe with indirect DM searches.

## 5 Conclusions

In this paper we have presented state-of-the-art profile likelihood global fits to three constrained versions of the minimal supersymmetric standard model, using GAMBIT. We have incorporated updated experimental data, additional observables and improved calculations for many quantities compared to previous global fits. We have also fully explored the parameter space in which the models are not excluded by any experimental measurements, specifically including areas

where the neutralino only constitutes a fraction of the dark matter in the Universe.

In the CMSSM, we show that the stau co-annihilation region is finally ruled out at more than 95% CL. This comes about due to Run II LHC constraints, difficulty in fitting the Higgs mass in this region, and an overall lifting of the isosymmetry contours defining the boundaries of this region, brought about by our improved sampling in this paper and resulting discovery of what is effectively a better best-fit than in previous works. The NUHM1 and NUHM2 allow more freedom, permitting lighter staus and a re-appearance of the stau co-annihilation region as a source of equally good fits as other mechanisms for depleting the relic density. Those include stop co-annihilation, chargino co-annihilation and resonant annihilation through the $A/H$ funnel. We find that the chargino co-annihilation region also widens substantially in the NUHM1 and NUHM2 compared to the CMSSM, extending to arbitrarily low values of $m_0$.

Current constraints from the LHC push superpartner masses towards the multi-TeV regime, even if one does not demand that the lightest neutralino is the only DM species. The important exceptions are the lightest neutralinos and charginos, which can still have masses as low as ~100 GeV without violating any experimental constraints, the lightest stau, which can be as light as ~200 GeV, and the lightest stop, which can be as light as ~500 GeV.

Despite very heavy spectra in many parts of the parameter space, future direct detection experiments will fully explore the chargino co-annihilation region, encompassing the so-called ‘focus point’.

We find a region of good fits at large negative trilinear coupling, where the neutralino relic abundance is set by co-annihilation with the lightest stop. The trilinear couplings in this region raise questions about colour- and charge-breaking vacua, but our tests indicate that large parts of this region remain unaffected by such considerations. More detailed investigation would however be interesting. This region has been properly seen only in very recent fits performed contemporaneously with this one [115]. Models in this region feature quite light stops (~500 GeV), making it very appealing from the point of view of electroweak naturalness. However, the stop co-annihilation mechanism requires the neutralino-stop mass difference to be quite small, which may constitute a fine-tuning in itself. A more detailed analysis of naturalness considerations, including a full Bayesian treatment of the fit, would be illuminating. Models in this region will be challenging to discover at the LHC, and next to impossible at direct detection experiments, but are promising targets for a future linear collider.

We began this study mainly intending to validate the new generic beyond-the-Standard-Model global fitting framework GAMBIT. In the end however, we have found quite a few genuinely new and interesting results. This serves to
illustrate the utility of a modern and adaptive global fitter such as **GAMBIT**, where the impacts of different searches on different models can be easily examined and compared whilst retaining a consistent treatment of theoretical assumptions, systematics, nuisances, scanning algorithms, statistical approaches, experimental analyses and external code interfaces.

All input files, samples and best-fit benchmarks produced for this paper are publicly accessible from Zenodo [163]. The GAMBIT software is available from gambit.hepforge.org.

**Acknowledgements** We thank Andrew Fowlie, Tomás Gonzalo, Julia Harz, Sebastian Hoof, Felix Kahlhoefer, James McKay, Roberto Trotta and Sebastian Wild for useful discussions, and Lucien Boland, Sean Crosby and Goncalo Borges of the Australian Centre of Excellence for Particle Physics at the Terascale for computing assistance and resources. We warmly thank the Casa Matemáticas Oaxaca, affiliated with the Banff International Research Station, for hospitality whilst part of this work was completed, and the staff at Cyfronet, for their always helpful supercomputing support. GAMBIT has been supported by STFC (UK; ST/K00414X/1, ST/P000762/1), the Royal Society (UK; UF110191), Glasgow University (UK; Leadership Fellowship), the Research Council of Norway (FRIPRO 230546/F20), NOTUR (Norway; NN9284K), the Knut and Alice Wallenberg Foundation (Sweden; Wallenberg Academy Fellowship), the Swedish Research Council (621-2014-5772), the Australian Research Council (CE110001004, FT130100018, FT140100244, FT160100274), The University of Sydney (Australia; IRCA-G162448), PLGrid Infrastructure (Poland), Red Española de Supercomputación (Spain; FI-2016-1-0021), Polish National Science Center (Sonata UMO-2015/17/D/ST2/03532), the Swiss National Science Foundation (PP00P2-144674), the European Commission Horizon 2020 Marie Skłodowska-Curie actions (H2020-MSCA-RISE-2015-691164), the ERA-CAN+ Twinning Program (EU & Canada), the Netherlands Organisation for Scientific Research (NWO-Vidi 680-47-532), the National Science Foundation (USA; DGE-1339067), the FRQNT (Québec) and NSERC/The Canadian Tri-Agencies Research Councils (BPDF-424460-2012).

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.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³.

## References

1. S. Profumo, T. Stefaniak, Alignment without decoupling: the portal to light dark matter in the MSSM. Phys. Rev. D **94**, 095020 (2016). arXiv:1608.06945

2. D.P. Roy, SUSY Dark Matter in Universal and Nonuniversal Gaugino Mass Models. arXiv:1608.03062

3. J. Harz, B. Herrmann, M. Klasen, K. Kovarik, P. Steppeler, Theoretical uncertainty of the supersymmetric dark matter relic density from scheme and scale variations. Phys. Rev. D **93**, 114023 (2016). arXiv:1602.08103

4. Q. Riffard, F. Mayet, G. Bélanger, M.H. Genest, D. Santos, Extracting constraints from direct detection searches of supersymmetric dark matter in the light of null results from the LHC in the squark sector. Phys. Rev. D **93**, 035022 (2016). arXiv:1602.01030

5. M. van Beekveld, W. Beenakker, S. Caron, R. Ruiz de Austri, The case for 100 GeV bino dark matter: A dedicated LHC tri-lepton search. JHEP **04**, 154 (2016). arXiv:1602.00590

6. M. Beneke, A. Bharucha et al., Relic density of wino-like dark matter in the MSSM. JHEP **03**, 119 (2016). arXiv:1601.04718

7. IceCube Collaboration: M.G. Aartsen et al., Improved limits on dark matter annihilation in the Sun with the 79-string IceCube detector and implications for supersymmetry. JCAP **04**, 022 (2016). arXiv:1601.00653

8. A. Berlin, D.S. Robertson, M.P. Solon, K.M. Zurek, Bino variations: Effective field theory methods for dark matter direct detection. Phys. Rev. D **93**, 095008 (2016). arXiv:1511.05964

9. J. Cao, Y. He, L. Shang, W. Su, Y. Zhang, Testing the light dark matter scenario of the MSSM at the LHC. JHEP **03**, 207 (2016). arXiv:1511.05386

10. G. Grilli di Cortona, Searching SUSY from below. PoS PLANCK2015, 054 (2015). arXiv:1510.07616

11. H. Baer, (Mainly) axion dark matter. arXiv:1510.07501 [AIP Conf. Proc.1743,050002(2016)]

12. J. Harz, B. Herrmann, M. Klasen, K. Kovarik, P. Steppeler, Precise prediction of the dark matter relic density within the MSSM. PoS EPS-HEP2015, 410 (2015). arXiv:1510.06295

13. K. Hamaguchi, K. Ishikawa, Prospects for Higgs- and Z-resonant Neutralino Dark Matter. Phys. Rev. D **93**, 055009 (2016). arXiv:1510.05378

14. A.L. Erickcek, K. Sinha, S. Watson, Bringing isolated dark matter out of isolation: late-time reheating and indirect detection. Phys. Rev. D **94**, 063502 (2016). arXiv:1510.04291

15. J. Ellis, J.L. Evans, F. Luo, K.A. Olive, Scenarios for gluino coannihilation. JHEP **02**, 071 (2016). arXiv:1510.03498

16. J. Bramante, N. Desai et al., Towards the final word on neutralino dark matter. Phys. Rev. D **93**, 063525 (2016). arXiv:1510.03460

17. M. Baumgart, V. Vaidya, Semi-inclusive wino and higgsino annihilation to LL'. JHEP **03**, 213 (2016). arXiv:1510.02470

18. T. Bringmann, A.J. Galea, P. Walia, Leading QCD corrections for indirect dark matter searches: A fresh look. Phys. Rev. D **93**, 043529 (2016). arXiv:1510.02473

19. H. Eberl, V.C. Spanos, GravitinoPack and decays of supersymmetric metastable particles. Comput. Phys. Commun. **202**, 310-325 (2016). arXiv:1509.09159

20. K. Freese, A. Lopez, N.R. Shah, B. Shakya, MSSM A-funnel and the galactic center excess: prospects for the LHC and direct detection experiments. JHEP **04**, 059 (2016). arXiv:1509.05076

21. M. Berggren, A. Cakir et al., Non-simplified SUSY: $\tilde{\tau}$-coannihilation at LHC and ILC. Eur. Phys. J. C **76**, 183 (2016). arXiv:1508.04383

22. S. Raby, Long-lived gluinos and stable axinos. Phys. Rev. Lett. **115**, 231801 (2015). arXiv:1508.04373

23. L. Roszkowski, S. Trojanowski, K. Turzynski, Axino dark matter with low reheating temperature. JHEP **11**, 139 (2015). arXiv:1507.06164

24. J. Bramante, P.J. Fox et al., Relic neutralino surface at a 100 TeV collider. Phys. Rev. D **91**, 054015 (2015). arXiv:1412.4789

25. G. Arcadi, L. Covi, M. Nardecchia, Gravitino dark matter and low-scale baryogenesis. Phys. Rev. D **92**, 115006 (2015). arXiv:1507.05584

26. A. Achterberg, M. van Beekveld, W. Beenakker, S. Caron, L. Hendriks, Comparing galactic center MSSM dark matter solutions to the Reticulum II gamma-ray data. JCAP **1512**, 013 (2015). arXiv:1507.04644

27. A. Arbey, M. Battaglia, L. Covi, J. Hasenkamp, F. Mahmoudi, LHC constraints on gravitino dark matter. Phys. Rev. D **92**, 115008 (2015). arXiv:1505.04595
28. A. Crivellin, M. Hoferichter, M. Procura, L.C. Tunstall, Dark matter: connecting LHC searches to direct detection, In *Gravitation: 100 years after GR. Proceedings of the 50th Rencontres de Moriond, La Thuile, Italy, March 21–28, 2015* (2015), pp. 237–242. arXiv:1505.02314

29. M. Ibe, S. Matsumoto, S. Shirai, T.T. Yanagida, Wino dark matter in light of the AMS-02 2015 data. Phys. Rev. D **91**, 111701 (2015). arXiv:1504.05554

30. A. Arbey, M. Battaglia, F. Mahmoudi, The Higgs boson, Super-symmetry and dark matter: relations and perspectives. Ann. Phys. **528**, 179–186 (2016). arXiv:1504.05091

31. J. Hisano, K. Ishiwata, N. Nagata, QCD effects on direct detection of Wino dark matter. JHEP **06**, 097 (2015). arXiv:1504.00915

32. N. Nagata, H. Otono, S. Shirai, Probing bino–gluino coannihilation at the LHC. Phys. Lett. B **748**, 24–29 (2015). arXiv:1504.00504

33. J. Ellis, F. Luo, K.A. Olive, Gluino coannihilation revisited. JHEP **09**, 127 (2015). arXiv:1503.07142

34. A. Crivellin, M. Hoferichter, M. Procura, L.C. Tunstall, Light stops, blind spots, and isospin violation in the MSSM. JHEP **07**, 129 (2015). arXiv:1503.03478

35. M.E. Cabrera-Catalan, S. Ando, C. Weniger, F. Zandanel, Indirect and direct detection prospect for TeV dark matter in the nine parameter MSSM. Phys. Rev. D **92**, 035018 (2015). arXiv:1503.00599

36. A. Berlin, S. Gori, T. Lin, L.-T. Wang, Pseudoscalar portal dark matter. Phys. Rev. D **92**, 015005 (2015). arXiv:1502.06000

37. A. Achterberg, S. Amoroso et al., A description of the Galactic Center excess in the minimal supersymmetric standard model. JCAP **1508**, 006 (2015). arXiv:1502.05703

38. L. Aparicio, M. Cicoli et al., Non-thermal CMSSM with a 125 GeV Higgs. JHEP **05**, 098 (2015). arXiv:1502.05672

39. G.L. Kane, P. Kumar, B.D. Nelson, B. Zheng, Dark matter production mechanisms with a nonthermal cosmological history: A classification. Phys. Rev. D **93**, 063527 (2016). arXiv:1502.05406

40. M. Baumgart, I.Z. Rothstein, V. Vaidya, Constraints on Galactic Wino densities from gamma ray lines. JHEP **04**, 106 (2015). arXiv:1412.8698

41. S. Liebler, S. Profumo, T. Stefaniak, Light stop mass limits from Higgs rate measurements in the MSSM: Is MSSM electroweak baryogenesis still alive after all? JHEP **04**, 143 (2016). arXiv:1512.09172

42. M. Dhuria, C. Hati, U. Sarkar, Moduli induced cogenesis of baryon asymmetry and dark matter. Phys. Lett. B **756**, 376–383 (2016). arXiv:1508.04144

43. G. Barenboim, W.-I. Park, Peccei–Quinn field for inflation, baryogenesis, dark matter, and much more. Phys. Lett. B **756**, 317–322 (2016). arXiv:1508.00011

44. P. Athron, J.-H. Park, T. Steudtner, D. Stöckinger, A. Voigt, Precise Higgs mass calculations in (non-)minimal supersymmetry at both high and low scales. JHEP **01**, 079 (2017). arXiv:1609.00371

45. R.K. Barman, B. Bhattachjerjee et al., Status of MSSM Higgs sector after ICHEP 2016. arXiv:1608.02573

46. H. Bahl, W. Hollik, Precise prediction for the light MSSM Higgs boson mass combining effective field theory and fixed-order calculations. Eur. Phys. J. C **76**, 499 (2016). arXiv:1608.01880

47. P. Bechtle, H.E. Haber et al., The light and heavy Higgs interpretation of the MSSM. Eur. Phys. J. C **77**, 67 (2017). arXiv:1608.00638

48. P. Draper, H. Rzehak, A review of Higgs mass calculations in supersymmetric models. Phys. Rep. **619**, 1–24 (2016). arXiv:1601.01890

49. M. Carena, J. Ellis, J.S. Lee, A. Pilaftsis, C.E.M. Wagner, CP violation in heavy MSSM Higgs scenarios. JHEP **02**, 123 (2016). arXiv:1512.00437

50. B. Bhattachjerjee, M. Chakraborti, A. Chakraborty, U. Chattopadhyay, D.K. Ghosh, Status of the 98–125 GeV Higgs bosons scenario with updated LHC-8 data. Phys. Rev. D **93**, 075004 (2016). arXiv:1511.08461

51. J. Baglio, A. Djouadi, J. Quevillon, Prospects for Higgs physics at energies up to 100 TeV. Rep. Prog. Phys. **79**, 116201 (2016). arXiv:1511.07853

52. S. Heinemeyer, C. Schappacher, Neutral Higgs boson production at $e^+e^-$ colliders in the complex MSSM: a full one-loop analysis. Eur. Phys. J. C **76**, 220 (2016). arXiv:1511.06002

53. G. Tsamis, Reduction of couplings in the MSSM. PoS ORFU2014, 068 (2015). arXiv:1507.04469

54. E. Arganda, J. Guasch, W. Hollik, S. Penaranda, Discriminating between SUSY and non-SUSY Higgs sectors through the ratio $H \to b\bar{b}/H \to \tau^+\tau^-$ with a 125 GeV Higgs boson. Eur. Phys. J. C **76**, 286 (2016). arXiv:1506.08462

55. A. Djouadi, Higgs physics. PoS ORFU2014, 018 (2015). arXiv:1505.01059

56. L. Wu, J.M. Yang, C.-P. Yuan, M. Zhang, *Higgs self-coupling in the MSSM and NMSSM after the LHC Run 1*. Phys. Lett. B **747**, 378–389 (2015). arXiv:1504.06932

57. H. Mantler, M. Wiesemann, Hadronic Higgs production through NLO + PS in the SM, the 2HDM and the MSSM. Eur. Phys. J. C **75**, 257 (2015). arXiv:1504.06625

58. J. Pardo Vega, G. Villadoro, SusyHD: Higgs mass Determination in Supersymmetry. JHEP **07** 159 (2015). arXiv:1504.05200

59. B. Bhattachjerjee, A. Chakraborty, A. Choudhury, Status of the MSSM Higgs sector using global analysis and direct search bounds, and future prospects at the High Luminosity LHC. Phys. Rev. D **92**, 093007 (2015). arXiv:1504.04308

60. A. Djouadi, L. Maiani, A. Polosa, J. Quevillon, V. Riquer, Fully covering the MSSM Higgs sector at the LHC. JHEP **06**, 168 (2015). arXiv:1502.05653

61. A. Hebbar, Q. Shafi, C.S. Un, Light Higgsinos, Heavy Gluino and $b-\tau$ Quasi-Yukawa Unification: Will the LHC find the Gluino?. arXiv:1702.05431

62. T. Fukuyama, N. Okada, H.M. Tran, Sparticle spectroscopy of the minimal SO(10) model. Phys. Lett. B **767**, 295–302 (2017). arXiv:1611.08341

63. E. Bagnaschi et al., Likelihood analysis of supersymmetric SU(5) GUTs. Eur. Phys. J. C **77**, 104 (2017). arXiv:1610.10084

64. J. Ellis, J.L. Evans, A. Mustafayev, N. Nagata, K.A. Olive, The super-GUT CMSSM revisited. Eur. Phys. J. C **76**, 592 (2016). arXiv:1608.05370

65. M. Crispim Romão, A. Karozas, S.F. King, G.K. Leontaris, A.K. Meadowcroft, MSSM from F-theory SU(5) with Klein monodromy. Phys. Rev. D **93**, 126007 (2016). arXiv:1512.09148

66. M. Cannoni, J. Ellis, M.E. Gómez, S. Lola, R. Ruiz de Ausasi, Supersymmetry searches in GUT models with non-universal scalar masses. JCAP **1603**, 041 (2016). arXiv:1511.06205

67. J. Ellis, J.L. Evans et al., Beyond the CMSSM without an accelerator: Proton decay and direct dark matter detection. Eur. Phys. J. C **76**, 8 (2016). arXiv:1509.08838

68. D. Bourilkov, Strong coupling running, gauge coupling unification and the scale of new physics. JHEP **11**, 117 (2015). arXiv:1508.04176

69. G. Bhattacharyya, T.T. Yanagida, N. Yokozaki, Focus point gauge mediation with incomplete adjoint messengers and gauge coupling unification. Phys. Lett. B **749**, 82–87 (2015). arXiv:1506.05962

70. C.S. Aulakh, New minimal SO(10) GUT: A theory for all Epochs. Pramana **86**, 207–221 (2016). arXiv:1506.05850

71. Z. Berezhiani, M. Chianese, G. Miele, S. Morisi, Chances for SUSY-GUT in the LHC Epoch. JHEP **08**, 083 (2015). arXiv:1505.04950
72. K.S. Babu, B. Bajc, V. Susi, A minimal supersymmetric $E_6$ unified theory. JHEP 05, 108 (2015). [arXiv:1504.00904]

73. F. Wang, W. Wang, J.M. Yang, Reconcile muon g-2 anomaly with LHC data in SUGRA with generalized gravity mediation. JHEP 06, 079 (2015). [arXiv:1504.00505]

74. S. Fichet, B. Herrmann, Y. Stoll, Tasting the SU(5) nature of supersymmetry at the LHC. JHEP 05, 091 (2015). [arXiv:1501.05307]

75. F. Wang, W. Wang, J.M. Yang, A split SUSY model from SUSY GUT. JHEP 03, 050 (2015). [arXiv:1501.02906]

76. S. Heinemeyer, M. Mondragon, G. Zoupanos, The LHC Higgs boson discovery: Implications for finite unified theories. Int. J. Mod. Phys. A 29, 1430032 (2014). [arXiv:1412.5766]

77. W.G. Hollik, A new view on vacuum stability in the MSSM. JHEP 08, 126 (2016). [arXiv:1606.08356]

78. E. Bagnaschi, F. Brümmer, W. Buchmüller, A. Voigt, G. Weiglein, Vacuum stability and supersymmetry at high scales with two Higgs doublets. JHEP 03, 158 (2016). [arXiv:1512.07761]

79. N. Blinov, D.E. Morrissey, Vacuum stability and the MSSM Higgs mass. JHEP 03, 106 (2014). [arXiv:1310.4174]

80. M. Carena, S. Gori, I. Low, N.R. Shah, C.E.M. Wagner, Vacuum stability and higgs diphon decay in the MSSM. JHEP 02, 114 (2013). [arXiv:1211.6136]

81. S. Ferrara, D. Roest, General sGoldstino inflation. JCAP 1610, 038 (2016). [arXiv:1608.03709]

82. J. Ellis, M.A.G. Garcia, D.V. Nanopoulos, K.A. Olive, Phenomenological aspects of no-scale inflation models. JCAP 1510, 003 (2015). [arXiv:1503.08867]

83. F.S. Queiroz, K. Sinha, W. Wester, Rich tapestry: Supersymmetric axions, dark radiation, and inflationary reheating. Phys. Rev. D 90, 115009 (2014). [arXiv:1407.4110]

84. S. Kasuya, F. Takahashi, Flat direction inflation with running kinetic term and baryogenesis. Phys. Lett. B 736, 526-532 (2014). [arXiv:1405.4125]

85. C. Pallis, Linking Starobinsky-type inflation in no-scale supergravity to MSSM. JCAP 1404, 024 (2014). [arXiv:1312.3623]

86. W. Buchmuller, V. Domcke, K. Kamada, K. Schmitz, A minimal supersymmetric model of particle physics and the early universe. [arXiv:1309.7788]

87. M. Bose, M. Dine, P. Draper, Moduli or not. Phys. Rev. D 88, 023533 (2013). [arXiv:1305.1066]

88. C. Pallis, Q. Shafi, Update on minimal supersymmetric hybrid inflation in light of PLANCK. Phys. Lett. B 725, 327-333 (2013). [arXiv:1304.5202]

89. L. Wang, E. Pukartas, A. Mazumdar, Visible sector inflation and the right thermal history in light of Planck data. JCAP 1307, 019 (2013). [arXiv:1303.5351]

90. C. Boehm, J. Da Silva, A. Mazumdar, E. Pukartas, Probing the supersymmetric inflaton and dark matter link via the CMB, LHC and XENON1T experiments. Phys. Rev. D 87, 023529 (2013). [arXiv:1205.2815]

91. A. Kobakhidze, M. Talia, L. Wu, Probing the MSSM explanation of the muon g-2 anomaly in dark matter experiments and at a 100 TeV pp collider. Phys. Rev. D 95, 055023 (2017). [arXiv:1608.03641]

92. P. Athron, M. Bach et al., GM2Calc: precise MSSM prediction for (g-2) of the muon. Eur. Phys. J. C 76, 62 (2016). [arXiv:1510.08071]

93. J.-H. Park, M. Bach, D. Stöckinger, H. Stöckinger-Kim, Large muon ($g-2$) from TeV-scale MSSM with infinite $\tan\beta$. PoS PLANCK2015, 099 (2015). [arXiv:1510.04263]

94. J. Hisano, D. Kobayashi, W. Kuramoto, T. Kuwahara, Nucleon electric dipole moments in high-scale supersymmetric models. JHEP 11, 085 (2015). [arXiv:1507.05836]

95. K. Harigaya, T.T. Yanagida, N. Yokozaki, Muon $g-2$ in focus point SUSY. Phys. Rev. D 92, 035011 (2015). [arXiv:1505.01987]

96. M. Bach, J.-H. Park, D. Stöckinger, H. Stöckinger-Kim, Large muon ($g-2$) with TeV-scale SUSY masses for $\tan\beta \to \infty$. JHEP 10, 026 (2015). [arXiv:1504.05500]

97. J. Chakrabortty, A. Choudhury, S. Mondal, Non-universal Gaugino mass models under the lamppost of muon ($g-2$). JHEP 07, 038 (2015). [arXiv:1503.08703]

98. K. Kowalska, L. Roszkowski, E.M. Sessolo, A.J. Williams, GUT-inspired SUSY and the muon $g-2$ anomaly: prospects for LHC 14 TeV. JHEP 06, 020 (2015). [arXiv:1503.08219]

99. A. Aboubrahim, T. Ibrahim, P. Nath, Neutron electric dipole moment and probe of PeV scale physics. Phys. Rev. D 91, 095017 (2015). [arXiv:1503.06850]

100. K. De Causmaecker, B. Fuks et al., General squark flavour mixing: constraints, phenomenology and benchmarks. JHEP 11, 125 (2015). [arXiv:1509.05414]

101. J.A. Evans, D. Shih, A. Thalapillil, Chiral flavor violation from extended gauge mediation. JHEP 07, 040 (2015). [arXiv:1504.00930]

102. B. Dutta, Y. Mimura, Enhancement of $Br(B_d \to \mu^+\mu^-)/Br(B_s \to \mu^+\mu^-)$ in supersymmetric unified models. Phys. Rev. D 91, 095011 (2015). [arXiv:1501.02044]

103. A.H. Chamseddine, R.L. Arnowitt, P. Nath, Locally supersymmetric grand unification. Phys. Rev. Lett. 49, 970 (1982)

104. R. Barbieri, S. Ferrara, C.A. Savoy, Gauge models with spontaneously broken local supersymmetry. Phys. Lett. B 119, 343 (1982)

105. L.E. Ibanez, Locally supersymmetric SU(5) grand unification. Phys. Lett. B 118, 73-78 (1982)

106. L.J. Hall, J.D. Lykken, S. Weinberg, Supergravity as the messenger of supersymmetry breaking. Phys. Rev. D 27, 2359-2378 (1983)

107. J.R. Ellis, D.V. Nanopoulos, K. Tamvakis, Grand unification in simple supergravity. Phys. Lett. B 121, 123-129 (1983)

108. L. Alvarez-Gaume, J. Polchinski, M.B. Wise, Minimal low-energy supergravity. Nucl. Phys. B 221, 495 (1983)

109. H.P. Nilles, Supersymmetry, supergravity and particle physics. Phys. Rep. 110, 1-162 (1984)

110. D. Matalliotakis, H.P. Nilles, Implications of nonuniversality of soft terms in supersymmetric grand unified theories. Nucl. Phys. B 435, 115-128 (1995). [arXiv:hep-ph/9407251]

111. M. Olechowski, S. Pokorski, Electroweak symmetry breaking with nonuniversal scalar soft terms and large tan beta solutions. Phys. Lett. B 344, 201-210 (1995). [arXiv:hep-ph/9407404]

112. V. Berezinsky, A. Bottino et al., Neutralino dark matter in supersymmetric models with nonuniversal scalar mass terms. Astropart. Phys. 5, 1-26 (1996). [arXiv:hep-ph/9508249]

113. M. Drees, M.M. Nojiri, D.P. Roy, Y. Yamada, Light Higgsino dark matter. Phys. Rev. D 56, 276-290 (1997). [arXiv:hep-ph/9701219 [Erratum: Phys. Rev. D 64, 039901 (2001)]]

114. P. Nath, R.L. Arnowitt, Nonuniversal soft SUSY breaking and dark matter. Phys. Rev. D 56, 2820-2832 (1997). [arXiv:hep-ph/9701301]

115. C. Han, K.-I. Hikasa, L. Wu, J.M. Yang, Y. Zhang, Status of CMSSM in light of current LHC Run-2 and LUX data. [arXiv:1612.02296]

116. P. Bechtle et al., How alive is constrained SUSY really?, in Proceedings, 37th International Conference on High Energy Physics (ICHEP 2014): Valencia, Spain, July 2-9, 2014 (2016). [arXiv:1410.6035]

117. L. Roszkowski, E.M. Sessolo, A.J. Williams, What next for the CMSSM and the NUHM: Improved prospects for superpartner and dark matter detection. JHEP 08, 067 (2014). [arXiv:1405.4289]

118. A. Fowlie, M. Raidal, Prospects for constrained supersymmetry at $\sqrt{s} = 33$ TeV and $\sqrt{s} = 100$ TeV proton-proton super-colliders. Eur. Phys. J. C 74, 2948 (2014). [arXiv:1402.5419]
119. O. Buchmueller et al., The CMSSM and NUHM1 after LHC Run 1. Eur. Phys. J. C 74, 2922 (2014). [arXiv:1312.5250](https://arxiv.org/abs/1312.5250)

120. O. Buchmueller et al., Implications of improved Higgs mass calculations for supersymmetric models. Eur. Phys. J. C 74, 2809 (2014). [arXiv:1312.5233](https://arxiv.org/abs/1312.5233)

121. P. Bechtle et al., Constrained supersymmetry after the Higgs Boson discovery: A global analysis with Fittino. PoS EPS-HEP2013, 313 (2013). [arXiv:1310.3045](https://arxiv.org/abs/1310.3045)

122. S. Henrot-Versillé, R. Lafaye et al., Constraining supersymmetry using the relic density and the Higgs boson. Phys. Rev. D 89, 055017 (2014). [arXiv:1309.6958](https://arxiv.org/abs/1309.6958)

123. N. Bornhauser, M. Drees, Determination of the CMSSM parameters using neural networks. Phys. Rev. D 88, 075016 (2013). [arXiv:1307.3383](https://arxiv.org/abs/1307.3383)

124. S. Akula, P. Nath, Gluino-driven radiative breaking, Higgs boson mass, muon g-2, and the Higgs diphoton decay in supergravity unification. Phys. Rev. D 87, 115022 (2013). [arXiv:1304.5526](https://arxiv.org/abs/1304.5526)

125. M. Citron, J. Ellis et al., End of the CMSSM coannihilation strip is nigh. Phys. Rev. D 87, 036012 (2013). [arXiv:1212.2886](https://arxiv.org/abs/1212.2886)

126. C. Strege, G. Bertone et al., Global fits of the cMSSM and NUHM including the LHC Higgs discovery and new XENON100 constraints. JCAP 4, 13 (2013). [arXiv:1212.2636](https://arxiv.org/abs/1212.2636)

127. A.V. Gladyshev, D.I. Kazakov, Is (Low Energy) SUSY still alive?. in Proceedings, 2012 European School of High-Energy Physics (ESHEP 2012): La Pommeraye, Anjou, France, June 06-19, 2012 (2014), pp. 107–159. [arXiv:1212.2548](https://arxiv.org/abs/1212.2548)

128. K. Kowalska, S. Munir et al., Constrained next-to-minimal supersymmetric standard model with a 126 GeV Higgs boson: A global analysis. Phys. Rev. D 87, 115010 (2013). [arXiv:1211.1693](https://arxiv.org/abs/1211.1693)

129. O. Buchmueller, R. Cavanaugh et al., The CMSSM and NUHM1 in light of 7 TeV LHC, $B_s \to \mu^+\mu^-$ and XENON100 data. Eur. Phys. J. C 72, 2243 (2012). [arXiv:1207.7315](https://arxiv.org/abs/1207.7315)

130. S. Akula, P. Nath, G. Peim, Implications of the Higgs Boson discovery for mSUGRA. Phys. Lett. B 717, 188–192 (2012). [arXiv:1207.1839](https://arxiv.org/abs/1207.1839)

131. H. Baer, V. Barger, A. Lessa, X. Tata, Discovery potential for SUSY at a high luminosity upgrade of LHC14. Phys. Rev. D 86, 117701 (2012). [arXiv:1207.4846](https://arxiv.org/abs/1207.4846)

132. L. Roszkowski, E.M. Sessolo, Y.-L.S. Tsai, Bayesian implications of current LHC supersymmetry and dark matter detection searches for the constrained MSSM. Phys. Rev. D 86, 095005 (2012). [arXiv:1202.1503](https://arxiv.org/abs/1202.1503)

133. C. Strege, G. Bertone et al., Updated global fits of the cMSSM including the latest LHC SUSY and Higgs searches and XENON100 data. JCAP 3, 30 (2012). [arXiv:1112.4192](https://arxiv.org/abs/1112.4192)

134. P. Bechtle, T. Bringmann et al., Constrained supersymmetry after two years of LHC data: A global view with Fittino. JHEP 6, 98 (2012). [arXiv:1204.4199](https://arxiv.org/abs/1204.4199)

135. O. Buchmueller, R. Cavanaugh et al., Supersymmetry in light of 1/fb of LHC data. Eur. Phys. J. C 72, 1878 (2012). [arXiv:1110.3568](https://arxiv.org/abs/1110.3568)

136. A. Fowlie, A. Kalinowski, M. Kazana, L. Roszkowski, Y.L.S. Tsai, Bayesian implications of current LHC and XENON100 search limits for the constrained MSSM. Phys. Rev. D 85, 075012 (2012). [arXiv:1111.6098](https://arxiv.org/abs/1111.6098)

137. P. Bechtle, K. Desch, M. Uhlenbrock, P. Wienemann, Constraining SUSY models with Fittino using measurements before, with and beyond the LHC. Eur. Phys. J. C 66, 215–259 (2010). [arXiv:0907.2589](https://arxiv.org/abs/0907.2589)

138. R. Trotta, F. Feroz, M. Hobson, L. Roszkowski, R. Ruiz de Austri, The impact of priors and observables on parameter inferences in the constrained MSSM. JHEP 12, 24 (2008). [arXiv:0809.3792](https://arxiv.org/abs/0809.3792)

139. R. Ruiz de Austri, R. Trotta, L. Roszkowski, A Markov chain Monte Carlo analysis of CMSSM. JHEP 5, 2 (2006). [arXiv:hep-ph/0602028](https://arxiv.org/abs/0602.028)

140. B.C. Allanach, C.G. Lester, Multidimensional mSUGRA likelihood maps. Phys. Rev. D 73, 015013 (2006). [arXiv:hep-ph/0507283](https://arxiv.org/abs/0507.283)

141. P. Bechtle, K. Desch, P. Wienemann, Fittino, a program for determining MSSM parameters from collider observables using an iterative method. Comput. Phys. Commun. 174, 47–70 (2006). [arXiv:hep-ph/0412012](https://arxiv.org/abs/0412.012)

142. E.A. Baltz, P. Gondolo, Markov chain Monte Carlo exploration of minimal supergravity with implications for dark matter. JHEP 10, 52 (2004). [arXiv:hep-ph/0407039](https://arxiv.org/abs/0407.039)

143. R. Lafaye, T. Plehn, D. Zerwas, SFITTER: SUSY Parameter Analysis at LHC and LC. [arXiv:hep-ph/0404282](https://arxiv.org/abs/0404.282)

144. H. Baer, V. Barger, M. Savoy, Upper bounds on sparticle masses from naturalness or how to disprove weak scale supersymmetry. Phys. Rev. D 93, 035016 (2016). [arXiv:1509.02929](https://arxiv.org/abs/1509.02929)

145. E.A. Bagnaschi, O. Buchmueller et al., Supersymmetric dark matter after LHC run 1. Eur. Phys. J. C 75, 500 (2015). [arXiv:1508.01173](https://arxiv.org/abs/1508.01173)

146. H. Baer, V. Barger, M. Savoy, Supergravity gauge theories strike back: There is no crisis for SUSY but a new collider may be required for discovery. Phys. Scr. 90, 068003 (2015). [arXiv:1502.04127](https://arxiv.org/abs/1502.04127)

147. J.L. Evans, M. Ibe, K.A. Olive, T.T. Yanagida, Light Higgsinos in pure gravity mediation. Phys. Rev. D 91, 055008 (2015). [arXiv:1412.3403](https://arxiv.org/abs/1412.3403)

148. O. Buchmueller et al., The NUHM2 after LHC Run 1. Eur. Phys. J. C 74, 3212 (2014). [arXiv:1408.4060](https://arxiv.org/abs/1408.4060)

149. D.A. Demir, C.S. Ün, Stop on top: SUSY parameter regions, fine-tuning constraints. Phys. Rev. D 90, 095015 (2014). [arXiv:1407.1481](https://arxiv.org/abs/1407.1481)

150. H. Baer, V. Barger, D. Mickelson, M. Padelfeke-Kirkland, SUSY models under siege: LHC constraints and electroweak fine-tuning. Phys. Rev. D 89, 115019 (2014). [arXiv:1404.2277](https://arxiv.org/abs/1404.2277)

151. H. Baer, V. Barger, D. Mickelson, How conventional measures overestimate electroweak fine-tuning in supersymmetric theory. Phys. Rev. D 88, 095013 (2013). [arXiv:1309.2984](https://arxiv.org/abs/1309.2984)

152. S. Bhattacharya, S. Chakdar, K. Ghosh, S. Nandi, Non-universal SUGRA at LHC: Prospects and discovery potential. Phys. Rev. D 89, 015004 (2014). [arXiv:1309.0036](https://arxiv.org/abs/1309.0036)

153. H. Baer, J. List, Post LHC8 SUSY benchmark points for ILC physics. Phys. Rev. D 88, 055004 (2013). [arXiv:1307.0782](https://arxiv.org/abs/1307.0782)

154. G. Bhattacharyya, A. Kundu, T.S. Ray, Minimal supersymmetry confronts $R_b$, $A_{FB}^b$ and $m_h$. J. Phys. G 41, 035002 (2014). [arXiv:1306.0344](https://arxiv.org/abs/1306.0344)

155. GAMBIT Collaboration: P. Athron, C. Balazs et al., GAMBIT: The Global and Modular Beyond-the-Standard-Model Inference Tool. [arXiv:1705.07908](https://arxiv.org/abs/1705.07908)

156. GAMBIT Dark Matter Workgroup: T. Bringmann, J. Conrad et al. DarkBit: A GAMBIT module for computing dark matter observables and likelihoods. [arXiv:1705.07920](https://arxiv.org/abs/1705.07920)

157. GAMBIT Collider Workgroup: C. Balázs, A. Buckley et al., ColliderBit: a GAMBIT module for the calculation of high-energy collider observables and likelihoods. [arXiv:1705.07919](https://arxiv.org/abs/1705.07919)

158. GAMBIT Flavour Workgroup: F.U. Bernlochner, M. Chrzaszcz et al., FlavBit: A GAMBIT module for computing flavour observables and likelihoods. [arXiv:1705.07933](https://arxiv.org/abs/1705.07933)

159. GAMBIT Models Workgroup: P. Athron, C. Balázs et al., SpecBit, DecayBit and PrecisionBit: GAMBIT modules for computing mass spectra, particle decay rates and precision observables. [arXiv:1705.07936](https://arxiv.org/abs/1705.07936)

160. GAMBIT Scanner Workgroup: G.D. Martinez, J. McKay et al., Comparison of statistical sampling methods with ScannerBit, the GAMBIT scanning module. [arXiv:1705.07959](https://arxiv.org/abs/1705.07959)

161. GAMBIT Collaboration: P. Athron, C. Balázs et al., Status of the scalar singlet dark matter model. [arXiv:1705.07931](https://arxiv.org/abs/1705.07931)
162. GAMBIT Collaboration: P. Athron, C. Balázs et al., *A global fit of the MSSM with GAMBIT*. arXiv:1705.07917

163. GAMBIT Collaboration, Supplementary Data: Global fits of GUT-scale SUSY models with GAMBIT (2017). doi:https://doi.org/10.5281/zenodo.801642. arXiv:1705.07935

164. T. Falk, K.A. Olive, M. Srednicki, Heavy sneutrinos as dark matter. Phys. Lett. B 339, 248-251 (1994). arXiv:hep-ph/9409270

165. C. Arina, N. Fornengo, Sneutrino cold dark matter, a new analysis: Relic abundance and detection rates. JHEP 11, 029 (2007). arXiv:0709.4477

166. L. Roszkowski, R. Ruiz de Austri, K.-Y. Choi, Gravitino dark matter in the CMSSM and implications for leptogenesis and the LHC. JHEP 08, 080 (2005). arXiv:hep-ph/0408227

167. B.C. Allanach et al., SUSY Les Houches Accord 2. Comput. Phys. Commun. 180, 8-25 (2009). arXiv:0801.0045

168. ATLAS, CDF, CMS and D0 Collaborations, First combination of Tevatron and LHC measurements of the top-quark mass. arXiv:1403.4427

169. Particle Data Group: K.A. Olive et al., Review of particle physics. Chin. Phys. C 38, 090001 (2014)

170. J.F. Navarro, C.S. Frenk, S.D.M. White, The structure of cold dark matter halos. ApJ 462, 563-575 (1996). arXiv:astro-ph/9508025

171. R. Schoenrich, J. Binney, W. Dehnen, Local Kinematics and the Local Standard of Rest. MNRAS 403, 1829 (2010). arXiv:0912.3693

172. K. Freese, M. Lisanti, C. Savage, Annual modulation of dark matter: A review. Rev. Mod. Phys. 85, 1561-1581 (2013). arXiv:1209.3339

173. M.J. Reid et al., Trigonometric parallaxes of massive star forming regions: VI. Galactic structure, fundamental parameters and non-circular motions. ApJ 700, 137-148 (2009). arXiv:0902.3913

174. J. Bovy, D.W. Hogg, H.-W. Rix, Galactic masers and the Milky Way circular velocity. ApJ 704, 1704-1709 (2009). arXiv:0907.5423

175. M.C. Smith et al., The RAVE Survey: Constraining the local galactic escape speed. MNRAS 379, 755-772 (2007). arXiv:astro-ph/0611671

176. Y. Akrami, C. Savage, P. Scott, J. Conrad, J. Edsjö, How well will ton-scale dark matter direct detection experiments constrain minimal supersymmetry? JCAP 1104, 012 (2011). arXiv:1011.4318

177. R.D. Young, Strange quark content of the nucleon and dark matter searches. PoS LATTICE2012, 014 (2012). arXiv:1301.1765

178. B.C. Allanach, K. Cranmer, C.G. Lester, A.M. Weber, Natural priors, CMSSM fits and LHC weather forecasts. JHEP 08, 023 (2007). arXiv:0705.0487

179. M.E. Cabrera, J.A. Casas, R. Ruiz de Austri, Bayesian approach and Naturalness in MSSM analyses for the LHC. JHEP 03, 075 (2009). arXiv:0812.0536

180. F. Feroz, M.P. Hobson, M. Bridges, MULTINEST: an efficient and robust Bayesian inference tool for cosmology and particle physics. MNRAS 398, 1601-1614 (2009). arXiv:0809.3437

181. R. Storn, K. Price, Differential evolution: A simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim. 11, 341-359 (1997)

182. J. Brest, S. Greiner, B. Boskovic, M. Mernik, V. Zumer, Self-adapting control parameters in differential evolution: A comparative study on numerical benchmark problems. Evol. Comput. IEEE Trans. 10, 646-657 (2006)

183. J. Skilling, Nested sampling, in *American Institute of Physics Conference Series*, ed. by R. Fischer, R. Preuss, U.V. Toussaint, vol. 735 (2004), pp. 395-405

184. Y. Akrami, P. Scott, J. Edsjö, J. Conrad, L. Bergström, A profile likelihood analysis of the Constrained MSSM with genetic algorithms. JHEP 4, 57 (2010). arXiv:0910.3950

185. F. Feroz, K. Cranmer, M. Hobson, R. Ruiz de Austri, R. Trotta, Challenges of profile likelihood evaluation in multi-dimensional SUSY scans. JHEP 6, 42 (2011). arXiv:1101.3296

186. P. Scott, Pippi—painless parsing, post-processing and plotting of posterior and likelihood samples. Eur. Phys. J. Plus 127, 138 (2012). arXiv:1206.2245

187. R. Catena, P. Ullio, A novel determination of the local dark matter density. JCAP 1008, 004 (2010). arXiv:0907.0018

188. M. Pato, F. Iocco, G. Bertone, Dynamical constraints on the dark matter distribution in the Milky Way. JCAP 1512, 001 (2015). arXiv:1504.06324

189. J.I. Read, The local dark matter density. J. Phys. G 41, 063101 (2014). arXiv:1404.1938

190. H.-W. Lin, Lattice QCD for precision nucleon matrix elements. arXiv:1112.2435

191. P. Athron, J.-H. Park, D. Stöckinger, A. Voigt, Flexible SUSY—A spectrum generator generator for supersymmetric models. Comput. Phys. Commun. 190, 139-172 (2015). arXiv:1406.2319

192. F. Staub, SARAH. arXiv:0806.0538

193. F. Staub, Automatic calculation of supersymmetric renormalization group equations and self energies. Comput. Phys. Commun. 182, 808-833 (2011). arXiv:1002.0840

194. B.C. Allanach, SOFTSUSY: a program for calculating supersymmetric spectra. Comput. Phys. Commun. 143, 305-331 (2002). arXiv:hep-ph/0104145

195. B.C. Allanach, P. Athron, L.C. Tunstall, A. Voigt, A.G. Williams, Next-to-Minimal SOFTSUSY. Comput. Phys. Commun. 185, 2322-2339 (2014). arXiv:1311.7659

196. G. Degrassi, P. Slavich, F. Zwirner, On the neutral Higgs boson masses in the MSSM for arbitrary stop mixing. Nucl. Phys. B 611, 403-422 (2001). arXiv:hep-ph/0105096

197. A. Brignole, G. Degrassi, P. Slavich, F. Zwirner, On the O(alpha(t)**2) two loop corrections to the neutral Higgs boson masses in the MSSM. Nucl. Phys. B 631, 195-218 (2002). arXiv:hep-ph/0112177

198. A. Brignole, G. Degrassi, P. Slavich, F. Zwirner, On the two loop sbottom corrections to the neutral Higgs boson masses in the MSSM. Nucl. Phys. B 643, 79-92 (2002). arXiv:hep-ph/0206101

199. A. Dedes, G. Degrassi, P. Slavich, On the two loop Yukawa corrections to the MSSM Higgs boson masses at large tan beta. Nucl. Phys. B 672, 144-162 (2003). arXiv:hep-ph/0305127

200. P. Draper, G. Lee, C.E.M. Wagner, Precise estimates of the Higgs mass in heavy supersymmetry. Phys. Rev. D 89, 055023 (2014). arXiv:1312.5743

201. E. Bagnaschi, G.F. Giudice, P. Slavich, A. Strumia, Higgs mass and unnatural supersymmetry. JHEP 09, 092 (2014). arXiv:1407.4081

202. J.P. Vega, G. Villadoro, SusyHD: Higgs mass determination in supersymmetry. JHEP 07, 159 (2015). arXiv:1504.05200

203. G. Lee, C.E.M. Wagner, Higgs bosons in heavy supersymmetry with an intermediate $m_A$. Phys. Rev. D 92, 075032 (2015). arXiv:1508.00576

204. Planck Collaboration, P.A.R. Ade, et al., Planck 2015 results. XIII. Cosmological parameters. A&A 594, A13 (2016). arXiv:1502.01589

205. J. Edsjö, P. Gondolo, Neutralino relic density including coannihilations. Phys. Rev. D 56, 1879-1894 (1997). arXiv:hep-ph/9704361

206. J. Edsjö, M. Schelke, P. Ullio, P. Gondolo, Accurate relic densities with neutralino, chargino and sfermion coannihilations in mSUGRA. JCAP 0304, 001 (2003). arXiv:hep-ph/0301106

207. P. Gondolo, J. Edsjö et al., DarkSUSY: computing supersymmetric dark matter properties numerically. JCAP 7, 8 (2004). arXiv:astro-ph/0406204

208. J. Harz, B. Herrmann, M. Klasen, K. Kovařík, M. Meinecke, SUSY-QCD corrections to stop annihilation into electroweak final