Theoretical uncertainties in the calculation of supersymmetric dark matter observables

We estimate the current theoretical uncertainty in supersymmetric dark matter predictions by comparing several state-of-the-art calculations within the minimal supersymmetric standard model (MSSM). We consider standard neutralino dark matter scenarios — coannihilation, well-tempering, pseudoscalar resonance — and benchmark models both in the pMSSM framework and in frameworks with Grand Unified Theory (GUT)-scale unification of supersymmetric mass parameters. The pipelines we consider are constructed from the publicly available software packages SOFTSUSY, SPheno, FeynHiggs, SusyHD, micrOMEGAs, and DarkSUSY. We find that the theoretical uncertainty in the relic density as calculated by different pipelines, in general, far exceeds the statistical errors reported by the Planck collaboration. In GUT models, in particular, the relative discrepancies in the results reported by different pipelines can be as much as a few orders of magnitude. We find that these discrepancies are especially pronounced for cases where the dark matter physics relies critically on calculations related to electroweak symmetry breaking, which we investigate in detail, and for coannihilation models, where there is heightened sensitivity to the sparticle spectrum. The dark matter annihilation cross section today and the scattering cross section with nuclei also suffer appreciable theoretical uncertainties, which, as experiments reach the relevant sensitivities, could lead to uncertainty in conclusions regarding the viability or exclusion of particular models.


Introduction
The search for supersymmetry and its connection to dark matter physics have been prominent areas of research in particle phenomenology, both theoretical and experimental, over the last few decades. Experimental results have provided significant constraints on the Minimal Supersymmetric Standard Model (MSSM) via the discovery of the Higgs boson [1,2], as well as via null results from collider searches for new particles [3,4] and dark matter direct and indirect detection experiments. The LHC in particular has pushed the limits on squark masses to roughly the TeV range, ruling out much of the constrained MSSM (CMSSM) (eg. [5]). Nevertheless, within the full MSSM as well as more minimal frameworks, much parameter space remains in which the thermal relic abundance of the lightest neutralino explains the astrophysical cold dark matter and a Higgs boson consistent with that discovered at the LHC is predicted . Furthermore, it can be argued that, despite not yet having discovered any new supersymmetric (SUSY) partners, the verdict is still out on (even weak-scale) supersymmetry (eg. [34] and [35]).
As the mechanism of supersymmetry breaking is far from clear, models may be defined either in the UV, near the so-called Grand Unification (GUT) scale, or in the IR, for example within the well-studied phenomenological MSSM (pMSSM) framework. In either case, one would like to calculate the superpartner spectrum at the weak scale and the masses and couplings in the Higgs sector. Finally, once the Lagrangian is known at the weak scale, it can be used to calculate the dark matter observables, including the relic density, the current annihilation cross section, and the scattering rates, all of which can be compared with experimental results.
A plethora of public software packages have been developed to facilitate the analysis of various SUSY models and compare their predictions with experimental results (for example, [36][37][38][39][40][41][42][43]). Our confidence in the accuracy of the calculations done by any package is based on the continual improvements made by its authors and on the agreement of results from different packages. Comparative studies of spectrum generators and Higgs sector calculators have been undertaken before (eg. [44][45][46][47][48][49]). Differences in the renormalization group running and predictions for sparticle masses within the same supersymmetric model have been observed [44][45][46], and sensitivities of the Higgs sector have also been explored [47][48][49]. However, previous studies have focused primarily on models with relatively light (O(100) GeV) sparticles, and the most recent comparison of the full sparticle spectrum was undertaken more than a decade ago [46].
There are several publicly-available software packages that calculate quantities that can be observed at dark matter direct and indirect detection experiments as well as the relic abundance of dark matter within a particular model. micrOMEGAs [42] and DarkSUSY [43] are two examples. Studies have been carried out regarding the accuracy with which a single observable is calculated by an individual software package [50][51][52], though there are relatively few studies that compare the calculation of dark matter observables by different software packages (eg. [53]), and none of which we are aware that address LHC-era supersymmetric benchmarks.

JHEP05(2018)113
In this report, we embark on a comparison study designed with three advancements over previous studies: first, we compare calculations for the sparticle spectrum, the Higgs sector, and the dark matter observables, discussing, when possible, differences in the implementations of the underlying physics of the calculations in each case. Second, we incorporate the various calculators into comprehensive pipelines to study not only the effects of the choice of an individual calculator, but also all downstream effects of those choices on subsequent calculations. Finally, we analyze the above choices and observables in the context of several SUSY benchmark models chosen as representative of models that are interesting in the light of LHC Run-1 and null results from recent dark matter searches as described below.
This study is conducted in two parts: to begin, we investigate a set of pMSSM models from ref. [54]. The dark matter scenarios we consider are coannihilation (bino-stop and bino-squark), A-funnel, well-tempered neutralinos, and pure higgsinos. We will see that the spectrum generators can differ by up to 1 -2 % in their predicted masses for the stop and the first two generations of squarks, and by up to 20% in the gauge composition of the lightest neutralino, for a given pMSSM model. As for the dark matter observables, differences of up to a factor of ∼ 3 − 5 in the relic density and current annihilation cross section, and up to a factor of ∼ 10 in the predicted scattering cross section are possible for the different pipelines. The theoretical uncertainty in the relic density of neutralino dark matter already far exceeds the statistical errors reported by the Planck collaboration, while the dark matter annihilation cross section today and the scattering cross section with nuclei also suffer appreciable theoretical uncertainties, which, as experiments reach the relevant sensitivities, could lead to uncertainty in conclusions regarding the viability or exclusion of particular models.
In the second part of our study, we consider four benchmark models defined at the GUT scale -two CMSSM points and two points from models with non-universal Higgs masses (NUHM) [26,30]. For GUT-scale models, we will find that discrepancies among the various pipelines are often amplified by the renormalization group running. For our CMSSM and NUHM benchmarks, we will see that the spectrum generators can give low energy values of the higgsino mass parameter µ and the pseudoscalar Higgs mass m A that differ by up to 150% -200% (though the differences can be much greater at larger m 0 than the values considered here). This leads to dramatic differences in the annihilation and scattering cross sections computed by the dark matter calculators.
Before proceeding, we would like to reflect on whether a study like this is anachronistic at this juncture. With the LHC failing to find new physics yet, and supersymmetric WIMP searches yielding null results, one might ask whether it makes sense to go back to benchmark SUSY scenarios yet again. We remind the reader that the connection of supersymmetry to dark matter physics, while robust from a high level perspective due to the WIMP miracle, was always fragile at the model-building level, at least under the assumption of a standard cosmological history. The dark matter relic density is often obtained in fine-tuned regions of parameter space, which either exhibit compressed spectra, or have suppressed interactions with nuclei. Many of these scenarios are difficult to probe at colliders or direct detection experiments, and also have small annihilation rates in the current Universe. Moreover, their fine-tuned nature means that detailed predictions for physical quantities in these scenarios are particularly sensitive to the approximations -3 -

JHEP05(2018)113
used. It is not a surprise that theoretical uncertainties in these scenarios can substantially outweigh experimental uncertainties. Given the current precision of the measurement of the dark matter abundance and the dramatically improving sensitivities to dark matternucleon scattering in the era of ton-scale experiments, it could be argued that it is more important than ever to examine the precision and accuracy of the predicted values for observable quantities in supersymmetric models. We also note that our findings may be relevant for some more general (non-supersymmetric) models of dark matter, so long as they share particular characteristics with the benchmarks considered here, for example, models in which the relic abundance is achieved via resonant annihilations.
The paper is organized as follows. In section 2 we present the benchmark MSSM points considered here and discuss the calculators and pipelines we study. In section 3 we discuss relevant aspects of the physics of neutralino dark matter. In sections 4 and 5 we present our results for the pMSSM and GUT-scale benchmarks, respectively. Finally, in section 6 we summarize the conclusions of our study. Numerical results for all benchmarks are compiled in appendices A and B.

Methodology of the comparison
In this section, we present our benchmark points, the calculator pipelines we study, and the details of the calculations undertaken by each of the calculators.

