Precision predictions for dark matter with DM@NLO in the MSSM

We present DM@NLO , a Fortran 77 based program with a C++ interface dedicated to precision calculations of dark matter (DM) (co)annihilation cross-sections and elastic dark matter-nucleon scattering amplitudes in the Minimal Supersymmetric (SUSY) Standard Model (MSSM) at next-to-leading order (NLO) in perturbative (SUSY) QCD. If the annihilating initial particles carry an electric or colour charge, the Sommer-feld enhanced cross section is included as well and can be matched to the NLO cross section. We review these calculations including technical details relevant for us-ing the code. We illustrate their impact by applying DM@NLO to an example scenario in the constrained MSSM.


Introduction
The relic density of cold dark matter (CDM) in the cosmological standard (ΛCDM) model as determined from Planck data [1] is one of the most stringent constraints on the nature and properties of DM.In the canonical freeze-out mechanism, today's DM abundance may be predicted for a given particle physics framework with a suitable DM candidate through all possible (co)annihilation processes of the thermal relic into Standard Model (SM) particles, and the most important channels are those with two particles in the final state [2,3].
In addition to the increased precision, a further advantage of full loop calculations is that they allow for a systematical evaluation of the theoretical uncertainties from missing higher-order corrections through variations of the renormalisation scheme as well as the renormalisation scale [16,30].
In this context, we present in this paper the public release of the DM precision code DM@NLO ("Dark Matter at Next-to-Leading Order"), which allows to numerically calculate the total DM (co)annihilation as well as DM-nucleon scattering cross-sections at nextto-leading order (NLO) in the strong coupling constant for most (co)annihilation channels within the MSSM.In addition to the fixed-order one-loop calculation, the code includes several resummed corrections, such as the SUSY-QCD ∆m b resummation and the Sommerfeld effect from the exchange of gluons or photons between the incoming particles.
Additional uncertainties not addressed within the context of DM@NLO stem from possible non-standard cosmologies [55], the neglect of multi-body final states [56], uncertainties of the SM effective number of degrees of freedom [57], differences in the numerical computation of the physical mass spectrum of the BSM model under consideration [58], thermal corrections [59], or the early-kinetic decoupling effect [60,61].Moreover, on the side of the experimental analysis, the measured DM relic density relies on the the six-parameter cosmological concordance model.Including additional physical parameters [62] or relaxing the assumption of a constant comoving DM density [63] may increase the allowed range of the DM density by up to about 20%.
In the remainder of this manuscript, we begin by briefly summarising in Sec. 2 all relevant formulae for the calculation of the DM relic density in the standard freeze-out scenario and the spin-independent as well as spin-dependent elastic DM-nucleon scattering cross sections which are necessary to correctly interpret the DM@NLO output.The different components of a oneloop cross section as well as the Sommerfeld enhancement effect are described in Sec. 3. Sec. 4, Appendix A, and Appendix B are dedicated to the installation and running of the program.An illustrative MSSM scenario is defined and analysed in Sec. 5. We summarise this paper in Sec. 6.

Standard calculation of the relic density
The standard procedure [2,64,65] of calculating the freeze-out abundance of a single relic particle is usually based on the assumptions that i) the total (co)annihilation cross-section is governed by 2 → n annihilation processes of DM into SM particles, ii) all dark sector particles are in kinetic equilibrium with the SM thermal bath at the photon temperature T due to sufficiently large elastic scattering rates between both sectors, iii) these share the same chemical potential and iv) are highly non-relativistic so that in-medium as well as finite temperature effects are negligible.This means in particular that even for n > 2, SM particles are assumed to obey a Maxwell-Boltzmann distribution.In this case, the evolution of the DM number density n χ over time t can be described with the single Boltzmann equation with the equilibrium number density of a particle species i with mass m i and g i degrees of freedom for a vanishing chemical potential.The Hubble expansion rate H is usually parameterised through the effective number of SM energy and entropy degrees of freedom.K i denote the modified Bessel functions of order i.Particle physics enters through the thermally averaged annihilation cross section including all possible 2 → n (co)annihilation channels of dark sector particles into a set X of n SM particles.For a given initial-state configuration, the thermal average can be cast into a single integral over the centre-ofmomentum energy with p cm being the relative annihilation momentum.
From the present-day number density n 0 χ , the theoretically predicted dark matter relic density is obtained via with ρ c being the critical density.The task of numerically integrating Eq. ( 2) in conjunction with computing the thermal average in Eq. ( 4) for a model dependent annihilation cross section can be performed with a high accuracy using adequate numerical codes, e.g., MicrOMEGAs [66], SuperIso Relic [67], MadDM [68], DarkSUSY [69] or DRAKE [70], which can even go beyond the assumption of kinetic equilibrium.

