Robustness of dark matter constraints and interplay with collider searches for New Physics

We study the implications of dark matter searches, together with collider constraints, on the phenomenological MSSM with neutralino dark matter and focus on the consequences of the related uncertainties in some detail. We consider, inter alia, the latest results from AMS-02, Fermi-LAT and XENON1T. In particular, we examine the impact of the choice of the dark matter halo profile, as well as the propagation model for cosmic rays, for dark matter indirect detection and show that the constraints on the MSSM differ by one to two orders of magnitude depending on the astrophysical hypotheses. On the other hand, our limited knowledge of the local relic density in the vicinity of the Earth and the velocity of Earth in the dark matter halo leads to a factor 3 in the exclusion limits obtained by direct detection experiments. We identified the astrophysical models leading to the most conservative and the most stringent constraints and for each case studied the complementarities with the latest LHC measurements and limits from Higgs, SUSY and monojet searches. We show that combining all data from dark matter searches and colliders, a large fraction of our supersymmetric sample could be probed. Whereas the direct detection constraints are rather robust under the astrophysical assumptions, the uncertainties related to indirect detection can have an important impact on the number of the excluded points.

We study the implications of dark matter searches, together with collider constraints, on the phenomenological MSSM with neutralino dark matter and focus on the consequences of the related uncertainties in some detail. We consider, inter alia, the latest results from AMS-02, Fermi-LAT and XENON1T. In particular, we examine the impact of the choice of the dark matter halo profile, as well as the propagation model for cosmic rays, for dark matter indirect detection and show that the constraints on the MSSM differ by one to two orders of magnitude depending on the astrophysical hypotheses. On the other hand, our limited knowledge of the local relic density in the vicinity of the Earth and the velocity of Earth in the dark matter halo leads to a factor 3 in the exclusion limits obtained by direct detection experiments. We identified the astrophysical models leading to the most conservative and the most stringent constraints and for each case studied the complementarities with the latest LHC measurements and limits from Higgs, SUSY and monojet searches. We show that combining all data from dark matter searches and colliders, a large fraction of our supersymmetric sample could be probed. Whereas the direct detection constraints are rather robust under the astrophysical assumptions, the uncertainties related to indirect detection can have an important impact on the number of the excluded points.

I. INTRODUCTION
The Large Hadron Collider (LHC) is currently operating at 13 TeV, with the main purpose of searching for new particles and new phenomena, in order to establish the first direct proofs of New Physics Beyond the Standard Model (BSM). The lack of a stable, neutral, massive and very weakly interacting dark matter particle (WIMP) in the Standard Model (SM) constitutes one of the main motivations for search for New Physics. The dark matter (DM) paradigm emerged in astrophysics and cosmology, through the observation of galaxy rotation curves and galaxy clusters, giving birth to the concept of dark matter haloes around galaxies. Since then, numerous observations have indicated that dark matter is probably cold, i.e. has small velocities. In addition, collisions of clusters such as the Bullet Cluster reveal that dark matter can be separated from the baryons [1,2]. Furthermore, Cosmic Microwave Background (CMB) observations, and in particular those of the Planck Collaboration [3], provide precise measurements of the cold dark matter density, which can be distinguished from the subdominant baryon and hot dark matter densities.
One of the most common hypotheses is that cold dark matter is made of new particles still undiscovered at colliders. Two types of experiments have been designed in order to discover dark matter particles. First, direct detection experiments are based on the assumption that dark matter particles interact weakly with matter, and in particular with nucleons. In view of the dark matter density in galactic haloes, it is assumed that a large number of such particles would cross the Earth at any time. Therefore, the design of direct detection experiments is basically to gather a large amount of crystals or gases in big tanks in order to measure the recoil energy of the nuclei when dark matter scatters with the nucleons. Second, dark matter indirect detection experiments aim at detecting the annihilation or decay products of dark matter particles. In particular, the density of dark matter should be larger in the galactic center making the annihilation of dark matter particles more probable. Dark matter annihilations can produce gamma rays, which can be observed for example with Cherenkov telescopes, or other SM particles which can populate the antiproton spectrum, that is observed in particle detectors. So far, none of the dark matter detection experiments have been able to find a solid evidence for dark matter.
In parallel, the LHC continues the quest for New Physics. In Summer 2012, the announcement of the discovery of the Higgs boson with a mass of 125 GeV during the 8 TeV run [4,5] was a final step towards the completion of the Standard Model spectrum. Thereafter, the measurements of the properties of the Higgs boson have established its compatibility with the SM predictions. Further, no new particle has been discovered yet, in spite of the higher energy of 13 TeV available in run 2. Searches for supersymmetric or BSM particles, or exotic phenomena such as monojets, are ongoing, and strong limits on BSM models are obtained. Monojet searches are generally considered as dark matter searches at the LHC, since they aim at finding evidence for missing energy in the final states through the presence of an initial state radiated hard jet. This missing energy would correspond to invisible particles leaving no energy in the detectors, which may be the dark matter constituents.
Searches for supersymmetry (SUSY) are still the main focus of the ATLAS and CMS experiments. In particular, the minimal supersymmetric extension of the Standard Model (MSSM) constitutes an excellent playground to design new physics searches or study dark matter, with however a limitation for systematic studies due to the large number of free parameters. If R-parity is conserved, supersymmetric particles can only interact in pair with SM particles, so that the lightest supersymmetric particle (LSP) is stable and can constitute dark matter, provided it is neutral and weakly interacting. The lightest neutralino is generally considered to be an adequate dark matter candidate.
The phenomenological MSSM (pMSSM) [6], with its 19 parameters, is a good compromise as it is the most general R-parity and CP-conserving MSSM scenario respecting minimal flavour violation, and has a manageable number of parameters to allow for systematic studies. In this paper, we will study dark matter and collider constraints within the pMSSM with 19 independent parameters, assuming that the lightest neutralino (χ 0 1 , simply labelled χ in the following) is the LSP. This scenario is general enough so that our main conclusions can hold in other supersymmetric scenarios. There have been several studies combining the LHC limits and dark matter constraints in the pMSSM, which either aim at determining the excluded regions or perform global fits in order to find the preferred parameter regions (see e.g. [7][8][9][10][11][12][13][14] for some recent studies). Instead, the focus here will be on the astrophysical and cosmological uncertainties that can affect the interpretation of the dark matter experiment results, and to show explicitly and quantify the impact of such uncertainties in the dark matter limits that are set on the MSSM parameters.
The rest of this paper is organised as follows. In Section II, we will review the theoretical framework of relic dark matter density, indirect and direct de-tections, and study the astrophysical and cosmological uncertainties that can affect them. In Section III, we describe the methods used for our analysis of the pMSSM, for the collider constraints, as well as for direct and indirect detections, and evaluate the general consequences of the choice of astrophysical and cosmological assumptions. In Section IV, we show the results in the pMSSM, considering the dark matter observables and their uncertainties, collider constraints, and the complementarity between dark matter and collider results. We will also briefly discuss the prospects for dark matter experiments. Finally, the conclusions will be given in Section V.

II. DARK MATTER SEARCHES AND UNCERTAINTIES
A. Relic density The dark matter abundance has been measured in the framework of the standard cosmological model, and the Planck Collaboration has provided a precise evaluation of the cold dark matter density [3]: Constraints on new physics scenarios which propose dark matter candidates can therefore be obtained by comparing the computed dark matter density to the Planck value. The standard assumption to compute the dark matter density is to consider that dark matter particles are thermal relics, i.e. were in thermal equilibrium in the early Universe and we observe today only the surviving part. A second assumption is that there is a single thermal relic candidate contributing to the dark matter density, which is generally the case in BSM scenarios where dark matter particles have to be stable, electrically neutral and very weakly interacting.
With these assumptions, the relic density can be obtained by considering that all the new physics particles were originally in thermal equilibrium. Then the expansion of the Universe, which lowers the temperature, eventually breaks this equilibrium, and the evolution of the number densities of all the new particles can be obtained using Boltzmann equations, in which the expansion of the Universe introduces a friction-like term and the collision terms include annihilations and co-annihilations of these new particles into SM particles. When the dark matter density is diluted enough so that the interactions become negligible, the relic density is frozen and becomes only diluted by the expansion of the Universe. A detailed description of the calculation can be found in [15,16].
Comparing the obtained relic density to the very precise dark matter measurement can lead to very strong constraints on new physics parameters.
In the MSSM with R-parity conservation, neutralino and gravitino constitute good dark matter candidates, provided they are the lightest supersymmetric particles. The gravitino is however produced nonthermally, and was considered recently in [17][18][19]. We therefore focus on the case of the lightest neutralino. The co-annihilations have in this case a very important role. If the lightest neutralino is mainly bino, it interacts weakly and the annihilation cross section is small, leading to a large relic density. To obtain the Planck limit, it is necessary to have coannihilations of the neutralino with SUSY particles which are close in mass in order to increase the effective (co-)annihilation cross section. Similarly, if the neutralino is mostly wino or Higgsino, it is accompanied by a chargino which is very close in mass, making the co-annihilations possible. Considering the Sommerfeld enhancement, a wino-like neutralino can have the correct relic density naturally for a mass of about 2.8 TeV in absence of other co-annihilations, and 1 TeV for a Higgsino [20,21]. Careful studies about the consequences of the Sommerfeld enhancement in the context of relic density and indirect detection in the MSSM can be found in [22][23][24][25][26][27].
Several assumptions can nevertheless limit the constraining power of the relic density constraint as will be discussed in the following.

