colorlinks=true,linkcolor=RedViolet,citecolor=RedViolet,filecolor=RedViolet,urlcolor=RedViolet

<EMAIL_ADDRESS><EMAIL_ADDRESS>e3michael.klasen@uni-
muenster.de<EMAIL_ADDRESS><EMAIL_ADDRESS>
11institutetext: PRISMA+ Cluster of Excellence & Mainz Institute for
Theoretical Physics, Johannes Gutenberg University, 55099 Mainz, Germany
22institutetext: LAPTh, Université Savoie Mont Blanc, CNRS, 74000 Annecy,
France 33institutetext: Institut für Theoretische Physik, Universität Münster,
Wilhelm-Klemm-Straße 9, 48149 Münster, Germany

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

Julia Harz e2,addr2  Björn Herrmanne1,addr1  Michael Klasen e3,addr3  Karol
Kovařík e4,addr3  Luca Paolo Wiggeringe5,addr3

(Received: date / Accepted: date)

###### Abstract

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 Sommerfeld 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 using the code. We
illustrate their impact by applying DM@NLO to an example scenario in the
constrained MSSM.

††journal: Eur. Phys. J. C

MITP-23-084, MS-TP-23-52, LAPTH-062/23

##  1 Introduction

The relic density of cold dark matter (CDM) in the cosmological standard
($\Lambda$CDM) model

$\displaystyle\Omega_{\rm CDM}h^{2}~{}=~{}0.120\pm 0.001$ (1)

as determined from Planck data Planck:2018vyg  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
Gondolo:1990dk ; Jungman:1995df .

Given the sub-percent accuracy of the Planck observation in Eq. (1), matching
the experimental and theoretical uncertainties clearly calls for the
evaluation of the associated cross sections beyond the tree-level
approximation. The impact of higher-order corrections to DM annihilation, both
strong and electroweak, has been discussed within many well-motivated and
intensely studied extensions of the SM such as the Minimal Supersymmetric
Standard Model (MSSM) Freitas:2007sa ; Herrmann:2007ku ; Baro:2007em ;
Baro:2008bg ; Baro:2009gn ; Herrmann:2009wk ; Herrmann:2009mp ; Harz:2012fz ;
Harz:2014tma ; Herrmann:2014kma ; Harz:2014gaa ; Schmiemann:2019czm ;
Branahl:2019yot ; Klasen:2022ptb , the Next-to-Minimal Supersymmetric Standard
Model (NMSSM) Belanger:2016tqb ; Belanger:2017rgu , or the Inert Doublet Model
Klasen:2013btp ; Banerjee:2019luv ; Banerjee:2021oxc . The relic density can
in addition potentially receive large corrections from non-perturbative
effects like the Sommerfeld enhancement coming from the exchange of light
mediators between the annihilating particles Iengo:2009ni ; Beneke:2016ync ;
Harz:2017dlj  or the formation of bound-states vonHarling:2014kha ;
Harz:2018csl ; Biondini:2018xor ; Harz:2019rro .

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
Harz:2016dql ; Branahl:2019yot .

The need for an increased theoretical precision extends to the calculation of
indirect detection signals from present-day DM annihilation processes
Hryczuk:2011vi ; Hryczuk:2014hpa ; Baumgart:2017nsr ; Baumgart:2018yed  as
well as to DM-nucleon interactions in the context of direct DM detection.
Calculations of higher-order corrections to the corresponding scattering cross
sections have been computed in many UV-complete models like the MSSM
Drees:1992rr ; Drees:1993bu ; Hisano:2004pv ; Berlin:2015njh ; Klasen:2016qyz
; Bisal:2023fgb , the Inert Doublet Model Klasen:2013btp ; Abe:2015rja ,
simple Higgs-portal models Azevedo:2018exj ; Ishiwata:2018sdi ;
Ghorbani:2018pjh ; Abe:2018emu , the Next-to-Minimal Two Higgs Doublet Model
Glaus:2022rdc , a vector DM model Glaus:2019itb , the Singlet-Extended Two
Higgs Doublet Model Biekotter:2022bxp , but also in simplified fermionic DM
models Berlin:2015ymu ; Ertas:2019dew ; Borschensky:2020olr , or frameworks
described through effective operators Haisch:2013uaa ; Crivellin:2014qxa .

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 next-to-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 $\Delta 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 Arbey:2008kv , the neglect of multi-body
final states Yaguna:2010hn , uncertainties of the SM effective number of
degrees of freedom Saikawa:2020swg , differences in the numerical computation
of the physical mass spectrum of the BSM model under consideration
Allanach:2003jw , thermal corrections Beneke:2014gla , or the early-kinetic
decoupling effect Binder:2017rgn ; Aboubrahim:2023yag . 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 Hamann:2006pf  or relaxing the assumption of a constant comoving DM
density Bringmann:2018jpr  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 one-loop
cross section as well as the Sommerfeld enhancement effect are described in
Sec. 3. Sec. 4, A, and 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.

##  2 Dark matter theory

###  2.1 Standard calculation of the relic density

The standard procedure PhysRevLett.39.165 ; Gondolo:1990dk ; Edsjo:1997bg  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\to 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_{\chi}$ over time $t$ can be
described with the single Boltzmann equation

$\derivative{n_{\chi}}{t}+3Hn_{\chi}~{}=~{}\langle\sigma_{\mathrm{ann}}v\rangle\Big{(}(n_{\chi}^{\mathrm{eq}})^{2}-n_{\chi}^{2}\Big{)}$
(2)

with the equilibrium number density

$\displaystyle
n_{i}^{\mathrm{eq}}~{}=~{}\frac{g_{i}Tm_{i}^{2}}{(2\pi)^{2}}\,K_{2}(m_{i}/T)$
(3)

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

$\langle\sigma_{\mathrm{ann}}v\rangle~{}=~{}\sum_{a,b}\langle\sigma_{ab\to
X}v\rangle\frac{n_{a}^{\mathrm{eq}}}{n_{\chi}^{\mathrm{eq}}}\frac{n_{b}^{\mathrm{eq}}}{n_{\chi}^{\mathrm{eq}}}$
(4)

including all possible $2\to 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

$\displaystyle\begin{split}\langle\sigma_{ab\to
X}v\rangle&=\frac{1}{2\,T\,m_{a}^{2}m_{b}^{2}\,K_{2}(m_{a}/T)\,K_{2}(m_{b}/T)}\\\
&~{}~{}~{}~{}\times\\!\int\\!\differential{s}\sqrt{s}\,p^{2}_{\mathrm{cm}}\sigma_{ab\to
X}(p_{\mathrm{cm}})K_{1}(\sqrt{s}/T)\end{split}$ (5)

can be cast into a single integral over the centre-of-momentum energy
$\sqrt{s}=\sqrt{m_{a}^{2}+p_{\mathrm{cm}}^{2}}+\sqrt{m_{b}^{2}+p_{\mathrm{cm}}^{2}}$
with $p_{\mathrm{cm}}$ being the relative annihilation momentum. From the
present-day number density $n_{\chi}^{0}$, the theoretically predicted dark
matter relic density is obtained via

$\displaystyle\Omega_{\chi}~{}=~{}\frac{m_{\chi}n^{0}_{\chi}}{\rho_{c}}$ (6)