Calculation of the direct detection rate
Results of direct dark matter detection experiments are usually presented as exclusion limits on the spindependent (SD) and spin-independent (SI) DM-nucleon scattering cross-sections, σ SD N and σ SI N , as a function of the DM mass.However, as the typical energies in a direct detection experiment are much smaller than the heavy particle masses of the microscopic theory mediating the interaction between DM and the constituents of a nucleon, it is customary to perform the calculation in the language of an effective field theory [71][72][73][74], i.e. by integrating out those heavy mediators.
The spin-independent cross section is then expressed through the SI effective DM coupling to nucleons g SI N with µ N = m N m χ /(m N +m χ ) being the reduced mass of the DM-nucleon system.The effective coupling is computed as where the sum runs over all six quark flavors q and α SI q is the Wilson coefficient describing the SI interaction between quarks and the DM particle.The nuclear matrix element ⟨N | qq |N ⟩ can be qualitatively understood as the probability of finding the quark q inside the nucleon N and is commonly expressed through the scalar nuclear form factors with the quark mass m q and the nucleon mass m N .The scalar coefficients f N T q are determined from experiment and lattice QCD and are another source of theoretical uncertainties.To highlight the latter point, we show in Tab.2.2 the associated values that are hardcoded in DM@NLO and the two other DM packages Dark-SUSY and MicrOMEGAs.The heavy quark form factors f N T q are obtained from those related to light quarks via the relation [75] The SD scattering cross section for DM on a single nucleon is given by where the effective SD coupling g SD N between DM and nucleons reads 0.0663 0.0595 0.0682 with the SD Wilson coefficient α SD q describing the DMquark interaction.In contrast to the SI case, the sum runs only over the light quarks u, d and s, as these carry the largest fraction of the nucleon spin which in turn is quantified through the axial-vector form factors (∆q) N .The corresponding numerical values in DM@NLO are identified with those in MicrOMEGAs 5.3, given by For a detailed discussion of direct detection in NLO SUSY-QCD, we refer the reader to Ref. [39].

Dark matter annihilation beyond tree-level
At next-to-leading order (NLO), the tree-level DM (co)annihilation cross-section is extended by the contribution which contains virtual (dσ V ) and real (dσ R ) corrections, contributing at the same order in the coupling constant.In the present work, we focus on one-loop corrections in SUSY-QCD at order α s , including the emission of a real gluon.Note that, for certain couplings, we include the resummation of higher-order SUSY-QCD contributions as discussed in Refs.[5,9,10].

