# Fast and accurate AMS-02 antiproton likelihoods for global dark matter fits

Sowmiya Balan,11footnotetext: Corresponding authors Felix Kahlhoefer  Michael
Korsmeier  Silvia Manconi  and Kathrin Nippel

###### Abstract

The antiproton flux measurements from AMS-02 offer valuable information about
the nature of dark matter, but their interpretation is complicated by large
uncertainties in the modeling of cosmic ray propagation. In this work we
present a novel framework to efficiently marginalise over propagation
uncertainties in order to obtain robust AMS-02 likelihoods for arbitrary dark
matter models. The three central ingredients of this framework are: the neural
emulator DarkRayNet, which provides highly flexible predictions of the
antiproton flux; the likelihood calculator pbarlike, which performs the
marginalisation, taking into account the effects of solar modulation and
correlations in AMS-02 data; and the global fitting framework GAMBIT, which
allows for the combination of the resulting likelihood with a wide range of
dark matter observables. We illustrate our approach by providing updated
constraints on the annihilation cross section of WIMP dark matter into bottom
quarks and by performing a state-of-the-art global fit of the scalar singlet
dark matter model, including also recent results from direct detection and the
LHC.

##  1 Introduction

The search for dark matter (DM) is a global effort, both in the geographical
sense, i.e. in terms of the sheer number of experiments and groups involved,
and in the methodological one, i.e. in terms of the variety of relevant
constraints that need to be considered. The scope of this endeavour calls for
computational tools that provide fast and accurate theory predictions of
relevant observables, efficient routines for likelihood calculations, and
methods for their statistical interpretation. Of particular importance are
tools that can perform DM relic density calculations and combine the results
with a variety of constraints from direct and indirect detection experiments,
such as micrOmegas [1], DarkSUSY [2], MadDM [3], SuperIso Relic [4] and the
DarkBit [5] module of the GAMBIT [6] global fitting framework.

An excellent example for the need for tools are satellite measurements of
cosmic rays (CRs), which are sensitive to charged antiparticles produced in DM
annihilations [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. In
particular the most recent measurements of the antiproton flux by AMS-02 [22]
have become so precise that we can expect strong constraints on, or potential
evidence for, many different DM models, provided the background from secondary
antiprotons produced in astrophysical systems can be modelled with sufficient
accuracy [23, 24].

Unfortunately, a detailed modeling of CR propagation in the Galaxy needed for
this purpose is computationally expensive. The state-of-the-art modeling for
the CR journey in our Galaxy solves such transport problem numerically. Many
astrophysical ingredients enter the computation, such as CR injection and
transport properties. This requires a large number of parameters to be
included in the analysis, some of them poorly constrained. While it is
possible to consider these uncertainties in the calculation of AMS-02
antiproton constraints [25], including these constraints in global fits of DM
models together with direct detection and collider constraints is very
challenging without simplifying assumptions (see Ref. [26]).

In the present work, we address this gap in the toolbox of DM phenomenologists
by providing three separate codes that together form a full analysis chain for
AMS-02 antiproton data in the context of analysing DM models:

  * •

DarkRayNet.v2 is an updated version of the tool first presented in Ref. [27].
It provides fully trained neural networks that predict the primary (i.e. DM-
produced) and secondary antiproton fluxes as a function of the DM properties
and a wide range of propagation parameters. Compared to the first release,
DarkRayNet.v2 contains more flexible propagation models and additional
training data.

  * •

pbarlike takes predictions for the antiproton flux in order to calculate the
likelihood of the AMS-02 measurements, including a state-of-the-art treatment
of correlations. Moreover, it provides the marginalised likelihood when the
propagation parameters and the effect of solar modulation are treated as
nuisance parameters.

  * •

GAMBIT 2.4 is the most recent version of the GAMBIT [6] global fitting
framework. It includes the necessary interface to calculate these likelihoods
in large-scale parameter scans in order to perform global fits of many
different DM models.222For additional features of the new GAMBIT release, we
refer to Refs. [28, 29].

We illustrate the usefulness of these codes in two applications. First, we
provide updated model-independent constraints on the annihilation cross
section of DM into bottom quarks. In agreement with recent studies [30, 31,
32, 24], we show that the previously observed excess around DM masses of 100
GeV [33, 34, 35] is reduced to negligible significance when correlations and
uncertainties are properly included. Moreover, the nature of the excess
depends on the specific propagation model being used, highlighting the need
for a sufficiently flexible propagation model.

As a second application, we provide the most up-to-date global analysis of the
scalar singlet DM model [36, 37, 38], which has been the subject of a number
of global analyses in recent years [39, 40, 41]. In addition to the AMS-02
antiproton data, we also include the most recent measurements of the Higgs
invisible width from the LHC [42, 43] and new constraints from the direct
detection experiments LZ [44] and PandaX-4T [45]. We find that with these
constraints, the high-mass region ($m_{s}>m_{h}$) is largely excluded and only
the resonance region ($m_{s}\approx m_{h}/2$) remains of interest.

The remainder of this work is structured as follows. In section 2 we present
the framework for predicting the CR flux and introduce the specific injection
and propagation models that we consider. The AMS-02 likelihood implementation
and the formalism for marginalisation of nuisance parameters, and results of
model-independent parameter scans are discussed in section 3. Finally, we
introduce the model of scalar singlet DM in section 4 and provide our results
from a global analysis of this model. Additional details on the numerical
tools are provided in the appendix.

##  2 Injection and propagation of cosmic rays

The transport of CRs in the Galaxy can be described by the following equation
[46]:

$\displaystyle\frac{\partial\psi_{i}(\bm{x},p,t)}{\partial t}=q_{i}(\bm{x},p)$
$\displaystyle+$
$\displaystyle\bm{\nabla}\cdot\left(D_{xx}\bm{\nabla}\psi_{i}-\bm{V}\psi_{i}\right)$
$\displaystyle+$ $\displaystyle\frac{\partial}{\partial
p}p^{2}D_{pp}\frac{\partial}{\partial
p}\frac{1}{p^{2}}\psi_{i}-\frac{\partial}{\partial
p}\left(\frac{\mathrm{d}p}{\mathrm{d}t}\psi_{i}-\frac{p}{3}(\bm{\nabla\cdot
V})\psi_{i}\right)-\frac{1}{\tau_{f,i}}\psi_{i}-\frac{1}{\tau_{r,i}}\psi_{i}\;.$

where $\psi_{i}(\bm{x},p,t)$ is the CR number density per volume and absolute
momentum $p$, for each CR species $i$, evaluated at position $\bm{x}$ in
Galactocentric coordinates. The various terms in eq. (2) are carefully
explained in the remainder of this section.

We employ Galprop version 56 [47] and Galtoollibs 855333
https://galprop.stanford.edu to solve numerically the transport equations for
each species required in the following, such as for the training of
DarkRayNet, see section 3.2. A number of custom modifications to the public
Galprop code have been implemented as described in Ref. [31]. The Galprop code
permits state-of-the-art modeling of CR transport in the Galaxy and is widely
used to compute CR spectra at Earth’s position, as well as CR non-thermal
emission at different wavelengths. Other, fully numerical codes such as DRAGON
[48, 49], PICARD [50], as well as semi-analytical frameworks such as USINE
[51] are also used in the recent literature. CRs are assumed to be in steady
state (left-hand side in eq. (2) is zero). The equations are solved on a
3-dimensional grid, where two dimensions describe the CR’s spatial
distribution, and the remaining one their kinetic energy. The spatial grid is
defined in cylindrical coordinates, being $r$ the radial distance from the
Galactic center and $z$ the distance perpendicular to the plane. The kinetic
energy grid is logarithmically spaced, and the ratio between successive grid
points is set to 1.3. The values of the step sizes for the spatial grid follow
our previous work [27].

###  2.1 Primary and secondary cosmic rays

The term $q_{i}(\bm{x},p)$ in eq. (2) incorporates the injection of CRs at
each Galactic position $\bm{x}$, with a given energy dependence. Depending on
the considered CR species $i$, we include primary and/or secondary
contributions, which are briefly illustrated in the next subsections. We
define as primary CR sources both the standard astrophysical sources such as
supernova remnants (section 2.1.1), along with more exotic CR injections such
as the production of (anti)particle CRs in the annihilation or decay of DM
particles in the Galactic DM halo (section 2.1.2). Secondary CRs are instead
produced by the interaction of the primary CRs with the gas in the
interstellar medium (ISM) through fragmentation or decay processes (section
2.1.3).

####  2.1.1 Primary p, He

CRs are dominated by protons, which account for about 90% of the observed
flux, followed by 10% of helium (He) and even smaller percentages of heavier
nuclei, electrons and positrons [22]. In the standard paradigm, Galactic CR
such as protons, electrons, He and other nuclei are accelerated and then
released in primary sources such as supernova remnants [52]. Diffusive shock
acceleration is considered to be the main mechanism for promoting these
particles up to multi-TeV energies. Particle-in-cell simulations have
demonstrated that this mechanism produces power-law spectra for the
accelerated particles [53], which can be modified during their escape from the
remnant [54, 55]. In addition, a break in the injection spectrum at low
energies is required to explain data at about few GeV, which could be
connected to self-confinement of particles in supernova remnants [56].

We model the source term for primary $p$ and He by factorizing the spatial and
energetic dependence. The first is characterized by a smooth distribution of
supernova remnants in Galactocentric cylindrical coordinates ($r,z$, with
$r=\sqrt{x^{2}+y^{2}}$) of the form
$f(r,z)=\left(\frac{r}{r_{\odot}}\right)^{\alpha_{1}}\exp(-\alpha_{2}(r-r_{\odot})/r_{\odot}))\exp({-\mid
z\mid/z_{0}})$, where the parameters
$\alpha_{1}=1.09,\alpha_{2}=3.87,z_{0}=0.2$ are fixed according to default
prescriptions of Galprop, and the Earth’s distance from the Galactic Center is
set to $r_{\odot}=8.5$ kpc. The modeling of the spatial distribution of
sources of CR nuclei has been demonstrated to have a very minor impact on the
resulting fluxes at Earth [57]. The energetic dependence is modeled as a power
law or as a smoothly broken power law, depending on the chosen setup (see
section 2.3). This approximation is widely used [58, 59, 60, 61] and permits
to well describe the data in the rigidity range we consider:

$\displaystyle q_{R}(R)$ $\displaystyle=$
$\displaystyle\left(\frac{R}{R_{0}}\right)^{-\gamma_{1}}\left(\frac{R_{0}^{1/s}+R^{1/s}}{2\,R_{0}^{1/s}}\right)^{-s(\gamma_{2}-\gamma_{1})},$
(2.2)

which is given as a function of rigidity $R$ (momentum $p$ divided by the
absolute charge value), a more natural quantity than momentum when dealing
with charged particle acceleration. As illustrated by figure 1, when we
consider a simple power law, only one spectral index $\gamma_{2}$ is relevant.
When a break in the injection spectrum at the rigidity value $R_{0}$ is
considered, this is regulated by a smoothing parameter $s$, and requires two
spectral indices below ($\gamma_{1}$) and above ($\gamma_{2}$) the break. We
assume a universal injection spectrum for all primary nuclei, except for
protons for which different spectral indices $\gamma_{1,p},\gamma_{2,p}$ are
introduced, see discussion in [59] and references therein.

Figure 1:  Sketch of the propagation setups INJ.BRK (upper row) and DIFF.BRK
(lower row). The different panels depict the effect of the model parameter on
the spectrum $\phi$ or individual terms of the transport equation 2.

####  2.1.2 Primary antiprotons from dark matter

DM interactions in the diffusion halo of our Galaxy could inject further,
exotic primary CRs through annihilation or decay. These processes produce an
equal amount of matter and antimatter CRs, such as protons and antiprotons.
However, in the standard paradigm the antimatter CRs such as antiprotons,
antideuterons are produced only as secondaries (see next section), and their
production is suppressed with respect to primaries by 4-5 orders of magnitude.
Thus antimatter CR fluxes are sensitive targets to search for exotic spectral
components [62]. Further astrophysical mechanism have been explored in the
literature to produce primary antimatter, such as the acceleration in old SNRs
[63, 64, 60, 65].

In this work we focus on annihilation processes of the type DM DM $\rightarrow
f+\bar{f}$, where $f$ indicates a specific standard model particle final
state, for example the quark-antiquark mode $b\bar{b}$. The source term for
primary CR antiprotons from DM annihilation is factorized into a spatial term
and a term dependent on the antiproton kinetic energy $E_{\mathrm{kin}}$:

