# Impact of QCD uncertainties on antiproton spectra from dark-matter
annihilation

Adil Jueid  Jochem Kip  Roberto Ruiz de Austri  and Peter Skands

###### Abstract

Dark-matter particles that annihilate or decay can undergo complex sequences
of processes, including strong and electromagnetic radiation, hadronisation,
and hadron decays, before particles that are stable on astrophysical time
scales are produced. Antiprotons produced in this way may leave footprints in
experiments such as AMS–02. Several groups have reported an excess of events
in the antiproton flux in the rigidity range of $10$–$20$ GV. However, the
theoretical modeling of baryon production is not straightforward and relies in
part on phenomenological models in Monte Carlo event generators. In this work,
we assess the impact of QCD uncertainties on the spectra of antiprotons from
dark-matter annihilation. As a proof-of-principle, we show that for a two-
parameter model that depends only on the thermally-averaged annihilation cross
section ($\langle\sigma v\rangle$) and the dark-matter mass ($M_{\chi}$), QCD
uncertainties can affect the best-fit mass by up to $\sim 14\%$ (with large
uncertainties for large DM masses), depending on the choice of $M_{\chi}$ and
the annihilation channel ($b\bar{b}$ or $W^{+}W^{-}$), and $\langle\sigma
v\rangle$ by up to $\sim 10\%$. For comparison, changes to the underlying
diffusion parameters are found to be within $1\%$–$5\%$, and the results are
also quite resilient to the choice of cosmic-ray propagation model. These
findings indicate that QCD uncertainties need to be included in future DM
analyses. To facilitate full-fledged analyses, we provide the spectra in
tabulated form including QCD uncertainties and code snippets to perform mass
interpolations and quick DM fits. The code can be found in this github
repository.

##  1 Introduction

In the Cosmological Standard Model, dark matter (DM) constitutes about $85\%$
of the matter in the universe (see e.g. [1]). Analyses of the cosmic microwave
background prefer DM particles that are nonrelativistic in the era of galaxy
formation. This requirement is relatively straightforward to realize in
particle-physics models by extending the Standard Model (SM) with weakly
interacting massive particles (WIMPs). Such scenarios can reach excellent
agreement with the measured relic density, $\Omega_{\rm DM}h^{2}=0.1188\pm
0.0010$ [2], and typically predict that WIMPs can annihilate (or decay) into
SM states, which — via often complex sequences of decays, radiation, and
hadronisation — produce observable final-state particles such as photons,
anti-protons, antineutrinos, or positrons. Contributions of such DM
annihilation products to cosmic-ray (CR) fluxes may be observable in
experiments like the Fermi Large Area Telescope (Fermi–LAT), or the Alpha
Magnetic Spectrometer (AMS–02).
After the discovery of secondary CR antiprotons in 1979, an excess over the
astrophysical background was observed in the antiproton-to-proton ratio [3, 4]
which at the time was explained by a photino DM with mass of a few GeV [5].
Despite the unprecedented precision on the measurement of the antiproton flux
performed by the AMS–02 collaboration [6], the excess seems to be still alive
especially in the rigidity range $10$–$20$ GV (see e.g. [7, 8, 9, 10]).
Interestingly, DM interpretations of this excess can be consistent with those
of the so-called Galactic Center Excess (GCE) in gamma rays [11, 12, 13, 14,
15, 16, 17, 18, 19, 20]. In other words, a DM of mass around $50$–$200$ GeV
and thermal annihilation cross section of about $\sim 10^{-26}~{}{\rm
cm}^{2}~{}{\rm s}^{-1}$ can address both of these anomalies. The statistical
uncertainties on the antiproton flux are now considered very subleading,
making the treatment of systematic uncertainties and their correlations
crucial to further progress in DM analyses [21, 22, 23, 24, 25, 26, 27]. It is
in this context that we note that, for DM masses above a few GeV, antiprotons
are produced through Quantum Chromodynamics (QCD) jet fragmentation. Since the
problem of confinement remains fundamentally unsolved, it is addressed either
via phenomenological fits such as Fragmentation Functions (FFs), and/or
explicit dynamical models embodied by Monte Carlo (MC) event generators. Both
approaches involve free parameters (see, e.g., [28, 29, 30]). The associated
uncertainties on the particle fluxes in dark-matter annihilation have,
however, so far been neglected in the literature. As experimental systematics
at the few-percent level have caused important shifts in the antiproton
excess, it is relevant to examine whether QCD uncertainties could be of a
comparable magnitude to this.

In this article, we perform a dedicated analysis with state-of-the-art MC
tools which shows that QCD uncertainties can alter, and even dominate, the
particle-physics component of the overall error budget for antiproton DM fits.
This is quite different to the case of gamma rays for which QCD uncertainties
are below $10\%$ in the peak region [31] (see also [32, 33] for short
summaries). We expect this to have a significant impact on DM interpretations
of the AMS–02 excess. We first perform several fits of the fragmentation-
function parameters in Pythia 8.244 [34] to salient LEP measurements, starting
from the default “Monash” tune [29]. We then discuss the different sources of
QCD uncertainties on the antiproton production from DM annihilation and show
that, in our analysis they are of similar size as the experimental errors on
the AMS–02 data, especially in the low rigidity range. As a first step towards
quantifying the impact on DM fits, we study the effects on the best-fit points
of both the mass of the DM as well as the thermally-averaged annihilation
cross section for a few annihilation channels and for fitted diffusion
parameters of the CR propagation model defined in DRAGON2. We encourage
complementary measurements to be developed, e.g., at the Large Hadron Collider
(LHC), which could be used to further improve the quality of the theoretical
models for baryon production. The rest of this article is organised as
follows. In section 2, we discuss the modeling of antiproton spectra in Pythia
8 and briefly highlight differences in the theory predictions of the state-of-
art multi-purpose MC event generators: Herwig 7, Pythia 8 and Sherpa 2\. In
section 3, we present the results of a new tune of the Lund fragmentation-
function parameters in Pythia 8 and estimate the QCD uncertainties for a few
DM masses. In section 4, we discuss the impact on DM observables. Finally we
conclude in section 5.