Higher order corrections
The first uncertainties arise from the numerical calculations of the annihilation and co-annihilation crosssections. Whereas in the simplest cases the calculation of the relic density relies on a few decay channels, in the most compressed scenarios of the MSSM, more than 3000 channels can get involved, severely limiting the calculation speed of relic density. For this reason, the cross-sections are generally considered at tree-level. Yet, in individual channels, higher-order corrections can lead to 30% modification or more [28]. However, in most cases, the relic density calculated at tree-level differs by less than 10% from the one calculated at one-loop [29,30]. Therefore, in the general case, about 10% uncertainty can be associated to treelevel calculations of the relic density.

QCD equations of state
A second limitation comes from QCD equations of state. Indeed, computing the relic density requires the knowledge of the number of effective degrees of freedom of radiation, which lead to energy and entropy content of the Universe. While it was originally thought that the primordial plasma could be treated as an ideal gas above the QCD phase transition temperature, non-perturbative studies showed that at high temperature, the ideal gas approximation does not work, and different models for this plasma have been studied [31][32][33], leading to different sets of QCD equations of state. The consequences on the relic density are however rather mild and can modify it by a few percent.

Early Universe properties
In the usual calculation of relic density, the expansion of the Universe is considered to be dominated purely by the radiation density. This hypothesis can however be falsified in many extensions of the standard model of cosmology [34][35][36][37][38]. Similarly, entropy injection or non-thermal production of dark matter particles can modify the relic density [39][40][41][42][43]. These modifications of the standard model of cosmology can result in a change of the relic density by orders of magnitude, but are more likely to increase it. As a consequence, the uncertainties due to these effects are completely dominating the relic density calculation over the previous uncertainties.
To be conservative, we add to the Planck measurement error a theoretical uncertainty of 10%, and consider the 3.5σ interval 0.0772 < Ωh 2 < 0.1604 .
Moreover, since a modification of the cosmological standard scenario can result in a large increase of the relic density, the lower dark matter density limit can be disregarded.

B. Indirect detection
Dark matter particles hosted in galaxies are supposed to annihilate into SM particles to yield, after hadronisation and decay, nuclei, electrons, photons and neutrinos. The emissivity of one of these particles i injected by the annihilation of two DM particles is where σv is the thermal average annihilating cross section, B j the branching ratio of the annihilation channel j, dN j i /dE i the multiplicity of the particle i and η is equal to 1/2 (1/4) for a Majorana (Dirac) type particle. The density distribution ρ of dark matter particles is discussed in Sec. II B 1. Indirect detection experiments try to find an excess of those messengers on top of their astrophysical background. Even in absence of signals, these experiments provide useful information about the dark matter nature.
Antiparticle cosmic rays are regarded with great interest. Indeed, their astrophysical background is composed of secondary particles i.e. particles produced by the interaction of primary cosmic rays (mostly proton and helium nuclei) on the interstellar medium (mostly hydrogen and helium atoms). Hence, their background is feeble and relatively under control compared to other species. Antiprotons (p) are the most abundant antinuclei in cosmic rays that could be produced by dark matter and their spectral shape is distinguishable from the astrophysical background. For a dark matter mass larger than a few GeV, the flux of antiprotons features a cut off at the dark matter mass. The most accurate measurements of thep flux at the Earth was reported by the space-borne detectors PAMELA [44] and AMS-02 [45]. The discovery of an excess around 100 GeV was recently claimed [46,47].
However, both secondary and primary antiprotons suffer from theoretical uncertainties which make the significance of such an excess uncertain. On the one hand, the astrophysical background of secondary antiprotons is affected by the lack of knowledge of the antiproton production cross section from proton-proton and proton-helium interactions, leading to an uncertainty for the flux at the Earth of ∼ 50% (see for example [48,49]). On the other hand, the antiproton flux produced by DM is very sensitive to the DM profile, altering the primary antiproton flux at the Earth by up to a factor of 2-6 [50,51]. In addition, both secondary and primary antiprotons are sensitive to uncertainties related to their propagation throughout the Galaxy [52,53]. Astrophysical uncertainties on galactic properties as well as the production cross sections used for secondary cosmic rays are the main uncertainties for the determination of the propagation parameters. The total uncertainty for the secondary antiproton component was assessed in [54,55] to be up to a factor 3 at ∼ 100 GeV. Moreover, the total uncertainty for the DM signal was shown to be as large as a factor of about 20 in [50] and 50 in [51]. These results show how the constraints on the DM particle annihilation cross sections are sensitive to astrophysical and nuclear uncertainties. In the following, we reconsider these uncertainties using the most recent cosmic ray propagation results as well as the most recent galactic mass models, and study their consequences within the MSSM.
Cosmic ray positrons could also be produced by the annihilation of DM particles. Above a few GeV, the astrophysical background of positrons is not under control as for antiprotons. As a matter of fact, the positron excess reported by AMS-02 [56] could be explained by the presence of young and nearby pulsars. In addition, the lack of knowledge about these systems makes it difficult to distinguish this hypothesis from an exotic component to explain the data. Therefore, this channel is not much useful to derive constraints on dark matter properties when m DM is larger than a few GeV.
Compared to charged cosmic rays, gamma rays have the advantage of propagating straight ahead. This allows us to characterise the morphology of their sources and to observe regions where the dark matter particle density is expected to be large and to produce a sizeable flux. The Fermi-LAT space-borne telescope covers the GeV energy range whereas the ground based Cherenkov telescopes HESS [57], MAGIC [58], VER-ITAS [59] and HAWC [60] are sensitive to the TeV range. Since the density of dark matter particle is peaked in the center of the galaxy, the galactic center is one of the best targets to look for a dark matter signal. Nevertheless, this region hosts important astrophysical activities and it is difficult to estimate both the astrophysical background and foreground. Indeed, the gamma ray excess exhibited in the Fermi-LAT data [61][62][63][64][65][66] could be interpreted either by annihilating DM particles, or by the presence of millisecond pulsars or the remnants of the past activity of the supermassive black hole lying in the center of the galaxy [67,68]. On the other hand, dwarf spheroidal galaxies are considered as very interesting targets to look for a dark matter signal. Indeed, these systems are expected to i) be dominated in mass by a DM component ii) exhibit feeble stellar activities and then a low astrophysical background. Despite the fact that the dark matter distribution and concentration inside these objects is still under debate, they provide one of the best bounds on the average annihilating cross section σv .

Dark matter halo profiles
Dark matter particles are assumed to be isotropically distributed in a spherical halo around the galactic center.
The radial density profile of dark matter arising from cosmological simulations were parametrised by Navarro, Frenk and White  [72], for Einasto [73] and for Burkert [74] profiles.
(NFW) [69] as where r s is the radius at which the logarithmic slope of the profile is −2 and ρ s the dark mater density normalisation.
The Einasto profile on the other hand is defined as and provides a better agreement with the latest simulations [70] and does not suffer from the central divergence of (4). The star activity occurring in the inner galaxy could sweep dark matter particles from the inner region, resulting in a core profile as observed in many galaxies. Such profiles were introduced by Burkert et al. [71] with the parametrisation The parameters r s and ρ s as well as the distance of the Solar system to the galactic center are determined by dynamical observations of the Galaxy. We have used the values determined by [72] for NFW, by [73] for Einasto and by [74] for Burkert as reported in Table I.