$\displaystyle
q_{\bar{p}}^{(\mathrm{DM})}(\bm{x},E_{\mathrm{kin}})=\frac{1}{2}\left(\frac{\rho(\bm{x})}{m_{\mathrm{DM}}}\right)^{2}\sum_{f}\left\langle\sigma
v\right\rangle_{f}\frac{\mathrm{d}N^{f}_{\bar{p}}}{\mathrm{d}E_{\mathrm{kin}}}\;,$
(2.3)

where the factor $1/2$ is for Majorana fermion DM, and $m_{\mathrm{DM}}$ is
the DM particle mass. The $\rho(\bm{x})$ term in eq. (2.3) indicates the DM
spatial density profile. As done in previous works, we assume that the DM halo
of our Galaxy is described by a simple NFW radial profile [66] with scale
radius $r_{h}=20$ kpc. We assign the characteristic halo density $\rho_{h}$ in
order to reproduce a local DM density at solar position $r_{\odot}$ of $0.43$
GeV/cm3 [67]. We note that the DM density is consistently rescaled among all
the considered observables when computing the constraints on DM models in
section 4. The sum over individual final states $f$ in eq. (2.3) includes the
thermally averaged annihilation cross section $\left\langle\sigma
v\right\rangle_{f}$ for each individual final state $f$, and
${\mathrm{d}N^{f}_{\bar{p}}}/{\mathrm{d}E_{\mathrm{kin}}}$ which is the
antiproton energy spectrum for a single DM annihilation. Here we follow [27]
and fix the annihilation cross section independently of $f$, and assign
branching fractions into different final states. We consider both a single
annihilation final state ($b\bar{b}$) as well as the branching fraction
structure of the scalar singlet DM model (see section 4 and [27]). For both
cases, we take the antiproton energy spectrum for each final state $f$ from
the widely used tabulated results of [68]. These results include electroweak
corrections [69]. For the annihilation channel into a pair of W and Z bosons,
we have extended these tables to include the contribution from the off-shell
production of one W or Z boson following Ref. [70].

####  2.1.3 Secondaries

The propagation of primary CRs in the interstellar medium (ISM) leads to the
production of secondary CRs by fragmentation reactions. We recall that the ISM
of our Galaxy is mainly composed of gas, with a small fraction (0.5-1%) of
dust. The interstellar gas is dominantly hydrogen (about $90$%), and a small
fraction of helium (about $10$%). When primary CRs such as $p$, He interact
with the hydrogen and helium in the ISM, secondary CR particles are produced.
Since some CR species, such as Boron are thought to be exclusively produced
through these processes, secondary-over-primary ratios of CR fluxes are
powerful probes of the transport properties of our Galaxy [71]. In the
standard picture for CRs, also antiprotons are exclusively produced as
secondaries [72].

In order to compute the source of secondary CRs we need to convolute the flux
of primary CRs ($\phi_{i}(E_{\mathrm{kin},i})$) to the ISM density
($n_{\mathrm{ISM}}$) with the energy-differential production cross section
($\frac{\mathrm{d}\sigma_{ij\rightarrow\mathrm{sec}}}{\mathrm{d}E_{\mathrm{kin},\mathrm{sec}}}$):

$\displaystyle q_{\mathrm{sec}}({\bm{x}},E_{\mathrm{kin},\mathrm{sec}})$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\sum\limits_{j\in\\{\mathrm{H},\mathrm{He}\\}}\\!\\!\\!\\!4\pi\,n_{\mathrm{ISM},j}({\bm{x}})\sum\limits_{i}\int\mathrm{d}E_{\mathrm{kin},i}\,\phi_{i}(E_{\mathrm{kin},i})\,\frac{\mathrm{d}\sigma_{ij\rightarrow\mathrm{sec}}}{\mathrm{d}E_{\mathrm{kin},\mathrm{sec}}}(E_{\mathrm{kin},i},E_{\mathrm{kin},\mathrm{sec}})\,.$
(2.4)

The source term in eq. 2.4 is evaluated for each CR species $q_{\mathrm{sec}}$
by considering the relevant CR primaries and the corresponding cross sections.
The ISM density as a function of the position in the Galaxy follows the
default Galprop model [47]. As for the production cross section, we make
different choices depending on the CR species. The secondary $p$ and He cross
sections are modeled following the default Galprop implementation. As for the
secondary antiprotons, we instead implement in Galprop the updated analytic
parametrisation of the Lorentz invariant cross section as obtained in Ref.
[73] as detailed in Ref. [31]. Being tuned to all available cross section data
recorded by colliders at low energies, this parametrisation implements a more
reliable treatment of the production cross section for antiproton energies
below about $10$ GeV. We note that for the fragmentation of primary CRs we
always assume
$\mathrm{d}\sigma_{ij}/\mathrm{d}E_{sec}=\sigma_{ij}\,\delta(E_{kin,i}-E_{kin,sec})$.

Secondary CRs such as antiprotons may further scatter inelastically with the
ISM and consequently lose energy. We consider this contribution, which is
suppressed with respect to the secondaries, and is usually referred to as
tertiary CRs [74].

###  2.2 Propagation

As outlined above, the source term of CR nuclei significantly differs between
primary and secondary CRs as well as between DM and astrophysical sources. In
contrast, the mechanism and, therefore, the description of CR propagation are
the same for all species. The dominant process for CR nuclei, especially at
high energies, is the scattering on the turbulent magnetic fields in our
Galaxy. Effectively, this leads to diffusion of CRs in a halo extending a few
kpc above and below the Galactic plane. From linear perturbation theory, the
diffusion coefficient is expected to be antiproportional to the spectrum of
magnetic wave turbulence [75]. The spectrum of magnetic turbulence depends on
the exact turbulence model; the typical benchmarks are Kolmogorov or
Kraichnan, both leading to power-law dependence of the diffusion coefficient
as a function of rigidity, but with different spectral indices of
$\delta=0.33$ and 0.5 [76], respectively. We use a phenomenological approach
and model the diffusion coefficient as a double-broken power law in rigidity

$\displaystyle D_{xx}\sim\beta
R^{\delta_{l}}\cdot\left(1+\left(\frac{R}{R_{D,0}}\right)^{1/s_{D,0}}\right)^{s_{D,0}\,(\delta-\delta_{l})}\cdot\left(1+\left(\frac{R}{R_{D,1}}\right)^{1/s_{D,1}}\right)^{s_{D,1}\,(\delta_{h}-\delta)}\,,$
(2.5)

which is normalized to $D_{xx}(R=4\;\mathrm{GV})=D_{0}$. The spectral indices
below, between, and above the two breaks at $R_{D,0}$ and $R_{D,1}$ are
labeled $\delta_{l}$, $\delta$, and $\delta_{h}$. At the positions of the
breaks, the power law is smoothed by the parameters $s_{D,0}$ and $s_{D,1}$,
respectively. Indeed, a smooth transition is expected in the rigidity
dependence of the diffusion coefficient, when the turbulence changes between
two regimes [77, 78]. The second break at about 300 GV is directly observed in
the AMS-02 data [22]. The fact that the break is more pronounced in
secondaries than in primaries [79] clearly points to a change in the diffusion
coefficient rather than in the injection spectrum [80]. A natural explanation
for this scenario is for example provided by self-generated turbulence, as
explored in Refs. [78, 77]. In contrast, a break at lower rigidities ($\sim
10$ GV) is more speculative because also other processes like convection and
reacceleration influence the CR spectra [71, 23]. Therefore, we will explore
two different models (as in Ref. [81]), one with the low-energy break and one
without. The models will be detailed further below.

Convective winds with the velocity $\bm{V}_{\rm c}$ can transport the CRs away
from the Galactic plane and induce adiabatic energy losses. We employ a
constant convection velocity that is perpendicular to the Galactic plane,
$\bm{V}_{\rm c}=v_{0,\rm c}\,{\rm sign}(z)\,\bm{e}_{z}$. Increasing the
convection velocity decreases the CR flux below $\sim 10$ GV. In contrast,
reacceleration is due to scattering of CRs on Alfvén magnetic waves and has
the opposite effect on the CR spectra. In a head-on scattering with the Alfvén
waves a CR gains energy while it loses energy in back-on collisions. However,
head-on collisions are more likely such that statistically the CRs gain
energy. In eq. (2), reacceleration is modeled as diffusion in momentum space
with the coefficient $D_{pp}$. Larger magnetic turbulence makes reacceleration
more efficient, such that $D_{pp}\sim 1/D_{xx}$ [82]. The amount of
reacceleration further depends on the Alfvén velocity, $D_{pp}\sim v_{\rm A}$.

The term $\partial/\partial p(\mathrm{d}p/\mathrm{d}t\psi_{i})$ represents
continuous energy losses from ionization and Coulomb collisions [83], while
the last two terms of eq. (2) describe catastrophic losses by fragmentation
and decay with the characteristic time scales $\tau_{f,i}$ and $\tau_{r,i}$,
respectively.

Finally, when the CRs enter the heliosphere, they are deflected and
decelerated by the solar winds. The effect is known as solar modulation and
varies in a 22-year cycle. The effect is most prominent at low energies. In
principle, solar modulation can be described by a diffusion equation, similar
to equation (2), but with the geometry, magnetic field and turbulence adjusted
to the heliosphere. There are semi-analytical [84] or fully numerical codes
[85, 86] solving this equation numerically. We use a simplified approach
instead, and treat solar modulation in the force-field approximation [87].
This approximates well enough the effect of solar winds for nuclei observed
above about few GeV.

###  2.3 Choice of models

We explore two distinct frameworks for CR propagation labelled INJ.BRK and
DIFF.BRK. The assumptions and free parameters of each model are detailed
below.

The INJ.BRK (injection break) model has been explored extensively in
literature [83, 88, 89, 90, 91, 92, 93]. It also matches the model that we
studied in the previous work on antiproton constraints [27]. The assumptions
and free parameters are sketched in figure 1 (upper row). We employ a broken
power law for the injection spectra of the primary (astrophysical) CRs with
slopes $\gamma_{1}$ and $\gamma_{2}$ below and above the break position at
$R_{0}$ with a smoothing $s$. We use different slopes for $p$ and He in the
fit to account for the observed difference of the slopes in the CR fluxes of
the two primaries. Understanding this difference is subject to current
research. Possible explanations are, for example, different source populations
[94, 95, 96] or a $Z/A$-dependence of efficiency of Fermi-shock acceleration
[97, 98]. The diffusion coefficient is modeled as a single broken power law
with the break $R_{D,1}$ around 300 GV (_i.e._ $\delta_{l}=\delta$ in eq.
(2.5)). Furthermore, we allow for reacceleration and convection444In the
recent update of Galprop [99] from version 56 to 57 the implementation of
convection was changed, correcting a false definition of the Crank–Nicolson
coefficients near the Galactic plane. However, this has negligible impact on
our results because our CR fits prefer very small (compatible with 0)
convection velocities. with $v_{\rm A}$ and $v_{0,c}$ as the two free
parameters, respectively. Finally, solar modulation is treated in the force-
field approximation. We allow for a slightly different solar modulation
potential for antiprotons because solar modulation is charge-sign dependent
[100].

The DIFF.BRK (diffusion break) model employs a single power law for the
injection spectrum of the primary CRs. This type of model has also been
studied in literature (see e.g. [83]), but only more recently it is tested
against AMS-02 data [101, 102, 23]. In contrast to the INJ.BRK model,
reacceleration is replaced by an additional break, $R_{D,0}$, in the diffusion
coefficient, see figure 1. It is conceivable that the interaction of CRs and
magnetic turbulence causes a damping of turbulence at low energies that
effectively leads to an increase of the diffusion coefficient [103].

For both setups, we fix the half-height of the diffusion halo to
$4\,\mathrm{kpc}$, which is roughly the lower bound compatible with the
beryllium data from AMS-02 [104, 105, 93, 23, 106]. This is different with
respect to the propagation setup used in our previous work [27], in which
$z_{h}$ was varied as free parameter. However, a precise determination of the
halo size is currently prevented by large systematic uncertainty in the
secondary fragmentation cross section [107, 104, 23]. Thus, at the moment,
larger values for the half-height are equally viable because of the well-known
degeneracy of $z_{h}$ with the normalization of the diffusion coefficient. We
chose a small value for the half-height because this corresponds to the most
conservative DM limit.