Benchmark points
We consider two sets of supersymmetric benchmark points, all of which assume that the lightest supersymmetric particle (LSP) is the lightest neutralino, which is therefore the dark matter candidate. For future reference, we describe our notation here. We choose the bino -wino -higgsino basis to write the neutralino mass matrix as where we follow the standard notation: s W = sin θ W , c W = cos θ W , s β = sin β, c β = cos β and tan β = v 2 /v 1 , with v 1,2 being the vacuum expectation values of the two Higgs fields H 1,2 . The mass matrix can be diagonalized by a unitary mixing matrix N , where the eigenvalues are the neutralino masses. The lightest neutralino mass eigenstate can be written asχ The first set of benchmark points consists of 5 pMSSM points from the Snowmass 2013 white paper ref. [54]. These points are representative of the pMSSM landscape and of the primary mechanisms by which the correct relic density of neutralino dark matter is achieved: sfermion coannihilation, rapid annihilation via a pseudo-scalar Higgs resonance, pure higgsino content, and the so-called well-tempered neutralino. We will discuss each of these in section 3. The spectrum for each point can be found at [55].
The second set consists of 4 points of interest defined at the GUT scale. The defining parameters of these 4 points may be found in table 1. Three of these points are based on the MasterCode Collaboration's post-LHC Run I best fit points from the CMSSM and NUHM 1 (hereafter, "NUHM 1" one will simply be referred to as "NUHM") [30]. The MasterCode analysis also includes constraints from dark matter direct detection experiments and the observed dark matter abundance. The MasterCode CMSSM best fit point will be denoted CMSSM 1. The final CMSSM point, denoted CMSSM 2, is inspired by theτ coannihilation benchmark point from ref. [26].
Both the MasterCode best fit NUHM point and the CMSSM 2 point are in coannihilation regions of parameter space, the former by virtue of having nearly pure higgsino dark matter, and are therefore extremely sensitive to variations in the RGE running. Significant variations in the running can occur between different spectrum calculators as well as between different versions of the same spectrum calculator. For example, a point that yields the correct dark matter abundance viaτ coannihilation may end up with aτ LSP if a different calculator or version is employed. Since the publication of ref. [30] and ref. [26], there have been several updates to SOFTSUSY, which was used to calculate the sparticle spectrum in both studies. Here, we consider two NUHM points inspired by the best fit point in ref. [30], denoted "NUHM A" and "NUHM B", chosen with the requirement that a valid relic density would be achieved by NUHM A via our SPheno pipelines and NUHM B via our SOFTSUSY pipelines. The original MasterCode NUHM point is included in table 1 for reference. Furthermore, the originalτ coannihilation benchmark from ref. [26], calculated with SOFTSUSY 3.3.7, yields aτ LSP in the more contemporary version SOFTSUSY 3.7.3. As such, we consider a similar point where m 0 has been increased by about 20 GeV over the original value to avoid aτ LSP.
As demonstrated in figure 1, each pipeline is composed of 3 calculators, one of each type -spectrum, Higgs, and dark matter. In this way, we consider 8 different pipelines, such as SOFTSUSY-FeynHiggs-micrOMEGAs or SOFTSUSY-SusyHD-DarkSUSY. The inclusion of the Higgs calculator is to ensure that details for the Higgs sector are achieved before computing dark matter observables. Files in SLHA [68] format are used to pass information between each calculator, with the input and output being retained at each stage. We note -6 -JHEP05(2018)113 that two separate input files are necessary for SOFTSUSY and SPheno, as there are minor differences in the expected format of the SLHA input files for the two calculators.
An important caveat in the passage of information between the programs is the handling of the branching ratios. While SLHA formatted files do include blocks for detailing particle decays, they are not universally utilized by all spectrum calculators. For example, SPheno does write the decay blocks for its SLHA output files, while SOFTSUSY does not. For SOFTSUSY pipelines, if FeynHiggs is used, Higgs decay widths will be written, but if SusyHD is used, since it only calculates the CP -even Higgs mass, no Higgs decay widths will be recorded. This means that there are no recorded widths in the SOFTSUSY-SusyHD pipelines, which can lead to discrepancies in the calculation of the dark matter abundance, for example, if dark matter annihilates primarily via the psuedoscalar resonance.
For our analysis, the versions of the calculators implemented (unless otherwise noted) are SOFTSUSY 3.7.3, SPheno 3.3.8, FeynHiggs 2.12.0, SusyHD 1.0.2, micrOMEGAs 4.3.1f, and DarkSUSY 5.1.2. Since all of the calculators studied here are continuously updated and improved, specifically since the publication of refs. [26,30,54], we also include the versions of the pipelines used in the Snowmass and MasterCode studies for proper comparisons with their results, as summarized in table 2. The Snowmass pipeline uses SOFTSUSY 3.1.7 and micrOMEGAs 2.4 and is denoted as "Snowmass", and the updated pipeline (still without FeynHiggs) is denoted as "Snowmass * ," i.e. SOFTSUSY 3.7.3 and micrOMEGAs 4.3.1f. Alternatively, the MasterCode pipeline utilized 1 SOFTSUSY 3.3.9, FeynHiggs 2.10.0, and micrOMEGAs 3.2. Since the updated MasterCode pipeline (MasterCode * ) is identical to that of our SOFTSUSY-FeynHiggs-micrOMEGAs pipeline we do not denote it separately from here forward.
Furthermore we denote pipelines with a dagger (Snowmass † /MasterCode † ) to indicate when we have reproduced the calculation of the original pipeline; otherwise the result is quoted as published. That said, not all versions of micrOMEGAs are currently available, in which case we use the closest available version.
As will be discussed below, there can be substantial variation in results and calculational techniques between different versions of the same software package. Indeed, the version numbers are critical to the interpretation of the results presented here. In the remainder of this paper, however, for the sake of brevity, we will suppress the version numbers for the packages that compose the pipelines unless otherwise specified, and refer the reader to figure 1 and table 2.

Details of the calculations
Here we discuss the details of the calculations performed by each software package. In particular we focus on the contrasting choices underlying the differences between packages, taking each tier of the pipeline in turn. Unless otherwise specified, we take the default settings for each calculator throughout the following analysis.
SOFTSUSY's optional modes are detailed in parentheses: a "high order mode" for 2-loop radiative corrections to the squark and gluino pole masses [70] and a "high accuracy mode" for 3-loop RGEs (requires CLN and GiNaC interfaces) [71]. The default modes are used in our analysis.

Spectrum calculators
For the evaluation of the sparticle mass spectrum, we consider SOFTSUSY and SPheno. First, they evaluate the gauge and Yukawa couplings at the electroweak scale before running them to the high scale and applying the soft SUSY breaking boundary conditions. After running back down to the electroweak scale, m Z , the initial tree-level values of the sparticle and Higgs masses are calculated. These masses are used as input for the iterative loop that comprises the calculation. In the iterative step, the current mass spectrum is evolved to a high scale M x , defined for GUT models to be the scale at which g 1 (M x ) = g 2 (M x ) and defined for the pMSSM to be some low scale near M SUSY = √ mt 1 mt 2 , where the soft SUSY breaking parameters are set from the specified boundary conditions. GUT models are then evolved down to M SUSY , where the two cases proceed in the same manner. At M SUSY , electroweak boundary conditions are applied and the sparticle and Higgs pole masses are evaluated at the loop level. These are now input for the next iteration of the loopbeginning with a new, and more accurate, calculation of the gauge and Yukawa couplings. When a stable solution of a given accuracy is reached, the iteration terminates and the spectrum is run down to m Z . The programs, however, do differ in the details of their calculations, as summarized in table 3 (see also table 1 from reference [45]). It is worth reminding that, even when the quoted loop level is the same, the scheme in which the calculation is handled can lead to important differences. This manifests in the choice between M S and DR schemes, where the latter amounts to a higher order correction to the former. While SOFTSUSY and SPheno both employ running DR masses in their calculations, their methods of calculating the DR corrections are different.

JHEP05(2018)113
This important difference enters in determination of the Yukawa couplings, which are calculated from the quark masses at scale Q: In running to the high scale, the quark masses must be shifted from M S to DR, and both programs ultimately follow reference [72] for the 2-loop QCD corrections and reference [73] for the 1-loop SUSY contributions. For the bottom mass the DR value is arrived at in both SOFTSUSY and SPheno by and then resummed with the SUSY corrections via .

(2.6)
For the top mass, SOFTSUSY employs a similar correction, with the 2-loop QCD corrections as where L = ln(m 2 t (m Z )/m 2 Z ) for the top mass. But SPheno uses a modified α 2 s term according to the large quark mass expansion in ref. [74], which results in an α 2 s coefficient of

Higgs calculators
The two Higgs mass calculators studied here are FeynHiggs and SusyHD. Prior to FeynHiggs 2.11.3, FeynHiggs consistently yielded SM Higgs masses ∼ 2 − 4 GeV above the value yielded by SusyHD [41]. We found that the differences between the two Higgs sector calculators were enough to present small but noticeable differences in the dark matter observables, particularly for A-funnel points. However, as of FeynHiggs 2.12.0, the differences in the results from FeynHiggs and SusyHD are far better understood, and it is possible to choose flags in FeynHiggs such that the numerical discrepancies are dramatically reduced. 2 Most notable are two changes that yielded large shifts [112,113]. The first change is the inclusion of next-to-next-to-leading order (NNLO) m M S t , which induces a downward shift in m h by as much as ∼ 2 GeV relative to when the NLO m M S t is used. The second change is the inclusion of electroweak contributions in evaluating m M S t , accounting for a downward shift of about 1 GeV. 2 We note that the hybrid approach employed by FeynHiggs has also been analytically compared to results from pure effective field theory, as employed by SusyHD, in ref. [114], which sheds light on the differences between the two approaches. For the versions employed in our primary pipelines (displayed in figure 1), FeynHiggs and SusyHD are in close agreement when FeynHiggs includes the resummation of large logs at the 2-loop level; note that this is not the default mode of the calculation, but it is used in this study. For this reason, we will focus on pipelines that include FeynHiggs in the remainder of our analysis, though results from the SusyHD pipelines are included in all tables in the appendix. We stress that shifts in the value for m h are significant not just for the calculation of the dark matter observables, but also because they introduce important caveats to previous analyses where the Higgs mass is germane.