Renormalisation
The virtual corrections are plagued by ultraviolet (UV) divergences whose removal requires a (numerically well behaved) renormalisation scheme, coming along with a suitable regularisation prescription.For the regularisation, the SUSY-preserving dimensional reduction (DR) scheme [78,79] is employed, i.e. the relevant loop integrals are evaluated in D = 4−2ϵ space-time dimensions.When it comes to the renormalisation of the MSSM, the squark masses have to be renormalised carefully since the stop and sbottom sectors have to be treated simultaneously due to the fact that the up-and downtype squarks share the common soft breaking parameter M q as a result of the SU(2) L gauge symmetry.The squark mass matrix can be diagonalised, with the two physical masses m 2 q1 and m 2 q2 being the eigenvalues of the non-diagonal mass matrix with the entries Out of the eleven parameters , θ t and θ b only five are completely independent.As the renormalisation scheme should be applicable to all (co)annihilation channels with squarks in a leading role, we replace the soft SUSY-breaking masses M Q, M Ũ , M D as input parameters by the physical on-shell masses m b1 , m b2 and m t1 .The three aforementioned soft parameters are then fixed through the requirement that Eq. ( 15) holds even at the one-loop order which, by inverting the corresponding eigenvalue equations, results in two possible solutions, for the diagonal entries of the mass matrix.Consequently, there are two possible values for M Q, M Ũ , and M D .However, not both of them may yield a numerically stable renormalisation scheme, the reason being that the diagonalisation may not correctly reproduce the mass of the lighter stop used as an input value, and, more importantly, the counterterm belonging to the heavier stop mass δm t2 ∼ (U q 21 U q 12 ) −1 may become singular for vanishing off-diagonal elements of the squark mixing matrix.The same problem might occur in the counterterm related to the squark mixing angle δθ q ∼ (U q 11 U q 22 + U q 12 U q 21 ) −1 .To avoid these issues, a scheme is defined as numerically stable if the following three conditions are fulfilled: Otherwise a scheme is declared as invalid (unstable).Given that both solutions are compatible, by default the solution is chosen where the dependent stop mass m t2 is closer to the corresponding physical value.
In a series of previous analyses [11-13, 15, 17], the following three renormalisation schemes, adapted to the present situation of DM (co)annihilation, have been introduced: The hybrid on-shell/DR scheme 1, which resembles the RS2 scheme presented in Ref. [80], is the recommended option, whereas the other two schemes are well suited for the estimation of theoretical uncertainties from scheme variations.The integration of an automated selection of the best renormalisation scheme as, e.g., discussed in Ref. [81] is left for a future update.

Infrared treatment
The real corrections on the other hand suffer from infrared (IR) divergent terms occurring in the soft or collinear phase-space regions which cancel against those singularities appearing in the one-loop diagrams of the virtual corrections.To make the analytic cancellation manifest and allow for the numerical integration over the real phase-space, we choose within DM@NLO the Catani-Seymour dipole subtraction method [82,83] for massive initial-state particles [84] over phase-space slicing methods [85][86][87].This subtraction technique is based on the introduction of a local counterterm dσ A that cancels the singularities in the real emission matrix element pointwise and is at the same time simple enough such that the integrals over the singular region can be performed analytically and the IR divergences appear as poles of the form ε −1 and ε −2 .The whole procedure can be schematically captured in the equation and has the advantage over the slicing approach that no introduction of an arbitrary cutoff value on the energy or emission angle of the emitted gluon is required.This makes subtraction methods in general more numerically stable and therefore better suited for parameter space scans 1 .

Intermediate on-shell resonance subtraction
Within the real emission contribution to the process χ0 n ti → bW + , the internal top propagator can become on-shell if the collisional energy √ s exceeds the top mass m t .To cure the singularity, we follow the "Prospino scheme" defined in Refs.[88,89] and substitute the top propagator with the Breit-Wigner form, according to in the resonant part M r of the total amplitude M tot = M r + M nr , whereas the non-resonant piece M nr remains unchanged.Since the corresponding process χ0 n ti → tg is already accounted for in the calculation of the neutralino relic density, the contribution from the leading order on-shell production of a top with the subsequent decay into a bottom quark and a W -boson is removed locally through the replacement with the physical top width Γ t .This procedure has the advantage that it retains the interference Re(M * r M nr ) containing only one on-shell propagator and thus finite principal-value integrals.However, to stabilise the numerical integration, we use a small artificial top width Γ t = 10 −3 • m t in the interference part instead of the physical width.