##  2 Antiprotons from DM annihilation

For a generic production of a final state $f$111$f$ can be any parton-level
final state, typically pairs such as
$f\in[\gamma\gamma,\ell\ell,q\bar{q},t\bar{t},gg,WW,ZZ,hh,\ldots]$. Since this
work focuses on antiprotons, we disregard leptonic final states and QED-like
processes connected to them. in DM annihilation or decay, the general
expression of the antiproton flux at the position of Earth can be written as
[35]

$\displaystyle\frac{{\rm d}\Phi_{\bar{p}}}{{\rm d}E_{\rm kin}}(E_{\rm
kin},r_{\odot})\equiv\frac{v_{\bar{p}}}{4\pi}\bigg{(}\frac{\rho_{\odot}}{M_{\chi}}\bigg{)}^{k}\mathcal{R}(E_{\rm
kin})\sum_{f}k^{-1}\mathcal{P}_{i\to f}(g_{\chi};M_{\chi})\bigg{(}\frac{{\rm
d}N}{{\rm d}E_{\rm kin}}\bigg{)}_{f\to\bar{p}},$ (2.1)

where $k=2~{}(1)$ for annihilation (decay) process, $\rho_{\odot}$ sets the
normalisation of the DM density in terms of its reference value at the Sun,
$E_{\rm kin}$ is the kinetic energy of the antiproton and $v_{\bar{p}}$ is its
velocity. $\mathcal{R}(E_{\rm kin})$ encodes all the astrophysical factors on
both the production and the propagation of antiprotons including the shape of
the DM distribution. Particle physics enters in the last two terms in the
right hand side of equation (2.1): (i) The model-dependent term is encoded in
$\mathcal{P}_{i\to f}(g_{\chi};M_{\chi})$ that gives the transition rate for
DM into the final state $f$ which can be either the thermal annihilation cross
section today ($\langle\sigma v\rangle_{\chi\chi\to f}$) for the annihilation
or the partial decay width ($\Gamma_{\chi\to f}$) for the decay of DM and (ii)
The model-independent part which represents the differential distribution of
the antiproton on the kinetic energy. The model-independent contribution to
the antiproton flux is usually computed using MC event generators. Coloured
particles produced in DM processes will undergo QCD bremsstrahlung wherein
additional coloured particles are produced and whose multiplicity depend on
how far the parent particles are from their production threshold. The rate of
QCD processes is controlled by the effective value of the strong coupling
constant evaluated a scale proportional to the shower evolution variable 222We
note that the value of $\alpha_{S}(M_{Z})$ in Pythia 8 is larger than
$\alpha_{S}(M_{Z})^{\overline{{\rm MS}}}$ by about $20\%$ [36, 29].. At a
scale $Q_{\rm IR}\simeq\mathcal{O}(1)~{}{\rm GeV}\ll Q_{\rm UV}$, any coloured
particle must hadronise to produce a set of colourless hadrons. This process,
called fragmentation is modeled within Pythia 8 with the Lund string model
[37, 38, 39] whose longitudinal properties (in lightcone coordinates) are
encoded in the following fragmentation function:

$\displaystyle f_{\mathrm{Lund}}(z,m_{\perp h})$ $\displaystyle\propto$
$\displaystyle\frac{(1-z)^{a}}{z}\exp\left(\frac{-bm_{\perp
h}^{2}}{z}\right)~{},$ (2.2)

with $m_{\perp h}\equiv\sqrt{m_{h}^{2}+p_{\perp h}^{2}}$ being the “transverse
mass” of the hadron $h$ with mass $m_{h}$ and transverse momentum $p_{\perp
h}$ relative to the string axis. In the standard Lund model, the latter is
assumed to follow a 2D Gaussian distribution. For each iterative string
fragmentation process, $f_{\mathrm{Lund}}(z)$ gives the probability that a
hadron $h$ with $m_{\perp h}$ and energy fraction $z\in[0,1]$ is produced. For
mesons, it depends on three tunable parameters: the longitudinal $a$ and $b$
parameters in equation (2.2), and $\sigma_{\perp}\approx\sqrt{\langle p_{\perp
h}^{2}\rangle}/\sqrt{2}$ which controls the transverse-momentum distribution.
For baryons, a fourth parameter enters, with $a\to a+a_{\rm QQ}$. The $a$ and
$b$ parameters act to suppress the relative rates of very hard and very soft
hadrons, respectively, and in practice are highly correlated [31]. For this
reason a new parameterization of $f_{\mathrm{Lund}}(z)$ was introduced wherein
$b$ is replaced by $\langle z_{\rho}\rangle$, the average longitudinal
momentum fraction taken by a $\rho$ meson; $\langle
z_{\rho}\rangle=\int_{0}^{1}{\rm d}z\ zf(z,\langle m_{\perp\rho}\rangle)$.

Figure 1: Comparison between the theory predictions and the experimental
measurements of $p/\bar{p}$ scaled-momentum (left) and $\Lambda^{0}$ Log of
scaled momentum at LEP. Pythia 8 with the default Monash tune is shown in red,
Sherpa 2 is shown in blue for the cluster (solid) and Lund (dashed)
hadronisation models while the result of Herwig 7 is shown in green for dipole
(solid) and angular-based (dashed) parton-shower algorithms. Data is taken
from Opal [40] (left) and Aleph [41] (right).