Cosmic ray propagation
Following the work of [75] (and reference therein), we describe the galaxy using the two-zone model. The first zone, in which the interstellar medium is homogeneously distributed, represents the galactic disc of half-height h = 100 pc. Atomic densities are taken to be n H = 0.9 cm −3 and n He = 0.1 cm −3 . The disc is embedded inside a magnetic halo of half-height L lying between 1 and 15 kpc. Both zones share the radius R = 20 kpc. As cosmic rays travel across the galaxy, they are affected by many processes as a result of their interactions with the galactic magnetic field and the interstellar medium. The scattering of cosmic rays on the galactic magnetic field is modelled by a homogeneous and isotropic diffusion in space. The diffusion coefficient reads K(E) = βK 0 (R/1 GV) where β is the velocity of the particle and R the rigidity related to its momentum p and its charge q by R = p/q. Since the diffusion centers move with the Alfvèn waves velocity V a , the second-order Fermi mechanism applies and cosmic rays undergo a diffusive reacceleration. This process can be modelled by a diffusion in energy space with coefficient D(E) = (2/9)V 2 a E 2 β 4 /K(E). Moreover, cosmic rays can interact with the interstellar medium, leading to energy losses (including ionisation and Coulomb interaction) and their destruction at rates b and Γ, respectively. Finally, cosmic rays undergo the effect of the galactic wind produced by supernova remnant explosions in the galactic disc. We assume the galactic wind to be homogeneous and perpendicular to the galactic disc, with velocity V c = sign(z)V c e z . This process leads to the adiabatic cooling of cosmic rays, which enters as an additional term in the energy loss rate b.
Under a steady state and thin disc approximation, the density of cosmic rays per unit of space and energy ψ ≡ dN/d 3 xdE obeys the transport equation where Q represents the injection rate of cosmic rays in the galaxy.
The interstellar flux of cosmic rays at the Earth is given by Φ(E, ) = v/4π ψ(E, ). For more details on the resolution method of the transport equation, we refer the reader to [76,77]. In this way, a semi-analytical method was used in [77] to derive the benchmark Min, Med, and Max propagation models presented in Table II. The Med model corresponds to the best fit to the boron over carbon (B/C) ratio whereas the Min and Max sets of parameters define the lower and upper bounds for the primaryp flux, consistent with the B/C ratio. We emphasise that recent papers, based on synchrotron radio emission [78][79][80][81], on cosmic ray positrons [82] as well as on gamma rays [83], find that the thin halo predicted by Min is disfavoured. Furthermore, in Ref. [84] it was pointed out that secondary positrons can be used to improve the determination of the propagation parameters since they do not undergo exactly the same prop- agation processes as for nuclei. Following this idea, in Ref. [85] the pinching method was used to compute properly the flux of secondary positrons below 10 GeV and derive stringent constraints on the propagation parameters, which rule out the models with L < 4 kpc at the 3σ level, including the Min benchmark propagation model. As a result, the Med model provides a conservative lower bound to the dark matter antiproton signal. The recent B/C data reported by AMS-02 and their future studies would result in an improved determination of the parameters of the propagation models. Finally, cosmic rays have to penetrate the heliosphere where they interact with the Solar wind and the Solar magnetic field. In this work, we use the forced-field approximation [86] parametrised by the Fisk potential φ F , to predict the flux of cosmic rays on the top of the atmosphere where they are measured by the space-borne detectors.

C. Direct detection
Direct dark matter searches aim at directly detecting wimps via tiny energy deposits when they scatter off target atomic nuclei in ultra-sensitive, low background detectors. No convincing dark matter signal has been detected so far, however, limits on the wimpnucleon cross-section are set by comparing the measured differential recoil rate per unit detector mass to the theoretical rate given by: where ρ 0 is the local wimp density, m DM is the wimp mass, M N the target nucleus mass, and v min = E th m N 2µ 2 , with E th the recoil energy threshold and µ = m N m DM m N +m DM the reduced mass of the wimp-nucleus system. f (v) is the wimp velocity distribution in the Earth's rest frame.
Usually, the wimp-nucleus cross-section is decomposed in a spin-independent (SI) and a spin-dependent (SD) contributions in the zero momentum transfer limit: where the F functions are form factors describing how the wimp interferes with the nucleon structure of the nucleus and depend on the recoil energy. The σ χ−N are the wimp-nucleus cross sections at zero momentum transferred. The SI form factors are experimentally well known from the study of elastic electronic scattering on nuclei and are reasonably well approximated by Helm form factors [87], while the SD form factors are obtained from nuclear shell model calculations [88]. Those functions can be subject to uncertainties at high recoil energy, but their study is beyond the scope of this paper.
The SI wimp-nucleus cross-section is given by where Z and (A − Z) are the number of protons and neutrons in the nucleus and f p and f n are the effective SI wimp-proton and wimp-neutron couplings. In the calculation of experimental cross section limits, the approximation f p ≈ f n is commonly used. Under this assumption, the wimp-proton and wimp-neutron SI cross sections are equal (σ SI χ−p ≈ σ SI χ−n ) and the limits then concern a general wimp-nucleon cross section σ SI χ−nucleon . Since this cross section scales as A 2 , heavy target nuclei, such as argon and xenon, are favoured. In this context, the strongest limits on SI cross section, for m DM 10 GeV, are given by xenon target experiments, the leader being currently the XENON1T experiment [89]. 1 Argon target experiments, such as DarkSide-50 [91], give limits which are two orders of magnitude weaker.
The SD wimp-nucleus cross-section is given by where G F is the Fermi constant, J is the total spin of the nucleus, f p and f n are the effective SI wimpproton and wimp-neutron couplings and S p,n are the average spin contributions from the protons and neutrons in the nucleus. On the one hand, to set constraints on the SD wimp-neutron cross section, nuclei with an even number of protons and an odd-number of neutrons are needed. In xenon target experiments, for instance, the spin is carried by neutrons in neutronodd isotopes ( 129 Xe, 131 Xe). The best SD wimpneutron cross section limit is currently given by the LUX experiment [92]. On the other hand, to put constraints on the SD wimp-proton cross section, it is necessary to use nuclei with an odd number of protons. One of the strongest limits is given by the PICO-60 experiment, using C 3 F 8 target [93]. This limit is in competition with the one coming from the 79-string IceCube detector [94] that aims to detect a neutrino excess from the Sun. DM would be captured in the Sun through scattering on the hydrogen nuclei. The captured dark matter would then annihilate and produce neutrinos that would be detected by the IceCube experiment. Their limits on the cross section are calculated by assuming equilibrium between DM capture and its annihilation and depend strongly on the DM annihilation channel.
In the next few years, these limits are expected to be drastically lowered by experiments that will increase their total target mass and time of exposure, starting by XENONnT [95], LZ [96], and reaching the neutrino background with DARWIN [97] in ten years or so. Once the neutrino background is reached, if no dark matter particle is discovered by then, directional detection will be key to pursue the search for DM particles.

Global and local dark matter densities
All the experimental limits are calculated using the benchmark value ρ 0 = 0.3 GeV/cm 3 for the local DM density, but recent studies give a best fit value closer to 0.4 GeV/cm 3 [73,74,98]. The uncertainties on the local density value are still quite large, one of the main source residing in the knowledge of the baryon density in the galaxy. There may also be a discrepancy between the value calculated from the study of the motion of nearby stars and the one calculated from a global fit of stellar dynamics over the galaxy, assuming a spherical dark matter halo. In our study, we will consider that the local DM density lies between 0.2 and 0.6 GeV/cm 3 (see [99] for a complete review) and will choose three different values to test the impact of those uncertainties on the exclusions in our sample of points: ρ 0 = 0.2, 0.4 and 0.6 GeV/cm 3 .

Velocities
Customarily, an isotropic Maxwellian distribution is assumed for the WIMP velocity distribution f (v), with the galactic disk rotation velocity v rot being the most probable speed. It corresponds to the Standard Halo Model describing the dark matter halo as a nonrotating isothermal sphere [100,101]. The canonical value for v rot is 220 km/s but it is believed that it can range from 200 to 250 km/s [102][103][104]. This velocity distribution is truncated at the escape velocity v esc at which a wimp can escape the galaxy potential well. Its value is subject to large uncertainties, v esc = 500 − 600 km/s, with a benchmark value v esc = 544 km/s [105]. However, for wimp masses m DM > 10 GeV, v min is relatively low. The velocity distribution is then integrated over a large range of velocities and dR/dE R is not sensitive to the tail of the distribution. Thus, the uncertainties on v esc should not impact our analysis.
Other halo models have been proposed, such as the King Model which describes the finite size of the halo and the gravitational interaction with ordinary matter in a more realistic way [106,107] or such as triaxial halo models [108]. In this study, we will focus only on the uncertainties related to the Standard Halo Model, which is the most widely used in the literature.

A. MSSM Scans
We consider in this analysis the pMSSM, which is the most general R-parity and CP-conserving MSSM scenario with minimal flavour violation. It was shown in [109][110][111] that CP violation does not have important consequences on the dark matter sector after imposing the experimental constraints from the electric dipole moments and the Higgs coupling measurements, so that the results presented in the following will remain valid also for CP violating scenarios. We impose the lightest neutralino to be the lightest supersymmetric particle which constitutes dark matter, using the set-up presented in [112,113]. As the neutralino can be bino-like, wino-like, Higgsino-like or a mixed state, such a scenario allows for a large flexibility making our analysis of the astrophysical and cosmological uncertainties relatively general that can hold also in other dark matter models. We will not consider here the case of very light neutralinos that were studied in detail in [114][115][116][117].
The pMSSM points are generated with SOFTSUSY [118], with a flat random sampling using the ranges given in Table III for the 19 parameters. After checking the theoretical validity of each point, we impose to have a neutralino LSP as well as a light Higgs of mass between 122 and 128 GeV. We then apply different constraints from the dark matter and collider experiments, which are described below. We do not intend to perform a robust statistical analysis as performed for example by the GAMBIT Collaboration [13,119,120]. Instead, the constraints are imposed separately at the 2σ level, apart from for the Higgs sector where a likelihood analysis is used. As we do not aim at finding the preferred parameter regions, this choice will not affect the conclusions of our study.