Sommerfeld enhancement
The Sommerfeld enhancement is an elementary quantum mechanical effect that increases (decreases) annihilation cross sections for small relative velocities in the 1 For historical reasons, squark annihilation into electroweak final states is computed using the phase space slicing method.For stop annihilation into gluons and light quarks, the quarks of the first two generations are considered as effectively massless.Consequently, a consistent cancellation of infrared divergent terms requires the combination of the two processes ti t * j → gg and ti t * j → q q at the loop level.To avoid doublecounting, we therefore adopt the convention that the process with two gluons in the final state automatically includes the light quark contribution at both LO and NLO.For this reason, the code returns a zero cross section if the final state is set to a light quark-antiquark pair.
presence of an attractive (repulsive) long-range potential affecting the incoming particles.From a field theory point of view, this effect is described by ladder diagrams involving the exchange of light mediators with some coupling λ to the initial particles.More quantitatively, the Sommerfeld factor is obtained as a solution ϕ(⃗ r) evaluated at the origin ⃗ r = 0 of the stationary Schrödinger equation for the potential describing the interaction of the annihilating particles transforming under the representation R of the corresponding force carriers.For an s-wave dominated annihilation process, the Sommerfeld factor simply multiplies the perturbative tree-level cross section giving the Sommerfeld corrected cross section where the sum runs over all irreducible representations contained in the decomposition of the initial particle pair.In this convention for the Sommerfeld factor, the free wave-function is normalised to one |ϕ [R] 0 (0)| 2 = 1 to ensure that σ Som → σ Tree is fulfilled if the interaction governing the enhancement effect is turned off (λ → 0).When combining the Sommerfeld effect with the full O λ 2 correction, one has to be careful to not overcount the single mediator exchange contained in both calculations.Therefore, we make the decision to match both by removing the O λ 2 contribution from the Sommerfeld factor.
4 Installing and running DM@NLO

Installation
The DM@NLO package is a high-energy physics program whose source code is written in Fortran 77 with a C++ interface similar to the precision code Resummino [90,91].It is publicly available for download at https://dmnlo.hepforge.organd is licensed under the European Union Public Licence v1.1.
The code can be compiled with the GNU compiler collection (GCC) and CMake version 3.0 or higher.As external dependencies, the libraries SLHALib-2.2[92] and LoopTools-2.16[93] are required for reading particle spectra following the Supersymmetry Les Houches Accord 2 (SLHA 2) convention [94,95] and for evaluating one-loop integrals, respectively.The code ships directly with slightly modified versions of both libraries, as well as the CUBA-1.1 [96] library for performing multidimensional phase space integrals through the VEGAS Monte Carlo algorithm.
Once downloaded, the code can easily be unpacked and installed by running the following commands in a shell: tar xvf DMNLO-X.Y.Z.tar cd DMNLO-X.Y.Z mkdir build cd build cmake .. [options] make make install The last command is optional and places the DM@NLO binary dmnlo as well as the static library libdmnlo.a in the top-level source directory, which is the setup we assume in the following.Otherwise, the executable can be found in build/bin and the library in build/src.To install the code, e.g.system wide, the installation directory can be set with the cmake option -DCMAKE_INSTALL_PREFIX= Compilers different from the default C, C++ and Fortran compilers identified by CMake can be set with

-DCMAKE_<LANG>_COMPILER=
The path to an alternative LoopTools installation can be specified with -DLOOPTOOLS=PATH after setting -DBUILD_LOOPTOOLS (default: ON) to OFF if libraries and headers are installed in the same folder, or through LOOPTOOLS_INCLUDE_DIR and LOOPTOOLS_LIB_DIR if not.
After successful compilation, the local installation can be tested by running the commands ./dmnlo--help ./dmnloinput/DMNLO.in in a shell.The source files of the C++ interface to DM@NLO are located in src, whereas the processes themselves implemented in Fortran 77 are collected in the folder run_dmnlo.The name of each subfolder for every process supported by DM@NLO is summarised in Tab. 2, together with the key references documenting the corresponding calculational details.The directory external contains external dependencies like Loop-Tools or SLHALib.The folder input/demo provides for every process available in DM@NLO, sorted according the arXiv number of the corresponding publication, the associated example scenarios as SLHA 2 files as well as Python 3 plotting routines that partially use PySLHA [97] and allow to reproduce the most important cross section plots.