If we had performed the whole analysis with a different value of $z_{h}$, to a
first approximation, we would have inferred a different value of $D_{0}$. The
other propagation parameters would only change marginally, and also the
astrophysical fluxes of primary and secondary CRs would not be affected by
this. However, the DM flux would increase as a function of $z_{h}$. The
difference between astrophysical CRs and CRs from DM annihilation is the
spatial extent of the source term. The astrophysical CRs are produced in a
thin layer of the Galactic plane, while DM annihilates and produces
antiprotons in the entire diffusion halo [8]. Empirically, we found the
following enhancement factor for the flux of DM antiprotons, which is valid
for $z_{h}$ between 1 and 10 kpc:

$\displaystyle f_{\bar{p},{\rm DM}}(z_{h})$ $\displaystyle=$ $\displaystyle
1.0+3.7\times 10^{-1}\left(\frac{z_{h}-4\,{\rm kpc}}{\rm kpc}\right)+5.0\times
10^{-3}\left(\frac{z_{h}-4\,{\rm kpc}}{\rm kpc}\right)^{2}$
$\displaystyle-8.1\times 10^{-3}\left(\frac{z_{h}-4\,{\rm kpc}}{\rm
kpc}\right)^{3}+8.4\times 10^{-4}\left(\frac{z_{h}-4\,{\rm kpc}}{\rm
kpc}\right)^{4}\,.$

We note that this factor is normalized to one at our benchmark of $z_{h}=4$
kpc.

##  3 AMS-02 antiproton likelihood

In this work, the antiproton analysis pipeline implemented using pbarlike
obtains antiproton flux predictions from DarkRayNet to calculate the
marginalised AMS-02 antiproton likelihood while considering correlated errors.
Correlations for AMS-02 data have not been published and hence have to be
modeled. The prescription used for modeling of data correlations and relevant
theoretical uncertainties are discussed in section 3.1.
For the choice of propagation models DIFF.BRK and INJ.BRK, the CR propagation
parameter space compatible with recent CR data was identified using MultiNest
scans. The posterior samples from these scans were then used for training the
neural networks in DarkRayNet. They are also used to perform marginalization
over propagation and solar modulation parameters. Details about the MultiNest
scans and marginalisation are discussed in section 3.2 and section 3.3
respectively.
In section 3.4, we illustrate the effect of the new AMS-02 antiproton
likelihood on constraints from the benchmark case of DM annihilation into
bottom quarks.

###  3.1 AMS-02 data and correlations

The likelihood $\mathcal{L}_{\bar{p}}$ for the AMS-02 measurements of the
antiproton flux can be written as
$-2\log\mathcal{L}_{\bar{p}}=\chi^{2}_{\bar{p}}$, with

$\chi^{2}_{\bar{p}}(\bm{x}_{\text{DM}},\bm{\theta}_{\text{prop}})=\sum_{i,j}\left(\phi_{\bar{p},i}^{(\text{AMS})}-\phi_{\bar{p},i}(\bm{x}_{\text{DM}},\bm{\theta}_{\text{prop}})\right)V^{-1}_{ij}\left(\phi_{\bar{p},j}^{(\text{AMS})}-\phi_{\bar{p},j}(\bm{x}_{\text{DM}},\bm{\theta}_{\text{prop}})\right)\;,$
(3.1)

where $\mathbf{\phi}_{\bar{p}}^{(\text{AMS})}$ denotes the AMS-02 measurements
and $V$ is the covariance matrix. We estimate the covariance matrix following
the approach of Refs. [71, 32], i.e. we write

$V=\text{diag}(\sigma^{2}_{\text{stat},i})+\sum_{\alpha}V^{\alpha}+V_{\text{\rm
xs}}\;.$ (3.2)

The first two terms represent the uncertainty of the flux measurement by
AMS-02, where $\sigma_{\text{stat},i}$ denotes the (uncorrelated) statistical
uncertainty in bin $i$ and $V^{\alpha}$ denote additional systematic
uncertainties.555AMS-02 provides only the diagonal entries of the sum of first
two terms in eq. (3.2), i.e, it gives the uncorrelated statistical and
systematic uncertainties in flux measurements. A crucial contribution to the
correlated uncertainty comes from the antiproton absorption cross section in
the AMS-02 detector, as carefully modeled in Ref. [32]. All other
contributions can be approximated by

$\left(V^{\alpha}\right)_{ij}=\Delta_{i}^{\alpha}\Delta_{j}^{\alpha}\exp\left(-\frac{(\log
R_{i}/R_{j})^{2}}{2(l^{\alpha})^{2}}\right)\;$ (3.3)

with $R_{i}$ denoting the rigidity of bin $i$, $\Delta_{i}^{\alpha}$ denoting
the systematic error and $l^{\alpha}$ denoting the corresponding correlation
length. Following Ref. [32] we consider errors in the effective acceptance,
errors in the rigidity scale, unfolding errors, geomagnetic cut-off errors,
template shape errors and selection errors.

In contrast to the first two terms, the third term takes into account the
uncertainties in modeling the production of secondary antiprotons. Here we
follow the procedure form Refs. [30, 31]. We translate the correlated
uncertainties in the production cross section from Ref. [73] into the
covariance matrix $V_{\rm xs}$ that can be applied directly on the antiproton
flux. In more detail, we generate 2400 vectors of the cross section parameters
using the correlated uncertainties from Ref. [73] and then calculate the
antiproton flux with Galprop for each parameter vector. Finally, we have used
the sample covariance for the antiproton flux evaluated at rigidities $R_{i}$
and $R_{j}$:

$\left(V_{\rm
xs}\right)_{ij}=\frac{1}{N-1}\sum\limits^{N}_{k=1}\big{(}\phi_{\bar{p},i}^{k}-\bar{\phi}_{\bar{p},i}\big{)}\big{(}\phi_{\bar{p},j}^{k}-\bar{\phi}_{\bar{p},j}\big{)}\;,$
(3.4)

where the superscript $k$ denotes the parameter vector and $N=2400$ is the
sample size. For this procedure, the CR propagation parameters are fixed to
the best-fit values, meaning that $V_{\text{\rm xs}}$ is slightly different
for the two propagation setups DIFF.BRK and INJ.BRK.

###  3.2 Cosmic ray scans with MultiNest

We want to gain an understanding of the reasonable parameter space of the
parameters relevant for CR propagation. These will be our nuisance parameters
in the likelihoods for the DM models we consider. We use MultiNest [108] as a
tool to efficiently sample the parameter space and obtain the likelihoods of
each sampled point based on our model. MultiNest uses a nested sampling
algorithm that initially samples from the prior and then iteratively samples
and stores parameter points from ordered prior volumes according to their
respective likelihoods in the parameter space. We thus get an ensemble of
samples covering the entire parameter space, with the most dense sampling
within the regions of high likelihood.

For each propagation model described in section 2.2, we perform this fit
twice. Our initial scan includes all of the propagation parameters described
previously with the assumption of no DM signal contributing to the measured CR
spectra, i.e. our null hypothesis. This fit results in our most general
understanding of the parameter space. We choose the priors as uniform
distributions that are either constrained by observations of different CR
species or have been indicated by previous parameter inferences. We
simultaneously fit the most recent AMS-02 antiproton-over-proton ratio
together with the fluxes of protons and helium and the
${}^{3}\mathrm{He}/^{4}\mathrm{He}$ ratio [22, 109]. For proton and Helium we
also use Voyager data [110]. For the MultiNest scans, we do not consider
correlations in the data sets, i.e. we assume that the given error bars can be
interpreted as uncorrelated statistical uncertainties. The correlations in the
anti-proton data, modeled as discussed above, will be included when deriving
DM constraints, see section 3.3.

In the fit, we consider the propagation parameters going into eq. (2) for each
model. We consider the effect of solar modulation, which is modeled with the
force-field potential $\varphi_{\mathrm{AMS-02}}$ for proton and helium.
Antiprotons are allowed to have a different value of the force-field
potential, as described earlier. In the fit this is taken into account by
using the difference of the force-field potentials,
$\varphi_{\mathrm{AMS-02}}-\varphi_{\bar{p}}$, as a free parameter.
Furthermore, the normalization of the proton flux, $\mathrm{norm}_{p}$,666The
global free normalization of all primary and secondary fluxes (except the DM
flux) in Galprop is fixed by choosing the proton flux at 100 GeV. and the 4He
isotopic abundance relative to the proton abundance, ${\rm
Abd}_{{}^{4}\mathrm{He}}$, are free parameters. Finally, we allow for cross
section uncertainties in the 3He production cross section by introducing
nuisance parameters for the cross section normalization and slope,
$A_{\mathrm{XS,^{4}\mathrm{He}\rightarrow^{3}\mathrm{He}}}$ and
$\delta_{\mathrm{XS,^{4}\mathrm{He}\rightarrow^{3}\mathrm{He}}}$, as in Ref.
[59]. The complete setup allows for independent constraints on the propagation
and a fully Bayesian interpretation of all parameters.

The prior ranges and resulting best fit points of the scan (plus their
1$\sigma$ intervals) are shown in table 1. We show the result of the best fit
parameter point applied to the fluxes of $p,\overline{p}$, 3He and He for both
propagation models in figure. 2. While the best fit values for different
physical properties vary significantly for both models, the overall fit to the
data is very similar, in line with the findings of Ref. [59]. The combined
$\chi^{2}$ here is 159.1 (163.7) for the DIFF.BRK (INJ.BRK) model. The total
number of likelihood evaluations for the scan is $1.4\,(3.7)\cdot 10^{6}$.

Table 1:  Priors and posterior for the CR parameters obtained from MultiNest fit. Results are provided for both CR propagation setups, INJ.BRK and DIFF.BRK.  Parameters |  | Priors |  | DIFF.BRK |  | INJ.BRK
---|---|---|---|---|---|---
$\gamma_{1,p}$ |  | 1.2 | – | 2.1 |  | ${2.338}$ | $\\!\\!\\!\\!{}^{+0.008}_{-0.009}$ |  | ${1.72}$ | $\\!\\!\\!\\!{}^{+0.03}_{-0.02}$
$\gamma_{1}$ |  | 1.2 | – | 2.1 |  | ${2.280}$ | $\\!\\!\\!\\!{}^{+0.007}_{-0.008}$ |  | ${1.73}$ | $\\!\\!\\!\\!{}^{+0.03}_{-0.03}$
$\gamma_{2,p}$ |  | 2.1 | – | 2.6 |  | ${2.338}$ | $\\!\\!\\!\\!{}^{+0.008}_{-0.009}$ |  | ${2.448}$ | $\\!\\!\\!\\!{}^{+0.008}_{-0.007}$
$\gamma_{2}$ |  | 2.1 | – | 2.6 |  | ${2.280}$ | $\\!\\!\\!\\!{}^{+0.007}_{-0.008}$ |  | ${2.393}$ | $\\!\\!\\!\\!{}^{+0.007}_{-0.007}$
$R_{0}\,\mathrm{[GV]}$ |  | 1.0 | – | 20 |  | - |  |  | ${6.43}$ | $\\!\\!\\!\\!{}^{+0.38}_{-0.48}$
$s$ |  | 0.1 | – | 0.7 |  | - |  |  | ${0.33}$ | $\\!\\!\\!\\!{}^{+0.02}_{-0.03}$
$D_{0}\,\mathrm{[10^{28}\;cm^{2}/s]}$ |  | 0.5 | – | 10.0 |  | ${3.78}$ | $\\!\\!\\!\\!{}^{+0.10}_{-0.13}$ |  | ${4.10}$ | $\\!\\!\\!\\!{}^{+0.11}_{-0.08}$
$\delta_{l}$ |  | $-1.0$ | – | 0.5 |  | ${-0.66}$ | $\\!\\!\\!\\!{}^{+0.07}_{-0.06}$ |  | ${0.372}$ | $\\!\\!\\!\\!{}^{+0.007}_{-0.008}$
$\delta$ |  | 0.3 | – | 0.7 |  | ${0.516}$ | $\\!\\!\\!\\!{}^{+0.010}_{-0.008}$ |  | ${0.372}$ | $\\!\\!\\!\\!{}^{+0.007}_{-0.008}$
$\delta_{h}-\delta$ |  | $-0.2$ | – | 0.0 |  | ${-0.16}$ | $\\!\\!\\!\\!{}^{+0.02}_{-0.01}$ |  | ${-0.09}$ | $\\!\\!\\!\\!{}^{+0.02}_{-0.01}$
$R_{D,0}\,\mathrm{[GV]}$ |  | 1.0 | – | 20.0 |  | ${3.91}$ | $\\!\\!\\!\\!{}^{+0.24}_{-0.25}$ |  | - |
$s_{D}$ |  | 0.1 | – | 0.9 |  | ${0.41}$ | $\\!\\!\\!\\!{}^{+0.02}_{-0.02}$ |  | - |
$R_{D,1}\,[10^{3}]$ |  | 100 | – | 500 |  | ${222}$ | $\\!\\!\\!\\!{}^{+21}_{-21}$ |  | ${234}$ | $\\!\\!\\!\\!{}^{+24}_{-35}$
$v_{\mathrm{A}}\,\mathrm{[km/s]}$ |  | 0 | – | 30 |  | - |  |  | ${20.40}$ | $\\!\\!\\!\\!{}^{+1.10}_{-0.64}$
$v_{0,\mathrm{c}}\,\mathrm{[km/s]}$ |  | 0 | – | 60 |  | ${1.91}$ | $\\!\\!\\!\\!{}^{+0.38}_{-1.91}$ |  | ${0.64}$ | $\\!\\!\\!\\!{}^{+0.10}_{-0.64}$
$\varphi_{\mathrm{AMS-02}}\,\mathrm{[GV]}$ |  | 0.4 | – | 0.8 |  | ${0.43}$ | $\\!\\!\\!\\!{}^{+0.01}_{-0.01}$ |  | ${0.62}$ | $\\!\\!\\!\\!{}^{+0.01}_{-0.01}$
$(\varphi_{\bar{p}}-\varphi_{\mathrm{AMS-02}})\,\mathrm{[GV]}$ |  | $-0.2$ | – | 0.2 |  | ${0.181}$ | $\\!\\!\\!\\!{}^{+0.019}_{-0.004}$ |  | ${-0.165}$ | $\\!\\!\\!\\!{}^{+0.006}_{-0.035}$
$\mathrm{norm}_{p}\,[10^{-8}\,\mathrm{MeV^{-1}cm^{-2}s^{-1}sr^{-1}]}$ |  | 0.3 | – | 0.5 |  | ${0.430}$ | $\\!\\!\\!\\!{}^{+0.001}_{-0.001}$ |  | ${0.434}$ | $\\!\\!\\!\\!{}^{+0.001}_{-0.001}$
${\rm Abd}_{{}^{4}\mathrm{He}}\,[10^{5}]$ |  | 0.7 | – | 1.3 |  | ${1.063}$ | $\\!\\!\\!\\!{}^{+0.007}_{-0.007}$ |  | ${0.984}$ | $\\!\\!\\!\\!{}^{+0.006}_{-0.006}$
$A_{\mathrm{XS,^{4}\mathrm{He}\rightarrow^{3}\mathrm{He}}}$ |  | 0.8 | – | 1.2 |  | ${1.17}$ | $\\!\\!\\!\\!{}^{+0.01}_{-0.01}$ |  | ${1.177}$ | $\\!\\!\\!\\!{}^{+0.023}_{-0.006}$
$\delta_{\mathrm{XS,^{4}\mathrm{He}\rightarrow^{3}\mathrm{He}}}$ |  | -0.2 | – | 0.2 |  | ${0.01}$ | $\\!\\!\\!\\!{}^{+0.02}_{-0.02}$ |  | ${0.190}$ | $\\!\\!\\!\\!{}^{+0.010}_{-0.002}$