B. Collider constraints
Collider searches are very important to constrain the supersymmetric parameter space, but are also relevant for dark matter, through their correlations with the dark matter detection experiments.
To the set of points in our analysis we apply constraints from LEP and Tevatron, from flavour physics, as well as the LHC constraints from the Higgs sector and supersymmetry and monojet direct searches.

Particle Limits
Conditions

LEP and Tevatron constraints
LEP and Tevatron have provided very robust constraints on the mass of the supersymmetric particles [121]. We apply to our set of points the limits summarised in Table IV. One can note that LEP also provides limits for the lightest neutralino, but they can be evaded in specific cases and since our analysis is focussed on dark matter and the lightest neutralino, we prefer not to apply it. The neutralino mass will nevertheless be constrained by the light Higgs signal strength measurements, which can lead to stronger limits than LEP [122][123][124][125][126][127].

Flavour constraints
Flavour constraints are complementary to the dark matter and collider searches. We consider here the three major decays, namely B s → µ + µ − , B → X s γ and B u → τ ν which capture the main constraints in the MSSM. B s → µ + µ − is a rare transition which has a very strong constraining power. Indeed the scalar and pseudoscalar contributions lead to enhancements of its branching fraction proportional to tan 6 β/M 4 A , strongly constraining the large tan β and small M A parameter regions [128][129][130][131][132]. The inclusive decay B → X s γ receives contributions from charged Higgstop and chargino-stop loops, which also restrict the Observable Experiment SM prediction BR(Bs → µ + µ − ) × 10 9 3.0 ± 0.65 [136] 3.54 ± 0.27 BR(B → Xsγ) × 10 4 3.32 ± 0.15 [137] 3.34 ± 0.22 BR(Bu → τ ντ ) × 10 4 1.06 ± 0.19 [137] 0.82 ± 0.29 charged Higgs, stop and chargino masses in the large tan β regions. It is worth noting that the pseudoscalar and charged Higgs masses are connected at tree level by the relation so that the pseudoscalar masses are also restricted. The third transition, B u → τ ν is a tree-level leptonic decay which can be mediated by a W -boson or a charged Higgs. It also restricts the small M H + and large tan β region. The value of the branching ratios of the three transitions is computed with SuperIso v3.7 [133][134][135], and we apply the constraints shown in Table V.