Dark matter calculators
The two dark matter calculators we consider are DarkSUSY and micrOMEGAs. Both are sophisticated programs that analyze dark matter observables and relevant collider observables (eg. b → sγ). We confine our interest to three astrophysical observables: the neutralino relic density, Ωh 2 ; the annihilation rate today, σv ; and the spin-independent (SI) scattering cross section with nuclei, σ SI .
The relic density is calculated by micrOMEGAs according to the relation where Y 0 is the abundance of dark matter today. The same relation is used by DarkSUSY except that the numerical factor is 0.5% larger. The aim of both codes is thus to calculate the abundance of dark matter at the current temperature Y 0 ≡ Y (T 0 ), where the abundance is defined as the ratio of the number density and entropy density of dark matter Y = n/s. Both programs start with the Boltzmann equation [75] and follow reference [64] to write the differential equation as such that where the temperature has been swapped for the dimensionless quantity X = T /mχ0 1 and M P l is the Planck mass. Y eq is the thermal equilibrium abundance, and is expressed as where h eff is the number of effective degrees of freedom in the entropy density and K n is a Bessel function of the second kind. The parameter g * is related to the number of degrees of freedom of the system as where g eff is the number of effective degrees of freedom in the energy density. g eff and h eff are drawn from hard-coded tables in both programs. In micrOMEGAs, the tables come from -10 -

JHEP05(2018)113
Olive et al. [76,77] by default, but there is an option to use the tables from Hindmarsh & Philipsen [78], which is the default in DarkSUSY. The effective degrees of freedom are calculated from their respective extensive quantities, typically by assuming an ideal gas as was done in Olive et al. However, interactions will be significant at in the early universe where temperatures are high, and allowing for interactions requires a modification to the equation of state. As the weak corrections will be suppressed by the W and Z masses, QCD corrections to the effective degrees of freedom will be dominant. Hindmarsh & Philipsen found that incorporating QCD corrections allows for as much as a 3.5% modification to the relic density [78]. From here, the two programs diverge in their treatment of equation (2.9), as care must be taken due to the stiffness of the ODE [42,43]. To calculate Y (T 0 ), micrOMEGAs employs the freeze-out approximation 3 [81,82], writing Letting ∆Y (X f 1 ) = δ Y eq (X f 1 ) where δ is a small number chosen to be 1.5, micrOMEGAs solves for X f 1 . This point is used as the starting point for the numerical evaluation of equation (2.9) via the Runge-Kutta method, stopping at a point X f 2 . This latter point is chosen such that Y eq (X f 2 ) < 0.01Y (X f 2 ), and allows for the integration of equation (2.9) to solve for Y 0 : Because T 0 = 2.725 K , X 0 ∼ 10 14 and micrOMEGAs takes the upper bound to be effectively infinity. Alternatively, DarkSUSY chooses to solve equation (2.9) without applying approximations. Stiffness is still a concern, so DarkSUSY solves the problem by first discretizing the function with trapezoids and then numerically solving the differential equation with an adaptive step-size approach to Euler's method. For the computation of the (co-)annihilation of sparticles contributing to the relic density, both micrOMEGAs and DarkSUSY follow reference [65]. Furthermore, both codes include all 2-body processes 4 between neutralinos, charginos, sneutrinos, sleptons, and squarks. micrOMEGAs includes processes with gluons, as well. External programs are incorporated into the distributions of DarkSUSY and micrOMEGAs to calculate the relevant cross sections. micrOMEGAs includes CalcHEP [84] for the evaluation of relevant tree-level annihilation and coannihilation diagrams at run-time for a given model. However, some processes are JHEP05(2018)113 significantly suppressed so as to be insignificant. micrOMEGAs calculates the Boltzmann suppression factor for each channel and, by default, neglects channels where B f < 10 −6 , though this cut-off can be changed by the user [61].
Within DarkSUSY, on the other hand, the programs REDUCE [85] and FORM [86] are used in the evaluation of annihilation and coannihilation cross sections. REDUCE is used for the evaluation of helicity amplitudes for all processes between charginos and neutralinos. This allows for the analytical determination of one type of diagram only once with a numerical sum over all initial and final states performed for the contributing diagrams afterwards. All other processes involving sfermions have their scattering amplitudes evaluated by FORM.
One interesting difference between micrOMEGAs and DarkSUSY is in the inclusion of internal bremsstrahlung (IB). micrOMEGAs includes final states with two SM particles plus an additional photon for the evaluation of the annihilation cross section both in the early Universe and today. While DarkSUSY includes IB when considering the gamma ray signature from annihilation in our Galaxy's halo, it is not incorporated by default in their calculation of the relic density. As we will see below, this will lead to significant differences in the dark matter observables.
Next, we consider the differences in approaches used to calculate the SI scattering cross sections. Both programs utilize loop corrections to the scattering amplitudes but follow different frameworks. DarkSUSY follows the effective Lagrangian framework laid out in ref. [67], while micrOMEGAs utilizes the framework of ref. [87]. As discussed, for example, in ref. [87], the effective Lagrangian framework can miss crucial QCD effects, though, with modification, it is capable of reliably reproducing the 1-loop result for most cases.
Another important difference between DarkSUSY and micrOMEGAs is in the nucleon form factors used to calculate the expected neutralino-nucleon elastic scattering cross section. The form factors for each calculator are tabulated in table 4. For all quarks, DarkSUSY uses larger form factors than micrOMEGAs. As discussed below, this leads to a difference in the predicted SI neutralino-nucleon scattering cross sections, however, as we will demonstrate, these differences alone are not enough to fully explain the discrepancies in the predictions. It is clear that the details of the loop corrections also play an important role in the calculation, and can provide significant contributions to the SI cross sections.
A final difference exists in how each calculator determines the relevant particle widths used in the calculations. When available, micrOMEGAs reads in the decay blocks from the SLHA input file, but otherwise employs their own calculation to find any necessary widths. Alternatively, DarkSUSY does not currently read SLHA decay blocks and always performs their own evaluation of the relevant particle widths, though they note that future versions of their SLHA reader should include this functionality. 5 The importance of including accurate widths in the relic density calculation is particularly relevant for funnel points; previous studies have found O(10%) difference in the calculation of the A-funnel  Table 4. Spin-independent, scattering form-factors used by micrOMEGAs and DarkSUSY for protons (f p T ) and neutrons (f n T ).
between micrOMEGAs and DarkSUSY [26]. As we will demonstrate, this difference in how the width is included can lead to seemingly inconsistent results among pipelines that end with micrOMEGAs, while results from DarkSUSY pipelines seem more consistent but are less robust. Finally, we would like to mention the important effect of Sommerfeld enhancement. In short, this is an effect from non-relativistic quantum mechanics where scattering is modified by the presence of a potential interaction between the two scattering particles. This enhancement to the scattering cross section goes as 1/v, and has been found to provide substantial enhancements to dark matter annihilation when m χ 0 2 TeV (e.g. see references [88][89][90][91]). This is of particular interest to high mass wino dark matter, where the presence of a highly degenerate wino-like chargino will decrease the relic density [92,93]. Higgsino-like dark matter would also experience Sommerfeld enhancement, thanks to a triple mass degeneracy betweenχ 0 1 ,χ 0 2 , andχ ± 1 [94]. However, higgsino-like dark matter already exhibits efficient s-channel annihilation, reducing the importance of Sommerfeld enhancement [93]. Of the benchmark models we study herein, we expect the pMSSM pureh to be the only model where Sommerfeld enhancement would play a role, albeit a small one as mχ0 1 is barely above 1 TeV. For this reason, and since neither micrOMEGAs nor DarkSUSY include a calculation of the Sommerfeld enhancement, we do not consider it further. We note that for a heavier wino-like dark matter candidate, it would be interesting and important to explore the theoretical uncertainties in the calculation of dark matter observables with the inclusion of the Sommerfeld effect.

The physics of neutralino dark matter benchmarks
In this section, we provide a brief introduction to aspects of the physics of neutralino dark matter that will be relevant for our analysis of the benchmark points.
The challenge in obtaining the correct relic density of neutralino dark matter observed in the Universe is well-known. 6 Specifically, one generically obtains an annihilation cross section at freeze-out that is too small, leading to too much neutralino dark matter in the present epoch.
The general idea of WIMP dark matter is that, in order to produce the correct relic density, the annihilation cross section of dark matter in the early Universe should have been JHEP05(2018)113 around 1 pb. This is approximately the value one obtains if dark matter is a new particle with approximately weak-scale mass and with electroweak couplings, a happy accident called the "WIMP miracle".
However, as has been long appreciated, the details of actual SUSY models are somewhat less attractive than the idea sketched above. Although supersymmetry predicts that the annihilation cross section of two neutralinos should be in the neighborhood of 1 pb, the exact numerical value can span a range that covers orders of magnitude, depending on the composition of the neutralino and the mass spectrum of the other supersymmetric particles. Dark matter that is predominantly higgsino-like (h) annihilates to W and Z bosons with a cross section that involves the full strength of the SU(2) gauge coupling, and is moreover enhanced by the presence of spin-one particles in the final state. Higgsino and wino dark matter thus have cross sections that are too large for the observed relic density. On the other hand, binos (B) mainly annihilate to quark and lepton pairs, a process that suffers from helicity suppression. Binos therefore typically have a cross section for annihilation that is too low.
The regions of supersymmetric parameter space that are compatible with the dark matter relic density thus tend to be fine-tuned. In the following subsections, we will consider the most well-studied regions, characterized by how the observed dark matter relic density is achieved: coannihilation models, well-tempered dark matter, A-funnel annihilation, and pure higgsino composition.