Figure 2:  Result of the CR fits with MultiNest for the two propagation setups
DIFF.BRK (left) and INJ.BRK (right). We show the best fit for $p$, He, 3He,
and $\bar{p}$. Solid lines correspond to the CR flux including solar
modulation for the period of AMS-02 data-taking, while dashed lines show the
local interstellar flux as measured by Voyager. In the lower panel, we show
the residuals of the $\bar{p}$ data.

In a second step, we also include a tentative DM signal and add the DM mass
and annihilation cross-section as free parameters to the MultiNest scan. For
simplicity, we fix the annihilation channel to $b\overline{b}$, which is a
standard benchmark choice. Here, we are not primarily interested in the
likelihood of the final result of the scan, but rather in how the posteriors
of the propagation parameters shift in order to accommodate the additional DM
signal. This is important for one of the two use-cases of the MultiNest
samples for our results. In this work we consider more general DM models than
the one considered in this scan, assuming annihilation only into
$b\overline{b}$, but the overall shift in the posterior for the individual
propagation parameters can well be approximated with this assumption. We have
found in particular a shift for $\gamma_{1,p}$ and $\gamma_{1}$ at the level
of about $1.7\sigma$ towards lower values, and for $\delta$ towards a higher
value at about $1.5\sigma$ in the DIFF.BRK model. The remaining fit
parameters, as well as the INJ.BRK model were affected less by the additional
signal, i.e. the shift was at most $1\sigma$.

We use the combination of both scans to set up an extensive training set for
the DarkRayNet simulation tool, a deep neural network emulator that we set up
as in Ref. [27]. In the remainder of this work, this tool is used to speed up
the required simulations of CR-spectra.

For both scans in each propagation model, we store all simulated spectra with
a $\Delta\chi^{2}\leq$ 30 and combine them as a training set. We include model
independent DM signals by sampling four sets of DM masses and annihilation
branching fractions for each set of propagation parameters. We set a flat
prior on a logarithmic scale for the DM mass between 5 GeV and 5 TeV, and
consider annihilation into $q\overline{q},\ c\overline{c},\ b\overline{b},\
t\overline{t},\ W^{+}W^{-},\ ZZ,\ gg,\ hh$, as we expect non-negligible
contributions to a resulting antiproton flux from these channels, and sample
random branching fractions, always summing up to one. We don’t need to vary
the annihilation cross section in the emulator, as it can be inserted as a
normalization to the DM signal later. With the final training sample for each
propagation model, the fully trained network can emulate accurate CR-spectra
within an extended, relevant parameter space. With the addition of the second
scan, the shift in preferred regions of the propagation parameters due to that
additional signal are accounted for. More details on the DarkRayNet and the
involved neural networks can be found in appendix A.1.

The second use-case of the MultiNest scans is the generation of a sample of
propagation parameters, which can then be used to perform the marginalization
as described in more detail in the following section. For this, we extract the
equally weighted posterior samples from the scan.

###  3.3 Marginalisation over propagation parameters

The likelihood introduced above can be used to perform a simultaneous scan
over DM and propagation parameters. In the context of constraining DM models,
it is however often unnecessary to infer the preferred propagation parameters.
In this case, it is convenient to eliminate the dependence on the propagation
parameters and obtain a likelihood function that depends exclusively on the DM
parameters. Since fixing the propagation parameters to specific values may
lead to overly aggressive constraints, the two possible options are either
profiling or marginalisation. The former approach, i.e. maximising the
likelihood with respect to the propagation parameters for every choice of
$\bm{x}_{\text{DM}}$ is not only computationally expensive, but may be
sensitive to finely-tuned parameter regions where small excesses in the data
can be fitted. Moreover, when using DarkRayNet this approach may be
susceptible to inaccurate network predictions for parameter regions with
insufficient training (see the more detailed discussion in Ref. [27]).

The most robust approach is therefore to calculate the marginalised likelihood

$\bar{\mathcal{L}}_{\text{tot}}(\bm{x}_{\text{DM}})=\int\mathrm{d}\bm{\theta}_{\text{prop}}\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}},\bm{\theta}_{\text{prop}})\pi(\bm{\theta}_{\text{prop}})\;,$
(3.5)

where
$\mathcal{L}_{\text{tot}}=\mathcal{L}_{\bar{p}}\cdot\mathcal{L}_{p}\cdot\mathcal{L}_{\text{He}}$
and $\pi(\bm{\theta}_{\text{prop}})$ denotes the assumed prior probability for
the propagation parameters. In the following, we will adopt flat priors and
consider the same parameter regions used for the MultiNest scans described
above, i.e. we set
$\pi(\bm{\theta}_{\text{prop}})=\Pi_{i}(\theta_{i,\text{max}}-\theta_{i,\text{min}})^{-1}$.
In principle, the marginalised likelihood can be obtained through Monte Carlo
integration by drawing a random sample $\\{\bm{\theta}_{{\rm prop},i}\\}$ of
size $N$ from the prior probability distribution and writing

$\bar{\mathcal{L}}_{\text{tot}}(\bm{x}_{\text{DM}})=\frac{1}{N}\sum_{\bm{\theta}_{{\rm
prop},i}\sim\pi}\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm
prop},i})\;.$ (3.6)

In practice, a more accurate estimate777The gain in accuracy is both due to
the fact that more weight is given to parameter regions where the likelihood
is large, thus decreasing the uncertainty of the Monte Carlo integration, and
due to the fact that DarkRayNet is only evaluated in parameter regions where
sufficient training data is available to guarantee reliable predictions. is
obtained by drawing samples from the posterior probability for the propagation
parameters in the absence of a DM signal
$\mathcal{P}(\bm{\theta}_{\text{prop}})\propto\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{\text{prop}})\pi(\bm{\theta}_{\text{prop}})$.
In this case,

$\bar{\mathcal{L}}_{\text{tot}}(\bm{x}_{\text{DM}})\propto\frac{1}{N}\sum_{\bm{\theta}_{{\rm
prop},i}\sim\mathcal{P}}\frac{\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm
prop},i})}{\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{{\rm
prop},i})}\;.$ (3.7)

The constant of proportionality is independent of $\bm{x}_{\text{DM}}$ and
therefore drops out when calculating likelihood ratios.

Now we can make use of the fact that the contribution of DM to the local
proton and Helium flux is completely negligible. Hence, the proton and Helium
likelihoods $\mathcal{L}_{p,\text{He}}$ are independent of
$\bm{x}_{\text{DM}}$ and therefore cancel in the likelihood ratio. We can
therefore replace $\mathcal{L}_{\text{tot}}\to\mathcal{L}_{\bar{p}}$ in the
equation above. The marginalised $\chi^{2}$ is then given by

$\displaystyle\bar{\chi}^{2}(\bm{x}_{\text{DM}})$
$\displaystyle=-2\log\left(\sum_{\bm{\theta}_{{\rm
prop},i}\sim\mathcal{P}}\frac{\mathcal{L}_{\bar{p}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm
prop},i})}{\mathcal{L}_{\bar{p}}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{{\rm
prop},i})}\right)$ $\displaystyle=-2\log\left[\sum_{\bm{\theta}_{{\rm
prop},i}\sim\mathcal{P}}\exp\left(-\frac{\chi^{2}_{\bar{p}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm
prop},i})-\chi^{2}_{\bar{p}}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{{\rm
prop},i})}{2}\right)\right]$ (3.8)

with $\chi^{2}_{\bar{p}}$ defined in eq. (3.1). We note that by construction
$\bar{\chi}^{2}(\bm{x}_{\text{DM}}=0)=0$.

In general, drawing a random sample from the posterior probability is far from
easy. Fortunately, this sample only needs to be generated once and can then be
used for arbitrary DM parameters. Moreover, such a sample is generated
automatically as a by-product of the MultiNest scans that we have run to
generate the training data for DarkRayNet and therefore requires no additional
calculations. This sample of about $10^{4}$ sets of propagation parameters is
included in pbarlike and used by default for marginalisation.

To conclude this discussion, we note that it is in fact possible to use a
different definition of the antiproton likelihood for the marginalised
likelihood than the one that has been used to calculate the posterior
probability. In this case, eq. (3.8) simply becomes

$\bar{\mathcal{L}}_{\text{tot}}(\bm{x}_{\text{DM}})\propto\frac{1}{N}\sum_{\bm{\theta}_{{\rm
prop},i}\sim\mathcal{P}_{0}}\frac{\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm
prop},i})}{\mathcal{L}_{\text{tot},0}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{{\rm
prop},i})}$ (3.9)

with
$\mathcal{P}_{0}(\bm{\theta}_{\text{prop}})\propto\mathcal{L}_{\text{tot},0}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{\text{prop}})\pi(\bm{\theta}_{\text{prop}})$.
This means in particular that it is possible to modify the covariance matrix
without the need to rerun the computationally expensive MultiNest scan.

###  3.4 Constraints on DM annihilations into $b\bar{b}$

Assuming all of DM annihilates only into bottom quarks, we calculate the
marginalised log-likelihood ratios pbarlike, for each point on a 2D grid of DM
model parameters $\\{m_{DM},\langle\sigma v\rangle\\}$. We consider the DM
parameter ranges, $m_{DM}\in[10\ \mathrm{GeV},5\ \mathrm{TeV}]$ and
$\langle\sigma v\rangle\in[10^{-28},10^{-24}]\ \mathrm{cm^{3}s^{-1}}$. An
entire routine consisting of obtaining antiproton flux prediction from
DarkRayNet, solar modulation and marginalised log-likelihood ratio calculation
for a single point in DM model parameter space takes $\mathcal{O}(10)\
\mathrm{s}$. The evaluation of $\mathcal{O}(10^{4})$ antiproton fluxes (for
the entire sample of $\mathcal{O}(10^{4})$ propagation parameter vectors used
for marginalization) using DarkRayNet takes only $\mathcal{O}(1)\,\mathrm{s}$,
whereas Galprop would take $\mathcal{O}(10)$ s for each set of propagation
parameters. The most time-consuming part of the likelihood calculation turns
out to be the modulation of each of $\mathcal{O}(10^{4})$ antiproton fluxes
with appropriate solar modulation parameters.