The fragmentation-function parameters cannot be computed from first
principles. But exploiting the fact that hadronisation occurs at very long
distance scales as compared to high-energy scattering processes, the
factorization theorem holds, and measurements of hadron spectra at e.g. the
Large Electron Positron Collider (LEP) can be used to constrain them. The same
principle implies that those parameters can then also be used to make
predictions for hadronic DM annihilation processes (at least for DM masses
above a few GeV); this is sometimes referred to as jet universality. At LEP,
hadronic $Z$ decays produced (anti-)protons either directly from the
fragmentation of quarks/gluons (dubbed primary) or from decays of heavier
baryons (dubbed secondary). Secondary (anti)protons come mainly from decays of
five baryons: $\Lambda^{0},\Delta^{++},\Delta^{+},\Delta^{0}$ and
$\Sigma^{\pm}$, with the following fractions to the total number of produced
antiprotons: $R_{\Lambda^{0}\to p}\approx 22\%,R_{\Delta^{++}\to p}\approx
10\%,R_{\Delta^{+}\to p}\approx 8\%,R_{\Delta^{0}\to p}\approx 3\%$ and
$R_{\Sigma^{\pm}\to p}\approx 5\%$333This finding also holds true for
antiprotons from dark matter annihilation except that there is a significant
contribution from antineutron decays, which are stable at LEP time scales
while they are not in dark-matter indirect detection experiments.
Nevertheless, if one considers the antiprotons that are not produced from
antineutron decays (roughly $50\%$ of the total), the contribution from the
other sources is almost independent of the dark matter mass and annihilation
final state.. We close this section by discussing the agreement between the
theory predictions of three state-of-art MC event generators and the
experimental measurements of the baryon spectra at LEP (proton and
$\Lambda^{0}$). In figure 1, we show the comparison between the predictions of
(i) Pythia 8 version 307 [42] with the default Monash tune [29], (ii) Sherpa
version 2.2.12 [43] with the CSS shower model [44] using the cluster
hadronisation model [45] and the Lund string model based on Pythia 6 [46] and
(iii) Herwig version 7.2.3 [47] with two radiation models: the angular [48]
and the dipole-type [49, 50] algorithms for two key observables used in our
fits. We can see that the predictions of Herwig 7 and Sherpa based on the
cluster models disagree with both the Pythia 8 prediction as well as with the
experimental measurements. Since the disagreement can reach up to $50\%$, we
conclude that the envelope spanned by the differenc MC event generators cannot
define a faithful estimate of the QCD uncertainties on the baryon spectra.

Parameter |  Pythia8 setting | Monash | This work | $a_{\rm q}$ | $\langle z_{\rho}\rangle$ |  $\sigma_{\perp}$ [GeV] | $a_{QQ}$
---|---|---|---|---|---|---|---
$a_{\rm q}$ | `StringZ:aLund` | $0.68$ | $0.601_{-0.038}^{+0.038}$ | $1.000$ |  |  |
$\langle z_{\rho}\rangle$ | `StringZ:avgZLund` | $(0.55)$ | $0.540_{-0.004}^{+0.004}$ | $0.718$ | $1.000$ |  |
$\sigma_{\perp}$ [GeV] | `StringPT:Sigma` | $0.335$ | $0.307_{-0.002}^{+0.002}$ | $0.057$ | $-0.270$ | $1.000$ |
$a_{QQ}$ | `StringZ:aExtraDiquark` | $0.97$ | $1.671_{-0.196}^{+0.196}$ | $0.415$ | $0.816$ | $-0.204$ | $1.000$
$b$ | `StringZ:bLund` | $0.98$ | $(0.897)$ |  |  |  |
Table 1:  The fit results for the parameters of the fragmentation function in
Pythia 8\. The errors on the parameters correspond to the Migrad errors
estimated at the fit. For comparison we show the corresponding values of the
parameters in the baseline Monash tune. The last four columns show the
symmetric correlation matrix. The goodness-of-fit per degrees of freedom
($\chi^{2}/N_{\rm df}$) is $\chi^{2}/N_{\rm df}=676.69/852$ where $N_{\rm
df}=N_{\rm bins}-N_{\rm params}$.

##  3 QCD uncertainties on antiproton fluxes