Coannihilation ofB with light scalars
The calculation of the relic density for the coannihilation region is very sensitive to the relative masses of the dark matter candidate and the light scalar(s) that coannihilate with it, as detailed eg. in [75]. Even small discrepancies in the mass spectrum given by different spectrum generators will result in significantly different predictions for the relic density.
Regarding the indirect detection prospects for bino-scalar coannihilation models, the annihilation cross section for bino dark matter in the present Universe occurs mainly through t-channel exchange of light scalars. Bino-nucleon scattering generally proceeds via Higgs-or squark-exchange. The Higgs exchange diagram is suppressed for models in which bino-squark coannihilations dominate in the early Universe, due to the pure bino nature of the dark matter in these cases. Furthermore, if the first and second generation squarks are heavy, as is the case in our bino-stop coannihilation scenario, the squark-exchange diagram is suppressed, resulting in very small scattering cross sections with nuclei, as we will see.
In the following analysis, three coannihilation scenarios will be relevant:B −τ 1 ,Bt 1 ; andB-q, whereq denotes any first or second generation squark. In each case, the composition of the neutralino LSP is 99.9% bino.

Well-tempering of dark matter
Well-tempered dark matter has been extensively studied in the context of recent direct, indirect, and collider searches [97][98][99][100]. The annihilation cross section depends on the mixture of bino and higgsino states, or, equivalently, on the higgsino fraction |N 13 | 2 +|N 14 | 2 . The full expression for the annihilation cross section, as well as various interesting limits, are -14 -

JHEP05(2018)113
available for example in reference [101]. The predicted relic density and indirect detection signals in the current Universe are both sensitive to the higgsino fraction.
The dominant neutralino-nucleon scattering occurs via CP-even Higgs exchange, which is a bino-higgsino-Higgs coupling. Since a well-tempered neutralino has sizable higgsino and bino fractions, the scattering rates with nuclei can be relatively large.

A-funnel annihilation
A neutralino LSP can annihilate resonantly by exchanging a pseudoscalar Higgs boson A in the s-channel, provided m A ∼ 2mχ0 1 . Prospects of probing the A-funnel at colliders and the correlation with the observed Higgs mass have been extensively studied [102][103][104][105][106]. Prospects for direct [107,108] and indirect [109] detection have also been studied.
The annihilation cross section can be expressed as [110] σv ∼ 3 2π is the Mandelstam variable, Γ A = 4.6 GeV is the width of the pseudoscalar at our benchmark, m A = 2042 GeV is its mass, and the Yukawa couplings to the b-quarks and neutralinos are given by To obtain the annihilation cross section, one must take a thermal average of eq. (3.1) or its more general equivalent. The resonance occurs in the zero velocity limit in the current Universe, and the annihilation cross section today is given by σv| v→0 . In the early Universe, the non-zero velocity leads to thermal broadening of the resonance, as is clear from the velocity dependence of the Mandelstam variable, eq. (3.2). Small differences in the calculation of m A and v by different generators can affect the relic density significantly. Note that from eq. (3.3), it is clear that to obtain the correct relic density, the A-funnel benchmark requires non-zero values of N 11 and either N 13 or N 14 . Thus, though the dark matter is typically primarily bino, a higgsino component must be retained, i.e. the neutralinoχ 0 1 must have some bino-higgsino mixture. In this case, the SI scattering cross section with nuclei is mediated primarily by Higgs exchange. This is enabled by the non-zero higgsino component N 13 or N 14 required to obtain the correct relic density. However, the values of the higgsino fraction required to obtain the correct relic density are typically very small, corresponding to feeble scattering cross sections with nuclei. The direct detection prospects for the A-funnel scenario are thus rather challenging, as we will see. This should be contrasted to the case of the welltempered neutralino, where the large higgsino fraction drives both the relic density and the scattering cross section.  Table 5. Masses, in GeV, of the lightest and second-lightest neutralinos and the lighter chargino for the pMSSM benchmark points. The corresponding higgsino fraction is given in parentheses. The percent differences are given relative to the SOFTSUSY values.

Pure Higgsino (h) composition
Higgsinos with mass of ∼ 1 TeV can satisfy the relic density constraint, with the dominant mechanism being annihilation to gauge bosons. Coannihilation among charged (χ ± 1 ) and neutral (χ 0 1 andχ 0 2 ) higgsinos is also important in this case. We refer to [111] for expressions of the relic density in various limits.
The contribution of Higgs exchange diagrams to the scattering cross section with nuclei is suppressed due to the small gaugino-higgsino mixing. The contribution of squark exchange diagrams is also suppressed at our benchmark point since the squark masses are at several TeV. Thus, we generally expect small direct detection signals for this benchmark point. We refer to [97,98] for detailed calculations of the scattering cross section of pure higgsinos.

Results: pMSSM analysis
Here we present our pMSSM analysis. We consider the five pMSSM points from the Snowmass 2013 benchmarks [54] as discussed in section 2.1.
The spectra in the neutralino, squark, and Higgs sectors obtained using the different spectrum generators we consider are displayed in tables 5, 6, and 7, respectively. For tables 5 and 6, the first column lists the dark matter benchmark scenario, while the next columns display the spectra of relevant sparticles obtained from SOFTSUSY and SPheno.  figure 1. We pause to elaborate on the Higgs sector, before moving on to analyze the dark matter results. In tables 7a & 7c we show the masses in the Higgs sector for all the dark matter benchmark scenarios studied. The masses corresponding to table 7a & 7b are computed by FeynHiggs, while those corresponding to tables 7c & 7d are computed by SusyHD. As per convention, Snowmass * is the updated Snowmass pipeline which amounts to the SOFTSUSY spectrum prior to FeynHiggs, here. We note that SusyHD only corrects m h and, thus, m H , m A , and m H ± are identical to those calculated by the relevant spectrum calculator.
From these results, we see that the Higgs sector masses are not significantly affected by the choice of SUSY spectrum generator, with differences amounting to less that 0.01%. However, the inclusion of either FeynHiggs or SusyHD can provide a significant shift in m h from the Higgs mass calculated by the spectrum generator itself. The Snowmass points were constrained by the Higgs mass (126 ±1 GeV), but utilizing FeynHiggs moves the mass down by as much as 2 GeV (1.5%) which puts the benchmark points in strong tension with the measured value of the Higgs mass (125 ± 0.24 GeV).    in table 12). For comparison, we also include the Planck 3-sigma range [115] for the relic density, which is highlighted by the red band. The upper right panel of figure 2 shows the annihilation cross section today (tabulated in table 14), with the Fermi-LAT 6-year limits from dwarf spheroidal galaxies [116] included for comparison (limits on annihilation to τ + τ − in red and bb in blue). In the lower left panel of figure 2 we show the spin-independent neutralino-nucleon elastic scattering cross section, where we plot the per-nucleon cross section averaged for Xe, computed from the proton and neutron values tabulated in table 16 (FeynHiggs pipelines), as well as exclusion contours from LUX [117] in solid red, PandaX [118] in solid blue, and LZ (projected) [119] in dashed black lines, for comparison.
We now turn to a description of our results. Throughout our discussion, we will concentrate on the SOFTSUSY (magenta) and SPheno (cyan) pipelines. In the figures, we also plot the predictions from the Snowmass * and Snowmass † pipelines shown in solid green and solid black, respectively. Since the physics of these pipelines is very similar to the SOFTSUSY -micrOMEGAs pipeline shown by solid magenta shapes, we will not discuss them separately.

Bino-stop coannihilation
The masses relevant for coannihilation are the lightest neutralino and the stop. From table 5, we see that there is good agreement between SOFTSUSY and SPheno in the neutralino spectrum for these two benchmark scenarios. The dark matter has a mass of ∼ 754 GeV, with very good agreement between the two spectrum generators, and is almost completely bino. On the other hand, from table 6, we see that the stop mass differs by 2.7% between SOFTSUSY and SPheno, although it is generally in the range where coannihilation is operational. This has a significant effect on the relic density, which, in these cases, is exponentially sensitive to the mass difference between the dark matter and the relevant coannihilation partner.
The effect of the variation in the mass of the coannihilation partner can be seen in the upper left panel of figure 2. The solid magenta and cyan diamonds correspond to the relic density values computed by micrOMEGAs, for spectra coming from SOFTSUSY and SPheno, respectively. While SOFTSUSY yields a value of Ωh 2 = 0.094, SPheno yields a value of Ωh 2 = 0.035, and the difference stems entirely from the difference in stop masses computed by the two generators. Similarly, comparing the hollow magenta and cyan diamonds, which correspond to the relic density values computed by DarkSUSY, we see that while SOFTSUSY gives a value of Ωh 2 = 0.120, SPheno yields a value of Ωh 2 = 0.045.
It is also interesting to compare the values yielded by micrOMEGAs and DarkSUSY for the same spectrum. For example, selecting the spectrum from SOFTSUSY and comparing the solid and hollow magenta diamonds, we see that micrOMEGAs gives Ωh 2 = 0.094 while -19 -