Figure 3:  $\bar{\chi}^{2}$ for DM annihilating into $b\bar{b}$, obtained
using the two CR propagation models-DIFF.BRK (left) and INJ.BRK (right). The
black line represents the $95\%$ CL upper limit for $\langle{\sigma v}\rangle$
as a function of $m_{DM}$.

The results are shown in figure 3. The region in colour correspond to
$\bar{\chi}^{2}>0$ and hence prefers the case of antiprotons from a purely
secondary origin. The region in gray corresponds to $\bar{\chi}^{2}<0$ and
thus exhibits a small preference for the presence of a DM signal. The black
line shows the $95\%$ CL upper bound for $\langle\sigma v\rangle$ as a
function of DM mass $m_{DM}$ (corresponding to $\bar{\chi}^{2}=3.84$). For the
DIFF.BRK propagation model (left panel), preference for DM signal is seen in
$\mathrm{TeV}$ masses and for INJ.BRK (right panel) around
$100\,\mathrm{GeV}$. The bound obtained using INJ.BRK model is similar to ones
previously seen in literature [30, 31, 32, 24]. For the DIFF.BRK model,
however, the DM preference is pushed to $\mathrm{TeV}$ masses. This can be
understood with the help of residuals in figure 2. The additional break at low
rigidities in the diffusion coefficient in the DIFF.BRK model allows for a
good fit to low rigidity data, whereas it performs a poorer fit to data at
high rigidities compared to INJ.BRK model.

To conclude this section, we note that pbarlike also returns log-likelihood
ratios calculated using uncorrelated errors given by AMS-02. For more details,
we refer to appendix B.

##  4 Application to scalar singlet dark matter

We are now in the position to apply the AMS-02 antiproton likelihood to more
realistic DM models. As an example, we consider the scalar singlet DM model,
which is briefly reviewed in section 4.1. We discuss the additional
constraints that we have implemented in section 4.2 and present our results in
section 4.3.

###  4.1 Model details

In the scalar singlet DM model, the SM is extended by a gauge-singlet real
scalar boson $s$ stabilised by a $\mathbb{Z}_{2}$ symmetry [36, 37, 38]. The
resulting scalar potential can then be written as

$V(s^{2},H^{\dagger}H)=\lambda_{h}\left[\left(H^{\dagger}H\right)-\frac{v^{2}}{2}\right]^{2}+\frac{1}{2}\,\lambda_{hs}\,s^{2}\,H^{\dagger}H+\frac{1}{4}\,\lambda_{s}\,s^{4}+\frac{1}{2}\,m_{s,0}^{2}\,s^{2}$
(4.1)

with $v=246\>\text{GeV}$ denoting the vacuum expectation value of the SM Higgs
field. After electroweak symmetry breaking, we can replace
$H=(h+v,0)/\sqrt{2}$ with a real scalar $h$ to obtain the potential

$V(s^{2},h)=V(h)+\frac{1}{2}\,m_{s}^{2}\,s^{2}+\frac{1}{4}\,\lambda_{s}\,s^{4}+\frac{1}{2}\,\lambda_{hs}\,v\,h\,s^{2}+\frac{1}{4}\,\lambda_{hs}\,h^{2}\,s^{2}\;,$
(4.2)

where $m_{s}=\left[m_{s,0}^{2}+\lambda_{hs}\,v^{2}/2\right]^{1/2}$ denotes the
physical mass of the scalar singlet.

The quartic self-coupling $\lambda$ typically has no consequences for
phenomenology, and we will therefore set it to zero in the following. The
model is then fully characterised by the scalar singlet mass $m_{s}$ and the
Higgs portal coupling $\lambda_{hs}$.888In the analysis below, we also include
a number of nuisance parameters, most notably the local DM density
$\rho_{\odot}$, which affects direct and indirect detection constraints. See
Ref. [40] for details. The latter is responsible for all interactions of $s$
with SM particles, which determine the relic density of $s$ via the freeze-out
mechanism and the potentially observable signatures of $s$ in satellites and
laboratory experiments.

###  4.2 Implementation of new constraints

The phenomenology of the scalar singlet DM model has been discussed in great
detail in Refs. [39, 40, 41]. Here, we repeat the analysis procedure from Ref.
[40], i.e. we calculate the relic density $\Omega_{s}$ of scalar singlets and
require $f_{s}\equiv\Omega_{s}/\Omega_{\text{DM}}\leq 1$, where
$\Omega_{\text{DM}}=0.12$ [111]. Furthermore, we consider constraints from
direct detection experiments (rescaled by the fractional abundance $f_{s}$ of
scalar singlets), indirect detection experiments (rescaled by a factor
$f_{s}^{2}$) and the LHC as well as the perturbativity requirement
$\lambda_{s}<\sqrt{4\pi}$. In addition to previously considered constraints,
we include the antiproton likelihood from AMS-02 discussed above, as well as a
number of additional new likelihoods as detailed below.

####  4.2.1 Higgs invisible width

The partial width for the decay of a Higgs boson into a pair of scalar
singlets is given by

$\Gamma(h\rightarrow
ss)=\frac{\lambda_{hs}^{2}v^{2}}{32\pi\,m_{h}}\sqrt{1-\frac{4\,m_{s}^{2}}{m_{h}^{2}}}\;,$
(4.3)

from which the invisible branching ratio can be calculated as

$\text{BR}(h\rightarrow ss)=\frac{\Gamma(h\rightarrow ss)}{\Gamma(h\rightarrow
ss)+\Gamma_{\text{SM}}}\;.$ (4.4)

Both CMS [42] and ATLAS [43] have recently published new constraints on the
Higgs invisible branching ratio using the vector boson fusion channel.

CMS. The CMS constraint is provided in the form of a likelihood as a function
of $\text{BR}(h\to ss)$, see figure 12 of Ref. [42]. This likelihood can be
very well approximated by the quadratic function

$-2\log\mathcal{L}=a(\text{BR}(h\to ss)-b)^{2}$ (4.5)

with $a\approx 339$ and $b\approx 0.089$.

ATLAS. For ATLAS, no detailed likelihood function is available. However,
assuming that the ATLAS likelihood can also be written in the form of eq.
(4.5), we can infer the two parameters from the stated bounds $\text{BR}(h\to
ss)<0.127\ (0.145)$ at 90% (95%) confidence level, which yields $a\approx 303$
and $b\approx 0.032$.

The combined constraint from CMS and ATLAS gives $\text{BR}(h\to ss)<0.14$ at
95% confidence level. The best-fit value is $\text{BR}(h\to ss)=0.06$, which
is preferred over $\text{BR}=0$ with a significance of $1.2\sigma$.

####  4.2.2 Direct detection

The DM–nucleon scattering cross section at zero momentum transfer is given by

$\sigma_{\text{SI}}=\frac{\lambda_{hs}^{2}f_{N}^{2}}{4\pi}\frac{\mu^{2}\,m_{N}^{2}}{m_{h}^{4}\,m_{s}^{2}}\;,$
(4.6)

where $m_{N}$ is the nucleon mass, $\mu=(m_{s}\,m_{N})/(m_{s}+m_{N})$ is the
reduced mass and $f_{N}$ is the effective Higgs-nucleon coupling. Recent
constraints on $\sigma_{\text{SI}}$ as a function of $m_{s}$ have been
published by LZ and PandaX.

LZ. We consider events below the 90% quantile of the nuclear recoil band, such
that the effective exposure is $2.97\cdot 10^{5}\,\mathrm{kg\,days}$. At
$E_{\mathrm{R}}=15\,\mathrm{keV}$ we split the search region into two bins,
assuming an energy resolution of $\sigma=1.5\,\mathrm{keV}$. We then fix the
background expectation in each bin in such a way that the expected limit is
recovered when setting the observed number of events equal to the background
expectation. We find that this procedure yields 1 (7) expected background
event in the lower (upper) bin. The actual observation gave no events in the
lower and twelve events in the upper bin. This leads to an observed exclusion
limit somewhat stronger (slightly weaker) than the expected one for small
(large) DM masses, in agreement with the published LZ result.

PandaX-4T. We consider events below the mean of the nuclear recoil band,
leading to an effective exposure of $1.15\cdot 10^{5}\,\mathrm{kg\,days}$. In
this region, the expected background was 9.8 events, compared to 6 observed
events. To calculate the signal prediction, we take the efficiency from figure
2 of Ref. [45].

###  4.3 Results

In this section, we present updated constraints on the scalar singlet DM
model. To do this, one could perform a comprehensive global scan with all the
relevant available datasets and the resulting $\mathcal{O}(10)$ nuisance
parameters. Instead, for a judicious use of computational resources, we take
advantage of the publicly available results [112] (hereafter referred to as
GC17) from the extensive global analysis performed by the GAMBIT collaboration
in Ref. [40]. This analysis combined the then available relic density
measurements, direct detection limits, limits from invisible Higgs decays, and
indirect detection limits from dark annihilating into neutrinos and gamma-
rays. The results in GC17 contain parameter spaces allowed by the limits
considered and their corresponding combined likelihoods.

The large sample of $\mathcal{O}(10^{7})$ points in GC17 contains combined
results from multiple sampling runs. These scans included the DM model
parameters in the range $m_{s}\in[45\,\mathrm{GeV},10\,\mathrm{TeV}]$ and
$\lambda_{hs}\in[10^{-4},10]$. For details on the nuisance parameters and
their ranges, see table 2 of Ref. [40]. The scans identified the best-fit at
$\lambda_{hs}=6.5\times 10^{-4},\ m_{s}=62.51\,\mathrm{GeV}$ corresponding to
$f_{s}^{2}\langle\sigma v\rangle_{0}=8.65\times
10^{-28}\,\mathrm{cm^{3}s^{-1}}$.

Figure 4: Allowed parameter space of the scalar singlet model in the
$f_{s}^{2}\langle\sigma v\rangle_{0}$–$m_{s}$ plane. The plots show the
modification of parameter space, on updating constraints in [40] with the
addition of AMS-02 antiproton likelihood. The top (bottom) row shows
constraints obtained using DIFF.BRK (INJ.BRK) model.

For our analysis, we restrict the parameter space to
$m_{s}\in[45\,\mathrm{GeV},5\,\mathrm{TeV}]$ to remain within the parameter
region where DarkRayNet has been trained. We further reduce the GC17 sample to
$\mathcal{O}(10^{5})$ points. We then update the old combined likelihoods in
GC17 by adding the new likelihoods discussed above using the postprocessor
scanner available within the GAMBIT sampling module ScannerBit [113]. We
emphasize that this procedure is possible because there is no strong DM signal
in the AMS-02 data, and hence we do not expect the new likelihoods to open up
previously disfavoured parameter regions.

In a first step, we only add the new AMS-02 antiproton likelihood, to explore
its impact on the allowed parameter space. The resulting updated constraints
on the $f_{s}^{2}\langle\sigma v\rangle_{0}$–$m_{s}$ plane are shown in figure
4.999Plots were made using pippi [114]. Each panel shows the profiled
likelihood in the DM model parameter space, normalized by the best-fit
likelihood.101010Note that the total likelihood here is profiled over the
nuisance parameters in GC17 only; pbarlike already returns the marginalized
likelihood and hence, the CR propagation and solar modulation parameters are
not included in the GAMBIT postprocessing run as nuisance parameters. The
contours show the $1\sigma$ and $2\sigma$ confidence regions, and the star
indicates the best-fit point. The viable regions identified in previous
studies, i.e. the resonance region with $m_{s}\approx m_{h}/2$ and small cross
sections (left column) and the high-mass regions with $m_{s}\gg m_{h}$ and
large cross sections (right column), are still retained. For comparison, see
the right panel of figure 3 in Ref. [40].