Running DM@NLO from the command line
As indicated above, DM@NLO can be executed in a shell through the command ./dmnlo<dmnlo-input-file> where the mandatory argument <dmnlo-input-file> provides the path to a configuration file in plain text format specifying the process and corresponding input parameters.Details on all the available options in such an input configuration file are extensively documented in Appendix A. One example input file delivered with the code is input/DMNLO.in.Alternatively, the parameter values defined in the input file can also be passed through the command line interface (CLI), which then supersedes the value included in the text file.In the following, all possible command line options are described.A concise summary is also provided in Appendix B. The command line options follow the same naming convention as the variables in the configuration file, so that the transfer from the command line to the input file is straightforward.We start with the general options that are valid for both the relic density as well as the direct detection module.
The path to the SLHA file containing the numerical values of masses, mixing angles and decay widths has to be defined with --slha.The value of the renormalisation scale in GeV is fixed through the --muR option whereas the renormalisation scheme must be set to one of the three schemes defined in Sec.3.1 with the option --renscheme.The mixed DR-OS scheme no. 1 is here the recommended option.
The --choosesol option defines the solution in the heavy quark sector, with 0 being the recommended option, where the solution is chosen such that the dependent stop mass m t2 is closest to the corresponding onshell value from the SLHA file (see discussion above).The arguments 1 and 2 then correspond to the two solutions in Eqs.(19) and (20), respectively.If DM@NLO is used from the command line and the renormalisation scheme fails, the code simply stops after issuing a warning.
Also included is a legacy option which can only be turned on through the CLI by passing the flag --legacy.This mode defines the weak mixing angle θ W and the W -mass as in the default MSSM model file in Mi-crOMEGAs 2.4.1, i.e. sin θ W = 0.481 and m W = cos θ W m Z with m Z being the on-shell Z-mass.We include this option since this definition was adopted in DM@NLO before the public release and allows to reproduce old results.Starting with v1.0.0 the electroweak mixing angle is defined through the on-shell Zand Wmass from the SLHA 2 file Note that the legacy option should only be used for the reproduction of previously published results.

Process Folder
References Sommerfeld Lastly, the perturbative order of the calculation needs to be specified.This is only possible through the CLI via the arguments --lo for LO-accurate predictions and --nlo for NLO accuracy.For the calculation of (co)annihilation cross sections there are two more accuracy options.The flag --sommerfeld returns the Sommerfeld enhancement alone, whereas --full returns the NLO result matched to the Sommerfeld enhancement.Otherwise the highest order available is assumed.If no Sommerfeld enhancement is available, the --full option returns just the NLO cross section.
The initial and final particles of the (co)annihilation process are fixed according to the PDG numbering scheme [98].The two options --particleA and --particleB fix the initial state, whereas the two produced SM particles must be referred to by setting --particle1 and --particle2.The collisional energy √ s has to be defined with --pcm, which is the center-of-mass momentum p cm of the incoming particles.The option --result controls whether the output contains the total cross section σ or the cross section times velocity σv, both in units of GeV −2 , where the relative velocity is defined as b )/s with λ being the Källén function.
The direct detection module is enabled through the --DD option, which supersedes the specified (co)annihilation settings.The output contains then the SI and SD scattering cross sections of DM on protons and neutrons in cm 2 , respectively.The scalar nuclear form factors from Tab. 2.2 can be found (and modified) in DD/DD_Init.F.The --formfactor option followed by an integer number (0 for DM@NLO, 1 for DarkSUSY and 2 for MicrOMEGAs) allows the user to select a set of values from Tab. 2.2.