JHEP05(2018)113
DarkSUSY gives Ωh 2 = 0.120. We point out that these theoretical uncertainties exceed the current experimental uncertainty, which, here results in the hollow diamond lying within the Planck-allowed band, while the solid diamond does not. These discrepancies occur due to differences in the calculation of the effective cross section for each annihilation and coannihilation channel and the different relative weights of contributing final states in the coannihilation channels assigned by the calculators.
We note that both micrOMEGAs and DarkSUSY give values for the coannihilation cross sections using an effective tree-level calculation. These values can be significantly altered if higher-order SUSY-QCD corrections are taken into account. After including loop diagrams containing a gluon, a gluino, a four-squark vertex, and incorporating gluon radiation, the authors of reference [120] have found a ∼ 20% discrepancy with the relic density computed by micrOMEGAs in the bino-stop coannihilation region. We will not consider these loop corrections further in this paper, but note that global ∼ 20% theoretical uncertainties are expected in all cases.
The annihilation cross section of the bino-stop benchmark model in the current Universe is shown in the upper right panel of figure 2. There is remarkable agreement among the various pipelines. Since coannihilation channels are irrelevant in the present Universe, the sensitivity to the mass difference between the stop and the bino is absent. The annihilation proceeds mainly through the t-channel exchange of a stop, and the small difference in the stop mass reported by SOFTSUSY and SPheno does not affect this diagram as significantly. We find that both micrOMEGAs and DarkSUSY give similar values for the different channels inχ 0 1χ 0 1 annihilation, the only difference being that DarkSUSY ascribes ∼ 5% contribution to the annihilation cross section fromχ 0 1χ 0 1 → gg final state, while micrOMEGAs does not return this channel.
The SI scattering cross section is shown by the diamonds in the lower left panel of figure 2. For very heavy squarks, the leading scattering cross section is mediated by Higgs exchange, which is suppressed if the lightest neutralino is a pure higgsino or gauge state, as discussed in section 3. Thus, the combination of heavy squarks and pure bino eigenstate conspire to give suppressed scattering cross sections for the bino-stop coannihilation benchmark. The cross sections lie below the projected LZ limits, at approximately 10 −11 pb.
There is a factor of ∼ 5 discrepancy between the SI scattering cross section yielded by SOFTSUSY -micrOMEGAs (solid magenta diamond) relative to SPheno-micrOMEGAs (solid cyan diamond). From tables 5, 6, and 7, we can see that while the higgsino fraction and the Higgs mass match to a high approximation for both SOFTSUSY and SPheno in the binostop coannihilation benchmark, there can be up to a 4 GeV difference in the masses of the squarks. While it is unlikely that this small variation in squark masses can account for the observed variation in the SI scattering cross section, we cannot pinpoint the exact source for the discrepancy. Comparing the dark matter calculators, we see that the solid and hollow magenta diamonds overlap entirely, meaning that after receiving the spectrum from SOFTSUSY, both micrOMEGAs and DarkSUSY computed the same scattering cross section. The matching is not quite as exact for the spectrum coming from SPheno (solid and hollow cyan diamonds), but the values are quite close. While DarkSUSY implements an effective Lagrangian in the heavy squark limit following [110] (see reference [67] and references therein for details), micrOMEGAs implements the full one-loop Lagrangian following reference [87].    [115] is highlighted in red in the upper left panel, the limits on the annihilation cross section today from Fermi-LAT's 6-year analysis of dwarf spheroidal galaxies [116] for annihilation to τ + τ − (red) and bb (blue) are shown in the upper right panel, and exclusion limits from LUX (red) [117], PandaX (blue) [118], and LZ (projected; black) [119] are shown in the lower left panel.

Bino-squark coannihilation
The physics for this benchmark is similar to that of the bino-stop case. Although there is good agreement in the neutralino spectrum, we see from table 6 that squarks can differ by up to 1.51% between SOFTSUSY and SPheno. The effect on the relic density can be seen in the upper left panel of figure 2. The solid magenta and cyan circles correspond to the relic density values computed by micrOMEGAs for spectra coming from SOFTSUSY and SPheno, respectively. While SOFTSUSY yields a value of Ωh 2 = 0.087, SPheno yields a value of Ωh 2 = 0.110. The difference can be attributed to the difference in squark masses. The fact that SOFTSUSY produces a lower value for the relic density than SPheno is due to the fact that squark masses from SOFTSUSY are lighter than those from SPheno, giving a stronger coannihilation effect. We also note that for a given spectrum generator, say SOFTSUSY, DarkSUSY gives a larger value of the relic density than micrOMEGAs.

JHEP05(2018)113
The annihilation cross section in the current Universe is shown by the circles in the upper right panel of figure 2. It is evident that there is much more convergence of results among the various pipelines compared to the relic density computation. This can again be ascribed to the fact that coannihilation channels are absent in the current Universe and small differences in squark masses do not translate to large differences in the calculation of the t-channel squark exchange diagram.
Comparing the diamonds (stop coannihilation) and circles (squark coannihilation) in the upper right panel of figure 2, we can see that the bino-squark coannihilation points exhibit greater spread in their current annihilation cross sections. This is due to the fact that the cross sections are driven by t-channel stop or squark exchange diagrams in the present Universe, and the different pipelines give a greater spread in the squark masses than they do the stop mass for the bino-stop benchmark.
The SI scattering cross section is shown by the circles in the lower left panel of figure 2. As expected, the values are higher than for the bino-stop coannihilation case, due to the lighter squarks which make the leading contributions to the squark-exchange diagram. However, there is also greater disagreement between the different pipelines. Firstly, comparing the solid magenta circle with the solid cyan circle, we see that the cross section computed with the spectrum coming from SOFTSUSY is a factor ∼ 3 greater than that coming from SPheno. The same trend can be seen using DarkSUSY (comparing the hollow magenta circle with the hollow cyan circle), although the effect is smaller.
For a given spectrum generator, it is clear that micrOMEGAs is giving larger values of the scattering cross section than DarkSUSY, since the solid circles are above the hollow ones. This is partly due to differences in the form factors used by the two calculators (see table 4). Using the form factors of DarkSUSY in micrOMEGAs, we find that the scattering cross sections reported by micrOMEGAs reduces by ∼ 30%, bringing the two calculators to greater agreement with each other.

Pure higgsino
The results for the relic density, current annihilation cross section, and scattering cross section of the pure higgsino benchmark are shown by stars in the upper left, upper right, and lower left panels of figure 2, respectively. For the relic density there is good agreement between the different pipelines. The two spectrum calculators produce similar mass and higgsino fraction of the lightest neutralino for this benchmark, as is evident from table 5. Thus, the magenta and cyan stars for a fixed dark matter calculator overlap in the upper left panel of figure 2. There is some discrepancy between the results returned by micrOMEGAs and DarkSUSY due to differences in calculation of coannihilation processes between the neutral and charged higgsinos. Since coannihilation is unimportant in the current Universe, these discrepancies disappear and all stars align perfectly in the upper right panel of figure 2.
For the SI scattering cross section, we notice that the rates for a pure higgsino are suppressed due to the small gaugino-higgsino mixing, similar to the case of a pure bino. As for the relic density, the main difference between pipelines comes from the choice of the dark matter calculator and whether an effective Lagrangian is used or the full one-loop Lagrangian is considered.

Well-tempered neutralino
The results for the relic density, current annihilation cross section, and SI scattering cross section of the well-tempered neutralino benchmark are shown by triangles in the upper left, upper right, and lower left panels of figure 2, respectively.
We first discuss the relic density. The upper left panel of figure 2 shows that there are considerable differences between the various pipelines. The SPheno pipelines (cyan triangles) give a larger relic density than the SOFTSUSY pipelines (magenta triangles). This can be traced to the lower higgsino content of the lightest neutralino given by SPheno relative to SOFTSUSY. In fact, from table 5, the lightest neutralino is 30% higgsino when calculated by SPheno, but 36.4% higgsino when calculated by SOFTSUSY, although the masses agree to better than a percent.
On the other hand, for a given spectrum generator, there is almost no discrepancy in the relic abundance coming from the dark matter calculators. Thus, the solid and hollow triangles approximately overlap for a given color. What small discrepancy that is evident can be attributed to differences in the computation of the effective annihilation and coannihilation cross sections between the bino, neutral higgsino, and charged higgsinos, as well as the computation of the t-channel chargino exchange diagram. The annihilation cross section in the present Universe presents far fewer differences, since the coannihilation channels are absent. The triangles thus overlap in the upper right panel of figure 2.
For the SI scattering cross section with nuclei, we see that the well-tempered benchmark is the only one of our pMSSM benchmarks that is constrained by current experiments, regardless of the pipeline adopted. This is because of the non-negligible higgsino fraction. The different higgsino fractions reported by SOFTSUSY and SPheno affect the relative positions of the magenta and cyan triangles in the lower left panel of figure 2. Since SOFTSUSY reports the larger higgsino fraction, the scattering cross section for magenta triangles are larger than those for the cyan triangles corresponding to the SPheno pipelines.
For a given spectrum generator, say SOFTSUSY, there is some disagreement between the SI scattering cross sections computed by micrOMEGAs (solid magenta triangle) versus DarkSUSY (hollow magenta triangle). This can be attributed to the different form factors, as well as the way in which loop effects are incorporated. To check the effect of the different form factors, we used DarkSUSY's form factors from table 4 in the micrOMEGAs code. For most points, doing so brings micrOMEGAs's values for the scattering cross sections into better agreement with those of DarkSUSY. However, some differences remain, suggesting that other details of the calculation are also important. We also carried out the tree level calculation for SI scattering, following the discussion of appendix C in reference [121] 7 and using the values in table 4. Surprisingly, we found O(10%) differences between our treelevel calculation and that of the codes', implying substantial loop contributions that are different between the two codes.