The two rows of figure 4 correspond to the DIFF.BRK (top) and INJ.BRK (bottom)
models. For both propagation models, the likelihood of the best-fit point is
changed by $-2\Delta\log\mathcal{L}<0.1$. Neither propagation models shift the
best-fit point to different masses. For DIFF.BRK, the best-fit point is
however shifted to smaller couplings and a smaller log-likelihood, as this
model leads to strong bounds on the rescaled annihilation cross section for DM
masses in the 10–100$\,\mathrm{GeV}$ range (see also figure 3). For the
INJ.BRK model, the small excess seen in figure 3 leads to a weaker bound on
the rescaled cross section and a correspondingly larger coupling at the best-
fit point with a larger log-likelihood. In the TeV range, on the other hand,
the INJ.BRK model gives the stronger constraint than the DIFF.BRK model,
leading to a smaller allowed parameter region. This observation highlights the
relevant impact of the AMS-02 antiproton constraints on the scalar singlet DM
model.

Figure 5: Allowed parameter space of the scalar singlet model in the
$\lambda_{s}$–$m_{s}$ plane. In addition to the constraints considered in
figure 4, we also include the most recent constraints from direct detection
experiments (PandaX-4T and LZ) and the most up-to-date measurements of the
Higgs invisible branching ratio from ATLAS and CMS. In the top (bottom) row,
we consider the AMS-02 antiproton likelihood using the DIFF.BRK (INJ.BRK+vA)
propagation model.

In a second step, we now also include the new direct detection constraints
using DDCalc [5, 115] and the Higgs invisible branching ratio likelihood,
which has been directly implemented within GAMBIT. The resulting updated
constraints on the $\lambda_{hs}$–$m_{s}$ plane are shown in figure 5 for the
DIFF.BRK (top) and INJ.BRK (bottom) models. As before, the figure also shows a
zoom-in of the resonance region with mass in linear scale (left). The impact
of the antiproton constraints is less visible in this case, as the relevant
regions are already excluded by the new direct detection constraints. Hence,
the top and bottom panels look almost identical. For both the propagation
models, the likelihood of the best-fit is changed by
$-2\Delta\log\mathcal{L}\approx 2.5$, driven by the small excess seen in Higgs
invisible decays. Thus, the coupling of the best-fit is also shifted to larger
values ($\sim 0.02$). In the resonance region, the allowed parameter space is
constrained from the left by the improved direct detection constraints, which
also exclude a large region of the parameter space in the high-mass region.
The remaining parameter space will be further explored by future direct DM
detection experiments [116].

##  5 Conclusions

Observations of cosmic rays (CRs), in particular antiparticle fluxes, are
among the most sensitive probes of dark matter (DM) annihilations, provided
that the uncertainties related to injection and propagation are under control.
While various codes exist to simulate CR propagation, the computational cost
in combination with the complexity of the underlying propagation models make
it very difficult to comprehensively explore indirect detection constraints
for different DM models. In this work we have addressed this challenge for the
case of the AMS-02 antiproton data by combining two key ingredients: neural
networks capable of accurately predicting primary and secondary antiproton
fluxes, and efficient marginalisation of propagation uncertainties in the
calculation of AMS-02 likelihoods. The required software tools, DarkRayNet.v2
and pbarlike, are publicly available and can be readily interfaced with the
most recent release of the GAMBIT global fitting framework, which automates
the calculation of cross sections and branching ratios from Lagrangian-level
inputs [117].

Using importance sampling, the calculation of the marginalised likelihood for
a specific DM hypothesis requires calculating the primary and secondary
antiproton fluxes and the resulting AMS-02 likelihood for approximately
$10^{4}$ different combinations of propagation and solar modulation
parameters. In our approach, this calculation takes approximately
$10\,\mathrm{s}$ on a single cpu, which is sufficiently fast to include
antiproton constraints in global fits of DM models. As an illustration, we
have considered the scalar singlet DM model, updating all relevant constraints
to include the most recent experimental results. While in principle this model
can give sizeable signals in indirect detection experiments, the relevant
high-mass regions of parameter space are largely excluded by the most recent
direct detection experiments. Given these constraints, as well as a slight
excess in LHC measurements of the Higgs invisible branching ratio, the
resonance region remains as the parameter region of primary interest for this
model.

Our work also highlights the need for flexible propagation models and an
accurate treatment of experimental uncertainties. All our analyses have been
carried out using two different propagation models, which introduce a break in
the injection spectrum and in the diffusion coefficient, respectively. We show
that, for the case of annihilation into $b\bar{b}$, the former model gives
rise to a slight excess for DM masses around $100\,\mathrm{GeV}$ and tight
constraints on TeV-scale DM, while the latter yields an excess for DM masses
in the TeV range and very strong constraints on DM masses below
$100\,\mathrm{GeV}$. Once correlations in AMS-02 data are taken into account,
the local significance of any of these excesses is well below $3\sigma$.

At present there is hence no evidence for a DM signal in AMS-02 data and only
upper bounds on the interactions of DM particles can be obtained.
Nevertheless, our analysis clearly demonstrates the strong potential of using
the antiproton flux to search for DM. The software tools that we release
together with this work make it possible to fully exploit this potential and
efficiently reinterpret experimental data for a wide range of DM models. These
methods will become more and more important as our understanding of CR
propagation continues to improve, flux measurements for heavier antiparticles
are released (e.g. by the GAPS balloon mission [118]), and a new generation of
indirect detection experiments begin to take data [119, 120].

## Acknowledgments

We would like to thank Michael Krämer for discussions and Tomás Gonzalo for
help with the GAMBIT interface. SB acknowledges the support by the Doctoral
School ’Karlsruhe School of Elementary and Astroparticle Physics: Science and
Technology’. FK is partially funded by the Deutsche Forschungsgemeinschaft
(DFG) through the Emmy Noether Grant No. KA 4662/1-1. MK is supported by the
Swedish Research Council under contracts 2019-05135 and 2022-04283 and the
European Research Council under grant 742104. SM acknowledges the European
Union’s Horizon Europe research and innovation programme for support under the
Marie Sklodowska-Curie Action HE MSCA PF–2021, grant agreement No.10106280,
project VerSi. This project used computing resources from the Swedish National
Infrastructure for Computing (SNIC) under project Nos. 2021/3-42, 2021/6-326
and 2021-1-24 partially funded by the Swedish Research Council through grant
no. 2018-05973. The training of the neural networks was performed with
computing resources granted by RWTH Aachen University under project
‘rwth0754’.

##  Appendix A Details on new software tools

###  A.1 DarkRayNet.v2

The DarkRayNet is a neural emulation tool based on recurrent neural networks,
specifically LSTMs [121], that can quickly simulate antiprotons, protons, and
Helium CR spectra at Earth, based on Galprop simulations. It was first
published in Ref. [27] and is designed to predict measurable CR spectra for
different parameters in seconds, allowing quick scans for indirect DM
searches. The neural network architecture of the DarkRayNet is designed to
handle the large set of relevant parameters for CR propagation and the DM
model parameters individually, and convert the physical parameters into the
resulting cosmic-ray spectra, using multiple densely connected layers and
finally an LSTM layer. The network was trained using a large set of fluxes for
each particle type and origin (primary/secondary/DM annihilation). The samples
in the training data were chosen from MultiNest scans used to evaluate the
compatible parameter space with recent AMS-02 data for the parameters
describing CR transport (see also section 3.2). The network accuracy for
predictions within the trained parameter regions is always sufficiently high,
such that a difference between the emulated spectra and the corresponding
conventionally simulated spectra from Galprop will always be significantly
below the uncertainty in the most recent AMS-02 data. As a result, no
systematic will be introduced when using this emulator for likelihood
evaluations. The simulation of cosmic-ray spectra with DarkRayNet is about two
to three orders of magnitude faster than with the standard Galprop setup,
depending on the number of different cosmic-ray species considered.

In this paper, we have updated the model for CR propagation (INJ.BRK) with
respect to the setup in Ref. [27], and included a new model (DIFF.BRK), as
described in section 2.3. With respect to the original release, the new
version v2111111https://github.com/kathrinnp/DarkRayNet of DarkRayNet has thus
been extended with newly trained networks for the updated models.

###  A.2 pbarlike

pbarlike121212https://github.com/sowmiya-balan/pbarlike is an open-source
python code developed for indirect DM searches with antiprotons. One of the
main aims of pbarlike is state-of-the-art treatment of antiproton production
cross section uncertainties and data correlations; the other is to interface
the powerful CR emulator, DarkRayNet with the global fitting framework GAMBIT.
In addition, pbarlike also computes solar modulation of the antiproton flux.
pbarlike treats solar modulation using the Force Field Approximation. It
currently includes the 7-year AMS-02 antiproton dataset [22] and is easily
extendable to include more datasets. For likelihood calculations, one can
choose to either marginalize over all nuisance parameters from CR propagation
and solar modulation, or to profile over only the solar modulation parameters,
keeping the propagation parameters fixed. For marginalization, the MultiNest
sample of all nuisance parameters (described in sec. 3.2) is made available;
for profiling, custom CR propagation parameters need to be provided.

For a chosen propagation model (options are those available within DarkRayNet,
i.e. DIFF.BRK and INJ.BRK), and a given set of CR propagation and DM model
parameters, pbarlike obtains predictions for antiproton flux at the local
interstellar region from DarkRayNet. It then simulates solar modulation of
this flux and finally compares it against the given dataset to calculate
likelihoods. pbarlike can be easily installed and used as either a standalone
for antiproton constraints or in conjunction with GAMBIT for global fits. The
complete documentation can be found at https://pbarlike.readthedocs.io.

###  A.3 GAMBIT interface

The interface between pbarlike and GAMBIT (also called the _frontend_) makes
it possible to automatically evaluate the AMS-02 antiproton likelihood for any
DM model implemented in GAMBIT. For any such model, all relevant annihilation
cross sections are provided by the ProcessCatalog object [5], such that the
branching fractions for different final states and the total annihilation
cross section can be readily calculated. For new models, the Process Catalogue
may be automatically generated from Lagrangian-level input using CalcHEP [122]
via the GAMBIT Universal Model Machine [117].

##  Appendix B Comparison of correlated and uncorrelated likelihoods

Figure 6:  $\bar{\chi}^{2}$ for DM annihilating into $b\bar{b}$, obtained
using the two CR propagation models - DIFF.BRK (top), INJ.BRK (bottom) - and
using uncorrelated (left) and correlated (right) errors.

Figure 6 compares the result from section 3.4 to the ones obtained when using
uncorrelated errors. It can be observed that the inclusion of correlations
leads to a decrease in local significance from $3.9\sigma$ to $1.0\sigma$ for
the INJ.BRK model (assuming that $-\bar{\chi}^{2}_{\text{min}}$ follows a
$\chi^{2}$ distribution with two degrees of freedom). The same effect is not
observed in DIFF.BRK model where the local significance does not change by a
lot when moving from uncorrelated ($2.3\sigma$) to correlated ($2.2\sigma$)
errors, the reason being that errors are dominated by correlated systematic
uncertainties only in the low rigidity range.