Figure 2: Antiproton differential distribution in $x\equiv E_{\rm
kin}/M_{\chi}$ for dark matter annihilation into $q\bar{q}$ for
$M_{\chi}=10~{}{\rm GeV}$ (red), $M_{\chi}=100~{}{\rm GeV}$ (blue) and
$M_{\chi}=1000~{}{\rm GeV}$ (dark green). The three lower panels show the
ratio to the nominal prediction for each DM mass. The nominal prediction for
each case were estimated using the results in Table 1. The light (dark) error
bands correspond to uncertainties on the predictions from hadronisation
(parton showers). The hadronisation uncertainties correspond to the $2\sigma$
eigentunes (see the text). Here, we assume that the dark matter is
annihilating to all the SM quarks but the top quark with equal probabilities,
i.e. ${\rm BR}(\chi\chi\to q\bar{q})\equiv\frac{\sigma(\chi\chi\to
q\bar{q})}{\sum_{q}\sigma(\chi\chi\to q\bar{q}}=20\%$ for $q=u,d,s,c,b$.

To assess properly the QCD uncertainties on the antiproton spectra, we first
perform several retunings of Pythia 8 using a set of constraining data from
LEP. A complete global analysis would involve measurements not only of
proton/antiproton spectra but also those of hyperons and $\Delta$ baryons
which decay to protons. There are, however, no measurements of $\Delta$ or
$\Sigma$ baryon spectra with any significant constraining power, so for
antiproton production via decays we are forced to focus solely on the
$\Lambda$ component. To ensure that the fits remain consistent with overall
event properties as well, we also include spectra of light mesons, mean
charged-particle multiplicities and multiplicity distributions, and the thrust
and $C$-parameter event shapes, using all relevant measurements available in
Rivet 3.1.3 [51], specifically those by Aleph [52, 53, 41, 54], Delphi [55,
56, 57, 58, 40], L3 [59] and Opal [60, 61, 62, 63, 64]. The fit is performed
using Professor 2.3.3 [65], with goodness-of-fit measure

$\chi^{2}=\sum_{\mathcal{O}}\sum_{b\in\mathcal{O}}\bigg{(}\frac{{\rm
MC}_{(b)}(\\{p_{i}\\})-{\rm Data}_{(b)}}{\Delta_{b}}\bigg{)}^{2},$ (3.1)

where ${\rm MC}_{(b)}(\\{p_{i}\\})$ denotes the MC prediction for observable
$\mathcal{O}$ at a bin $b$, which depends on the parameter set
$\\{p_{i}\\}=\\{a_{\rm q},\langle z_{\rho}\rangle,\sigma_{\perp},a_{QQ}\\}$,
cast as a fourth-order polynomial interpolation (see e.g. [31] for details),
and $\Delta_{b}$ represents the total error which is estimated as the
quadratic sum of the MC statistical errors, experimental errors and a $5\%$
flat theory uncertainty; $\Delta_{b}=\sqrt{\Delta_{\rm MC,b}^{2}+\Delta_{\rm
Data,b}^{2}+(0.05\times{\rm MC}_{(b)})^{2}}$. The results of the tunes along
with the corresponding correlation matrix are shown in Table 1. We have
verified that the predictions of Pythia 8 at the best-fit point agree very
well with both the results from the baseline Monash tune and with the bulk of
the experimental data.

The QCD uncertainties on the antiproton spectra can be categorised into: (i)
perturbative uncertainties arising from parton-shower variations and (ii)
hadronisation uncertainties on the parameters of the fragmentation function
obtained at the best-fit point. The shower uncertainties are estimated using
the automated method developed in [66] and include a set of next-to-leading
order correction terms which act to reduce the renormalization-scale
variations. For the hadronisation uncertainties, we use the Professor toolkit
to estimate the Hessian errors which can be computed from the diagonalisation
of the $\chi^{2}$ covariance matrix near the minimum (details can be found in
e.g. [67]). One gets $2N_{\rm param}=8$ variations, known as eigentunes, and
one can associate $1\sigma$, $2\sigma$ and $3\sigma$ eigentunes to variations
satisfying $\Delta\chi^{2}\equiv\chi^{2}-\chi^{2}_{\rm min}\simeq N_{\rm df}$,
$4N_{\rm df}$, and $9N_{\rm df}$ respectively. Given that there are a priori
unknown parametric uncertainties on the theory side and that a “good”
uncertainty estimate should at least roughly cover the range of the
experimental uncertainties, we find that the $2\sigma$ eigentunes are well
suited to obtain reasonably conservative estimates, hence we will use those as
our baseline in the rest of this article. To illustrate the effects of the QCD
uncertainties, we calculate the antiproton spectra for DM annihilation into
$q\bar{q}$ choosing three DM masses of $10,100$ and $1000$ GeV. The results
are shown in figure 2 where we can see that hadronisation uncertainties can go
from $10$–$25\%$ in the peak regions to about $50$–$200\%$ in the high
$x$-regions. The shower uncertainties compete with hadronisation errors in the
peak and for $M_{\chi}\geqslant 100~{}{\rm GeV}$ while they are smaller for
$x\geqslant 0.3$.

##  4 Effects on dark matter observables

In this section we quantify to what extent the QCD uncertainties affect two DM
observables: the velocity-weighted cross section $\langle\sigma v\rangle$ and
the DM mass $M_{\chi}$. A total shift in the height of the spectrum, such that
it is still contained within the QCD uncertainties of the spectrum,
corresponds to an uncertainty on $\langle\sigma v\rangle$. The uncertainty on
the DM mass, $\Delta M_{\chi}$, for a DM particle of mass $M_{\chi}$ can be
determined by finding the masses of DM particles whose spectra including QCD
uncertainties match the spectrum of the DM particle with mass $M_{\chi}$. For
both DM uncertainties we determine the upper and lower bounds by demanding
that $\chi^{2}/N_{\text{df}}\approx 1$ holds. Thus the spectra of the
parameters corresponding to the upper and lower bounds of the DM uncertainties
can give the same spectra as the original $\langle\sigma v\rangle$ or
$M_{\chi}$ when the QCD uncertainties are accounted for. We note that the
uncertainty in $\langle\sigma v\rangle$ is equivalent to an uncertainty in the
DM density, $\rho_{\odot}$. The exact conversion between $\langle\sigma
v\rangle$ and $\rho_{\odot}$ however depends on the specific DM halo profile,
thus we shall only use $\langle\sigma v\rangle$ in the following.

Figure 3: The nominal value of two DM particles with a mass of 100 GeV
annihilating into antiprotons via a $W^{+}W^{-}$ pair is shown in blue. The
upper and lower limit of the DM masses whose antiproton spectra, again via
$W^{+}W^{-}$, encompass the nominal 100 GeV spectra in their uncertainties is
shown in green and orange respectively. The peak of the spectrum is magnified
for visibility purposes. Notably while the spectrum down to 1 MeV is shown,
the fit only considers energies higher than 1 GeV

In order to estimate the aforementioned uncertainties we have written
publicly-available code, obtainable at GitHub, that can estimate the DM
uncertainties for a wide variety of DM masses and channels. The DM antiproton
spectrum after Galactic propagation is computed by diffusing the DM
annihilation spectrum at each bin with pre-computed diffusion values. Thus by
diffusing all bins the entire propagated spectrum is obtained. We determine
the diffusion values for a specific energy bin with DRAGON 2 [68, 69] by using
a peaked DM annihilation spectrum at the energy of the bin. The settings for
DRAGON 2 are determined by fitting the proton ($p$), antiproton-over-proton
($\bar{p}/p$), and boron-over-carbon ratios ($B/C$) from AMS–02 [70] with an
artificial bee colony [71, 72]444The propagation parameters can be found in
the README of the code. In order to make the code flexible, the pre-computed
diffusion values can be supplied by the user, such that it can be tailored for
specific use. Additionally, the energy range which is considered by the fit
can be supplied, such that only the region of interest is fitted. We stress
here explicitly that our code should not be used in place of an actual CR
propagation code and should only be used in order to obtain the uncertainties
on the aforementioned DM observables. Notably, our fit of the AMS-02 spectra
resulted in a DM mass of $\mathcal{O}(100)~{}\text{GeV}$ and additionally is
only one parameter set. In order to assess the impact of different propagation
parameters on the DM uncertainties, we have performed cross checks with three
different propagation parameter sets: the aforementioned best-fit scenario,
the DRAGON 2 default settings, and one taken from [73]. We found that the
impact of different propagation parameters results in a variation on the
uncertainty of the DM observables of approximately 1-5%. It should be noted
that only a low-energy correction factor on the diffusion coefficient is
implemented, so solar modulation is not fully taken into account for the pre-
computed diffusion values. In the following all uncertainties are fitted by
disregarding spectra energies lower than 1 GeV. While solar modulation does
play a role at 1 GeV, the spectrum uncertainties are fairly uniform in this
region, thus the implementation of solar modulation is not expected to change
the results by a significant amount.
In figure 3 we showcase the diffused spectrum of a 100 GeV DM particle
annihilating into antiprotons via a $W^{+}W^{-}$ pair, and two spectra which
contain the 100 GeV spectrum in their uncertainty bounds, i.e. $\Delta
M_{\chi}$. One of the main constraints determining the DM mass uncertainty is
the relatively small uncertainty band at $x\approx 0.2$, which can most
clearly be seen in figure 2. While the diffusion of the spectra alters the
pre-diffused spectrum, the small uncertainty band remains a clear feature.
This region is especially important for the uncertainty on $\langle\sigma
v\rangle$.

$M$ | Channel |  $\langle\sigma v\rangle$ [%] |  $\Delta M_{\chi}$ [GeV]
---|---|---|---
100 | $b\bar{b}$ | $\phantom{}{}^{+9.2}_{-7.0}$ | $\phantom{}{}^{+12.8}_{-4.5}$
1000 | $b\bar{b}$ | $\phantom{}{}^{+7.9}_{-6.3}$ | $\phantom{}{}^{+143.3}_{-65.5}$
100 | $W^{+}W^{-}$ | $\phantom{}{}^{+6.7}_{-9.4}$ | $\phantom{}{}^{+10.0}_{-6.4}$
1000 | $W^{+}W^{-}$ | $\phantom{}{}^{+7.6}_{-9.3}$ | $\phantom{}{}^{+63.0}_{-56.9}$
Table 2: The uncertainties on the two different DM observables (see text), for
two different annihilation channels for both a 100 GeV and 1000 GeV DM
particle. The uncertainty on $\langle\sigma v\rangle$ is given % since its
absolute value only changes the normalization of the spectrum and is
additionally model specific. All fits again only fit the spectrum at energies
higher than 1 GeV.

In table 2 we show the uncertainties for two fairly commonplace annihilation
channels in DM studies; $b\bar{b}$ and $W^{+}W^{-}$, for two example masses,
100 GeV and 1000 GeV. Notably, the uncertainties for $b\bar{b}$ is fairly
asymmetrical, this is mainly due to the fact that for $b\bar{b}$ the lower QCD
uncertainties are larger than the upper uncertainties. Thus the specific
annihilation channel is highly relevant for the resulting uncertainties on the
DM observables, which can reach up to 14% in certain scenarios.

##  5 Conclusions

In this work, we have shown for the first time that QCD uncertainties on
antiproton spectra may have an important effect on the quality of DM fits.
Therefore, these uncertainties are relevant for DM interpretations of the
AMS–02 excess and should be used in future DM analyses555A. Jueid, J. Kip, R.
Ruiz de Austri, P. Skands, in preparation.. We analysed the effects of QCD
uncertainties from both parton showers and hadronisation on the best-fit point
and found that they can change the mass of DM by $\Delta
M_{\chi}=\phantom{}^{+12.8}_{-4.5}~{}(\phantom{}^{+143.3}_{-65.5})$ GeV for
$M_{\chi}=100~{}(1000)$ GeV in the $b\bar{b}$ annihilation channel. The
effects of the QCD uncertainties on the DM mass in the $W^{+}W^{-}$
annihilation channel are of similar size as in the $b\bar{b}$ channel.
Moreover, we found that the thermally-averaged annihilation cross section is
affected by $6.3\%$–$9.4\%$ depending on the DM mass and the annihilation
channel. The comparatively large size of the QCD uncertainties on antiproton
spectra (relative, e.g., to those for photon spectra [31]) are partly due to
tensions between different measurements of baryon spectra at LEP.
Complementary detailed measurements of baryon spectra at the LHC could
therefore potentially deliver additional constraints that would be
indispensable not only for DM analyses in indirect-detection experiments but
also for extrapolations of fragmentation functions needed for high-energy
cosmic rays. The final states at LHC are of course far more complex than those
at LEP, but relatively clean observables may still be constructed, e.g. by
using one of the recently developed quark/gluon jet discrimination observables
for jets with very high transverse momenta [74, 75, 76, 77]. To facilitate
further full-fledged analyses and help the user to quickly estimate the
effects of the QCD uncertainties on the fit results, we provide tabulated
spectra with uncertainties and the code to determine DM observable
uncertainties at github. The tables can also be found in the latest releases
of DarkSusy 6 [78], MicrOmegas 5 [79] and MadDM [80].

## Acknowledgments

The work of AJ is supported in part by the Institute for Basic Science (IBS)
under the project code, IBS-R018-D1. JK is supported by the NWO Physics Vrij
Programme “The Hidden Universe of Weakly Interacting Particles" with project
number 680.92.18.03 (NWO Vrije Programma), which is (partly) financed by the
Dutch Research Council (NWO). R. RdA acknowledges the Ministerio de Ciencia e
Innovación (PID2020-113644GB-I00). PS is funded by the Australian Research
Council via Discovery Project DP170100708 — “Emergent Phenomena in Quantum
Chromodynamics”. This work was also supported in part by the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie grant agreement No 722105 — MCnetITN3.

## References

  * [1] G. Bertone, D. Hooper, and J. Silk, Particle dark matter: Evidence, candidates and constraints, Phys.Rept. 405 (2005) 279–390, [hep-ph/0404175].
  * [2] Planck Collaboration, P. A. R. Ade et al., Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594 (2016) A13, [arXiv:1502.01589].
  * [3] R. L. Golden, S. Horan, B. G. Mauger, G. D. Badhwar, J. L. Lacy, S. A. Stephens, R. R. Daniel, and J. E. Zipse, EVIDENCE FOR THE EXISTENCE OF COSMIC RAY ANTI-PROTONS, Phys. Rev. Lett. 43 (1979) 1196–1199.
  * [4] A. Buffington, S. M. Schindler, and C. R. Pennypacker, A measurement of the cosmic-ray antiproton flux and a search for antihelium, Astrophys. J. 248 (1981) 1179–1193.
  * [5] J. Silk and M. Srednicki, Cosmic Ray anti-Protons as a Probe of a Photino Dominated Universe, Phys. Rev. Lett. 53 (1984) 624.
  * [6] AMS Collaboration, M. Aguilar et al., Antiproton Flux, Antiproton-to-Proton Flux Ratio, and Properties of Elementary Particle Fluxes in Primary Cosmic Rays Measured with the Alpha Magnetic Spectrometer on the International Space Station, Phys. Rev. Lett. 117 (2016), no. 9 091103.
  * [7] A. Cuoco, M. Kramer, and M. Korsmeier, Novel dark matter constraints from antiprotons in the light of AMS-02, arXiv:1610.03071.
  * [8] M.-Y. Cui, Q. Yuan, Y.-L. S. Tsai, and Y.-Z. Fan, Possible dark matter annihilation signal in the AMS-02 antiproton data, Phys. Rev. Lett. 118 (2017), no. 19 191101, [arXiv:1610.03840].
  * [9] A. Cuoco, J. Heisig, M. Korsmeier, and M. Krämer, Probing dark matter annihilation in the Galaxy with antiprotons and gamma rays, JCAP 10 (2017) 053, [arXiv:1704.08258].
  * [10] A. Reinert and M. W. Winkler, A Precision Search for WIMPs with Charged Cosmic Rays, JCAP 01 (2018) 055, [arXiv:1712.00002].
  * [11] L. Goodenough and D. Hooper, Possible Evidence For Dark Matter Annihilation In The Inner Milky Way From The Fermi Gamma Ray Space Telescope, arXiv:0910.2998.
  * [12] Fermi/LAT Collaboration, Indirect Search for Dark Matter from the center of the Milky Way with the Fermi-Large Area Telescope, arXiv:0912.3828.
  * [13] D. Hooper and L. Goodenough, Dark Matter Annihilation in The Galactic Center As Seen by the Fermi Gamma Ray Space Telescope, Phys. Rev. B697 (2011) 412–428, [arXiv:1010.2752].
  * [14] C. Gordon and O. Macias, Dark Matter and Pulsar Model Constraints from Galactic Center Fermi-LAT Gamma Ray Observations, Phys. Rev. D88 (2013) 083521, [arXiv:1306.5725].
  * [15] D. Hooper and T. Linden, On The Origin Of The Gamma Rays From The Galactic Center, Phys. Rev. D84 (2011) 123005, [arXiv:1110.0006].
  * [16] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden, S. K. N. Portillo, et al., The Characterization of the Gamma-Ray Signal from the Central Milky Way: A Compelling Case for Annihilating Dark Matter, arXiv:1402.6703.
  * [17] F. Calore, I. Cholis, and C. Weniger, Background model systematics for the Fermi GeV excess, JCAP 1503 (2015) 038, [arXiv:1409.0042].
  * [18] A. Achterberg, S. Amoroso, S. Caron, L. Hendriks, R. Ruiz de Austri, and C. Weniger, A description of the Galactic Center excess in the Minimal Supersymmetric Standard Model, JCAP 1508 (2015), no. 08 006, [arXiv:1502.05703].
  * [19] M. van Beekveld, W. Beenakker, S. Caron, and R. Ruiz de Austri, The case for 100 GeV bino dark matter: A dedicated LHC tri-lepton search, JHEP 04 (2016) 154, [arXiv:1602.00590].
  * [20] A. Achterberg, M. van Beekveld, S. Caron, G. A. Gómez-Vargas, L. Hendriks, and R. R. de Austri, Implications of the fermi-lat pass 8 galactic center excess on supersymmetric dark matter, 2018.
  * [21] M. di Mauro, F. Donato, A. Goudelis, and P. D. Serpico, New evaluation of the antiproton production cross section for cosmic ray studies, Phys. Rev. D 90 (2014), no. 8 085017, [arXiv:1408.0288]. [Erratum: Phys.Rev.D 98, 049901 (2018)].
  * [22] R. Kappl and M. W. Winkler, The Cosmic Ray Antiproton Background for AMS-02, JCAP 09 (2014) 051, [arXiv:1408.0299].
  * [23] M. Kachelriess, I. V. Moskalenko, and S. S. Ostapchenko, New calculation of antiproton production by cosmic ray protons and nuclei, Astrophys. J. 803 (2015), no. 2 54, [arXiv:1502.04158].
  * [24] M. W. Winkler, Cosmic Ray Antiprotons at High Energies, JCAP 02 (2017) 048, [arXiv:1701.04866].
  * [25] M. Korsmeier, F. Donato, and M. Di Mauro, Production cross sections of cosmic antiprotons in the light of new data from the NA61 and LHCb experiments, Phys. Rev. D 97 (2018), no. 10 103019, [arXiv:1802.03030].
  * [26] J. Heisig, M. Korsmeier, and M. W. Winkler, Dark matter or correlated errors: Systematics of the AMS-02 antiproton excess, Phys. Rev. Res. 2 (2020), no. 4 043017, [arXiv:2005.04237].
  * [27] F. Calore, M. Cirelli, L. Derome, Y. Genolini, D. Maurin, P. Salati, and P. D. Serpico, AMS-02 $\overline{p}$ ’s and dark matter: Trimmed hints and robust bounds, arXiv:2202.03076.
  * [28] A. Buckley et al., General-purpose event generators for LHC physics, Phys. Rept. 504 (2011) 145–233, [arXiv:1101.2599].
  * [29] P. Skands, S. Carrazza, and J. Rojo, Tuning PYTHIA 8.1: the Monash 2013 Tune, Eur. Phys. J. C74 (2014), no. 8 3024, [arXiv:1404.5630].
  * [30] Particle Data Group Collaboration, P. A. Zyla et al., Review of Particle Physics, PTEP 2020 (2020), no. 8 083C01.
  * [31] S. Amoroso, S. Caron, A. Jueid, R. Ruiz de Austri, and P. Skands, Estimating QCD uncertainties in Monte Carlo event generators for gamma-ray dark matter searches, JCAP 05 (2019) 007, [arXiv:1812.07424].
  * [32] S. Amoroso, S. Caron, A. Jueid, R. Ruiz de Austri, and P. Skands, Particle spectra from dark matter annihilation: physics modelling and QCD uncertainties, PoS TOOLS2020 (2021) 028, [arXiv:2012.08901].
  * [33] A. Jueid, Studying QCD modeling of uncertainties in particle spectra from dark-matter annihilation into jets, in 17th International Conference on Topics in Astroparticle and Underground Physics, 10, 2021.  arXiv:2110.09747.
  * [34] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, An Introduction to PYTHIA 8.2, Comput. Phys. Commun. 191 (2015) 159–177, [arXiv:1410.3012].
  * [35] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, et al., PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection, JCAP 1103 (2011) 051, [arXiv:1012.4515].
  * [36] P. Z. Skands, Tuning Monte Carlo Generators: The Perugia Tunes, Phys. Rev. D82 (2010) 074018, [arXiv:1005.3457].
  * [37] B. Andersson, G. Gustafson, G. Ingelman, and T. Sjöstrand, Parton Fragmentation and String Dynamics, Phys. Rept. 97 (1983) 31–145.
  * [38] T. Sjöstrand, The Lund Monte Carlo for Jet Fragmentation, Comput. Phys. Commun. 27 (1982) 243.
  * [39] T. Sjöstrand, Jet Fragmentation of Nearby Partons, Nucl. Phys. B 248 (1984) 469–502.
  * [40] DELPHI Collaboration, P. Abreu et al., pi+-, K+-, p and anti-p production in Z0 —$>$ q anti-q, Z0 —$>$ b anti-b, Z0 —$>$ u anti-u, d anti-d, s anti-s, Eur. Phys. J. C 5 (1998) 585–620.
  * [41] ALEPH Collaboration, R. Barate et al., Inclusive production of pi0, eta, eta-prime (958), K0(S) and lambda in two jet and three jet events from hadronic Z decays, Eur. Phys. J. C 16 (2000) 613.
  * [42] C. Bierlich et al., A comprehensive guide to the physics and usage of PYTHIA 8.3, arXiv:2203.11601.
  * [43] T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert, and J. Winter, Event generation with SHERPA 1.1, JHEP 02 (2009) 007, [arXiv:0811.4622].
  * [44] S. Schumann and F. Krauss, A Parton shower algorithm based on Catani-Seymour dipole factorisation, JHEP 03 (2008) 038, [arXiv:0709.1027].
  * [45] J.-C. Winter, F. Krauss, and G. Soff, A Modified cluster hadronization model, Eur. Phys. J. C36 (2004) 381–395, [hep-ph/0311085].
  * [46] T. Sjöstrand, S. Mrenna, and P. Z. Skands, PYTHIA 6.4 Physics and Manual, JHEP 0605 (2006) 026, [hep-ph/0603175].
  * [47] J. Bellm et al., Herwig 7.0/Herwig++ 3.0 release note, Eur. Phys. J. C76 (2016), no. 4 196, [arXiv:1512.01178].
  * [48] S. Gieseke, P. Stephens, and B. Webber, New formalism for QCD parton showers, JHEP 12 (2003) 045, [hep-ph/0310083].
  * [49] S. Platzer and S. Gieseke, Coherent Parton Showers with Local Recoils, JHEP 01 (2011) 024, [arXiv:0909.5593].
  * [50] S. Platzer and S. Gieseke, Dipole Showers and Automated NLO Matching in Herwig++, Eur. Phys. J. C72 (2012) 2187, [arXiv:1109.6256].
  * [51] C. Bierlich et al., Robust Independent Validation of Experiment and Theory: Rivet version 3, SciPost Phys. 8 (2020) 026, [arXiv:1912.05451].
  * [52] ALEPH Collaboration, D. Buskulic et al., Inclusive pi+-, K+- and (p, anti-p) differential cross-sections at the Z resonance, Z. Phys. C66 (1995) 355–366.
  * [53] ALEPH Collaboration, R. Barate et al., Studies of quantum chromodynamics with the ALEPH detector, Phys. Rept. 294 (1998) 1–165.
  * [54] ALEPH Collaboration, A. Heister et al., Studies of QCD at e+ e- centre-of-mass energies between 91-GeV and 209-GeV, Eur. Phys. J. C35 (2004) 457–486.
  * [55] DELPHI Collaboration, P. Abreu et al., Charged particle multiplicity distributions for fixed number of jets in Z0 hadronic decays, Z. Phys. C 56 (1992) 63–76.
  * [56] DELPHI Collaboration, P. Abreu et al., Production of Lambda and Lambda anti-Lambda correlations in the hadronic decays of the Z0, Phys. Lett. B 318 (1993) 249–262.
  * [57] DELPHI Collaboration, W. Adam et al., Measurement of inclusive pi0 production in hadronic Z0 decays, Z. Phys. C69 (1996) 561–574.
  * [58] DELPHI Collaboration, P. Abreu et al., Inclusive measurements of the K+- and p / anti-p production in hadronic Z0 decays, Nucl. Phys. B 444 (1995) 3–26.
  * [59] L3 Collaboration, P. Achard et al., Studies of hadronic event structure in $e^{+}e^{-}$ annihilation from 30-GeV to 209-GeV with the L3 detector, Phys. Rept. 399 (2004) 71–174, [hep-ex/0406049].
  * [60] OPAL Collaboration, P. D. Acton et al., A Study of charged particle multiplicities in hadronic decays of the Z0, Z. Phys. C53 (1992) 539–554.
  * [61] OPAL Collaboration, R. Akers et al., Measurement of the production rates of charged hadrons in e+ e- annihilation at the Z0, Z. Phys. C63 (1994) 181–196.
  * [62] OPAL Collaboration, G. Alexander et al., Strange baryon production in hadronic Z0 decays, Z. Phys. C73 (1997) 569–586.
  * [63] OPAL Collaboration, K. Ackerstaff et al., Measurements of flavor dependent fragmentation functions in Z0 –> q anti-q events, Eur. Phys. J. C7 (1999) 369–381, [hep-ex/9807004].
  * [64] OPAL Collaboration, G. Abbiendi et al., Measurement of event shape distributions and moments in e+ e- —> hadrons at 91-GeV - 209-GeV and a determination of alpha(s), Eur. Phys. J. C40 (2005) 287–316, [hep-ex/0503051].
  * [65] A. Buckley, H. Hoeth, H. Lacker, H. Schulz, and J. E. von Seggern, Systematic event generator tuning for the LHC, Eur. Phys. J. C65 (2010) 331–357, [arXiv:0907.2973].
  * [66] S. Mrenna and P. Skands, Automated Parton-Shower Variations in Pythia 8, Phys. Rev. D94 (2016), no. 7 074005, [arXiv:1605.08352].
  * [67] J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai, and W. K. Tung, Uncertainties of predictions from parton distribution functions. 2. The Hessian method, Phys. Rev. D 65 (2001) 014013, [hep-ph/0101032].
  * [68] C. Evoli, D. Gaggero, A. Vittino, G. D. Bernardo, M. D. Mauro, A. Ligorini, P. Ullio, and D. Grasso, Cosmic-ray propagation with dragon2: I. numerical solver and astrophysical ingredients, Journal of Cosmology and Astroparticle Physics 2017 (Feb, 2017) 015–015.
  * [69] C. Evoli, D. Gaggero, A. Vittino, M. Di Mauro, D. Grasso, and M. N. Mazziotta, Cosmic-ray propagation with dragon2: Ii. nuclear interactions with the interstellar gas, Journal of Cosmology and Astroparticle Physics 2018 (Jul, 2018) 006–006.
  * [70] AMS Collaboration, M. Aguilar et al., The alpha magnetic spectrometer (ams) on the international space station: Part ii — results from the first seven years, Physics Reports 894 (2021) 1–116.
  * [71] D. Karaboga, An idea based on honey bee swarm for numerical optimization, technical report - tr06, Technical Report, Erciyes University (01, 2005).
  * [72] B. Akay and D. Karaboga, A modified artificial bee colony algorithm for real-parameter optimization, Information Sciences 192 (2012) 120–142. Swarm Intelligence and Its Applications.
  * [73] G. D. Bernardo, C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, Unified interpretation of cosmic ray nuclei and antiproton recent measurements, Astroparticle Physics 34 (dec, 2010) 274–283.
  * [74] J. Gallicchio and M. D. Schwartz, Quark and Gluon Tagging at the LHC, Phys. Rev. Lett. 107 (2011) 172001, [arXiv:1106.3076].
  * [75] P. T. Komiske, E. M. Metodiev, and M. D. Schwartz, Deep learning in color: towards automated quark/gluon jet discrimination, JHEP 01 (2017) 110, [arXiv:1612.01551].
  * [76] E. M. Metodiev and J. Thaler, Jet Topics: Disentangling Quarks and Gluons at Colliders, Phys. Rev. Lett. 120 (2018), no. 24 241602, [arXiv:1802.00008].
  * [77] P. Gras, S. Höche, D. Kar, A. Larkoski, L. Lönnblad, S. Plätzer, A. Siódmok, P. Skands, G. Soyez, and J. Thaler, Systematics of quark/gluon tagging, JHEP 07 (2017) 091, [arXiv:1704.03878].
  * [78] 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, [arXiv:1802.03399].
  * [79] G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov, and B. Zaldivar, micrOMEGAs5.0 : Freeze-in, Comput. Phys. Commun. 231 (2018) 173–186, [arXiv:1801.03509].
  * [80] F. Ambrogi, C. Arina, M. Backovic, J. Heisig, F. Maltoni, L. Mantani, O. Mattelaer, and G. Mohlabeng, MadDM v.3.0: a Comprehensive Tool for Dark Matter Studies, Phys. Dark Univ. 24 (2019) 100249, [arXiv:1804.00044].