A-funnel
The results for the relic density, current annihilation cross section, and scattering cross section of the A-funnel neutralino benchmark are shown by pentagons in the upper left, upper right, and lower left panels of figure 2, respectively.
We first discuss the relic density. From the upper left panel of figure 2 and table 12, it is evident that there is a large variation, of more than factor of two, among the different pipelines, though the variation is 2 if one neglects the Snowmass and Snowmass* pipelines. We note that for a given spectrum generator, say SOFTSUSY, there is some difference in the calculation performed by micrOMEGAs (solid magenta pentagon) versus DarkSUSY (hollow magenta pentagon). This difference amounts to an uncertainty of 69% and 37% for SPheno and SOFTSUSY, respectively, both of which far exceed the experimental uncertainty in the measurement of the dark matter abundance.
The resonance region is notoriously sensitive to the approximations used to compute the relic density, especially for sharp peaks. The sharpness parameter in our case (following the same notation as reference [75]) is Reference [75] compared various approximation schemes in the relic density calculation (such as Taylor expansion in v) to a full numerical computation for values of near this value. Depending on the approximation, the relic density can vary over several orders of magnitude. Even for a full numerical computation, the resonance region is sharp enough that a factor of ∼ 2 can easily appear, unless there is an exact matching of calculation. Finally, we note that there is a significant difference between the SOFTSUSY-FeynHiggs-micrOMEGAs and SOFTSUSY-SusyHD-micrOMEGAs pipelines for the relic density of the Afunnel benchmark point, as can be seen by comparing tables 12a and 12b. This discrepancy is due to the inclusion of the pseudoscalar width, as in eq. (3.1). FeynHiggs calculates the pseudoscalar width, which is read in by micrOMEGAs (but not DarkSUSY, as discussed below). However, SusyHD does not calculate the width, so any program further down the pipeline either takes the width approximation from the spectrum generator or calculates the width itself. SPheno does estimate a width with reasonable agreement between its width and that calculated by FeynHiggs, yielding decent agreement between the SPheno-FeynHiggs-micrOMEGAs and SPheno-SusyHD-micrOMEGAs pipelines. But SOFTSUSY does not report a width that can be used by micrOMEGAs, so for the SOFTSUSY-SusyHD-micrOMEGAs pipeline, micrOMEGAs calculates its own width, which differs by nearly a factor of 1.7 from that calculated by FeynHiggs. This leads to a discrepancy of nearly a factor of 2 in the relic abundances, as can be verified, for example, by evaluating eq. 41 of ref. [75]. We do not see a discrepancy in the relic densities of any of the DarkSUSY pipelines because DarkSUSY always calculates all relevant particle widths, since, as mentioned above, up to version 5.1.2 DarkSUSY does not read in SLHA decay blocks.
The annihilation cross section today is given by an even sharper peak in σv centered at m A = 2mχ0 1 . From the upper right panel of figure 2, and comparing tables 12a and 14a, -24 -

JHEP05(2018)113
we see that the values reported by the various pipelines are in agreement with what we expect from the upper left panel; that is, there is a factor of ∼ 2 spread in the annihilation cross sections, similar to the factor of ∼ 2 spread in the relic densities (with the largest annihilation cross section corresponding to the smallest relic abundance, and so forth). We note that the current annihilation cross sections in the upper right panel of figure 2 are a factor ∼ 5 reduced compared to their values in the early Universe, used to calculate the abundances in the upper left panel of figure 2. This happens due to the thermal broadening of the resonance region in the early Universe. The SI scattering cross sections reported in the lower left panel of figure 2 match very closely due to the close agreement in the relevant masses and higgsino fractions. For a given spectrum generator, say SOFTSUSY, there is some disagreement between the scattering cross section computed by micrOMEGAs (solid magenta pentagon) versus DarkSUSY (hollow magenta pentagon). As for other benchmarks, this can be attributed to the different form factors, as in table 4 and the way in which loop effects are incorporated.

Summary: broad trends in the pMSSM analysis
It is instructive to look back at our analysis and draw some broad conclusions. The dark matter models studied in this section are a well-known subset of pMSSM benchmarks that satisfy the observed dark matter relic density. These benchmarks have been used in numerous studies, typically relying on one of the pipelines described in our work. Furthermore, connections between supersymmetric model building and cosmology often concern regions of parameter space based around one of the (co)annihilation mechanisms that our benchmarks capture.
It should be clear from our work that the theoretical uncertainty in the relic density calculation using standard pipelines far exceeds the experimental uncertainty. From the upper left panel of figure 2, it is apparent that only a minority of pipeline choices for any given benchmark actually fall within the red band that delineates the Planck range. For the coannihilation and funnel models, this spread is especially broad, with theoretical calculations yielding results that can vary by as much as 300%. The spread is somewhat lower for the well-tempered neutralino and pure higgsino benchmarks, but even in those cases it far exceeds the experimental uncertainty. We note also that there can be significant spread in the relic density due to updates to the software packages even for the same sequence of calculators.
There are several underlying reasons for the discrepancies in the relic density and other dark matter observables among the pipeline choices: • Small variations in the spectrum. -Coannihilations and A-resonance models are critically dependent on the relative masses of the dark matter particle and other light supersymmetric states. Within a pMSSM framework, the spectrum generators can easily produce 1% -2% variation in the low energy spectrum of squarks, stops, or the pseudoscalar Higgs, leading to a large variation of the relic density, evident in the upper left panel of figure 2. However, the annihilation cross section in the present Universe is far less dependent on the masses of other light states, since coannihilation • Composition of the lightest neutralino. -The well-tempered neutralino framework depends critically on the higgsino fraction of the dark matter for its relic density. We see from table 5 that the spectrum generators can vary by up to 20% in their calculation of the higgsino fraction relative to each other. Furthermore, the higgsino fraction plays a crucial role in the SI scattering cross section, most evident in the well-tempered benchmark.
• LO and NLO calculation of coannihilation channels. -For a given spectrum, there is wide variation (∼ 50%) in the relic density reported by micrOMEGAs and DarkSUSY (discrepancies between solid and hollow points of the same shape and color in the upper left panel of figure 2), especially in scenarios where coannihilation channels become important. This stems from differences in the tree level computation implemented in these programs. Moreover, as we have pointed out, the incorporation of NLO SUSY-QCD will further change the relic density calculation, by up to as much as 20%.
• Differences in form factors. -The form factors employed by micrOMEGAs and DarkSUSY are displayed in table 4. These differences affect the SI scattering cross sections reported in the lower left panel of figure 2 for a given spectrum calculator (discrepancies between solid and hollow points of the same shape and color). For a given spectrum, using the form factors of DarkSUSY in micrOMEGAs, we find a definitive shift in the SI cross sections, bringing them into closer agreement.
• NLO effects in scattering cross section. -Even accounting for the differences in spectrum and form factors, we see that different dark matter calculators report different SI scattering cross sections, especially for very low cross sections. These differences are likely coming from the fact that DarkSUSY implements an effective Lagrangian in the heavy squark limit following reference [110] (see reference [67] and references therein for details) while micrOMEGAs implements the full one-loop Lagrangian following reference [87]. For example, the pure higgsino or pure bino benchmarks in the lower left panel of figure 2 show a lot of variation. The theory calculations for pure higgsino and wino scattering cross sections have only converged recently [97,98]. The discrepancies among the pipelines is likely to become a pressing issue in the future, when experimental sensitivity reaches the relevant cross sections.

Results: GUT analysis
In the previous sections, we have presented our analysis of several standard pMSSM benchmark scenarios. For electroweak-scale models like the pMSSM, the supersymmetric mass -26 -  Table 8. Masses, in GeV, of the lightest and second-lightest neutralinos and the lighter chargino for the GUT benchmark points. The corresponding higgsino fraction is given in parentheses.