The DM@NLO library
To facilitate the usage of DM@NLO from within other codes, the static library libdmnlo.aprovides the two functions double cs_dmnlo(order, na, nb, n1, n2, PcmIn, muR, &slha, rs, sol, &corrFlags) void dd_dmnlo(order, muR, &slha, rs, sol, ff, &cs) where the former returns the total (co)annihilation cross section and the latter writes the SI and SD DM-nucleon cross sections into the array cs.The renormalisation scale is set with muR, the SLHA 2 input file with slha, the renormalisation scheme through the flag rs and the associated solution for the three soft-breaking parameters with rs.The parameter sol corresponds to the choosesol option and ff in the argument set of the direct detection function to the formfactor option.
The integer order specifies the perturbative order of the calculation.Possible vales are 0 for the LO result, 1 for the NLO result.For the computation of the (co)annihilation cross section, two additional options are available for the order parameter, namely 2 for the full result (including NLO calculation and Sommerfeld enhancement) and 3 for the Sommerfeld enhanced cross section alone (without including the NLO calculation).
The parameters na and nb are needed to fix the incoming particles through their respective PDG numbers, while n1 and n2 are meant to specify the two particles in the final state.The center-of-mass momentum is set through PcmIn.Finally, the integer array corrFlags allows to turn certain processes on and off, which may be useful if the corresponding contribution to the relic density is known to be negligible.
The static library libdmnlo.aalso provides the two functions int canImprove_dmnlo(na, nb, n1, n2) int consistent_RS_dmnlo(rs, &slha, muR) The former allows to check whether a given process can be corrected with DM@NLO, while the latter verifies whether the particle spectrum contained in slha yields a stable renormalisation scheme.
Alternative to the manual decision what annihilation channels to include, the file minimal_example.cpplocated in external/micromegas_5.3.41/MSSMexemplifies the use of these functions with MicrOMEGAs in a way that only those channels contributing more than 2 % to the relic density are corrected.Before compiling the minimal example file through make main=minimal_example.cpp the file micromegas_5.3.41/include/modelMakefile has to be replaced with the associated modified version shipped with DM@NLO.This can be achieved by running tar xvfk micromegas_5.3.41.tgz in the external/ directory where the tar option -k (or --keep-old-files) ensures that our modified version of modelMakefile containing the paths to the required libraries according to the default installation of DM@NLO is retained.For different paths or an alternative MicrOMEGAs version, modelMakefile has to be adjusted accordingly by the user.After successful compilation, the MicrOMEGAs interface can be tested by running ./minimal_exampleScenario.spc in a shell from the MSSM folder.For more details on the usage of the corrFlags argument, we refer to the explanation given in minimal_example.cpp.