with $\rho_{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., MicrOMEGAsBelanger:2018ccd
, SuperIso Relic Arbey:2009gu , MadDM Ambrogi:2018jqj ,
DarkSUSYBringmann:2018lay  or DRAKE Binder:2021bmg , which can even go beyond
the assumption of kinetic equilibrium.

###  2.2 Calculation of the direct detection rate

Results of direct dark matter detection experiments are usually presented as
exclusion limits on the spin-dependent (SD) and spin-independent (SI) DM-
nucleon scattering cross-sections, $\sigma_{N}^{\mathrm{SD}}$ and
$\sigma_{N}^{\mathrm{SI}}$, 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 Fan:2010gt ;
Hill:2014yka ; Hill:2014yxa ; Hisano:2015bma , i.e. by integrating out those
heavy mediators.

The spin-independent cross section

$\sigma_{N}^{\mathrm{SI}}=\frac{\mu_{N}^{2}}{\pi}|g_{N}^{\mathrm{SI}}|^{2}$
(7)

is then expressed through the SI effective DM coupling to nucleons
$g_{N}^{\mathrm{SI}}$ with $\mu_{N}=m_{N}m_{\chi}/(m_{N}+m_{\chi})$ being the
reduced mass of the DM-nucleon system. The effective coupling is computed as

$g_{N}^{\mathrm{SI}}~{}=~{}\sum_{q}\bra{N}\bar{q}q\ket{N}\alpha^{\mathrm{SI}}_{q}\,,$
(8)

where the sum runs over all six quark flavors $q$ and
$\alpha^{\mathrm{SI}}_{q}$ is the Wilson coefficient describing the SI
interaction between quarks and the DM particle. The nuclear matrix element
$\bra{N}\bar{q}q\ket{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 $f^{N}_{Tq}$ as

$\bra{N}m_{q}\bar{q}q\ket{N}~{}=~{}f^{N}_{Tq}\,m_{N}$ (9)

with the quark mass $m_{q}$ and the nucleon mass $m_{N}$. The scalar
coefficients $f^{N}_{Tq}$ are determined from experiment and lattice QCD and
are another source of theoretical uncertainties. To highlight the latter
point, we show in Tab. 1 the associated values that are hardcoded in DM@NLO
and the two other DM packages DarkSUSY and MicrOMEGAs. The heavy quark form
factors $f^{N}_{Tq}$ are obtained from those related to light quarks via the
relation Shifman:1978zn

$f_{Tc}^{N}~{}=~{}f_{Tb}^{N}~{}=~{}f_{Tt}^{N}~{}=~{}\frac{1}{27}\bigg{(}1-\sum_{q=u,d,s}f_{Tq}^{N}\bigg{)}\,.$
(10)

The SD scattering cross section for DM on a single nucleon is given by

$\sigma_{N}^{\mathrm{SD}}~{}=~{}\frac{3\mu_{N}^{2}}{\pi}\big{|}g_{N}^{\mathrm{SD}}\big{|}^{2}\,,$
(11)

where the effective SD coupling $g_{N}^{\mathrm{SD}}$ between DM and nucleons
reads

$g_{N}^{\mathrm{SD}}~{}=~{}\sum_{q=u,d,s}(\Delta
q)_{N}\,\alpha^{\mathrm{SD}}_{q}$ (12)

with the SD Wilson coefficient $\alpha^{\mathrm{SD}}_{q}$ describing the DM-
quark 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
$(\Delta q)_{N}$. The corresponding numerical values in DM@NLO are identified
with those in MicrOMEGAs 5.3, given by

$\displaystyle\begin{split}(\Delta u)_{p}~{}&=~{}(\Delta
d)_{n}~{}=~{}0.842\,,\\\ (\Delta d)_{p}~{}&=~{}(\Delta
u)_{n}~{}=~{}-0.427\,,\\\ (\Delta s)_{p}~{}&=~{}(\Delta
s)_{n}~{}=~{}-0.085\,.\end{split}$ (13)

Table 1: Scalar nuclear form factors $f_{Tq}^{N}$ used in DM@NLO based on Ref. Hoferichter:2015dsa , DarkSUSY 6.4 based on Ref. Gasser:1990ce  and MicrOMEGAs 5.3 Belanger:2018ccd . Scalar coefficient | DM@NLO | DarkSUSY | MicrOMEGAs
---|---|---|---
$f_{Tu}^{p}$ | 0.0208 | 0.023 | 0.0153
$f_{Tu}^{n}$ | 0.0189 | 0.019 | 0.0110
$f_{Td}^{p}$ | 0.0411 | 0.034 | 0.0191
$f_{Td}^{n}$ | 0.0451 | 0.041 | 0.0273
$f_{Ts}^{p}=f_{Ts}^{n}$ | 0.043 | 0.14 | 0.0447
$f_{Tc}^{p}=f_{Tb}^{p}=f_{Tt}^{p}$ | 0.0663 | 0.0595 | 0.0682
$f_{Tc}^{n}=f_{Tb}^{n}=f_{Tt}^{n}$ | 0.0661 | 0.0592 | 0.0679

For a detailed discussion of direct detection in NLO SUSY-QCD, we refer the
reader to Ref. Klasen:2016qyz .

##  3 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

$\Delta\sigma^{\mathrm{NLO}}~{}=~{}\int_{n}\differential{\sigma^{\mathrm{V}}}+\int_{n+1}\differential{\sigma^{\mathrm{R}}}\,,$
(14)

which contains virtual ($\differential{\sigma^{\mathrm{V}}}$) and real
($\differential{\sigma^{\mathrm{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 $\alpha_{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. Herrmann:2007ku ;
Herrmann:2009wk ; Herrmann:2009mp .

###  3.1 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 Siegel:1979wq ;
Signer:2008va  is employed, i.e. the relevant loop integrals are evaluated in
$D=4-2\epsilon$ 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 down-type squarks share the
common soft breaking parameter $M_{\tilde{q}}$ as a result of the
$\operatorname{SU}(2)_{L}$ gauge symmetry. The squark mass matrix can be
diagonalised,

$\displaystyle U^{\tilde{q}}\begin{pmatrix}m^{2}_{LL}&m^{2}_{LR}\\\
m^{2}_{RL}&m^{2}_{RR}\end{pmatrix}(U^{\tilde{q}})^{\dagger}=\begin{pmatrix}m^{2}_{\tilde{q}_{1}}&0\\\
0&m^{2}_{\tilde{q}_{2}}\end{pmatrix}$ (15)

with the two physical masses $m^{2}_{\tilde{q}_{1}}$ and
$m^{2}_{\tilde{q}_{2}}$ being the eigenvalues of the non-diagonal mass matrix
with the entries

$\displaystyle m^{2}_{LL}=M^{2}_{\tilde{Q}}+(I^{3L}_{q}-e_{q}s_{W}^{2})\cos
2\beta m_{Z}^{2}+m_{q}^{2}\,,$ (16) $\displaystyle
m^{2}_{RR}=M^{2}_{\tilde{U},\tilde{D}}+e_{q}s_{W}^{2}\cos 2\beta
m_{Z}^{2}+m_{q}^{2}\,,$ (17) $\displaystyle
m^{2}_{LR}=m^{2}_{RL}=m_{q}(A_{q}-\mu(\tan\beta)^{-2I^{3L}_{q}})\,.$ (18)

Out of the eleven parameters $M_{\tilde{Q}}$, $M_{\tilde{U}}$,
$M_{\tilde{D}}$, $A_{t}$, $A_{b}$, $m_{\tilde{t}_{1}}$, $m_{\tilde{t}_{2}}$,
$m_{\tilde{b}_{1}}$, $m_{\tilde{b}_{2}}$, $\theta_{\tilde{t}}$ and
$\theta_{\tilde{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_{\tilde{Q}}$, $M_{\tilde{U}}$, $M_{\tilde{D}}$ as input parameters by the
physical on-shell masses $m_{\tilde{b}_{1}}$, $m_{\tilde{b}_{2}}$ and
$m_{\tilde{t}_{1}}$. 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,

$\displaystyle m^{2}_{LL}$
$\displaystyle=\frac{m^{2}_{\tilde{q}_{1}}+m^{2}_{\tilde{q}_{2}}}{2}\pm\frac{1}{2}\sqrt{(m_{\tilde{q}_{1}}-m_{\tilde{q}_{2}})^{2}-4m_{LR}^{4}}\,,$
(19) $\displaystyle m^{2}_{RR}$
$\displaystyle=\frac{m^{2}_{\tilde{q}_{1}}+m^{2}_{\tilde{q}_{2}}}{2}\mp\frac{1}{2}\sqrt{(m_{\tilde{q}_{1}}-m_{\tilde{q}_{2}})^{2}-4m_{LR}^{4}}\,,$
(20)

for the diagonal entries of the mass matrix. Consequently, there are two
possible values for $M_{\tilde{Q}}$, $M_{\tilde{U}}$, and $M_{\tilde{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 $\delta
m_{\tilde{t}_{2}}\sim(U_{21}^{\tilde{q}}U_{12}^{\tilde{q}})^{-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
$\delta\theta_{\tilde{q}}\sim(U_{11}^{\tilde{q}}U_{22}^{\tilde{q}}+U_{12}^{\tilde{q}}U_{21}^{\tilde{q}})^{-1}$.
To avoid these issues, a scheme is defined as numerically stable if the
following three conditions are fulfilled:

  * $\bullet$

$|m^{\rm out}_{\tilde{t}_{1}}-m^{\rm in}_{\tilde{t}_{1}}|/m^{\rm
in}_{\tilde{t}_{1}}<10^{-5}$ ,

  * $\bullet$

$|\mathrm{Re}\,U_{11}^{\tilde{q}}U_{21}^{\tilde{q}}|>10^{-4}$ ,

  * $\bullet$

$|\mathrm{Re}\,(U_{11}^{\tilde{t}}U_{22}^{\tilde{t}}+U_{12}^{\tilde{t}}U_{21}^{\tilde{t}})|>10^{-4}$
.

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_{{\tilde{t}_{2}}}$ is closer to the corresponding
physical value.

In a series of previous analyses Harz:2012fz ; Herrmann:2014kma ; Harz:2014tma
; Schmiemann:2019czm ; Klasen:2022ptb , the following three renormalisation
schemes, adapted to the present situation of DM (co)annihilation, have been
introduced:

  * 0:

$m_{b}$, $m_{t}$, $m_{\tilde{f}}$, $\theta_{\tilde{f}}$, $A_{f}$ are all
$\smash{\overline{\mathrm{DR}}}$ parameters.

  * 1:

$m_{b}$, $A_{b}$ and $A_{t}$ are $\smash{\overline{\mathrm{DR}}}$ input
parameters whereas $m_{t}$, $m_{\tilde{t}_{1}}$ $m_{\tilde{b}_{1}}$ and
$m_{\tilde{b}_{2}}$ are on-shell (OS) masses. $\theta_{\tilde{t}}$,
$\theta_{\tilde{b}}$ and $m_{{\tilde{t}_{2}}}$ are then dependent quantities.

  * 2:

$m_{t}$, $m_{b}$, $A_{b}$ and $A_{t}$ are $\smash{\overline{\mathrm{DR}}}$
input parameters and $m_{\tilde{t}_{1}}$, $m_{\tilde{b}_{1}}$ and
$m_{\tilde{b}_{2}}$ are OS masses. $\theta_{\tilde{t}}$, $\theta_{\tilde{b}}$
and $m_{{\tilde{t}_{2}}}$ are then dependent quantities.

The hybrid on-shell/$\smash{\overline{\mathrm{DR}}}$ scheme 1, which resembles
the RS2 scheme presented in Ref. Heinemeyer:2010mm , 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. Heinemeyer:2023pcc  is left for a future update.

###  3.2 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 Catani:1996vz ; Catani:2002hc  for
massive initial-state particles Harz:2022ipe  over phase-space slicing methods
Fabricius:1981sx ; Giele:1991vf ; Harris:2001sx . This subtraction technique
is based on the introduction of a local counterterm
$\differential{\sigma^{\mathrm{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 $\varepsilon^{-1}$ and
$\varepsilon^{-2}$. The whole procedure can be schematically captured in the
equation

$\displaystyle\begin{split}\Delta\sigma^{\mathrm{NLO}}~{}&=~{}\int_{n}\left[\differential{\sigma^{\mathrm{V}}}+\int_{1}\differential{\sigma^{\mathrm{A}}}\right]_{\varepsilon=0}\\\
&~{}~{}~{}~{}~{}~{}+\int_{n+1}\Big{[}\differential{\sigma_{\varepsilon=0}^{\mathrm{R}}}-\differential{\sigma_{\varepsilon=0}^{\mathrm{A}}}\Big{]}\end{split}$
(21)

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 scans111For historical
reasons, squark annihilation into electroweak final states is computed using
the phase space slicing method. For stop annihilaton 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
$\tilde{t}_{i}\tilde{t}^{\ast}_{j}\to gg$ and
$\tilde{t}_{i}\tilde{t}^{\ast}_{j}\to q\bar{q}$ at the loop level. To avoid
double-counting, 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..

###  3.3 Intermediate on-shell resonance subtraction

Within the real emission contribution to the process
$\tilde{\chi}^{0}_{n}\tilde{t}_{i}\to bW^{+}$, the internal top propagator can
become on-shell if the collisional energy $\sqrt{s}$ exceeds the top mass
$m_{t}$. To cure the singularity, we follow the "Prospino scheme" defined in
Refs. Beenakker:1996ch ; Binoth:2011xi  and substitute the top propagator with
the Breit-Wigner form, according to

$\frac{1}{p^{2}-m^{2}_{t}}~{}\to~{}\frac{1}{p^{2}-m^{2}_{t}+im_{t}\Gamma_{t}}\,,$
(22)

in the resonant part $\mathcal{M}_{\rm r}$ of the total amplitude
$\mathcal{M}_{\rm tot}=\mathcal{M}_{\rm r}+\mathcal{M}_{\rm nr}$, whereas the
non-resonant piece $\mathcal{M}_{\rm nr}$ remains unchanged. Since the
corresponding process $\tilde{\chi}^{0}_{n}\tilde{t}_{i}\to 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

$\big{|}\mathcal{M}_{\rm r}\big{|}^{2}~{}\to~{}\big{|}\mathcal{M}_{\rm
r}\big{|}^{2}-\frac{m^{2}_{t}\Gamma^{2}_{t}}{(p^{2}-m^{2}_{t})^{2}+m^{2}_{t}\Gamma^{2}_{t}}\big{|}\mathcal{M}_{\rm
r}\big{|}^{2}_{p_{t}^{2}=m_{t}^{2}}$ (23)

with the physical top width $\Gamma_{t}$. This procedure has the advantage
that it retains the interference $\mathrm{Re}(\mathcal{M}^{\ast}_{\rm
r}\mathcal{M}_{\rm 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 $\Gamma_{t}=10^{-3}\cdot
m_{t}$ in the interference part instead of the physical width.

###  3.4 Sommerfeld enhancement

The Sommerfeld enhancement is an elementary quantum mechanical effect that
increases (decreases) annihilation cross sections for small relative
velocities in the 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 $\lambda$ to the initial particles. More
quantitatively, the Sommerfeld factor

$S_{0}^{[\mathbf{R}]}~{}=~{}\big{|}\phi^{[\mathbf{R}]}(0)\big{|}^{2}$ (24)

is obtained as a solution $\phi(\vec{r})$ evaluated at the origin $\vec{r}=0$
of the stationary Schrödinger equation for the potential describing the
interaction of the annihilating particles transforming under the
representation $\mathbf{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

$\sigma^{\mathrm{Som}}=\sum_{\mathbf{R}}S_{0}^{[\mathbf{R}]}\sigma^{\mathrm{Tree}}_{\mathbf{R}}\,,$
(25)

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
$|\phi_{0}^{[\mathbf{R}]}(0)|^{2}=1$ to ensure that
$\sigma^{\mathrm{Som}}\to\sigma^{\mathrm{Tree}}$ is fulfilled if the
interaction governing the enhancement effect is turned off ($\lambda\to 0$).
When combining the Sommerfeld effect with the full $\order{\lambda^{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 $\order{\lambda^{2}}$ contribution from the
Sommerfeld factor.

##  4 Installing and running DM@NLO

###  4.1 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 Fuks:2013vua ; Fiaschi:2023tkq . It is publicly available for
download at https://dmnlo.hepforge.org and 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
Hahn:2006nq  and LoopTools-2.16 Hahn:1998yk  are required for reading particle
spectra following the Supersymmetry Les Houches Accord 2 (SLHA 2) convention
Skands:2003cj ; Allanach:2008qq  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 Hahn:2004fe  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
        ./dmnlo input/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 LoopTools 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 Buckley:2013jua  and allow to reproduce the most important cross
section plots.

Table 2: List of (co)annihilation and elastic DM-nucleon scattering processes included in DM@NLO, given together with the location of the corresponding source code in run_dmnlo, the references to the orignal publication and whether the Sommerfeld enhancement is included. Here, $\phi=\\{h^{0},H^{0},A^{0},H^{\pm}\\}$, $V=\\{Z^{0},W^{\pm},\gamma\\}$, $\bar{V}=V\setminus\\{\gamma\\}$, and $\ell$ ($\bar{\ell}$) can be any (anti)lepton. The indices can take the values $\\{m,n\\}=\\{1,2,3,4\\}$, $\\{i,j\\}=\\{1,2\\}$. Process | Folder | References | Sommerfeld
---|---|---|---
$\tilde{\chi}^{0}_{m}\tilde{\chi}^{0}_{n},\tilde{\chi}^{\pm}_{i}\tilde{\chi}^{\pm}_{j},\tilde{\chi}^{0}_{n}\tilde{\chi}^{\pm}_{i}\to q\bar{q},q\bar{q}^{\prime}$ | ChiChi2QQ | Herrmann:2007ku ; Herrmann:2009wk ; Herrmann:2009mp ; Herrmann:2014kma  | ✗
$\tilde{\chi}^{0}_{n}\tilde{q}_{i}\to q^{\prime}\phi,q^{\prime}\bar{V},q^{\prime}g$ with $q,q^{\prime}\in\\{t,b\\}$ | NeuQ2qx | Harz:2012fz ; Harz:2014tma ; Harz:2022ipe  | ✗
$\tilde{t}_{1}\tilde{t}^{\ast}_{1}\to VV,V\phi,\phi\phi,\ell\bar{\ell}$ | QQ2xx | Harz:2014gaa  | ✓
$\tilde{q}_{i}\tilde{q}^{\prime}_{j}\to qq^{\prime}$ with $q,q^{\prime}\in\\{t,b\\}$ | stst2QQ | Schmiemann:2019czm  | ✓
$\tilde{\tau}_{1}\tilde{\tau}^{\ast}_{1}\to t\bar{t}$ | staustau2QQ | Branahl:2019yot  | ✓
$\tilde{t}_{i}\tilde{t}^{\ast}_{j}\to gg,q\bar{q}$ with $q\in\\{u,d,c,s\\}$ | stsT2xx | Klasen:2022ptb  | ✓
$\tilde{\chi}^{0}_{1}N\to\tilde{\chi}^{0}_{1}N$ | DD | Klasen:2016qyz  | –

###  4.2 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 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 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 $\mathrm{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
$\smash{\overline{\mathrm{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_{{\tilde{t}_{2}}}$ is closest to the corresponding on-
shell 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
$\theta_{W}$ and the $W$-mass as in the default MSSM model file in MicrOMEGAs
2.4.1, i.e. $\sin\theta_{W}=0.481$ and $m_{W}=\cos\theta_{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 $Z$\- and $W$-mass from the SLHA 2 file

$\sin^{2}\theta_{W}~{}=~{}1-\frac{m^{2}_{W}}{m^{2}_{Z}}\,.$ (26)

Note that the legacy option should only be used for the reproduction of
previously published results.

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 ParticleDataGroup:2022pth . 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 $\sqrt{s}$ has to be defined with \--pcm,
which is the center-of-mass momentum $p_{\mathrm{cm}}$ of the incoming
particles. The option \--result controls whether the output contains the total
cross section $\sigma$ or the cross section times velocity $\sigma v$, both in
units of $\mathrm{G}\mathrm{eV}^{-2}$, where the relative velocity is defined
as $v=2\lambda^{1/2}(s,m_{a}^{2},m_{b}^{2})/s$ with $\lambda$ 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
$\mathrm{c}\mathrm{m}^{2}$, respectively. The scalar nuclear form factors from
Tab. 1 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. 1.

###  4.3 The DM@NLO library

To facilitate the usage of DM@NLO from within other codes, the static library
libdmnlo.a provides 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.a also 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.cpp located in external/micromegas_5.3.41/MSSM
exemplifies the use of these functions with MicrOMEGAs in a way that only
those channels contributing more than $2\text{\,}\mathrm{\char 37\relax}$ 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_example Scenario.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.

##  5 Illustrative example calculations

Table 3: Example scenario in the cMSSM with a positive Higgs supersymmetric mixing parameter $\mu$ where stop (co)annihliation is the dominant dark matter mechanism. All dimensionful quantities are in $\mathrm{GeV}$. $m_{0}$ | $m_{1/2}$ | $\tan\beta$ | $A_{0}$ | $m_{\tilde{\chi}^{0}_{1}}$ | $m_{\tilde{t}_{1}}$
---|---|---|---|---|---
3000 | 1400 | 20 | 12000 | 606.3 | 648.3

[width=0.49]figures/cMSSM_GG.pdf [width=0.49]figures/cMSSM_Gt.pdf
[width=0.49]figures/cMSSM_Wb.pdf [width=0.49]figures/cMSSM_tt.pdf
[width=0.49]figures/cMSSM_Zt.pdf [width=0.49]figures/cMSSM_ht.pdf

Figure 1:  Tree-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 $\mu_{R}$ by a factor of two around
the central scale as shaded bands. The upper part of each plot shows the
absolute value of $\sigma v$ together with the thermal velocity distribution
(in arbitrary units), whereas the lower part shows the corresponding relative
shift (second item in the legend).

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 of the neutral components of the two Higgs
doublets $\tan\beta$, the universal trilinear coupling $A_{0}$, and the sign
of the Higgs mixing parameter $\mu$.

Table 4:  Dominant annihilation channels contributing to $\langle\sigma_{\rm ann}v\rangle$ for the cMSSM scenario in Tab. 3. Further contributions below 2% are omitted. Channel |  Contribution
---|---
$\tilde{t}_{1}\ \tilde{t}_{1}^{\ast}\rightarrow g\ g$ |  $36\text{\,}\mathrm{\char 37\relax}$
$\tilde{\chi}^{0}_{1}\ \tilde{t}_{1}\rightarrow t\ g$ |  $29\text{\,}\mathrm{\char 37\relax}$
$\tilde{\chi}^{0}_{1}\ \tilde{t}_{1}\rightarrow W^{+}\ b$ |  $8\text{\,}\mathrm{\char 37\relax}$
$\tilde{t}_{1}\ \tilde{t}_{1}\rightarrow t\ t$ |  $6\text{\,}\mathrm{\char 37\relax}$
$\tilde{t}_{1}\ \tilde{t}_{1}^{\ast}\rightarrow\gamma\ g$ |  $4\text{\,}\mathrm{\char 37\relax}$
$\tilde{\chi}^{0}_{1}\ \tilde{t}_{1}\rightarrow Z^{0}\ t$ |  $4\text{\,}\mathrm{\char 37\relax}$
$\tilde{\chi}^{0}_{1}\ \tilde{t}_{1}\rightarrow h^{0}\ t$ |  $4\text{\,}\mathrm{\char 37\relax}$
DM@NLO total |  $87\text{\,}\mathrm{\char 37\relax}$

In the following, we use, inspired by the recent search for non-excluded
regions in the cMSSM parameter space Ellis:2022emx , the parameter point given
in Tab. 3, where stop (co)annihilation is the dominant dark-matter
annihilation mechanism222Note that numerical differences in the physical mass
spectrum occur with respect to Ref. Ellis:2022emx  since SPheno 4.0.5 is used
as spectrum generator in this work whereas Ref. Ellis:2022emx  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. Ellis:2022emx .. The
(co)annihilation channels contributing most to $\langle\sigma_{\rm
ann}v\rangle$ are listed in Tab. 4 with stop-antistop annihilation into gluons
having the largest relative contribution. For each channel that can be
corrected with DM@NLO, we show our tree-level (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_{\rm cm}$ region the cross-section contributes to the
total annihilation cross-section $\langle\sigma_{\rm ann}v\rangle$. 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\bar{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
$\mu_{R}=\sqrt{m_{\tilde{t}_{1}}m_{\tilde{t}_{2}}}$, which for the particular
scenario in Tab. 3 yields $\mu_{R}=1368.2$ GeV, whereas MicrOMEGAs 5.3.41 uses
the scale $Q=2m_{\tilde{\chi}^{0}_{1}}/3$ for the evaluation of the strong
coupling stored in the global variable GGscale333This scale choice is
different from MicrOMEGAs 2.4.1, where $Q=2m_{\tilde{\chi}^{0}_{1}}$ is used..
We also show the uncertainties from variations of the renormalisation scale
$\mu_{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 $\Omega_{\chi}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 ($\Omega^{\rm MO}_{\chi}$) indicates the region
consistent with the observed value $\Omega_{\rm 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 ($\Omega^{\rm
Tree}_{\chi}$) corresponds to the prediction where the DM@NLO tree-level cross
sections replace the CalcHEP result, and the yellow band ($\Omega^{\rm
Full}_{\chi}$) is based on our full calculation.

The width of the three bands reflects the experimental $2\sigma$ 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 tree-level 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
$\Omega_{\rm 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.

[width=.49]figures/cMSSM_m12_m0.pdf [width=.49]figures/cMSSM_mneut_mstop.pdf

Figure 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 $(\Omega_{\chi}^{\rm Tree}-\Omega_{\chi}^{\rm
Full})/\Omega_{\chi}^{\rm Tree}$ in the relic density compared to our tree-
level result.

[width=0.49]figures/cMSSM_m0_m12_Oh2_SI_p_QCDoff.pdf
[width=0.49]figures/cMSSM_m0_m12_Oh2_SI_n_QCDoff.pdf
[width=0.49]figures/cMSSM_m0_m12_Oh2_SD_p_QCDoff.pdf
[width=0.49]figures/cMSSM_m0_m12_Oh2_SD_n_QCDoff.pdf

Figure 3: Spin-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.

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 $\order{\alpha_{s}}$ corrections to the
dominant effective operators (blue solid line), MicrOMEGAs (orange solid line)
and the corresponding analytic tree-level calculation444More precisely, the
orange solid line corresponds to the MicrOMEGAs function nucleonAmplitudes,
which includes additionally the box corrections calculated in Ref.
Drees:1993bu , while the green dotted line corresponds to the output of the
function MSSMDDtest(loop=0,…) with the additional QCD corrections turned off
QCDcorrections=0..

We also show all three values of the resulting relic density ($\Omega^{\rm
MO}_{\chi}$, $\Omega^{\rm Tree}_{\chi}$, $\Omega^{\rm Full}_{\chi}$) 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.

##  6 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 cross-sections 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.

###### Acknowledgements.

We would like to thank all colleagues who contributed to DM@NLO over the
years, i.e. (in alphabetical order), J. Branahl, Q. Le Boulc’h, M. Meinecke,
S. Schmiemann, P. Steppeler, and M. Younes Sassi. We would also like to thank
A. Pukhov for useful discussions throughout the development of DM@NLO, and A.
Puck Neuwirth for useful technical discussions on the build system. The
research of MK and LPW was supported by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) through the Research Training Group 2149
“Strong and Weak Interactions - from Hadrons to Dark matter”. This work was
also supported by _Investissements d’avenir_ , Labex ENIGMASS, contrat
ANR-11-LABX-0012. JH acknowledges support from the Emmy Noether grant
“Baryogenesis, Dark Matter and Neutrinos: Comprehensive analyses and accurate
methods in particle cosmology” (HA 8555/1-1, Project No. 400234416) funded by
the DFG and the Cluster of Excellence “Precision Physics, Fundamental
Interactions, and Structure of Matter” (PRISMA+ EXC 2118/1), funded by the DFG
within the German Excellence Strategy (Project No. 390831469).

##  Appendix A The DM@NLO configuration file

In this appendix, the options available in a DM@NLO configuration file are
described. Such an input file consists out of a series of keywords that can be
set to user defined values. These keywords are given by:

  * •

slha = <string>: path to the SLHA input file defining the SUSY scenario to
investigate (mass spectrum, mixing matrices, decay widths, etc.).

  * •

muR = <double>: renormalisation scale $\mu_{R}$ in $\mathrm{GeV}$.

  * •

renscheme = <int>: renormalisation scheme according to the numbering scheme
introduced in Sec. 3.1. The mixed $\smash{\overline{\mathrm{DR}}}$-OS scheme
no. 1 is the recommended option.

  * •

choosesol = <int>: solutions for $M_{\tilde{Q}}$, $M_{\tilde{U}}$,
$M_{\tilde{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.

  * •

pcm = <double>: centre-of-mass momentum in $\mathrm{GeV}$.

  * •

result = <string>: defines whether the output should contain the total cross
$\sigma$ corresponding to the value s or the total cross section times
relative velocity $\sigma v$ defined through the value sv.

  * •

formfactor = <int>: sets the scalar nuclear form factors $f^{N}_{Tq}$ to one
of the sets of values in Tab. 1.

##  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 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 $\mathrm{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_{\tilde{Q}}$, $M_{\tilde{U}}$, $M_{\tilde{D}}$ as explained in Sec. 3.1.

  * •

\--legacy, defines the weak mixing angle $\theta_{W}$ and the $W$-mass as in
the default MSSM model file in MicrOMEGAs 2.4.1.

  * •

\--lo, returns the result at LO accuracy.

  * •

\--nlo, returns the result at NLO accuracy.

  * •

\--sommerfeld, returns the Sommerfeld enhanced cross section.

  * •

\--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_{\rm cm}$ in $\mathrm{GeV}$.

  * •

\--result, followed by a string corresponding to s for the total cross
$\sigma$ or sv for the total cross section times the relative velocity $\sigma
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_{Tq}^{N}$ to one of the value sets shown in
Tab. 1 with zero for DM@NLO and two for MicrOMEGAs.

## References

  * (1) Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * (2) P. Gondolo and G. Gelmini, _Cosmic abundances of stable particles: Improved analysis_ , _Nucl. Phys. B_ 360 (1991) 145–179.
  * (3) G. Jungman, M. Kamionkowski and K. Griest, _Supersymmetric dark matter_ , _Phys. Rept._ 267 (1996) 195–373, [hep-ph/9506380].
  * (4) A. Freitas, _Radiative corrections to co-annihilation processes_ , _Phys. Lett. B_ 652 (2007) 280–284, [0705.4027].
  * (5) B. Herrmann and M. Klasen, _SUSY-QCD Corrections to Dark Matter Annihilation in the Higgs Funnel_ , _Phys. Rev. D_ 76 (2007) 117704, [0709.0043].
  * (6) N. Baro, F. Boudjema and A. Semenov, _Full one-loop corrections to the relic density in the MSSM: A Few examples_ , _Phys. Lett. B_ 660 (2008) 550–560, [0710.1821].
  * (7) N. Baro, F. Boudjema and A. Semenov, _Automatised full one-loop renormalisation of the MSSM. I. The Higgs sector, the issue of tan(beta) and gauge invariance_ , _Phys. Rev. D_ 78 (2008) 115003, [0807.4668].
  * (8) N. Baro and F. Boudjema, _Automatised full one-loop renormalisation of the MSSM II: The chargino-neutralino sector, the sfermion sector and some applications_ , _Phys. Rev. D_ 80 (2009) 076010, [0906.1665].
  * (9) B. Herrmann, M. Klasen and K. Kovarik, _Neutralino Annihilation into Massive Quarks with SUSY-QCD Corrections_ , _Phys. Rev. D_ 79 (2009) 061701, [0901.0481].
  * (10) B. Herrmann, M. Klasen and K. Kovarik, _SUSY-QCD effects on neutralino dark matter annihilation beyond scalar or gaugino mass unification_ , _Phys. Rev. D_ 80 (2009) 085025, [0907.0030].
  * (11) J. Harz, B. Herrmann, M. Klasen, K. Kovarik and Q. L. Boulc’h, _Neutralino-stop coannihilation into electroweak gauge and Higgs bosons at one loop_ , _Phys. Rev. D_ 87 (2013) 054031, [1212.5241].
  * (12) J. Harz, B. Herrmann, M. Klasen and K. Kovarik, _One-loop corrections to neutralino-stop coannihilation revisited_ , _Phys. Rev. D_ 91 (2015) 034028, [1409.2898].
  * (13) B. Herrmann, M. Klasen, K. Kovarik, M. Meinecke and P. Steppeler, _One-loop corrections to gaugino (co)annihilation into quarks in the MSSM_ , _Phys. Rev. D_ 89 (2014) 114012, [1404.2931].
  * (14) J. Harz, B. Herrmann, M. Klasen, K. Kovařík and M. Meinecke, _SUSY-QCD corrections to stop annihilation into electroweak final states including Coulomb enhancement effects_ , _Phys. Rev. D_ 91 (2015) 034012, [1410.8063].
  * (15) S. Schmiemann, J. Harz, B. Herrmann, M. Klasen and K. Kovařík, _Squark-pair annihilation into quarks at next-to-leading order_ , _Phys. Rev. D_ 99 (2019) 095015, [1903.10998].
  * (16) J. Branahl, J. Harz, B. Herrmann, M. Klasen, K. Kovařík and S. Schmiemann, _SUSY-QCD corrected and Sommerfeld enhanced stau annihilation into heavy quarks with scheme and scale uncertainties_ , _Phys. Rev. D_ 100 (2019) 115003, [1909.09527].
  * (17) M. Klasen, K. Kovařík and L. P. Wiggering, _Radiative corrections to stop-antistop annihilation into gluons and light quarks_ , _Phys. Rev. D_ 106 (2022) 115032, [2210.05260].
  * (18) G. Bélanger, V. Bizouard, F. Boudjema and G. Chalons, _One-loop renormalization of the NMSSM in SloopS: The neutralino-chargino and sfermion sectors_ , _Phys. Rev. D_ 93 (2016) 115031, [1602.05495].
  * (19) G. Bélanger, V. Bizouard, F. Boudjema and G. Chalons, _One-loop renormalization of the NMSSM in SloopS. II. The Higgs sector_ , _Phys. Rev. D_ 96 (2017) 015040, [1705.02209].
  * (20) M. Klasen, C. E. Yaguna and J. D. Ruiz-Alvarez, _Electroweak corrections to the direct detection cross section of inert higgs dark matter_ , _Phys. Rev. D_ 87 (2013) 075025, [1302.1657].
  * (21) S. Banerjee, F. Boudjema, N. Chakrabarty, G. Chalons and H. Sun, _Relic density of dark matter in the inert doublet model beyond leading order: The heavy mass case_ , _Phys. Rev. D_ 100 (2019) 095024, [1906.11269].
  * (22) S. Banerjee, F. Boudjema, N. Chakrabarty and H. Sun, _Relic density of dark matter in the inert doublet model beyond leading order for the low mass region: 1. Renormalisation and constraints_ , _Phys. Rev. D_ 104 (2021) 075002, [2101.02165].
  * (23) R. Iengo, _Sommerfeld enhancement: General results from field theory diagrams_ , _JHEP_ 05 (2009) 024, [0902.0688].
  * (24) M. Beneke, A. Bharucha, F. Dighera, C. Hellmann, A. Hryczuk, S. Recksiegel et al., _Relic density of wino-like dark matter in the MSSM_ , _JHEP_ 03 (2016) 119, [1601.04718].
  * (25) J. Harz and K. Petraki, _Higgs Enhancement for the Dark Matter Relic Density_ , _Phys. Rev. D_ 97 (2018) 075041, [1711.03552].
  * (26) B. von Harling and K. Petraki, _Bound-state formation for thermal relic dark matter and unitarity_ , _JCAP_ 12 (2014) 033, [1407.7874].
  * (27) J. Harz and K. Petraki, _Radiative bound-state formation in unbroken perturbative non-Abelian theories and implications for dark matter_ , _JHEP_ 07 (2018) 096, [1805.01200].
  * (28) S. Biondini, _Bound-state effects for dark matter with Higgs-like mediators_ , _JHEP_ 06 (2018) 104, [1805.00353].
  * (29) J. Harz and K. Petraki, _Higgs-mediated bound states in dark-matter models_ , _JHEP_ 04 (2019) 130, [1901.10030].
  * (30) J. Harz, B. Herrmann, M. Klasen, K. Kovarik and P. Steppeler, _Theoretical uncertainty of the supersymmetric dark matter relic density from scheme and scale variations_ , _Phys. Rev. D_ 93 (2016) 114023, [1602.08103].
  * (31) A. Hryczuk and R. Iengo, _The one-loop and Sommerfeld electroweak corrections to the Wino dark matter annihilation_ , _JHEP_ 01 (2012) 163, [1111.2916].
  * (32) A. Hryczuk, I. Cholis, R. Iengo, M. Tavakoli and P. Ullio, _Indirect Detection Analysis: Wino Dark Matter Case Study_ , _JCAP_ 07 (2014) 031, [1401.6212].
  * (33) M. Baumgart, T. Cohen, I. Moult, N. L. Rodd, T. R. Slatyer, M. P. Solon et al., _Resummed Photon Spectra for WIMP Annihilation_ , _JHEP_ 03 (2018) 117, [1712.07656].
  * (34) M. Baumgart, T. Cohen, E. Moulin, I. Moult, L. Rinchiuso, N. L. Rodd et al., _Precision Photon Spectra for Wino Annihilation_ , _JHEP_ 01 (2019) 036, [1808.08956].
  * (35) M. Drees and M. M. Nojiri, _New contributions to coherent neutralino - nucleus scattering_ , _Phys. Rev. D_ 47 (1993) 4226–4232, [hep-ph/9210272].
  * (36) M. Drees and M. Nojiri, _Neutralino - nucleon scattering revisited_ , _Phys. Rev. D_ 48 (1993) 3483–3501, [hep-ph/9307208].
  * (37) J. Hisano, S. Matsumoto, M. M. Nojiri and O. Saito, _Direct detection of the Wino and Higgsino-like neutralino dark matters at one-loop level_ , _Phys. Rev. D_ 71 (2005) 015007, [hep-ph/0407168].
  * (38) A. Berlin, D. S. Robertson, M. P. Solon and K. M. Zurek, _Bino variations: Effective field theory methods for dark matter direct detection_ , _Phys. Rev. D_ 93 (2016) 095008, [1511.05964].
  * (39) M. Klasen, K. Kovarik and P. Steppeler, _SUSY-QCD corrections for direct detection of neutralino dark matter and correlations with relic density_ , _Phys. Rev. D_ 94 (2016) 095002, [1607.06396].
  * (40) S. Bisal, A. Chatterjee, D. Das and S. A. Pasha, _Radiative Corrections to Aid the Direct Detection of the Higgsino-like Neutralino Dark Matter: Spin-Independent Interactions_ , 2311.09937.
  * (41) T. Abe and R. Sato, _Quantum corrections to the spin-independent cross section of the inert doublet dark matter_ , _JHEP_ 03 (2015) 109, [1501.04161].
  * (42) D. Azevedo, M. Duch, B. Grzadkowski, D. Huang, M. Iglicki and R. Santos, _One-loop contribution to dark-matter-nucleon scattering in the pseudo-scalar dark matter model_ , _JHEP_ 01 (2019) 138, [1810.06105].
  * (43) K. Ishiwata and T. Toma, _Probing pseudo Nambu-Goldstone boson dark matter at loop level_ , _JHEP_ 12 (2018) 089, [1810.08139].
  * (44) K. Ghorbani and P. H. Ghorbani, _Leading Loop Effects in Pseudoscalar-Higgs Portal Dark Matter_ , _JHEP_ 05 (2019) 096, [1812.04092].
  * (45) T. Abe, M. Fujiwara and J. Hisano, _Loop corrections to dark matter direct detection in a pseudoscalar mediator dark matter model_ , _JHEP_ 02 (2019) 028, [1810.01039].
  * (46) S. Glaus, M. Mühlleitner, J. Müller, S. Patel and R. Santos, _Electroweak corrections to dark matter direct detection in the dark singlet phase of the N2HDM_ , _Phys. Lett. B_ 833 (2022) 137342, [2204.13145].
  * (47) S. Glaus, M. Mühlleitner, J. Müller, S. Patel and R. Santos, _Electroweak Corrections to Dark Matter Direct Detection in a Vector Dark Matter Model_ , _JHEP_ 10 (2019) 152, [1908.09249].
  * (48) T. Biekötter, P. Gabriel, M. O. Olea-Romacho and R. Santos, _Direct detection of pseudo-Nambu-Goldstone dark matter in a two Higgs doublet plus singlet extension of the SM_ , _JHEP_ 10 (2022) 126, [2207.04973].
  * (49) A. Berlin, D. Hooper and S. D. McDermott, _Dark matter elastic scattering through Higgs loops_ , _Phys. Rev. D_ 92 (2015) 123531, [1508.05390].
  * (50) F. Ertas and F. Kahlhoefer, _Loop-induced direct detection signatures from CP-violating scalar mediators_ , _JHEP_ 06 (2019) 052, [1902.11070].
  * (51) C. Borschensky, G. Coniglio, B. Jäger, J. Jochum and V. Schipperges, _Direct detection of dark matter: Precision predictions in a simplified model framework_ , _Eur. Phys. J. C_ 81 (2021) 44, [2008.04253].
  * (52) U. Haisch and F. Kahlhoefer, _On the importance of loop-induced spin-independent interactions for dark matter direct detection_ , _JCAP_ 04 (2013) 050, [1302.4454].
  * (53) A. Crivellin, F. D’Eramo and M. Procura, _New Constraints on Dark Matter Effective Theories from Standard Model Loops_ , _Phys. Rev. Lett._ 112 (2014) 191304, [1402.1173].
  * (54) A. Arbey and F. Mahmoudi, _SUSY constraints from relic density: High sensitivity to pre-BBN expansion rate_ , _Phys. Lett. B_ 669 (2008) 46–51, [0803.0741].
  * (55) C. E. Yaguna, _Large contributions to dark matter annihilation from three-body final states_ , _Phys. Rev. D_ 81 (2010) 075024, [1003.2730].
  * (56) K. Saikawa and S. Shirai, _Precise WIMP Dark Matter Abundance and Standard Model Thermodynamics_ , _JCAP_ 08 (2020) 011, [2005.03544].
  * (57) B. C. Allanach, S. Kraml and W. Porod, _Theoretical uncertainties in sparticle mass predictions from computational tools_ , _JHEP_ 03 (2003) 016, [hep-ph/0302102].
  * (58) M. Beneke, F. Dighera and A. Hryczuk, _Relic density computations at NLO: infrared finiteness and thermal correction_ , _JHEP_ 10 (2014) 045, [1409.3049].
  * (59) T. Binder, T. Bringmann, M. Gustafsson and A. Hryczuk, _Early kinetic decoupling of dark matter: when the standard way of calculating the thermal relic density fails_ , _Phys. Rev. D_ 96 (2017) 115010, [1706.07433].
  * (60) A. Aboubrahim, M. Klasen and L. P. Wiggering, _Forbidden dark matter annihilation into leptons with full collision terms_ , _JCAP_ 08 (2023) 075, [2306.07753].
  * (61) J. Hamann, S. Hannestad, M. S. Sloth and Y. Y. Y. Wong, _How robust are inflation model and dark matter constraints from cosmological data?_ , _Phys. Rev. D_ 75 (2007) 023522, [astro-ph/0611582].
  * (62) T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg and P. Walia, _Converting nonrelativistic dark matter to radiation_ , _Phys. Rev. D_ 98 (2018) 023543, [1803.03644].
  * (63) B. W. Lee and S. Weinberg, _Cosmological lower bound on heavy-neutrino masses_ , _Phys. Rev. Lett._ 39 (Jul, 1977) 165–168.
  * (64) J. Edsjo and P. Gondolo, _Neutralino relic density including coannihilations_ , _Phys. Rev. D_ 56 (1997) 1879–1894, [hep-ph/9704361].
  * (65) G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov and B. Zaldivar, _micrOMEGAs5.0 : Freeze-in_ , _Comput. Phys. Commun._ 231 (2018) 173–186, [1801.03509].
  * (66) A. Arbey and F. Mahmoudi, _SuperIso Relic: A Program for calculating relic density and flavor physics observables in Supersymmetry_ , _Comput. Phys. Commun._ 181 (2010) 1277–1292, [0906.0369].
  * (67) F. Ambrogi, C. Arina, M. Backovic, J. Heisig, F. Maltoni, L. Mantani et al., _MadDM v.3.0: a Comprehensive Tool for Dark Matter Studies_ , _Phys. Dark Univ._ 24 (2019) 100249, [1804.00044].
  * (68) T. Bringmann, J. Edsjö, P. Gondolo, P. Ullio and L. Bergström, _DarkSUSY 6 : An Advanced Tool to Compute Dark Matter Properties Numerically_ , _JCAP_ 07 (2018) 033, [1802.03399].
  * (69) T. Binder, T. Bringmann, M. Gustafsson and A. Hryczuk, _Dark matter relic abundance beyond kinetic equilibrium_ , _Eur. Phys. J. C_ 81 (2021) 577, [2103.01944].
  * (70) J. Fan, M. Reece and L.-T. Wang, _Non-relativistic effective theory of dark matter direct detection_ , _JCAP_ 11 (2010) 042, [1008.1591].
  * (71) R. J. Hill and M. P. Solon, _Standard Model anatomy of WIMP dark matter direct detection I: weak-scale matching_ , _Phys. Rev. D_ 91 (2015) 043504, [1401.3339].
  * (72) R. J. Hill and M. P. Solon, _Standard Model anatomy of WIMP dark matter direct detection II: QCD analysis and hadronic matrix elements_ , _Phys. Rev. D_ 91 (2015) 043505, [1409.8290].
  * (73) J. Hisano, R. Nagai and N. Nagata, _Effective Theories for Dark Matter Nucleon Scattering_ , _JHEP_ 05 (2015) 037, [1502.02244].
  * (74) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, _Remarks on Higgs Boson Interactions with Nucleons_ , _Phys. Lett. B_ 78 (1978) 443–446.
  * (75) M. Hoferichter, J. Ruiz de Elvira, B. Kubis and U.-G. Meißner, _High-Precision Determination of the Pion-Nucleon $\sigma$ Term from Roy-Steiner Equations_, _Phys. Rev. Lett._ 115 (2015) 092301, [1506.04142].
  * (76) J. Gasser, H. Leutwyler and M. E. Sainio, _Sigma term update_ , _Phys. Lett. B_ 253 (1991) 252–259.
  * (77) W. Siegel, _Supersymmetric Dimensional Regularization via Dimensional Reduction_ , _Phys. Lett. B_ 84 (1979) 193–196.
  * (78) A. Signer and D. Stockinger, _Using Dimensional Reduction for Hadronic Collisions_ , _Nucl. Phys. B_ 808 (2009) 88–120, [0807.4424].
  * (79) S. Heinemeyer, H. Rzehak and C. Schappacher, _Proposals for Bottom Quark/Squark Renormalization in the Complex MSSM_ , _Phys. Rev. D_ 82 (2010) 075010, [1007.0689].
  * (80) S. Heinemeyer and F. von der Pahlen, _Automated choice for the best renormalization scheme in BSM models_ , _Eur. Phys. J. C_ 83 (2023) 865, [2302.12187].
  * (81) S. Catani and M. H. Seymour, _A General algorithm for calculating jet cross-sections in NLO QCD_ , _Nucl. Phys. B_ 485 (1997) 291–419, [hep-ph/9605323].
  * (82) S. Catani, S. Dittmaier, M. H. Seymour and Z. Trocsanyi, _The Dipole formalism for next-to-leading order QCD calculations with massive partons_ , _Nucl. Phys. B_ 627 (2002) 189–265, [hep-ph/0201036].
  * (83) J. Harz, M. Klasen, M. Y. Sassi and L. P. Wiggering, _Dipole formalism for massive initial-state particles and its application to dark matter calculations_ , _Phys. Rev. D_ 107 (2023) 056020, [2210.03409].
  * (84) K. Fabricius, I. Schmitt, G. Kramer and G. Schierholz, _Higher Order Perturbative QCD Calculation of Jet Cross-Sections in e+ e- Annihilation_ , _Z. Phys. C_ 11 (1981) 315\.
  * (85) W. T. Giele and E. W. N. Glover, _Higher order corrections to jet cross-sections in e+ e- annihilation_ , _Phys. Rev. D_ 46 (1992) 1980–2010.
  * (86) B. W. Harris and J. F. Owens, _The Two cutoff phase space slicing method_ , _Phys. Rev. D_ 65 (2002) 094032, [hep-ph/0102128].
  * (87) W. Beenakker, R. Hopker, M. Spira and P. M. Zerwas, _Squark and gluino production at hadron colliders_ , _Nucl. Phys. B_ 492 (1997) 51–103, [hep-ph/9610490].
  * (88) T. Binoth, D. Goncalves Netto, D. Lopez-Val, K. Mawatari, T. Plehn and I. Wigmore, _Automized Squark-Neutralino Production to Next-to-Leading Order_ , _Phys. Rev. D_ 84 (2011) 075005, [1108.1250].
  * (89) B. Fuks, M. Klasen, D. R. Lamprea and M. Rothering, _Precision predictions for electroweak superpartner production at hadron colliders with Resummino_ , _Eur. Phys. J. C_ 73 (2013) 2480, [1304.0790].
  * (90) J. Fiaschi, B. Fuks, M. Klasen and A. Neuwirth, _Electroweak superpartner production at 13.6 Tev with Resummino_ , _Eur. Phys. J. C_ 83 (2023) 707, [2304.11915].
  * (91) T. Hahn, _SUSY Les Houches Accord 2 I/O made easy_ , _Comput. Phys. Commun._ 180 (2009) 1681–1693, [hep-ph/0605049].
  * (92) T. Hahn and M. Perez-Victoria, _Automatized one loop calculations in four-dimensions and D-dimensions_ , _Comput. Phys. Commun._ 118 (1999) 153–165, [hep-ph/9807565].
  * (93) P. Z. Skands et al., _SUSY Les Houches accord: Interfacing SUSY spectrum calculators, decay packages, and event generators_ , _JHEP_ 07 (2004) 036, [hep-ph/0311123].
  * (94) B. C. Allanach et al., _SUSY Les Houches Accord 2_ , _Comput. Phys. Commun._ 180 (2009) 8–25, [0801.0045].
  * (95) T. Hahn, _CUBA: A Library for multidimensional numerical integration_ , _Comput. Phys. Commun._ 168 (2005) 78–95, [hep-ph/0404043].
  * (96) A. Buckley, _PySLHA: a Pythonic interface to SUSY Les Houches Accord data_ , _Eur. Phys. J. C_ 75 (2015) 467, [1305.4194].
  * (97) Particle Data Group collaboration, R. L. Workman et al., _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * (98) J. Ellis, K. A. Olive, V. C. Spanos and I. D. Stamou, _The CMSSM survives Planck, the LHC, LUX-ZEPLIN, Fermi-LAT, H.E.S.S. and IceCube_ , _Eur. Phys. J. C_ 83 (2023) 246, [2210.16337].


# Raspberry Pi Bee Health Monitoring Device

1st Jakub Nevláčil, 2nd Šimon Bilík, 3rd Karel Horák  Department of Control
and Instrumentation
Faculty of Electrical Engineering and Communication
Brno University of Technology
Brno, Czech Republic
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A declining honeybee population could pose a threat to a food resources of the
whole world one of the latest trend in beekeeping is an effort to monitor a
health of the honeybees using various sensors and devices. This paper
participates on a development on one of these devices. The aim of this paper
is to make an upgrades and improvement of an in-development bee health
monitoring device and propose a remote data logging solution for a continual
monitoring of a beehive.

###### Index Terms:

Appis melifera, Remote monitoring, Raspberry Pi, IoT, MQTT

##  I Introduction

Honeybees are crucial for our ecosystem and also for food production around
the globe. One of the main purposes of the honey bees is to pollinate flowers
and plants and produce honey. The pollination of the plants plays a key role
in the agriculture, because without the bees the plants could not reproduce
and grow its fruits or other crops.[1] Beekeeping is a thousands of years old
craft and as times goes the bees appear to be more and more endangered.
Whether it is endangering by pests and diseases or today especially by
pesticides from nearby agricultural fields.[2]

This is the reason why it is important for the beekeepers to be aware of the
honeybees health and activity in real time and in case of an irregular bee
behavior react upon this danger as soon as possible. Therefore this project
aims to modify a device that collects useful data from a beehive to be able to
log this data to a database where the data could be visualized using a
graphical interface. Later the data will be used to determine the health and
the activity of the bees using machine learning algorithms.

##  II Data acquisition device

###  II-A Original data acquisition device

The device for data acquisition is based on a previous work of [4] and [3].
The original concept of the data acquisition device can be seen in the Fig. 1.
This device is based on a Raspberry Pi platform with a Grove Base Hat for
connecting sensors to capture weather and atmospheric conditions inside and
outside of the beehive. Furthermore this device also monitors the traffic of
the incoming and outcoming bees using the built-in camera module and takes
sound samples from the inside of the beehive with a microphone.

Figure 1: Original data acquisition device

###  II-B Modified data acquisition device

Although the data acquisition from the sensors and microphone usually occurs
every 15 minutes (the interval can be adjusted, but there is most likely no
use in collecting the data more frequently), the SGP30 sensor that measures
CO2 concentration and TVOC (total volatile organic compounds) inside the
beehive needs to be activated every second in order to provide correct
readings. This task can be unnecessarily demanding on computing power, which
can be later utilized for later implementation of the machine learning
algorithms for determining health and activity of the bees.

For this purpose a Raspberry Pico was incorporated into the data acquisition
chain. It was connected to the Raspberry Pi over USB port and serial
communication between the two devices was established. This allows to transfer
all the weather and atmospheric data from the sensors to the Raspberry Pico
using the corresponding sensor Hat for this device. The Raspberry Pico is
programmed in such a way that it acquires data from all the weather and
atmospheric sensors every second. Then once the Raspberry Pi asks for the
measured data, the Raspberry Pico calculates an average value over the inter-
request period of the data from every sensor and sends them to the Raspberry
Pi. The outcome of this modification is that the Raspberry Pi does not have to
dedicate an entire computational thread to call the SGP30 sensor every second
and can use this thread for other purposes. It should also increase
reliability of the camera capturing.

Figure 2: Modified data acquisition device

In the process of modifying, the DHT10 sensor which measures the temperature
and humidity inside of the hive was updated to a newer version DHT20. This
modification should bring a higher measurement accuracy and a greater
measurement range. In addition to the weather and atmospheric sensors was
added a HX711 amplifier with 4 load cells in full bridge configuration. This
sensor will later be used for measuring the weight changes of individual
beehives. A scheme of our modified data acquisition device in its final
configuration can be seen from Fig. 2. In the Fig. 3 we can then see a real
implementation of the data acquisition device inside a box with a tunnels for
the passage of the bees.

Figure 3: Modified data acquisition device inside a box

##  III Logging data to a database

###  III-A Internet connection

After acquiring the useful data from the hive, we want to send them over the
internet to a database where we can look at the and assess the state of the
hive from anywhere we like. This of course requires an internet connection.
Due to the fact that beehives are usually placed on remote locations, where
Wi-Fi or wired internet connection is not an option, this can be a problem.
Luckily the place of deployment of our devices should have sufficient Wi-Fi
coverage. In situations where a Wi-Fi connection cannot be established,
solutions as GSM module or LoRaWAN could be used easily.

###  III-B MQTT messages

The device currently stores all the acquired data from the hive on an USB
flash drive. This includes sound samples from a hive, pictures of the bees and
data from the weather and atmospheric sensors. For the purpose of assessing
the state of health of the beehive, only the data from the weather sensors
needs to be logged to the database. In future this will also include the
outputs from a machine learning models that should quantify the health of the
beehive and activity of the bees. Due to this fact that only a few variables
needs to be transmitted to a database, the MQTT messages were chosen for this
purpose.

MQTT is a protocol used in IoT. One of the main advantages of the MQTT
protocol is its lightweight design. Usually this protocol is used on a
microcontrollers. This means that the MQTT protocol requires minimal
computational resources and also the internet connection requirements are
relatively low.[5] The MQTT space consists of two kind of devices, MQTT
clients and a MQTT broker. MQTT clients can either subscribe or publish
messages. MQTT broker controls the flow of the messages by receiving messages
from publishing clients and forwarding them to the clients that are subscribed
to a particular topic. For an easier understanding an example of a MQTT
network can be seen in Fig. 4. The MQTT broker can be hosted on a cloud or it
can be installed on a local device.

In our particular case the MQTT broker will be cloud based. This solution has
been assess as more suitable for our particular problem. The upside of this
solution is that it doesn’t require a public IP address or a VPN server on the
side of the subscribed client. In conclusion our setup will consist of a
Raspberry Pi data acquisition device which will be configured as publishing
client, then cloud based MQTT broker and one Raspberry Pi which will be
configured as a subscribing client. The subscribing Raspberry Pi will also act
as a database and will run a web application for displaying the logged data.

###  III-C InfluxDB database

As mentioned above, the MQTT messages will be sent from the data acquisition
device to an MQTT broker. All of the data from the sensors representing a
specific period of time will be sent in a single message stored as a JSON
object. The MQTT broker will the distribute every message to the subscribed
clients. In our case the subscribed client will be a Raspberry Pi device where
the messages will be received a processed by Telegraf, which is a data
collection agent. Telegraf the takes the message sorts the data from
individual sensors and saves them into an InfluxDB database. The reason the
InfluxDB database is used instead of other popular databases like MySQL is
because The Influx DB database is a time series database. This means that this
database is designed to store a data which are associated with time. This
makes it the best choice for our application of plotting the measurements from
sensors over time.

Figure 4: MQTT architercture[6]

###  III-D Grafana

Once the measurements from the sensors are saved in a database, we can them
plot the values from individual sensors over time. For this purpose Grafana
will be used. Grafana is a web application for visualizing data and metrics.
This application will also run on the Raspberry Pi next to the InfluxDB
database and Telegraf and will read metrics from the sensirs directly from the
database and plot them into a graphs. The Grafana application can be accessed
from a web browser on a device that is connected to the same local network as
the device it is installed on. Grafana can also send notifications on Slack
and email. This feature could be useful once the machine learning models for
recognition of bee activity and health will be implemented to warn the
beekeepers about emergency situations in the hive.

Figure 5: Grafana dashboard

In the Fig. 5 we can see a Grafana dashboard with two sets of plots. This is a
data from two prototypes of the data acquisition device that was described
above. This data was collected during testing phase before deployment of the
devices on a real beehives.

##  IV Conclusion

In this paper we modified previous work of [4] to use a Raspberry Pico for
collecting data from weather and atmospheric sensors to release some
computational resources of a Raspberry Pi. An upgrade of a sensor was made for
a newer version with better measuring accuracy and also a new weight sensor
for measuring a weight changes of a beehive.

Then for an easier remote monitoring of the beehives a data logging solution
was implemented. This solution starts at MQTT messages that are being sent
from the data acquisition device to a cloud based MQTT broker from where the
messages travel to a server with Influx DB database installed. As the last
step the data is visualized using a web application Grafana, from where the
visualized data can be accessed through a web browser.

In the future work an analysis of a sensor data fusion collected during a
deployment on real beehives is planned. This analysis should include
especially cluster analysis of a bee buzz audio samples, but also other
machine learning techniques applied to both audio data and also sensor data.

## Acknowledgment

The completion of this paper was made possible by the grant No. FEKT-S-23-8451
- ”Research on advanced methods and technologies in cybernetics, robotics,
artificial intelligence, automation and measurement” financially supported by
the Internal science fund of Brno University of Technology.

## References

  * [1] Why are Bees Important?. Friend of the Earth [online]. [cit. 2023-03-12]. Available at: https://foe.org/blog/why-are-bees-important/
  * [2] MEIKLE, William G., Guy MERCADIER, Fatiha GUERMACHE a Marie-Claude BON. Pseudomonas contamination of a fungus-based biopesticide: Implications for honey bee (Hymenoptera. Biological Control [online]. 2012, 60(3), 312-320 [cit. 2023-03-12]. ISSN 10499644. Available at: https://doi.org/10.1016/j.biocontrol.2011.12.004
  * [3] MVG-BVI: Bee Visual Inspector [online]. In: . s. 2 [cit. 2023-03-12]. Available at: http://vision.uamt.feec.vutbr.cz/data/MVG_2021_BVI_tech_spec.pdf
  * [4] Bee-Health-Monitor. GitHub[online]. 2022 [cit. 2022-12-15].
Available at: https://github.com/boortel/Bee-Health-Monitor.

  * [5] MQTT. MQTT [online]. [cit. 2023-03-09]. Available at: https://mqtt.org/
  * [6] IT Explained: MQTT. In: Paessler [online]. [cit. 2023-03-10]. Available at: https://www.paessler.com/it-explained/mqtt