JHEP05(2018)113
spectrum is used as an input for the spectrum generators at low energies, and the final sparticle spectrum is used as an input for the dark matter calculators. Thus, there is very little running of the sparticle masses and the results reported by different spectrum generators are generally in good agreement. In this section, we will analyze dark matter benchmark points in the context of supersymmetric models with boundary conditions for soft terms specified at the GUT scale. In particular, we will study two models: CMSSM and NUHM. In this case, the effects of RG running performed by the two spectrum calculators are expected to become more important, and can substantially affect the dark matter observables. We begin this section by first discussing the sparticle spectra of the GUT benchmark models. The GUT models are defined in table 1. The low energy spectrum of neutralinos and squarks generated by SOFTSUSY and SPheno are shown in tables 8 and 9, and Higgs masses are shown in table 10. Masses presented in tables 10a & 10b are computed by FeynHiggs while those in tables 10c & 10d SusyHD. The higgsino mass parameter µ is shown in table 11. We also display in tables 8-11 the values obtained via the MasterCode † pipeline, as described in table 2.
From table 8, we see that for the CMSSM points, the lightest neutralino is mainly bino, while for the NUHM points, it is mainly higgsino. The higgsino fraction calculated by the different spectrum generators has appreciable differences in the cases of the CMSSM 1 and NUHM B, which will significantly affect the dark matter observables reported for the two pipelines, as we shall see. Even in the case of the NUHM A point, the small difference in higgsino fraction (< 1% difference between SOFTSUSY and SPheno) will be important.
There is good agreement between SOFTSUSY and SPheno for the mass of the dark matter in the CMSSM cases. However, there is a vast disagreement in the mass of dark matter -  Table 9. Masses, in GeV, of the squarks for the GUT benchmark points.
between SOFTSUSY and SPheno for both NUHM points. There are also large discrepancies in the mass of the second lightest neutralino and the charginos in all cases except CMSSM 2. These discrepancies all stem from differences in the value of µ, as is evident from table 11. From table 10, we see that there are also large discrepancies in the mass of the pseudoscalar Higgs A, which is critical for the dark matter relic density computation in the A-funnel region of parameter space. On the other hand, the squark spectrum agrees among different generators quite well, as is evident from table 9, although there can be variations of up to ∼ 2 % in the calculation of squark masses. Even these ∼ 2 % discrepancies will be important in the following analysis.
We thus see that the largest discrepancies seem to occur in the calculation of the electroweak symmetry breaking (EWSB) sector between different spectrum generators. We now turn to a more detailed study of these differences: as a prelude to our investigation of dark matter observables, we discuss in detail the differences between the calculation of the higgsino mass parameter µ and the pseudoscalar Higgs mass m A by the different spectrum generators. Finally, we discuss the dark matter observables obtained from the different pipelines.

Comparison of EWSB sectors
Here we explore the EWSB calculations performed by the different spectrum generators. The differences in the EWSB calculations are particularly exacerbated at large values of tan β and m 0 , which is the regime where our GUT benchmark models CMSSM 1 and NUHM A/B lie. We discuss these issues in this section.   Table 11. Higgsino mass parameter, µ, after FeynHiggs.

JHEP05(2018)113
The RGE's for the higgsino mass parameters in the MSSM are, following the notation of reference [126], In the above equations, in standard notation. The higgsino mass parameter is given, in the large tan β limit, by where all quantities are defined at m Z . The pseudoscalar Higgs mass is given at tree level by with all quantities defined at the scale M SU SY , and the quantities on the right hand side of eq. (5.4) related to those at m Z by radiative corrections. Obviously, the calculated value of m A depends on the computation of µ, m 2 Hu , and m 2 H d . From eq. (5.1), eq. (5.3), and eq. (5.4), it is clear that the values of µ and m A reported by the programs will be greatly affected by their calculation of the top and bottom Yukawas y t and y b , as well as the calculation of squark and stop masses that enter into X t and X b . In table 9, there are variations of ∼ 2% in the stop masses between SOFTSUSY and SPheno. In fact, SOFTSUSY consistently reports higher values of the stop masses than SPheno across benchmark models. There are also variations in y t and y b between the programs, as studied in [45]. These factors result in vastly different values of µ and m A , especially for large values of m 0 where the squark and Yukawa calculations differ substantially. Similarly, we find that τ masses are reported ∼ 30% higher by SOFTSUSY than SPheno, following from a ∼ 20% difference in the calculation of the τ Yukawa couplings.
In  Hu with increasing m 0 near the benchmark point, for both SOFTSUSY and SPheno. We have also verified this using the approximate relations for the renormalization group running given in reference [127].
On the other hand, we have found that near the benchmark, the effect of increasing m 0 is to increase m 2 H d . This is because while increasing m 0 increases the slope of the renormalization group running from eq. (5.1), this increase is suppressed compared to the m 2 Hu case by the small value of the bottom Yukawa. The cumulative effect is that the increase in the boundary value of m 2 H d for increasing m 0 dominates, so m 2 H d increases with increasing m 0 . We see from eq. (5.3) that m 2 Hu and m 2 H d contribute with opposite signs to µ. However, the dominant contribution is from m 2 Hu , since m 2 H d is suppressed by the large value of tan β. Thus, following the behavior of m 2 Hu , µ too decreases with increasing m 0 . The values of µ (solid curves) from SOFTSUSY and SPheno start to diverge radically for m 0 > 2 TeV in the left panel of figure 3. This is the regime where differences in the squark masses and the top Yukawa calculated by the two spectrum generators start to become important in determining µ. From table 9, we see that SOFTSUSY produces heavier stop masses than SPheno for the same CMSSM model point. Thus, µ runs to smaller values faster in SOFTSUSY compared to SPheno. The values are µ = 1129 GeV for SOFTSUSY and µ = 1552 GeV for SPheno at the CMSSM 1 benchmark.
From the right panel of figure 3, we see that for both programs the mass of the lightest CP-even Higgs increases with increasing m 0 . This is expected, due to the usual loop corrections to the Higgs mass. We also note that SPheno reports a slightly larger Higgs mass than SOFTSUSY due to a combination of the low energy values of the stop mass and the trilinear coupling. Finally, we point out that the Higgs mass calculation is relatively robust to uncertainties in the EWSB calculations, since m h is sensitive only to m A (not µ independently) at tree level, which is reflected in the behavior of m h at very large m 0 near where µ becomes unphysical. We note that the uncertainties in the EWSB calculations tend to cancel each other in the calculations of m h , while, in contrast, they do not cancel each other in the calculation of m H , which tracks m A quite closely.
We now move on to a discussion of the pseudoscalar Higgs mass m A (dashed curves) in the left panel of figure 3. From the tree level relation for m A in eq. (5.4), we expect that the value of m A reported will depend on the relative magnitudes of µ, m 2 Hu , and m Finally, the CMSSM 2 point, which is aτ coannihilation model, is largely insensitive to uncertainties in the EWSB sector, with the exception of those associated with the τ Yukawa coupling. The dark matter physics of the CMSSM 2 benchmark primarily concerns the LSP, which is strongly bino-like with good agreement among the pipelines, and the lighterτ , with masses of 538.29 GeV and 390.01 GeV from SOFTSUSY and SPheno, respectively. The large difference in theτ masses comes primarily from the difference in the calculation of the τ Yukawa coupling, with SOFTSUSY reporting a value ∼ 20% larger than SPheno. Since the lightest neutralino mass for the CMSSM 2 is ∼ 386 GeV for both SOFTSUSY and SPheno, this means that only SPheno pipelines represent true coannihilation models, while SOFTSUSY pipelines do not coannihilate efficiently enough to achieve the correct relic abundance.

Dark matter observables
In this section, we study the dark matter observables for the GUT benchmarks. For each of the CMSSM and NUHM benchmark models, we discuss the relic density, the annihilation cross section today, and the predicted scattering cross section, all of which are plotted in figure 5. As in figure 2, for comparison, the Planck 3-sigma range for the dark matter abundance [115] is highlighted in red in the upper left panel, the limits on the annihilation cross section today from Fermi-LAT's 6-year analysis of dwarf spheroidal galaxies [116] for annihilation to τ + τ − (red) and bb (blue) are shown in the upper right panel, and exclusion limits from LUX (red) [117], PandaX (blue) [118], and LZ (projected; black) [119] are shown in the lower left panel. Here, we introduce a new set of unique shapes to denote the various pipelines, but follow the same colour scheme as in figure 2 -black is used again for our implementation of the original pipeline (MasterCode † , rather than Snowmass † ). 8

CMSSM benchmarks
We first examine the relic density as plotted in the upper left panel of figure 5 and tabulated in table 20.
For the CMSSM 1 point, the mass of the lightest neutralino from table 8 is mχ0 1 = 936.10 GeV from SOFTSUSY, and mχ0 1 = 938.0 GeV from SPheno. The dark matter is mostly bino in both cases. It is, however, the mass of the second lightest neutralino that differs radically between the two programs. For SOFTSUSY, we have mχ0 2 = 1147.50 GeV, while for SPheno, we have mχ0 2 = 1581.69 GeV. The second lightest neutralino is mostly higgsino. The huge discrepancy in masses is due to the very different values of µ reported by the two programs, as discussed above and shown in figure 3. The large difference in µ is also re- 8 There are no green markers to indicate an updated MasterCode pipeline, since the updated MasterCode pipeline is the same as our SOFTSUSY-FeynHiggs-micrOMEGAs pipeline. 10 -8   Figure 5. GUT model results. The upper left, upper right, and lower left panels show the neutralino relic density, annihilation cross section today, and the SI neutralino-nucleon elastic scattering cross section, respectively, all as functions of the dark matter mass, with a common legend in the lower right panel. We note that relic density for some pipelines yields values far larger than those plotted here, and for this reason the CMSSM 1 SPheno points and the CMSSM 2 SOFTSUSY points do not appear in the upper left panel. For comparison, the Planck 3-sigma range for the dark matter abundance [115] is highlighted in red in the upper left panel, the limits on the annihilation cross section today from Fermi-LAT's 6-year analysis of dwarf spheroidal galaxies [116] for annihilation to τ + τ − (red) and bb (blue) are shown in the upper right panel, and exclusion limits from LUX (red) [117], PandaX (blue) [118], and LZ (projected; black) [119] are shown in the lower left panel.
flected in the different higgsino fractions of the lightest neutralino. From SOFTSUSY, the higgsino fraction inχ 0 1 is around ∼ 3%, while from SPheno, the higgsino fraction is only 0.3%. We note another large discrepancy among the results from the different pipelines is that the values of m A obtained from SOFTSUSY and SPheno are 1470 GeV and 2594 GeV, respectively.
These differences in the spectra have a profound effect on the relic density. From the upper left panel of figure 5, we find that the SOFTSUSY-micrOMEGAs pipeline (solid magenta triangles) and SOFTSUSY-DarkSUSY pipeline (hollow magenta triangles) give values for the relic density Ωh 2 = 0.037 and Ωh 2 = 0.648, respectively, while the pipelines that involve SPheno as the spectrum calculator give relic densities that are Ωh 2 10, with the value returned from the DarkSUSY pipeline larger than that from the micrOMEGAs pipeline by a factor of 2.