Illustrative example calculations
To illustrate the usage of DM@NLO, we present example calculations in the constrained minimal supersymmetric extension of the Standard Model (cMSSM) which contains the simplifying assumption that the soft supersymmetry-breaking parameters unify at the gauge coupling unification scale of about 10 16 GeV.This setup is entirely characterised through the universal scalar mass parameter m 0 , the universal gaugino mass parameter m 1/2 , the ratio of the vacuum expectation values Channel Contribution of the neutral components of the two Higgs doublets tan β, the universal trilinear coupling A 0 , and the sign of the Higgs mixing parameter µ.
In the following, we use, inspired by the recent search for non-excluded regions in the cMSSM parameter space [99], the parameter point given in Tab. 3, where stop (co)annihilation is the dominant dark-matter annihilation mechanism 2 .The (co)annihilation channels contributing most to ⟨σ ann v⟩ are listed in Tab. 4 with stopantistop annihilation into gluons having the largest relative contribution.For each channel that can be corrected with DM@NLO, we show in Fig. 1 our treelevel (black dashed line), the one-loop (blue solid line), and the full cross-section (red solid line), containing in addition to the NLO result the Sommerfeld enhancement, if available.For reference, the cross-section produced with the default MicrOMEGAs setup (orange solid line) is also shown.The grey shaded area depicts (in arbitrary units) the thermal velocity distribution, in order to demonstrate in which p cm region the cross-section contributes to the total annihilation cross-section ⟨σ ann v⟩.In the lower part, the corresponding relative shifts of the different cross-section values (second item in the legend) are shown.Note that the difference between the MicrOMEGAs prediction and our tree-level result, which is particularly large for the gg (q q), tg and tt final states, is mainly due to a different choice of the renormalisation scale.We define the renormalisation scale through the tree-level stop masses as µ R = √ m t1 m t2 , which for the particular scenario in Tab. 3 yields µ R = 1368.2GeV, whereas Mi-crOMEGAs 5.3.41uses the scale Q = 2m χ0 1 /3 for the evaluation of the strong coupling stored in the global 2 Note that numerical differences in the physical mass spectrum occur with respect to Ref. [99] since SPheno 4.0.5 is used as spectrum generator in this work whereas Ref. [99] makes use of a private code.This is also the reason why A 0 = 4m 0 is chosen in the example scenario in Tab. 3 versus A 0 = 3m 0 in Ref. [99].variable GGscale3 .We also show the uncertainties from variations of the renormalisation scale µ R by a factor of two around the central scale as shaded bands.
In Fig. 2 we show the impact of our corrections on the neutralino relic density Ω χ h 2 by performing a scan in the m 1/2 -m 0 plane around the reference scenario of Tab. 3, which is indicated by a red star.The orange band (Ω MO χ ) indicates the region consistent with the observed value Ω CDM h 2 from Eq. ( 1) purely based on MicrOMEGAs and under the assumption that the lightest neutralino solely accounts for all of the observed DM, the blue band (Ω Tree χ ) corresponds to the prediction where the DM@NLO tree-level cross sections replace the CalcHEP result, and the yellow band (Ω Full χ ) is based on our full calculation.
The width of the three bands reflects the experimental 2σ uncertainty shown in Eq. ( 1).One can observe a clear separation between all three bands everywhere across the shown m 1/2 -m 0 plane.The black contour lines quantify the relative difference between our treelevel and our full calculation of the neutralino relic density.The increase amounts to roughly 16% to 18% in the regions consistent with the observed relic density Ω CDM .Let us mention that DM@NLO allows to correct a large portion of the different contributions to the relic density which is between 85% and 90% in the relevant regions.
Finally, we illustrate the application of DM@NLO to direct detection, i.e. we discuss the neutralino-nucleon cross-sections for different values of m 1/2 around the example scenario.In the upper panels of Fig. 3, the SI proton (left) and neutron (right) cross sections are shown, whereas the corresponding SD quantities are presented in the two lower panels.All quantities have been calculated with our DM@NLO code at tree level (black solid line), including the full O(α s ) corrections to the dominant effective operators (blue solid line), Mi-crOMEGAs (orange solid line) and the corresponding analytic tree-level calculation4 .
We also show all three values of the resulting relic density (Ω MO χ , Ω Tree χ , Ω Full χ ) with the same colour coding as in Fig. 2, as well as the Planck compatible value through a grey band.Note that the three curves increase as expected with the neutralino mass.