Higgs constraints
Higgs searches provide strong constraints on the MSSM and the dark matter sector. On the one side, the measurements of the 125 GeV Higgs signal strengths provide constraints on the pMSSM Higgs sector parameters. Indeed, the Higgs mass is given by where only the leading terms are given. Here M S = Mt 1 Mt 2 and X t = A t −µ cot β are the most relevant parameters to achieve a mass of 125 GeV.
The combined measurements of the Higgs mass by ATLAS and CMS from Run 1 gives [138] M hSM = 125.09 ± 0.21(stat.) ± 0.11(syst.) GeV . (14) However, the calculation of the Higgs mass in the MSSM is still subject to larger uncertainties (see for example [139]), and we adopt the constraint: The measurement of the Higgs couplings is also particularly constraining. Indeed, the couplings of the  Higgs bosons depend on the MSSM parameters as given at tree level in Table VI. α is the CP-even Higgs mixing angle: In the decoupling limit, which corresponds to M A M Z , the light Higgs couplings become SM-like. However, the couplings can receive higher-order corrections from the presence of supersymmetric particles, which can also lead to constraints on other MSSM parameters. LHC experiments have measured the signal strengths of different channels of the light Higgs boson, i.e. the product of the production cross sections times branching ratios. We use these measurements in our analyses, as given in Table VII. The decays h → W W, ZZ, bb, τ τ provide direct constraints on the couplings given in Table VI. On the other hand, h → γγ is a loop-level decay, in which the main contributions arise from top, stop, sbottom, chargino and charged Higgs loops [140]. Its measurement is therefore particularly important to constrain the MSSM. Additionally, limits on the Higgs decays into supersymmetric particles can be set, which become particularly relevant if the dark matter particle is lighter than half of the Higgs mass. Since these decays to new particles participate to the total decay width of the light Higgs boson, they lower the branching ratios to SM particles. The Higgs decay branching ratios and widths are computed using HDECAY v6.51 [141]. The production cross sections are calculated using Sushi 1.5.0 [142], VV2H v1.10 and V2HV v1.10 [143]. The constraints are obtained through a likelihood analysis using the experimental and theoretical correlations from [138] and [144], respectively. Constraints are applied at the 95% C.L.
Other relevant searches in the context of dark matter are searches for heavier Higgs bosons [122,[166][167][168][169][170]. As can be seen from Table VI, in the limit when M A is large, the light Higgs couplings are SM-like, and compatible with the current data. The heavier states are therefore expected to be heavy. Nevertheless, the couplings of the H/A to the b quarks and τ leptons are enhanced by tan β, so that it is possible to set strong limits in the small M A and large tan β region when searching for (pp)bb → H/A → τ τ . We use the results of CMS with 12.9 fb −1 [171], and assess the exclusion by comparing the calculated cross section times branching ratio with the published tables. We note that it is sensitive to the same region which is probed by the branching ratio of B s → µ + µ − .

LHC direct search constraints
Direct searches from supersymmetric particles at the LHC provide amongst the most important constraints on the MSSM parameter space. We consider in our study the LHC searches presented in Table VIII. Even if this list in not exhaustive, the most relevant searches for our study are considered, i.e. the channels with the highest sensitivity which are rather uncorrelated.
The SUSY direct searches correspond to final states with at least two SM particles and a large missing energy, carried by the invisible neutralinos. To assess the sensitivity of the LHC searches at 8 and 13 TeV, we generate inclusive samples of SUSY events with PYTHIA 8.150 [172,173], using the CTEQ6L1 parton distribution functions [174]. Delphes 3.0 [175] is then used to simulate the detector response and obtain the physics objects of the signal events. For each of the analyses, the signal selection cuts are applied to the simulated events, and the SM background events are taken from experimental publications. The CLs method [176] is used to obtain the 95% confidence level (C.L.) exclusion in presence of background only.
Monojet and mono-W, Z searches on the other hand have been designed in order to detect invisible particles in the final states through the detection of a hard jet emitted by the initial states. The basic idea is to search for a jet with high p T associated to a large missing E T . The main background for monojet searches stems from Z or W -boson and a jet, with the Z-boson decaying to neutrinos and the W -boson decaying to leptons which are missed by the detector.  Considering models in which a single mediator relates the dark matter particles to the SM particles reveals that the LHC can have a competitive or even superior reach compared to the dark matter detection experiments [177][178][179][180]. However, the simple description of dark matter production at the LHC based on a single mediator is not realistic with regard to concrete models such as the pMSSM, in which co-annihilations are favoured by the relic density constraints. Indeed, SUSY particles such as squarks or gluinos can be close in mass to the lightest neutralino, so that the production of two squarks or gluinos associated to a hard jet can still be seen as a monojet, because the jets produced in their decays would be soft enough to remain undetected [181][182][183][184][185]. In addition, several mediators can be involved. As a consequence, the single mediator limits cannot be recast in the pMSSM in a simple way.
To study the exclusion by the monojet and mono-W, Z searches, we use MadGraph 5 [186] to compute the full 2 → 3 matrix elements for all the combinations of pp →q/g +q/g + j/W/Z, pp →˜ +˜ + j/W/Z and pp →χ+χ+j/W/Z, whereq refers to a squark of any type and generation,g to the gluino,˜ to any type of sleptons,χ to any electroweakino. j corresponds to a hard jet as required for the monojet searches, and W/Z for mono-W, Z searches. As for the SUSY searches, we adopt the CTEQ6L1 parton distribution functions, hadronisation is performed using PYTHIA 8.150, and detector simulation with DELPHES 3.0. The cuts, selection efficiencies, acceptances and backgrounds for the 8 and 13 TeV runs are taken from the experimental publications cited in Table VIII. In addition, as the systematic uncertainties can have an important effect on these limits [184,187,188], we account for them by adding a 30% uncertainty to the cross sections.

Indirect detection
The annihilation cross sections necessary for the interpretation of indirect detection data are calculated with MicrOMEGAs [189][190][191], and PPPC4DMID [192] is used for the antiproton and gamma spectra.
a. Antiprotons We derived constraints on the dark matter annihilation cross section σv from the cosmic ray antiproton flux measured by PAMELA [44] as well as AMS-02 [45], following the procedure described in [75]. Secondary antiprotons constitute the astrophysical background. They are mostly produced by the interaction of primary proton and helium cosmic ray nuclei on the hydrogen and helium atoms lying in the interstellar medium, and their production rate is given by where T i is the kinetic energy of the nucleon i. The differential cross section dσ ij→pX /dTp is computed from the one for proton-proton interactions taken from [48] and the threshold T 0 p of this reaction is taken to be 7m p . The factor N IS takes into account the production of antineutrons decaying into antiprotons. Note that the bulk of the antiprotons are produced in proton-proton reactions. Using the retropropagation technique, we computed the fluxes Φ i (x) everywhere in the Galaxy from the fluxes Φ i ( ) measured at the Earth by PAMELA [193] or AMS-02 [194,195] when deriving constraints from the PAMELA or AMS-02 antiprotons flux, respectively. We further renormalise the production rate (17) by the factor A which takes into account the energy dependent uncertainties on the production cross section as well as on the antineutron yield based on the analysis of [48].
The production rate Q Ī p of primary antiprotons produced by the annihilation of two dark matter particles into the channel j is given by the expression (3) where the energy distribution of antiprotons per annihilation dN j p /dTp is taken from PPPC4DMID. The propagation of primary and secondary antiprotons is computed with the semi-analytical scheme described in Sec. II B 2 for the Med and Max sets of propagation parameters. The inelastic but non-annihilating interactions of antiprotons with the interstellar medium (tertiary component) are treated as in [75]. Practically, we apply the same procedure described in [75] to derive the 95% C.L. upper limit on the annihilation cross section σv , considering the secondary antiproton production uncertainty factor A and the Fisk potential φ F as profiling parameters.
The 95% C.L. upper limits on the annihilation cross section derived from the AMS-02 data are plotted in Fig. 1 with respect to the dark matter mass for the W , b, t and Z annihilation channels. These limits were computed for the Med (dashed) and Max (solid) propagation models and for the Einasto (red), NFW (green) and Burkert (blue) galactic mass models. The primary antiproton flux is maximised (minimised) by the Max (Med) model and the upper bound value of σv is thus minimised (maximised). In addition, the DM density is much larger in the Galactic center for the cuspy Einasto and NFW profiles than for the Burkert one. As a result, the annihilation rate integrated out over the Galaxy is higher for the former, leading to more stringent constraints (even if the local DM density is larger for the latter). As expected, the limits derived using the Med model characterised by a smaller value of the halo size L, is less sensitive to the shape of the halo profile since L is the "propagation horizon". In any case, the theoretical uncertainties coming from the poor knowledge of the propagation parameters is larger (up to a factor 4 on σv ) than the one arising from the choice of the DM profile (up to a factor 2). We found that for all annihilation channels and for the whole energy range we consider, the upper limit of σv is maximised for the Burkert-Med couple and minimised for the Einasto-Max one, providing an assessment of the theoretical uncertainties of our limits. Note that the constraints for the b quark channel become very stringent when m DM falls below 50 GeV, excluding the thermal relic cross section up to 3 orders of magnitude at 10 GeV.
For the sake of consistency, we performed the same analysis using the PAMELA proton, helium and antiproton data. The comparison between the results for the W boson channel obtained with AMS-02 and PAMELA data are plotted in Fig. 2 for the Burkert-Med and Einasto-Max cases. For m DM 1 TeV, the constraints derived from PAMELA data are more stringent than the AMS-02 ones. This can be understood by the fact that below ∼1 TeV, the proton flux measured by PAMELA is larger than the one reported by AMS-02 by a factor of up to 10%, leading to a larger yield of secondary antiprotons and thus a smaller room left for the primary component. The proton fluxes reported by the two experiments become similar above ∼1 TeV and in this regime, the experimental errors of AMS are much smaller than the PAMELA ones, leading to slightly stronger constraints for the more recent experiment. In the following of this paper, we consider only the results obtained using the AMS-02 data since they are more recent and they provide globally more conservative results.
As pointed out in [54] and in this analysis, one of the main uncertainties on the constraints derived from the cosmic ray antiprotons arise from the lack of knowledge of the propagation parameters. Hopefully, this uncertainty will substantially shrink since the AMS-02 Collaboration have recently released the long-awaited boron to carbon ratio [196], crucial measurement for the determination of the galactic transport properties.
b. Gamma rays We now turn to a combined analysis of the 19 confirmed dwarf spheroidal galaxies (dSphs) recently observed by Fermi-LAT [197].
For each dSphs and for each model point, we calculated the expected gamma ray flux per energy bin from dark matter annihilation Φ(∆Ω, E min , E max ) = annihilation channels 1 4π where dNγ dEγ channel is the gamma ray spectrum produced from dark matter annihilation, which depends on the dark matter mass and its annihilation channel and is obtained by interpolating the spectra tabulated in the PPPC4DMID [192,198]. The energy bins are those indicated in [197]. We were able to compute a delta-log likelihood for each of the points using the tabulated bin-by-bin likelihoods released by the Fermi-LAT Collaboration for each target [199] and we excluded points at the 95% C.L. We include statistical uncertainties on the J-factors of each dwarf spheroidal galaxy by adding an additional J-factor likelihood term, as prescribed by the Fermi-LAT Collaboration in their study. Those J-factors were calculated assuming a NFW profile, but previous work showed that the limits calculated with other halo profiles differed only by ∼30%, the strongest difference being for Burkert halo profile [200]. One of the largest uncertainties on these limits seems to reside in the choice of the dSphs sample used in the analysis. As pointed out in [197], adding galaxies with low-significance excesses, such as Reticulum II and Tucana III, can weaken significantly these limits.
Assessing the effects of such uncertainties seems to be very delicate and we will use these limits only as comparison with the constraints coming from antiprotons.
In addition, we considered the limits given by the HESS Collaboration [57]. As they do not use the same set of parameter values for the DM halo profiles as ours, we renormalised their limits following the J-Factors calculated for our different halo profiles NFW, Einasto and Burkert to be consistent with the rest of our study. The strongest limit being obtained for the NFW profile, we noticed that it barely reaches the distribution of our points without excluding any.

Direct detection
We calculated the WIMP-nucleon effective couplings f p and f n with MicrOMEGAs [191,201]. In our sample of points, the approximation f p ≈ f n used in the calculation of experimental SI cross section limits is reasonable for Higgsino-like neutralinos, but is not necessarily correct for other neutralino types. There is a simple way to cope with this problem. The experimental WIMP-nucleon limits can be described more generally in terms of WIMP-nucleus limits, averaged over all the target isotopes, and renormalised to a WIMP-nucleon limit in the case f p = f n . Consequently, for a given point, the appropriate quantity to be compared to the experimental limit is: where the subscript i stands for the various isotopes present in the experiment and η i is their corresponding abundance. These quantities depend on the target nucleus and are, a priori, different for xenon and argon. However, in our sample of points, we noticed that the was quite small, verifying δ 10% (δ 1% for the great majority of the points). The limits coming from xenon and argon experiments can then be easily compared, the XENON1T limit being the strongest one for our points.
Concerning the SD cross section limits, such problems do not exist. For the WIMP-neutron cross section, we will apply the limit given by the LUX experiment on our sample of points and for the WIMPneutron cross section, we will use the one given by the PICO-60 experiment. We also tested the limits given by IceCube, using the W + W − channel which is dominant for our wino-like and Higgsino-like neutralinos, and verified that the points excluded by the IceCube limit were already excluded by XENON1T or PICO-60.
In addition, we examine how the uncertainties on the local dark matter density and on the disc rotation velocity impact these limits. By rescaling the cross section coordinates, we obtain the limits for the three local density values 0.2, 0.4 and 0.6 GeV/cm 3 . In order to test the impact of v rot uncertainties, we proceeded to a variable substitution in the integral of the velocity distribution appearing in the calculation of the differential recoil rate per unit detector mass. To perform such a calculation, it was necessary to consider that v esc ≈ ∞. This approximation induces errors only for low WIMP masses that are not concerned by our study. We were then able to rescale the upper limits originally calculated with v rot = 220 km/s for two other values v rot = 200 and 250 km/s. Basically, taking smaller values for v rot shifts the limit to the right relative to m DM , and taking larger values shifts it to the left. The impact of ρ 0 and v rot uncertainties on the XENON1T 90% C.L. upper limit is shown in Fig. 3. The uncertainties on v rot within the considered values have a small impact compared to the local density uncertainties. Moreover, the uncertainties on v rot have a mild influence on the neutralino type of the excluded points and change the fraction of excluded points by less than 1%. For these reasons, and for sake of simplicity, we keep, in the rest of this study, the benchmark value v rot = 220 km/s and vary only the local dark matter density value.

IV. RESULTS
We consider only model-points which have the lightest neutralino as LSP and dark matter particle, as described in section III A. More than 20 million pMSSM points have been generated randomly with parameters in the ranges given in Table III. Prior to studying the effects of the different constraints, we impose the light Higgs mass to comply with the interval given in Eq. (15).
In the following, the neutralino 1 (denoted χ) will be said to be bino-/wino-/Higgsino-like if it is composed of 90% of bino-/wino-/Higgsino component, respectively, or mixed state otherwise. In Fig. 4, the composition of our sample of pMSSM points after imposing the light Higgs mass interval is shown. Binolike χ are the most represented points in our sample, followed by the winos and Higgsinos, with an almost equal share of each component. The fraction of mixed states is negligible.

A. Relic density constraints
We first consider the relic density constraint. The value of the neutralino relic density is computed with SuperIso Relic v3.4 [202,203]. In Fig. 5, the relic density is shown as a function of the neutralino 1 mass, for the different types. Bino-like neutralinos 1 have in general large relic densities, above the Planck measurement. This can be explained by the smaller couplings of the binos with SM particles, which leads to smaller annihilation cross sections and therefore larger relic densities. On the other hand, the Higgsino-like χ give smaller relic densities which are close to the Planck measurements for χ masses around 1.3 TeV. The wino-like χ tend to have even smaller relic densities, and the Planck line is naturally reached for a mass of 2.7 TeV. The line at about 90 GeV in the figure corresponds to cross section enhancements through a Z-boson resonance, which lower the relic density.
Imposing both the upper and lower relic density bounds generally leads to a selection of scenarios with co-annihilations, for which the mass splitting of the neutralino 1 with the next-to-lightest supersymmetric particle is small, or of scenarios where χ annihilations are enhanced through a resonance of the Zboson or one of the neutral Higgs bosons. This is demonstrated in Fig. 6. The valid points require in general small mass splitting, apart from some spread binos with larger mass splittings, which have a heavy Higgs boson or Z-boson resonance. For the case of winos, the small mass splitting is due to a chargino with a mass very close to the χ mass. For the Higgsino case, both the chargino 1 and the neutralino 2 have masses close to the neutralino 1 mass.
As discussed in Section II A, we consider only the upper bound of the Planck dark matter density interval, which favours light wino-and Higgsino-like χ, and bino-like χ with strong co-annihilations.

Constraints from AMS-02 and Fermi-LAT
We consider the constraints from AMS-02 antiproton and Fermi-LAT gamma ray data, which probe specific dark matter annihilation channels. For both sets of constraints, the most important parameters are the χ annihilation cross sections into specific channels. Annihilations to W W and bb are particularly interesting in the context of the pMSSM.
In Fig. 7, the total annihilation cross section times velocity σv tot is shown as a function of the neutralino 1 mass, for the different χ types. σv tot is the sum of all the σv of the different channels. The wino-and Higgsino-like neutralino 1 regions form two separate strips. The different types of neutralinos 1 have specific main decay channels: binos annihilate mainly into tt, bb, and in a lesser extent into W h, Zh and τ τ , Higgsinos into W W and ZZ, and winos into W W , when the decay channels are open. When the above-mentioned channels are closed because of a small neutralino 1 mass, the χ mostly decays to bb and τ τ , and less frequently into cc and ss, independently from their type. As seen earlier, winos more strongly annihilate than the other χ types, followed by the Higgsinos. The binos, apart from the case of a resonant annihilation, are more weakly annihilating and are far below the experimental limits.
In Fig. 8, the exclusion by Fermi-LAT and AMS-02 is shown in the σv tot vs. neutralino 1 mass parameter plane. In order to quantify the uncertain- ties related to indirect detection, we consider separately the most conservative limits, i.e. obtained using Burkert profile and Med propagation model, and the most stringent ones, i.e. using Einasto profile and Max propagation model. The conservative limits lead to the exclusion of neutralinos 1 with masses between 90 and 550 GeV, which are mainly wino-like. The stringent limits exclude points with χ masses between 0 and 850 GeV. In the small mass region, as well as for masses above 90 GeV, the stringent exclusion limit is strengthened by one order of magnitude in comparison to the conservative case. The stringent case excludes large zones of the wino strip, and of the Higgsino one in a lesser extent. AMS-02 alone brings very strong constraints in the stringent case, beyond the Fermi-LAT limits.

Connections with relic density
Indirect detection constraints may be considered to be redundant with the relic density constraint. This is generally true for simplified dark matter models [180], because the relic density is directly related to the annihilation cross sections. However, in a complete model such as the MSSM, the value of the relic density is often led by the co-annihilations, especially when both the upper and lower bounds of the Planck dark matter density measurements are applied. This was already shown in Fig. 6.
Yet, there is a strong complementarity between indirect detection and relic density, as shown in Fig. 9. Considering the gray points, we see an anti-correlated region where the relic density increases when the annihilation cross section decreases, which is due to the relation between the relic density and the annihilation cross sections. This region is largely excluded by the upper Planck bound. The points with small relic density have in general efficient co-annihilations, which reduces the relic density. While these points are far from being excluded by the Planck upper bound, they can be probed by the stringent AMS-02 limits obtained using the Einasto profile and Max propagation model. This clearly shows the complementarity between indirect detection and relic density constraints.
In Section II A, we discussed how the relic density constraint can be falsified. One of the possibilities is that the dark matter density measured by Planck is made only in part of neutralinos, the rest being made of other types of particles or more exotic objects. In such a case, galactic haloes would also be composed of different types of dark matters. Assuming that the mixture of dark matters is in the same proportion in galaxies as in the large scale Universe, the neutralino relic density is smaller than the measured dark matter density, and the dark matter density in galactic haloes has to be rescaled by the ratio of the neutralino relic density over the dark matter density, hence impacting the indirect detection limits. This is done in Fig. 10, in the total annihilation cross section vs. neutralino 1 mass parameter plane, where the total annihilation cross section is rescaled by the neutralino relic density over the measured dark matter density. Such a rescaling strongly weakens the indirect detection limits. Indeed, even using the most stringent AMS-02 constraints, only a very few points in the low mass region are still excluded, mostly in the bb channel. The large negative impact of the rescaling is due to the fact that the constraints from indirect detection scale as the squared density, leading to a strong loss of sensitivity.

Constraints from XENON1T, LUX and PICO-60
Contrary to relic density and indirect detection, which mainly depend on the annihilation and coannihilation cross sections, direct detection relies on the scattering cross section of neutralino 1 with nucleons. Direct detection is therefore complementary to indirect detection and relic density.
In Fig. 11, the generalised spin-independent WIMP-nucleon cross section -which roughly corresponds to the χ-xenon scattering cross section normalised to one nucleon, and which applies to xenonbased experiments -is shown as a function of the neutralino mass, for the different neutralino 1 types. Higgsinos are in general more strongly interacting than the winos, leading to larger cross sections. In order to assess the consequences of the uncertainties on the obtained constraints, the recent limits of the XENON1T experiment are superimposed, for three values of the local dark matter density, namely ρ 0 = 0.2, 0.4 and 0.6 GeV/cm 3 . Between the conservative line corresponding to ρ 0 = 0.2 GeV/cm 3 and the most stringent limit obtained for ρ 0 = 0.6 GeV/cm 3 , there is at most a factor 3 difference. While this is a large factor, in the context of pMSSM it does not change much the excluded region, which contains mainly Higgsino-like neutralinos 1.
In Fig. 12 the exclusion by the XENON1T data with ρ 0 = 0.4 GeV/cm 3 is shown in the (M A , tan β) parameter plane. For each bin, the fraction of excluded points is presented. This parameter plane is of interest since the neutral Higgs bosons can mediate the scattering, with couplings proportional to tan β. About 100% of the points are excluded in a triangle region starting from the origin of the plot and up to tan β = 60 and M A = 600 GeV. A large fraction of the points with larger M A can also be excluded. For comparison, the exclusion line from the CMS heavy Higgs searches for H/A → τ τ is also shown [171]. While the CMS limit extends beyond the 100% exclusion triangle and constitutes a well-defined and robust exclusion in this parameter plane, direct detection still adds complementary constraints for larger M A and smaller tan β values.
LUX and PICO-60 also provide important constraints on the spin-dependent scattering cross section with protons and neutrons. This is shown in Fig. 13, for ρ 0 between 0.2 and 0.6 GeV/cm 3 . The distribution of the points is different for the proton and neutron scatterings, because the wino-neutralino 1 mixing term in the neutralino-quark-squark coupling is proportional to the isospin. In both cases however, only the most strongly interacting Higgsinos are excluded, and the value of ρ 0 does not affect much the results. Practically, LUX and PICO-60 spin-dependent constraints are redundant, since both exclude the same points. The spin-independent XENON1T results give quite stringent constraints, which exclude most of the points probed by LUX and PICO-60. After imposing the XENON1T constraints, the spin-dependent results exclude about 0.5% of the remaining points, with dominantly Higgsino-like χ.

Connections with relic density
Direct detection constraints are not related to the relic density through annihilation cross sections, as for indirect detection. They are nevertheless complementary, since they provide constraints on different pMSSM parameters.
The same paradigm as for indirect detection can apply: if the relic density is smaller than the observed dark matter density, it may be because the neutralino is not the sole component of dark matter, thus the local dark matter density has to be rescaled accordingly to obtain the local neutralino density. As a consequence, the limits become less constraining, since the effective scattering cross sections are lowered by a factor proportional to the relic density. In comparison with indirect detection, the impact of the rescaling is less pronounced, because the rescaling is proportional to the dark matter density for direct detection, whereas it is proportional to the density squared for indirect detection.
In Table IX, the fractions of excluded points are given for the different neutralino 1 types, with rescaling and without rescaling, for ρ 0 = 0.2, 0.4, 0.6 GeV/cm 3 , after the upper limit of the relic density is applied. First, in absence of rescaling, even in the most conservative case corresponding to ρ 0 = 0.2 GeV/cm 3 , direct detection imposes strong limits, and one third of the points are excluded. The Higgsinos are the most affected, followed by the binos and winos. The mixed states are almost completely excluded by direct detection, but their number is too small to draw statistically significant conclusions. When increasing the density to ρ 0 = 0.6 GeV/cm 3 , the sensitivity is enhanced, with about half of the points excluded, and 70% of the Higgsinos. With the relic density rescaling, the exclusion power decreases strongly, as only 15% of the points remain excluded in the most favourable case. The exclusion hierarchy is also modified in presence of rescaling, with the binos being the most excluded neutralinos 1.

Combined dark matter constraints
Dark matter observables can lead to very strong constraints. The relic density, if compared to both upper and lower bounds of the measured dark matter density, leads to a very strong exclusion and a selection of the points with small mass splittings or resonant annihilations. However, the relic density constraint suffers from uncertain hypotheses about the Early Universe, which prevents us from considering the lower bound. Applying only the upper bound, most of the bino-like χ which are not accompanied by other SUSY particles close in mass are excluded. Indirect detection brings complementary constraints, and probes wino-like and Higgsino-like neutralinos. We showed that the different assumptions on the galactic halo profiles and propagation models for the cosmic rays can however modify the limits on the cross sections by a few orders of magnitude and strongly lower the constraining power of indirect detection limits. Direct detection on the other hand sets strong constraints on the pMSSM, with an exclusion of 25-40% of our scan points, depending on the local dark matter density. The main uncertainty is due to the local dark matter density, but even if the conservative choice lowers the exclusion power, it does not affect much the excluded parameter region.
It is however important to recall that in the case the neutralino 1 is not responsible for the whole dark matter density and a rescaling of the dark matter density is necessary for direct and indirect detections, dark matter observables severely lose their constraining power. In the following, we focus on the case where neutralinos constitute the whole dark matter and do not consider the scaling possibility any further. We will now quantitatively study the interplay between the different dark matter constraints. We define three cases: • CONSERVATIVE: ρ 0 = 0.2 GeV/cm 3 for direct detection, Burkert dark matter profile and cosmic ray Med propagation model using AMS-02 data for indirect detection.
• STANDARD: ρ 0 = 0.4 GeV/cm 3 for direct detection, NFW dark matter profile using the combined analysis of the 19 confirmed dwarf spheroidal galaxies observed by Fermi-LAT for indirect detection.
• STRINGENT: ρ 0 = 0.6 GeV/cm 3 for direct detection, Einasto dark matter profile and cosmic ray Max propagation model using AMS-02 data for indirect detection.
In Fig. 14, the fraction of pMSSM points initially satisfying the light Higgs mass constraint, which are excluded by the upper bound of the dark matter density, direct detection and indirect detection constraints, is shown for the three cases of astrophysical assumptions. The & symbol corresponds to the exclusive "and". Points excluded simultaneously by the relic density and indirect detection constraints represent less than 1% of the total number of points, and are not shown.
The relic density constraint excludes about 36% of the points. As already seen, direct detection constraints are relatively insensitive to the choice of the local density of dark matter, and direct detection excludes 25% of the points in the conservative case and 35% in the stringent case. Indirect detection is more sensitive to the choice of profile and propagation model and excludes less than 20% of the points in the conservative case and 30% in the stringent one. In all cases, the simultaneous application of the dark matter constraints is very important, and allows us to strongly reduce the number of valid points, even in the most conservative case.
In Fig. 15, the same analysis is performed for the different neutralino types separately.
First, the bino-like neutralinos 1 have in general weaker couplings, leading to large relic densities and small annihilation and scattering cross sections. Thus, the bino-like points are strongly excluded by the relic density, slightly probed by direct detection, and negligibly by indirect detection. Therefore, the choice of the conservative or stringent constraints has a negligible effect, since the exclusion is dominated by the relic density. Second, wino-like neutralinos 1 are dominantly excluded by indirect detection, followed by direct detection. After these constraints, relic density only affects a negligible fraction of points, which is why the exclusion by relic density does not appear in the figure. For the winos, the choice of the conservative or stringent cases strongly affects the results, leaving 50% of the points valid in the conservative case, and 28% in the stringent case. Again, the standard case leads to results similar to the stringent case. Third, the Higgsino-like neutralinos 1 are mainly excluded by direct detection, which mildly depends on the astrophysical hypotheses. Indirect detection also excludes a number of points, even if a large fraction of them is already excluded by direct detection. As for the winos, relic density only excludes a negligible fraction of points after the direct and indirect detection constraints. At the end, 40% of the Higgsinos remain valid in the conservative case, and 20% in the stringent case. Finally, the mixed-state neutralinos 1 are completely excluded independently from the astrophysical hypotheses, and predominantly by direct detection.
To summarise this section, dark matter constraints set strong constraints on the pMSSM parameter space. However, while direct detection leads to relatively robust constraints, indirect detection is more sensitive to the choice of galaxy halo profiles and cosmic ray propagation models.

D. Collider and Dark Matter constraints
In this section, the complementarity of collider and dark matter constraints will be studied.
Whereas dark matter limits are subject to astrophysical and cosmological uncertainties or hypotheses, collider constraints are obtained in environments under control, which therefore lead to relatively hypothesis-free limits.
As explained in Section III B, using the LHC results requires the computation of numerous cross sections, generation of events and detector simulation, which are computationally heavy and CPU-time consuming. In order to gain CPU time, we perform the event generation and detector simulation only for model points which respect the light Higgs mass constraint, flavour physics, and LEP and Tevatron constraints, as well as the upper bound of the relic density. The points satisfying these constraints will be referred to as "Accepted points" in the following.
In Fig 16, we present the type of neutralinos 1 for the accepted points. A comparison with Fig. 4 showing the type of the points satisfying only the light Higgs mass limit, reveals that most of the binos have been excluded, but that the fraction of winos in comparison with the Higgsinos is unchanged. This is mainly due to the upper bound of the relic density, as explained in Section IV A. The LEP and flavour constraints do not probe directly the neutralino 1, but can affect scenarios with light wino-like and Higgsino-like χ through the constraints on the charginos and heavier neutralinos. Nevertheless, the exclusion power of these constraints is limited in comparison to the relic density one.
In the Higgs sector, the light Higgs mass constraint favours the decoupling limit where the heavy Higgs bosons are heavy, and heavy stop masses with maximal mixing [122,[204][205][206]. Measurements of the light Higgs production and decay channels also point towards large heavy Higgs masses. In particular, the diphoton channel favours heavy charginos, stops and charged Higgs bosons [122,123,[207][208][209]. In addition, light Higgs decays into supersymmetric particles are rather limited [122,[124][125][126][127]. These important limits provide strong constraints in the (µ, M 2 ) parameter plane. Indeed, both parameters are important for the neutralino and chargino mixings, and µ is also important for the third generation squark mixings. The limits obtained from the measurements of the light Higgs couplings are complemented by the electroweakino direct searches at LEP and the LHC. This is illustrated in Fig. 17, where the small µ values are excluded. The complementarity with dark matter constraints is rather clear. Direct detection excludes points spread over the plane. Indirect constraint severely excludes points with M 2 600 GeV and |µ| 150 GeV. One should however note that due to the multi-dimensional parameter space, there could be points below the coloured regions that still survive the dark matter and collider constraints.
The heavy Higgs searches, and in particular H/A → τ τ searches, impose strong constraints in the (M A , tan β) parameter plane which is also relevant for direct detection as seen in Fig. 12. In Fig. 18, we superimpose over the points in agreement with the LHC constraints those which are excluded by direct detection. Similarly to direct detection, H/A → τ τ searches probe the large tan β and small M A region (corresponding to the empty region in the upper right part in the figure). We can see from the figure that the exclusion by direct detection is not well defined and spread. Comparing with Fig. 12 reveals that the strongest and well defined exclusion by direct detection in this plane occurs below the H/A → τ τ limit. Both constraints are nevertheless complementary and allow us to exclude points beyond the large tan β and small M A region.
As a hadron collider, LHC is more sensitive to strongly interacting particles. In particular, gluinos and squarks of the first and second generations are amongst the most actively searched particles, and LHC can probe masses as large as a few TeV in the most favourable scenarios. In Fig. 19, the accepted pMSSM points are plotted in the minimum mass amongst the gluino and first and second generation squark masses vs. neutralino 1 mass plane. We note that gluinos or squarks as light as a few hundred GeV can still escape LHC searches in a general scenario as the pMSSM. These points correspond mainly to compressed scenarios [210][211][212][213], where one or more supersymmetric particles have masses close-by, leading to decays with particles or jets in the final state which can leave the detectors undetected because of their small transverse energies. Dark matter searches the squark and gluino masses. We also see that after the LHC constraints, light squarks or gluinos of a few hundred GeV in compressed or complicated scenarios are still allowed, but after the dark matter constraints, they are less numerous and the surviving points correspond to very small squark/gluino-neutralino 1 mass splittings, and in the stringent case the squark and gluino masses are pushed beyond 450 GeV. So the complementarity is obvious, as dark matter experiments can probe parameter regions which are not accessible at the LHC, and vice versa. Similar result for the lightest stop is presented in Fig 20. As for the gluino and squark case, light stops are still allowed by collider constraints in compressed scenarios, which can still be probed by dark matter detection experiments. Light stop scenarios which escape LHC detection are still allowed, but the stop 1 mass is pushed beyond 500 GeV in the conservative case and 600 GeV in the stringent case, after imposing the direct and indirect detection limits.
Finally, the interplay of the LHC and dark matter constraints is presented in a quantitative way in Fig. 21. It can be seen that the LHC has the major role in probing the pMSSM parameter space, but dark matter detection constraints further probe the parameter space. The combination of all constraints leads to an exclusion of between 85% and 92% of our sample. Fig. 22 presents a more detailed view of the exclusion for the different neutralino 1 types. In particular, it reveals that LHC excludes more than 65% of the points independent of the neutralino 1 type. The role of dark matter constraints on the contrary is more type-dependent. As we showed earlier, binos, Higgsinos and mixed states are more strongly probed by direct detection, while indirect detection rather excludes winos. And whereas direct detection is mildly sensitive to the choice of the astrophysical parameters, indirect detection is more sensitive to it.
Finally, in Fig. 23, the fraction of neutralino 1 types after imposing all the constraints is shown. This figure is to be compared with Fig. 21, where only LEP, flavour and relic density constraints were applied. We can see that the final fractions are still similar after applying all constraints, with a larger proportion of winos, followed by a large proportion of Higgsino, and a small amount of binos. This shows that the relic density constraint is the most type-selecting constraint. However, we note that the proportion of winos is much larger in the conservative dark matter case than in the stringent case.
An important caveat here is in order. The fraction of points has no real statistical meaning, but rather shows the tendency of the constraints to select certain types. To illustrate this, we show in Fig. 24 the fraction of the types after applying all the constraints, including the Planck lower bound. In this case, the Higgsinos are now the dominant surviving species, followed by the binos, and the winos survive only in small proportion. It is interesting to note that in this case, the choice of conservative or stringent astrophysical hypotheses does not affect much the results.
Before concluding, it is worth mentioning that great improvements in the sensitivity of the direct and indirect detection experiments are expected in the coming years. Concerning direct detection, in the next few years XENONnT [95] and LZ [96] will push the XENON1T limit by two orders of magnitude, and within ten years DARWIN [97] will allow us to gain one extra order of magnitude. This is illustrated in Fig. 25. For comparison, the gray points correspond to a sample of our points which are in agreement with the current LHC 8 TeV and 13 TeV limits. Practically, XENONnT/LZ will exclude most of the Higgsino points, and DARWIN will be able to probe a large part of the wino region. In addition, we have shown that the constraining power of direct detection is only mildly affected by the choice of the astrophysical assumptions, thus these limits will provide relatively robust constraints on the pMSSM parameter space. The DARWIN limit will however be close to the neutrino background, which constitutes a large obstacle to further improvements. Nevertheless, the remaining points after DARWIN will have mainly winolike neutralinos 1, which will be probed by indirect detection.
For indirect detection, the Cherenkov Telescope Array (CTA) [214], dedicated to gamma rays, will use a Cherenkov imagery technique similar to HESS, VER-ITAS or MAGIC, and will be able to probe an energy range between a few tenths of GeV to above 100 TeV. Before 2030, CTA will also further push the in-direct detection limits by observing gamma rays at the center of the Milky Way, as shown in Fig. 26. It is important to remark however, that contrary to the Fermi-LAT limits, which are obtained from the observations of spheroidal dwarves and which are therefore less affected by the dark matter profile, since CTA will focus on the galaxy center, it is subject to strong uncertainties from the dark matter profile. Since the question of the existence of cuspy profiles is unresolved [215], dark matter density distributions such as NFW or Einasto which incorporate cuspy profiles, will lead to fundamentally different exclusion limits than a Burkert profile with a core. This is illustrated in the figure, a Burkert profile will lead to limits which are two orders of magnitude less constraining than the NFW or Einasto profile. Therefore, CTA will be even more subject to astrophysical uncertainties, even if we can hope for an improvement of our knowledge of the galactic center within the next decade.

V. SUMMARY
In this paper, we studied the impact of dark matter direct and indirect detections, in conjunction with relic density and collider constraints, on the phenomenological MSSM with neutralino dark matter and addressed in some detail the consequences of the related uncertainties.
First, the calculation of the relic density is based on the assumption of radiation dominating the Early Universe properties. Any deviation from this hypoth-esis can modify the relic density by orders of magnitude. In addition, dark matter may be made of different components, the neutralino being only one of them. These considerations justify the usage of solely the upper relic density bound. In the pMSSM, applying both relic density bounds selects compressed scenarios with co-annihilations or with annihilation through a Z-boson or Higgs boson resonance, and favours bino-like neutralinos 1. Disregarding the lower bound strongly changes this picture, and rather selects Higgsino-and wino-like neutralinos 1.
We then reviewed the calculation of indirect detection constraints, which relies on the choice of a dark matter halo profile for gamma rays, as well as a propagation model for cosmic rays. We showed that between the more conservative case, corresponding to the Burkert halo and the Med propagation model, and the most stringent case, Einasto profile and Max propagation model, the limits from AMS-02 differ by one order of magnitude. In the context of pMSSM, indirect detection excludes more strongly the wino-like neutralinos 1, followed by the Higgsino-like ones.
Turning to direct detection, we showed that the constraints are affected by the local relic density in the vicinity of the Earth, whose evaluations can vary within a factor 3, and to a lesser extent by the velocity of Earth in the dark matter halo. As a consequence, the exclusion limits obtained by direct detection experiments can vary by a factor 3. Even if this uncertainty is rather large, in the context of pMSSM, it does not strongly affect the excluding power of direct detection. Independent of the choice of the local density, direct detection is particularly efficient in probing scenarios with Higgsino-like neutralino 1.
An interesting aspect of the connection of direct and indirect detections with relic density comes from the possibility that dark matter could be made of several components. In such a case, the neutralino relic density is smaller than the measured dark matter density. Therefore, the local and galactic neutralino dark matter densities are smaller than the measured ones, and have to be rescaled by the ratio of the relic density over the dark matter density measured by Planck. Such a rescaling strongly alleviates the indirect detection constraints since they are proportional to the density squared, and decreases to a lesser extent the direct detection limits, which are proportional to the density.
Apart from this specific case, direct detection, indirect detection and (the lower bound of) relic density are very efficient in the pMSSM, excluding a large part of our sample. Even with the most conservative choice of astrophysical uncertainties, 70% of our sample is excluded by the dark matter constraints, and 85% in the most stringent case. As expected, constraints from indirect detection are the ones which are the most affected by the astrophysical assumptions, but the complementarity between the dark matter constraints is still of major importance when studying BSM scenarios.
We then studied the interplay of dark matter constraints with collider constraints from LEP, Tevatron, B-factories and LHC searches. The LHC is by design more apt to probe the strong interaction sector than the weak interaction one. On the contrary, dark matter searches probe the weak interaction sector. As a consequence, there is an interesting complementarity between the two kinds of searches. The Higgs boson however is also related to the weak sector. In addition, flavour physics and heavy Higgs searches allow us to explore the chargino and heavy Higgs sectors.
There exists two parameter planes where the complementarity between collider and dark matter constraints is obvious: the (M A , tan β) plane, which is probed by BR(B s → µ + µ − ) in flavour physics, H/A → τ τ searches at the LHC and dark matter direct detection. We showed that, while heavy Higgs searches provide very robust constraints in this parameter plane, dark matter direct detection can provide less strict constraints which spread beyond those of the LHC. The second plane is (M 2 , µ), which is very important for the light Higgs and chargino sectors. LEP and LHC searches for electroweakinos and Higgs results allowed us to exclude small |µ| and M 2 , but direct and indirect detections further probe this parameter region up to M 2 600 GeV and |µ| 150 GeV.
Furthermore, the LHC excludes strongly gluino and squarks of intermediate masses, but is less sensitive to scenarios with compressed spectra, which lead to invisible final states. As a consequence, light squarks or gluinos of a few hundreds of GeV can still escape detection. We showed that dark matter constraints ex-clude these light gluinos and squarks, including light stops. This highlights the importance of the complementarity between dark matter and collider constraints.
Finally, the latest collider constraints alone exclude 70% of our sample of points, and dark matter constraints alone between 55% and 80% depending on the astrophysical assumptions. Altogether, the exclusion reaches between 85% and 93% of our scan points, showing again the complementary of the collider and dark matter experiments, regardless of the astrophysical hypotheses.
In the future, there will be interesting prospects for dark matter direct and indirect detections. In particular, in the coming years XENONnT and LZ will improve the current limits by two orders of magnitude, and DARWIN within ten years by one extra order of magnitude. Similarly, the gamma ray telescope CTA will strongly improve indirect detection limits by 2030. Yet, astrophysical uncertainties constitutes a limitation. We can however hope for improvements in our knowledge of the halo profiles, local dark matter density and cosmic ray propagation models, which would lead to more robust constraints from direct and indirect detection experiments.