JHEP05(2018)113
In fact, the CMSSM 1 benchmark model is not even a cosmologically-favored point (at least within a thermal history) if one uses SPheno as the spectrum calculator. The lightest neutralino in that case is an almost pure bino, far away from the A-resonance, and without any possible contributions from coannihilation channels. The pipelines involving SPheno give a relic density that is 2 orders of magnitude larger than those given by SOFTSUSY pipelines.
For pipelines involving SOFTSUSY, the CMSSM 1 benchmark falls approximately into the category of well-tempered dark matter. 9 With the SOFTSUSY spectrum, we see a factor of ∼ 20 difference in Ωh 2 between micrOMEGAs (solid magenta triangles) and DarkSUSY (hollow magenta triangles), with neither giving a value within the limits of experimental uncertainty. As in the pMSSM, for well-tempered dark matter, DarkSUSY tends to give a relic density value that is larger than the one given by micrOMEGAs for the same spectrum. This could arise due to differences in the way the effective annihilation cross section between the dark matter and the neutral and charged higgsinos is computed by the two programs, as well as differences in the calculation of the t-channel chargino exchange diagram. While in the case of the pMSSM the difference was small, in the case of the GUT model benchmark the difference is enormous, possibly due to the proximity to the A-funnel region. With micrOMEGAs we obtain a dark matter candidate that annihilates too efficiently in the early Universe, while with DarkSUSY we obtain a candidate that does not annihilate efficiently enough.
Regarding the relic density for the CMSSM 2 point, as mentioned above, only SPheno pipelines represent true coannihilation models, while SOFTSUSY pipelines do not coannihilate efficiently enough to achieve the correct relic abundance. Thus, the relic abundance from SOFTSUSY pipelines appears at large values of Ωh 2 beyond the range shown in the upper left panel of figure 5.
We now turn to the annihilation cross section in the current Universe for the CMSSM benchmarks. The results are plotted in the upper right panel of figure 5 and tabulated in table 21. We see that the enormous difference in the relic density computation between spectra coming from SOFTSUSY and SPheno continues to persist in the computation of the current annihilation cross section for the CMSSM 1. For this benchmark we see the largest discrepancies in the calucation of the annihilation cross section today, nearly three orders of magnitude. Since the charged higgsinos are much heavier in the spectrum generated by SPheno, the t−channel chargino exchange diagram is suppressed in this case. This leads to a much smaller annihilation cross section, ∼ 10 −28 cm 3 s −1 . For a given spectrum, the difference between micrOMEGAs and DarkSUSY persists, with DarkSUSY giving smaller annihilation cross section as before. For the CMSSM 2 benchmark, neutralino-stau coannihilations play no role in the annihilation today, so we see reasonably good agreement among the pipelines, albeit with a very low annihilation cross section.
Finally, we consider the neutralino-nucleon scattering cross sections, which are presented in the lower left panel of figure 5 (as per-nucleon scattering cross sections) and JHEP05(2018)113 tabulated in table 22 and 23, where the subtables are organized by proton then neutron scattering cross sections, first for the spin-independent then for spin-dependent scattering.
From the lower left panel of figure 5, we see that the scattering cross sections for the CMSSM 1 spectrum coming from SOFTSUSY are much larger than those for the corresponding spectrum coming from SPheno for both DarkSUSY and micrOMEGAs. In fact, the SOFTSUSY model is already being constrained by current experiments. This is due to the larger higgsino content of the dark matter in the SOFTSUSY case, which leads to an enhancement of the Higgs exchange diagram. Moreover, as was observed in the pMSSM cases, there are discrepancies between the scattering cross sections reported by micrOMEGAs versus DarkSUSY, even for the same spectrum, from the differences in form factors employed by each code (see table 4) and the way in which loop effects are incorporated. For the CMSSM 2, micrOMEGAs gives the same SI scattering cross section no matter which spectrum generator is employed, while DarkSUSY yields somewhat larger SI scattering cross sections that do depend somewhat on the details of the spectrum that differ. Since the LSP for the CMSSM 2 benchmark is nearly pure bino and the SI cross sections are strongly suppressed, differences in the SI scattering cross sections from the SPheno-DarkSUSY versus SOFTSUSY-DarkSUSY pipelines likely come from loop corrections.

NUHM benchmarks
As discussed in section 2.1, since the original MasterCode NUHM best fit point has a nearly pure higgsino LSP and is therefore very sensitive to the details of the spectrum, our NUHM A and B benchmarks were chosen with the requirements that a valid relic density would be achieved by NUHM A via the SPheno pipelines and by NUHM B via the SOFTSUSY pipelines.
From table 8, the mass of the lightest neutralino obtained by SOFTSUSY is far smaller than that obtained by SPheno for both NUHM A and B. The dark matter is mostly higgsino in all cases, and the radically different masses for the LSP (and other light-inos) are due to the very different values of µ reported by the two programs, as discussed previously in section 5.1.
The analyses of the dark matter observables for both NUHM benchmarks follow the trends of the pure higgsino case discussed previously for the pMSSM. For both benchmarks, we expect that the relic density for pipelines involving SOFTSUSY should be lower than that given by pipelines involving SPheno, due to the smaller higgsino mass in the former case. This expectation is borne out in the upper left panel of figure 5. For a given spectrum, the relic densities computed by micrOMEGAs and DarkSUSY match quite well, as evidenced by the fact that solid and hollow NUHM markers more or less overlap.
For the annihilation cross section in the current Universe, we expect that the lower higgsino mass of the SOFTSUSY pipelines should correspond to a larger annihilation cross section than the SPheno pipelines. This is borne out by the relative positions of the magenta and cyan NUHM markers in the upper right panel of figure 5, though the difference is less pronounced for micrOMEGAs pipelines for the NUHM B benchmark.
The scattering cross sections with nuclei are shown in the lower left panel of figure 5. We first note that the SI scattering cross sections for higgsino dark matter in this case are within a factor of a few of ∼ 10 −9 pb, which is much larger than the cross section for the -36 -

JHEP05(2018)113
pure higgsino case in the pMSSM analysis. This can be attributed to the fact that the dark matter in the NUHM benchmarks is less pure higgsino than in the "pure higgsino" pMSSM benchmark, as can be seen by comparing tables 8 and 5. The more pure the LSP, the smaller the SI scattering cross section. We see also that DarkSUSY gives larger scattering cross sections than micrOMEGAs, consistent with results discussed above. As before, we can attribute this to the difference in form factors between the two programs.

Conclusions
In this paper, we have performed comparative studies of the physics of supersymmetric dark matter calculated with a sequence of spectrum generators (SOFTSUSY and SPheno), Higgs sector calculators (FeynHiggs and SusyHD) and to dark matter observable calculators (micrOMEGAs and DarkSUSY). We placed our study in the context of several SUSY benchmark models that are interesting in light of LHC Run-1 and null results from recent dark matter searches as studied previously by [26,30,54]. We have compared calculations for the sparticle spectra, the Higgs sectors, and the dark matter observables for each benchmark, and we have incorporated the various generators and calculators into comprehensive pipelines to study not only the effects of the choice of an individual package, but also all downstream effects of those choices on subsequent calculations.
This study was conducted in two parts. In the first part, we investigated a set of pMSSM models from ref. [54]. The dark matter scenarios we considered were coannihilation (bino-stop and bino-squark), A-funnel, well-tempered neutralinos, and pure higgsinos. We discovered that the spectrum generators can differ by up to 1 -2 % in their predicted masses for the stop and the first two generations of squarks, and by up to 20% in the gauge composition of the lightest neutralino, for a given pMSSM model. As for the dark matter observables, differences of up to a factor of ∼ 3 − 5 in the relic density and current annihilation cross section, and up to a factor of ∼ 10 in the predicted scattering cross section, were found to exist between the different pipelines. These discrepancies are already pressing in the case of the relic abundance of dark matter, for which the uncertainty in the experimental value is small compared to the theoretical uncertainty in pMSSM predictions. For the annihilation cross section today and the SI scattering cross section, the discrepancies will become important if/when future dark matter indirect and direct detection experimental sensitivities reach the predicted levels.
In the second part of our study, we considered four interesting benchmark models defined at the GUT scale -two CMSSM points and two NUHM points. For GUT-scale models, we found that discrepancies among the various pipelines are often amplified by the renormalization group running. For our CMSSM and NUHM benchmarks, we found that the spectrum generators can give low energy values of the higgsino mass parameter µ and the pseudoscalar Higgs mass m A that differ by up to 150% -200% (though the differences can be much greater at larger m 0 ). This can lead to large differences in the annihilation and scattering cross sections computed by the dark matter calculators.