## References

  * [1] 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].
  * [2] 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].
  * [3] F. Ambrogi, C. Arina, M. Backovic, J. Heisig, F. Maltoni, et al., MadDM v.3.0: a Comprehensive Tool for Dark Matter Studies,  Phys. Dark Univ. 24 (2019) 100249, [1804.00044].
  * [4] A. Arbey, F. Mahmoudi, and G. Robbins, SuperIso Relic v4: A program for calculating dark matter and flavour physics observables in Supersymmetry,  Comput. Phys. Commun. 239 (2019) 238–264, [1806.11489].
  * [5] GAMBIT Dark Matter Workgroup Collaboration, T. Bringmann et al., DarkBit: A GAMBIT module for computing dark matter observables and likelihoods,  Eur. Phys. J. C 77 (2017), no. 12 831, [1705.07920].
  * [6] GAMBIT Collaboration, P. Athron et al., GAMBIT: The Global and Modular Beyond-the-Standard-Model Inference Tool,  Eur. Phys. J. C 77 (2017), no. 11 784, [1705.07908]. [Addendum: Eur.Phys.J.C 78, 98 (2018)].
  * [7] L. Bergstrom, J. Edsjo, and P. Ullio, Cosmic anti-protons as a probe for supersymmetric dark matter?,  Astrophys. J. 526 (1999) 215–235, [astro-ph/9902012].
  * [8] F. Donato, N. Fornengo, D. Maurin, and P. Salati, Antiprotons in cosmic rays from neutralino annihilation,  Phys. Rev. D69 (2004) 063501, [astro-ph/0306207].
  * [9] T. Bringmann and P. Salati, The galactic antiproton spectrum at high energies: Background expectation vs. exotic contributions,  Phys. Rev. D75 (2007) 083006, [astro-ph/0612514].
  * [10] F. Donato, D. Maurin, P. Brun, T. Delahaye, and P. Salati, Constraints on WIMP Dark Matter from the High Energy PAMELA $\bar{p}/p$ data,  Phys. Rev. Lett. 102 (2009) 071301, [0810.5292].
  * [11] N. Fornengo, L. Maccione, and A. Vittino, Constraints on particle dark matter from cosmic-ray antiprotons,  JCAP 1404 (2014), no. 04 003, [1312.3579].
  * [12] C. Evoli, I. Cholis, D. Grasso, L. Maccione, and P. Ullio, Antiprotons from dark matter annihilation in the Galaxy: astrophysical uncertainties,  Phys. Rev. D85 (2012) 123511, [1108.0664].
  * [13] T. Bringmann, M. Vollmann, and C. Weniger, Updated cosmic-ray and radio constraints on light dark matter: Implications for the GeV gamma-ray excess at the Galactic center,  Phys. Rev. D90 (2014), no. 12 123001, [1406.6027].
  * [14] M. Cirelli, D. Gaggero, G. Giesen, M. Taoso, and A. Urbano, Antiproton constraints on the GeV gamma-ray excess: a comprehensive analysis,  JCAP 1412 (2014), no. 12 045, [1407.2173].
  * [15] J. A. R. Cembranos, V. Gammaldi, and A. L. Maroto, Antiproton signatures from astrophysical and dark matter sources at the galactic center,  JCAP 1503 (2015), no. 03 041, [1410.6689].
  * [16] D. Hooper, T. Linden, and P. Mertsch, What Does The PAMELA Antiproton Spectrum Tell Us About Dark Matter?,  JCAP 1503 (2015), no. 03 021, [1410.1527].
  * [17] M. Boudaud, M. Cirelli, G. Giesen, and P. Salati, A fussy revisitation of antiprotons as a tool for Dark Matter searches,  JCAP 1505 (2015), no. 05 013, [1412.5696].
  * [18] G. Giesen, M. Boudaud, Y. Genolini, V. Poulin, M. Cirelli, et al., AMS-02 antiprotons, at last! Secondary astrophysical component and immediate implications for Dark Matter,  JCAP 1509 (2015), no. 09 023, [1504.04276].
  * [19] A. Reinert and M. W. Winkler, A Precision Search for WIMPs with Charged Cosmic Rays,  JCAP 1801 (2018), no. 01 055, [1712.00002].
  * [20] C. Evoli, D. Gaggero, and D. Grasso, Secondary antiprotons as a Galactic Dark Matter probe,  JCAP 1512 (2015), no. 12 039, [1504.05175].
  * [21] P. D. L. T. Luque, Combined analyses of the antiproton production from cosmic-ray interactions and its possible dark matter origin,  JCAP 11 (2021) 018, [2107.06863].
  * [22] AMS Collaboration, M. Aguilar et al., The Alpha Magnetic Spectrometer (AMS) on the international space station: Part II — Results from the first seven years,  Phys. Rept. 894 (2021) 1–116.
  * [23] M. Korsmeier and A. Cuoco, Implications of Lithium to Oxygen AMS-02 spectra on our understanding of cosmic-ray diffusion,  Phys. Rev. D 103 (2021), no. 10 103016, [2103.09824].
  * [24] F. Calore, M. Cirelli, L. Derome, Y. Genolini, D. Maurin, et al., AMS-02 antiprotons and dark matter: Trimmed hints and robust bounds,  SciPost Phys. 12 (2022), no. 5 163, [2202.03076].
  * [25] A. Cuoco, J. Heisig, M. Korsmeier, and M. Krämer, Constraining heavy dark matter with cosmic-ray antiprotons,  JCAP 04 (2018) 004, [1711.05274].
  * [26] C. Arina, A. Beniwal, C. Degrande, J. Heisig, and A. Scaffidi, Global fit of pseudo-Nambu-Goldstone Dark Matter,  JHEP 04 (2020) 015, [1912.04008].
  * [27] F. Kahlhoefer, M. Korsmeier, M. Krämer, S. Manconi, and K. Nippel, Constraining dark matter annihilation with cosmic ray antiprotons using neural networks,  JCAP 12 (2021), no. 12 037, [2107.12395].
  * [28] C. Chang, P. Scott, T. E. Gonzalo, F. Kahlhoefer, and M. White, “Global fits of simplified models for dark matter with GAMBIT: II. Vector dark matter with an $s$-channel vector mediator.”  In preparation.
  * [29] GAMBIT Collaboration, V. Ananyev et al., “Collider constraints on electroweakinos in the presence of a light gravitino.”  In preparation.
  * [30] A. Reinert and M. W. Winkler, A Precision Search for WIMPs with Charged Cosmic Rays,  JCAP 01 (2018) 055, [1712.00002].
  * [31] A. Cuoco, J. Heisig, L. Klamt, M. Korsmeier, and M. Krämer, Scrutinizing the evidence for dark matter in cosmic-ray antiprotons,  Phys. Rev. D 99 (2019), no. 10 103014, [1903.01472].
  * [32] 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, [2005.04237].
  * [33] A. Cuoco, M. Krämer, and M. Korsmeier, Novel Dark Matter Constraints from Antiprotons in Light of AMS-02,  Phys. Rev. Lett. 118 (2017), no. 19 191102, [1610.03071].
  * [34] 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, [1610.03840].
  * [35] I. Cholis, T. Linden, and D. Hooper, A Robust Excess in the Cosmic-Ray Antiproton Spectrum: Implications for Annihilating Dark Matter,  Phys. Rev. D 99 (2019), no. 10 103026, [1903.02549].
  * [36] V. Silveira and A. Zee, SCALAR PHANTOMS,  Phys. Lett. B 161 (1985) 136–140.
  * [37] J. McDonald, Gauge singlet scalars as cold dark matter,  Phys. Rev. D 50 (1994) 3637–3649, [hep-ph/0702143].
  * [38] C. P. Burgess, M. Pospelov, and T. ter Veldhuis, The Minimal model of nonbaryonic dark matter: A Singlet scalar,  Nucl. Phys. B 619 (2001) 709–728, [hep-ph/0011335].
  * [39] J. M. Cline, K. Kainulainen, P. Scott, and C. Weniger, Update on scalar singlet dark matter,  Phys. Rev. D 88 (2013) 055025, [1306.4710]. [Erratum: Phys.Rev.D 92, 039906 (2015)].
  * [40] GAMBIT Collaboration, P. Athron et al., Status of the scalar singlet dark matter model,  Eur. Phys. J. C 77 (2017), no. 8 568, [1705.07931].
  * [41] P. Athron, J. M. Cornell, F. Kahlhoefer, J. Mckay, P. Scott, et al., Impact of vacuum stability, perturbativity and XENON1T on global fits of $\mathbb{Z}_{2}$ and $\mathbb{Z}_{3}$ scalar singlet dark matter,  Eur. Phys. J. C 78 (2018), no. 10 830, [1806.11281].
  * [42] CMS Collaboration, A. Tumasyan et al., Search for invisible decays of the Higgs boson produced via vector boson fusion in proton-proton collisions at s=13 TeV,  Phys. Rev. D 105 (2022), no. 9 092007, [2201.11585].
  * [43] ATLAS Collaboration, G. Aad et al., Search for invisible Higgs-boson decays in events with vector-boson fusion signatures using 139 fb-1 of proton-proton data recorded by the ATLAS experiment,  JHEP 08 (2022) 104, [2202.07953].
  * [44] LZ Collaboration, J. Aalbers et al., First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment, 2207.03764.
  * [45] PandaX-4T Collaboration, Y. Meng et al., Dark Matter Search Results from the PandaX-4T Commissioning Run,  Phys. Rev. Lett. 127 (2021), no. 26 261802, [2107.13438].
  * [46] A. W. Strong, I. V. Moskalenko, and V. S. Ptuskin, Cosmic-ray propagation and interactions in the Galaxy,  Ann. Rev. Nucl. Part. Sci. 57 (2007) 285–327, [astro-ph/0701517].
  * [47] A. W. Strong, I. V. Moskalenko, and O. Reimer, Diffuse continuum gamma-rays from the galaxy,  Astrophys. J. 537 (2000) 763–784, [astro-ph/9811296]. [Erratum: Astrophys. J.541,1109(2000)].
  * [48] C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, Cosmic-Ray Nuclei, Antiprotons and Gamma-rays in the Galaxy: a New Diffusion Model,  JCAP 10 (2008) 018, [0807.4730]. [Erratum: JCAP 04, E01 (2016)].
  * [49] C. Evoli, D. Gaggero, A. Vittino, M. Di Mauro, D. Grasso, et al., Cosmic-ray propagation with DRAGON2: II. Nuclear interactions with the interstellar gas,  JCAP 07 (2018) 006, [1711.09616].
  * [50] R. Kissmann, PICARD: A novel code for the Galactic Cosmic Ray propagation problem,  Astropart. Phys. 55 (2014) 37–50, [1401.4035].
  * [51] D. Maurin, USINE: semi-analytical models for Galactic cosmic-ray propagation,  Comput. Phys. Commun. 247 (2020) 106942, [1807.02968].
  * [52] S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, et al., The origin of Galactic cosmic rays: challenges to the standard paradigm,  Int. J. Mod. Phys. D 28 (2019), no. 15 1930022, [1903.11584].
  * [53] J. Park, D. Caprioli, and A. Spitkovsky, Simultaneous Acceleration of Protons and Electrons at Nonrelativistic Quasiparallel Collisionless Shocks,  Phys. Rev. Lett. 114 (2015), no. 8 085003, [1412.0672].
  * [54] S. Celli, G. Morlino, S. Gabici, and F. A. Aharonian, Exploring particle escape in supernova remnants through gamma rays,  Mon. Not. Roy. Astron. Soc. 490 (2019), no. 3 4317–4333, [1906.09454].
  * [55] G. Morlino and S. Celli, Cosmic ray electrons released by supernova remnants,  Mon. Not. Roy. Astron. Soc. 508 (2021), no. 4 6142–6154, [2106.06488].
  * [56] H. Jacobs, P. Mertsch, and V. H. M. Phan, Self-confinement of low-energy cosmic rays around supernova remnants,  JCAP 05 (2022), no. 05 024, [2112.09708].
  * [57] M. Korsmeier and A. Cuoco, Galactic cosmic-ray propagation in the light of AMS-02: Analysis of protons, helium, and antiprotons,  Phys. Rev. D 94 (2016), no. 12 123019, [1607.06093].
  * [58] M. J. Boschini et al., Solution of heliospheric propagation: unveiling the local interstellar spectra of cosmic ray species,  Astrophys. J. 840 (2017), no. 2 115, [1704.06337].
  * [59] M. Korsmeier and A. Cuoco, Testing the universality of cosmic-ray nuclei from protons to oxygen with AMS-02,  Phys. Rev. D 105 (2022), no. 10 103033, [2112.08381].
  * [60] P. Mertsch, A. Vittino, and S. Sarkar, Explaining cosmic ray antimatter with secondaries from old supernova remnants,  Phys. Rev. D 104 (2021), no. 10 103029, [2012.12853].
  * [61] P. De La Torre Luque, M. N. Mazziotta, F. Loparco, F. Gargano, and D. Serini, Implications of current nuclear cross sections on secondary cosmic rays with the upcoming DRAGON2 code,  JCAP 03 (2021) 099, [2101.01547].
  * [62] P. Salati, F. Donato, and N. Fornengo, Indirect dark matter detection with cosmic antimatter, in Particle Dark Matter : Observations, Models and Searches (G. Bertone, ed.), p. 521.  2010\.
  * [63] P. Blasi and P. D. Serpico, High-energy antiprotons from old supernova remnants,  Phys. Rev. Lett. 103 (2009) 081103, [0904.0871].
  * [64] P. Mertsch and S. Sarkar, AMS-02 data confront acceleration of cosmic ray secondaries in nearby sources,  Phys. Rev. D90 (2014) 061301, [1402.0855].
  * [65] K. Kohri, K. Ioka, Y. Fujita, and R. Yamazaki, Can we explain AMS-02 antiproton and positron excesses simultaneously by nearby supernovae without pulsars or dark matter?,  PTEP 2016 (2016), no. 2 021E01, [1505.01236].
  * [66] J. F. Navarro, C. S. Frenk, and S. D. M. White, The Structure of cold dark matter halos,  Astrophys. J. 462 (1996) 563–575, [astro-ph/9508025].
  * [67] P. F. de Salas and A. Widmark, Dark matter local density determination: recent observations and future prospects,  Rept. Prog. Phys. 84 (2021), no. 10 104901, [2012.11477].
  * [68] 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, [1012.4515]. [Erratum: JCAP1210,E01(2012)].
  * [69] P. Ciafaloni, D. Comelli, A. Riotto, F. Sala, A. Strumia, et al., Weak Corrections are Relevant for Dark Matter Indirect Detection,  JCAP 03 (2011) 019, [1009.0224].
  * [70] A. Cuoco, J. Heisig, M. Korsmeier, and M. Krämer, Probing dark matter annihilation in the Galaxy with antiprotons and gamma rays,  JCAP 1710 (2017), no. 10 053, [1704.08258].
  * [71] L. Derome, D. Maurin, P. Salati, M. Boudaud, Y. Génolini, et al., Fitting B/C cosmic-ray data in the AMS-02 era: A cookbook,  Astron. Astrophys. 627 (2019) A158, [1904.08210].
  * [72] M. Boudaud, Y. Génolini, L. Derome, J. Lavalle, D. Maurin, et al., AMS-02 antiprotons’ consistency with a secondary astrophysical origin,  Phys. Rev. Res. 2 (2020), no. 2 023022, [1906.07119].
  * [73] 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, [1802.03030].
  * [74] I. V. Moskalenko, A. W. Strong, J. F. Ormes, and M. S. Potgieter, Secondary anti-protons and propagation of cosmic rays in the galaxy and heliosphere,  Astrophys. J. 565 (2002) 280–296, [astro-ph/0106567].
  * [75] R. Aloisio and P. Blasi, Propagation of galactic cosmic rays in the presence of self-generated turbulence,  JCAP 07 (2013) 001, [1306.2018].
  * [76] J. A. Miller and D. A. Roberts, Stochastic Proton Acceleration by Cascading Alfven Waves in Impulsive Solar Flares,  ApJ 452 (Oct., 1995) 912.
  * [77] C. Evoli, P. Blasi, G. Morlino, and R. Aloisio, Origin of the Cosmic Ray Galactic Halo Driven by Advected Turbulence and Self-Generated Waves,  Phys. Rev. Lett. 121 (2018), no. 2 021102, [1806.04153].
  * [78] P. Blasi, E. Amato, and P. D. Serpico, Spectral breaks as a signature of cosmic ray induced turbulence in the Galaxy,  Phys. Rev. Lett. 109 (2012) 061101, [1207.3706].
  * [79] AMS Collaboration, M. Aguilar et al., Observation of New Properties of Secondary Cosmic Rays Lithium, Beryllium, and Boron by the Alpha Magnetic Spectrometer on the International Space Station,  Phys. Rev. Lett. 120 (2018), no. 2 021101.
  * [80] Y. Génolini et al., Indications for a high-rigidity break in the cosmic-ray diffusion coefficient,  Phys. Rev. Lett. 119 (2017), no. 24 241101, [1706.09812].
  * [81] M. Korsmeier and A. Cuoco, Testing the universality of cosmic-ray nuclei from protons to oxygen with AMS-02,  Phys. Rev. D 105 (2022), no. 10 103033, [2112.08381].
  * [82] E. S. Seo and V. S. Ptuskin, Stochastic Reacceleration of Cosmic Rays in the Interstellar Medium,  ApJ 431 (Aug., 1994) 705.
  * [83] A. W. Strong and I. V. Moskalenko, Propagation of cosmic-ray nucleons in the galaxy,  Astrophys. J. 509 (1998) 212–228, [astro-ph/9807150].
  * [84] M. Kuhlen and P. Mertsch, Time-dependent AMS-02 electron-positron fluxes in an extended force-field model,  Phys. Rev. Lett. 123 (2019), no. 25 251104, [1909.01154].
  * [85] R. Kappl, SOLARPROP: Charge-sign Dependent Solar Modulation for Everyone,  Comput. Phys. Commun. 207 (2016) 386–399, [1511.07875].
  * [86] M. J. Boschini, S. Della Torre, M. Gervasi, G. La Vacca, and P. G. Rancoita, Propagation of cosmic rays in heliosphere: The HELMOD model,  Adv. Space Res. 62 (2018) 2859–2879, [1704.03733].
  * [87] L. A. Fisk, Solar Modulation and a Galactic Origin for the Anomalous Component Observed in Low-Energy Cosmic Rays,  Astrophys. J. 206 (1976) 333–341.
  * [88] D. Maurin, F. Donato, R. Taillet, and P. Salati, Cosmic rays below z=30 in a diffusion model: new constraints on propagation parameters,  Astrophys. J. 555 (2001) 585–596, [astro-ph/0101231].
  * [89] G. Di Bernardo, C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, Unified interpretation of cosmic-ray nuclei and antiproton recent measurements,  Astropart. Phys. 34 (2010) 274–283, [0909.4548].
  * [90] G. Jóhannesson et al., Bayesian analysis of cosmic-ray propagation: evidence against homogeneous diffusion,  Astrophys. J. 824 (2016), no. 1 16, [1602.02243].
  * [91] M. Korsmeier and A. Cuoco, Galactic cosmic-ray propagation in the light of AMS-02: Analysis of protons, helium, and antiprotons,  Phys. Rev. D 94 (2016), no. 12 123019, [1607.06093].
  * [92] M. J. Boschini et al., Inference of the Local Interstellar Spectra of Cosmic-Ray Nuclei Z $\leq$ 28 with the GalProp–HelMod Framework,  Astrophys. J. Suppl. 250 (2020), no. 2 27, [2006.01337].
  * [93] P. D. L. T. Luque, M. N. Mazziotta, F. Loparco, F. Gargano, and D. Serini, Markov chain Monte Carlo analyses of the flux ratios of B, Be and Li with the DRAGON2 code,  JCAP 07 (2021) 010, [2102.13238].
  * [94] V. Ptuskin, V. Zirakashvili, and E.-S. Seo, Spectra of Cosmic Ray Protons and Helium Produced in Supernova Remnants,  Astrophys. J. 763 (2013) 47, [1212.0381].
  * [95] N. Tomassetti, Testing universality of cosmic-ray acceleration with proton/helium data from AMS and Voyager-1,  Adv. Space Res. 60 (2017) 815–825, [1610.06187].
  * [96] A. A. Lagutin, N. V. Volkov, R. I. Raikin, and A. G. Tyumentsev, Origin of hardening and universality of cosmic rays spectra in GV–PV rigidity region,  J. Phys. Conf. Ser. 1181 (2019), no. 1 012023, [1905.06699].
  * [97] D. Caprioli, D. T. Yi, and A. Spitkovsky, Chemical Enhancements in Shock-Accelerated Particles: Ab initio Simulations,  Phys. Rev. Lett. 119 (2017), no. 17 171101, [1704.08252].
  * [98] A. Hanusch, T. Liseykina, and M. Malkov, Acceleration of Cosmic Rays in Supernova Shocks: elemental selectivity of the injection mechanism,  Astrophys. J. 872 (2019), no. 1 108, [1803.00428].
  * [99] T. A. Porter, G. Johannesson, and I. V. Moskalenko, The GALPROP Cosmic-ray Propagation and Nonthermal Emissions Framework: Release v57,  Astrophys. J. Supp. 262 (2022), no. 1 30, [2112.12745].
  * [100] AMS Collaboration, M. Aguilar et al., Observation of Complex Time Structures in the Cosmic-Ray Electron and Positron Fluxes with the Alpha Magnetic Spectrometer on the International Space Station,  Phys. Rev. Lett. 121 (2018), no. 5 051102.
  * [101] A. Vittino, P. Mertsch, H. Gast, and S. Schael, Breaks in interstellar spectra of positrons and electrons derived from time-dependent AMS data,  Phys. Rev. D 100 (2019), no. 4 043007, [1904.05899].
  * [102] N. Weinrich, Y. Génolini, M. Boudaud, L. Derome, and D. Maurin, Combined analysis of AMS-02 (Li,Be,B)/C, N/O, 3He, and 4He data,  Astron. Astrophys. 639 (2020) A131, [2002.11406].
  * [103] V. S. Ptuskin, I. V. Moskalenko, F. C. Jones, A. W. Strong, and V. N. Zirakashvili, Dissipation of magnetohydrodynamic waves on energetic particles: impact on interstellar turbulence and cosmic ray transport,  Astrophys. J. 642 (2006) 902–916, [astro-ph/0510335].
  * [104] C. Evoli, G. Morlino, P. Blasi, and R. Aloisio, AMS-02 beryllium data and its implication for cosmic ray transport,  Phys. Rev. D 101 (2020), no. 2 023013, [1910.04113].
  * [105] N. Weinrich, M. Boudaud, L. Derome, Y. Genolini, J. Lavalle, et al., Galactic halo size in the light of recent AMS-02 data,  Astron. Astrophys. 639 (2020) A74, [2004.00441].
  * [106] D. Maurin, E. Ferronato Bueno, and L. Derome, A simple determination of the halo size from 10Be/9Be data,  Astron. Astrophys. 667 (2022) A25, [2203.07265].
  * [107] Y. Genolini, D. Maurin, I. V. Moskalenko, and M. Unger, Current status and desired precision of the isotopic production cross sections relevant to astrophysics of cosmic rays: Li, Be, B, C, and N,  Phys. Rev. C 98 (2018), no. 3 034611, [1803.04686].
  * [108] F. Feroz, M. P. Hobson, and M. Bridges, MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics,  Mon. Not. Roy. Astron. Soc. 398 (2009) 1601–1614, [0809.3437].
  * [109] AMS Collaboration, M. Aguilar et al., Properties of Cosmic Helium Isotopes Measured by the Alpha Magnetic Spectrometer,  Phys. Rev. Lett. 123 (2019), no. 18 181102.
  * [110] A. C. Cummings, E. C. Stone, B. C. Heikkila, N. Lal, W. R. Webber, et al., Galactic Cosmic Rays in the Local Interstellar Medium: Voyager 1 Observations and Model Results,  Astrophys. J. 831 (2016), no. 1 18.
  * [111] Planck Collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters,  Astron. Astrophys. 641 (2020) A6, [1807.06209]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [112] GAMBIT Collaboration, P. Athron et al., Supplementary Data: Status of the scalar singlet dark matter model (arXiv:1705.07931), .
  * [113] T. G. S. Workgroup, G. D. Martinez, J. McKay, B. Farmer, P. Scott, et al., Comparison of statistical sampling methods with ScannerBit, the GAMBIT scanning module,  The European Physical Journal C 77 (Nov., 2017).
  * [114] P. Scott, Pippi - painless parsing, post-processing and plotting of posterior and likelihood samples,  The European Physical Journal Plus 127 (Nov., 2012). arXiv:1206.2245.
  * [115] GAMBIT Collaboration, P. Athron et al., Global analyses of Higgs portal singlet dark matter models using GAMBIT,  Eur. Phys. J. C 79 (2019), no. 1 38, [1808.10465].
  * [116] J. Aalbers et al., A next-generation liquid xenon observatory for dark matter and neutrino physics,  J. Phys. G 50 (2023), no. 1 013001, [2203.02309].
  * [117] S. Bloor, T. E. Gonzalo, P. Scott, C. Chang, A. Raklev, et al., The GAMBIT Universal Model Machine: from Lagrangians to likelihoods,  Eur. Phys. J. C 81 (2021), no. 12 1103, [2107.00030].
  * [118] GAPS Collaboration, T. Aramaki, C. J. Hailey, S. E. Boggs, P. von Doetinchem, H. Fuke, et al., Antideuteron Sensitivity for the GAPS Experiment,  Astropart. Phys. 74 (2016) 6–13, [1506.02513].
  * [119] R. Battiston et al., High precision particle astrophysics as a new window on the universe with an Antimatter Large Acceptance Detector In Orbit (ALADInO),  Exper. Astron. 51 (2021), no. 5 1299–1330. [Erratum: Exper.Astron. 51, 1331–1332 (2021)].
  * [120] S. Schael et al., AMS-100: The next generation magnetic spectrometer in space – An international science platform for physics and astrophysics at Lagrange point 2,  Nucl. Instrum. Meth. A 944 (2019) 162561, [1907.04168].
  * [121] S. Hochreiter and J. Schmidhuber, Long Short-Term Memory,  Neural Computation 9 (11, 1997) 1735–1780, [https://direct.mit.edu/neco/article-pdf/9/8/1735/813796/neco.1997.9.8.1735.pdf].
  * [122] A. Belyaev, N. D. Christensen, and A. Pukhov, CalcHEP 3.4 for collider physics within and beyond the Standard Model,  Comput. Phys. Commun. 184 (2013) 1729–1769, [1207.6082].