Summary
In this paper, we have presented the DM@NLO package, dedicated to precision calculation of dark matter (co)annihilation processes and direct detection in the MSSM.The program allows to compute total crosssections at leading order and next-to-leading order in perturbative SUSY-QCD including the Sommerfeld enhancement effect.
To illustrate the usage of the code and the impact of the higher-order corrections, various computations in a typical supersymmetric dark matter scenario were performed using DM@NLO.Apart from the benefit of having more precise predictions, we emphasise the possibility of estimating theory errors as a major advantage of using NLO cross sections and beyond.In addition, the general structure of the code is well-suited for the extension to non-supersymmetric models.choosesol = <int>: solutions for M Q, M Ũ , M D as explained in Sec.3.1.
-particleA = <int> and particleB = <int>: PDG numbers of the first and second particle in the initial state.
-particle1 = <int> and particle2 = <int>: PDG numbers of the first and second particle in the final state.
result = <string>: defines whether the output should contain the total cross σ corresponding to the value s or the total cross section times relative velocity σv defined through the value sv.formfactor = <int>: sets the scalar nuclear form factors f N T q to one of the sets of values in Tab.2.2.
Appendix B: Options available from the command line interface of DM@NLO The DM@NLO program can be run with several command line options by typing in a shell ./dmnlo <input-file> [options] The keyword <input-file> provides the path to a configuration file specifying the details of the computation to be achieved according to the standard defined in Appendix A. The following options are allowed: ---help, prints a help message to the screen, indicating how to execute the code.
---slha, followed by a string sets the path to the SLHA 2 parameter file containing the numerical values of masses, mixing angles, decay widths, etc.
---muR, followed by a double-precision number sets the value of the renormalisation scale in GeV.
---renscheme, followed by an integer number sets the renormalisation scheme according to the numbering scheme introduced in Sec.3.1.
---choosesol, followed by an integer number defines which solution to use for M Q, M Ũ , M D as explained in Sec.3.1.
---legacy, defines the weak mixing angle θ W and the W -mass as in the default MSSM model file in MicrOMEGAs 2.4.1.---lo, returns the result at LO accuracy.
---full, returns the NLO result matched to the Sommerfeld enhancement if the latter is available.Otherwise the output is identical to --nlo.
---particleA and --particleB, followed by integer numbers defines the nature of the two initial-state particles through their PDG numbers.
---particle1 and --particle2, followed by integer numbers defines the nature of the two final-state particles through their PDG numbers.
---pcm, followed by a double-precision number sets the centre-of-mass momentum p cm in GeV.
---result, followed by a string corresponding to s for the total cross σ or sv for the total cross section times the relative velocity σv.
---DD, enables the direct detection module.This option supersedes (co)annihilation settings.---formfactor, followed by an integer number ranging from zero to two sets the scalar nuclear form factors f N T q to one of the value sets shown in Tab.2.2 with zero for DM@NLO and two for MicrOMEGAs.
all DR parameters.1: m b , A b and A t are DR input parameters whereas m t , m t1 m b1 and m b2 are on-shell (OS) masses.θ t, θ b and m t2 are then dependent quantities.2: m t , m b , A b and A t are DR input parameters and m t1 , m b1 and m b2 are OS masses.θ t, θ b and m t2 are then dependent quantities.

− 2 ]σ− 2 ]σ− 2 ]σFig. 1
Fig.1Tree-level (black dashed line), one-loop (blue solid line), full (red solid line) if present and MicrOMEGAs (orange solid line) cross sections for the dominant (co)annihilation channels shown in Tab. 4 that can be corrected with DM@NLO including the corresponding uncertainties from variations of the renormalisation scale µ R by a factor of two around the central scale as shaded bands.The upper part of each plot shows the absolute value of σv together with the thermal velocity distribution (in arbitrary units), whereas the lower part shows the corresponding relative shift (second item in the legend).

Fig. 2
Fig. 2 Bands compatible with the Planck measurement in Eq. (1) in the m 1/2 -m 0 plane (left) and the plane spanned by the associated physical masses of the lightest neutralino and the lightest stop (right) surrounding the example scenario from Tab. 3 shown in form of a red star.The three bands correspond to the MicrOMEGAs calculation (orange), our tree-level (blue) and our full corrections (yellow).The black solid lines indicate the relative change (Ω Tree χ − Ω Full χ )/Ω Tree

10 −46 cm 2 ] 9 Ω χ h 2 Fig. 3
Fig.3Spin-independent (top) and spin-dependent (bottom) neutralino-nucleon cross sections for protons (left) and neutrons (right) in the example scenario in Tab. 3 for different values of the universal gaugino mass parameter as well as the corresponding neutralino relic density obtained with MicrOMEGAs (MO), our tree-level calculation (Tree) and with our full calculation (Full).The upper and lower limits imposed by Eq. (1) are indicated through the grey band.

Table 3
Example scenario in the cMSSM with a positive Higgs supersymmetric mixing parameter µ where stop (co)annihliation is the dominant dark matter mechanism.All dimensionful quantities are in GeV.

Table 4
Dominant annihilation channels contributing to ⟨σ ann v⟩ for the cMSSM scenario in Tab. 3. Further contributions below 2% are omitted.