#  CosmiXs: Cosmic messenger spectra for indirect dark matter searches

Chiara Arina  , Mattia Di Mauro† , Nicolao Fornengo  , Jan Heisig  ,
Adil Jueid†  and Roberto Ruiz de Austri†

###### Abstract

The energy spectra of particles produced from dark matter (DM) annihilation or
decay are one of the fundamental ingredients to calculate the predicted fluxes
of cosmic rays and radiation searched for in indirect DM detection. We revisit
the calculation of the source spectra for annihilating and decaying DM using
the Vincia shower algorithm in Pythia to include QED and QCD final state
radiation and diagrams for the EW corrections with massive bosons, not present
in the default Pythia shower model. We take into account the spin information
of the particles during the entire EW shower and the off-shell contributions
from massive gauge bosons. Furthermore, we perform a dedicated tuning of the
Vincia and Pythia parameters to LEP data on the production of pions, photons,
and hyperons at the $Z$ resonance and discuss the underlying uncertainties. To
enable the use of our results in DM studies, we provide the tabulated source
spectra for the most relevant cosmic messenger particles, namely antiprotons,
positrons, $\gamma$ rays and the three neutrino flavors, for all the fermionic
and bosonic channels and DM masses between 5 GeV and 100 TeV, on github.

††footnotetext: † Contact authors

##  1 Introduction

The particle origin of dark matter (DM) remains one of the most puzzling
mysteries in Physics. Different strategies are pursued to search for a
particle physics signal originated by DM interactions beyond the gravitational
one (see, e.g., [1] for a review). Indirect detection searches seek to find
excesses in the fluxes of cosmic messenger particles, like positrons
($e^{+}$), antiprotons ($\bar{p}$), $\gamma$ rays, neutrinos ($\nu$) and
antinuclei [2], that stem from DM annihilation into SM particles and their
subsequent decays in locally overdense regions like the Galactic center. They
constitute an important pillar for testing the mechanism behind DM genesis in
the early Universe because they directly probe the annihilating nature of DM
required by thermal freeze-out.

A theoretically well-motivated DM candidate that allows for thermal freeze-out
is a Weakly Interacting Massive Particle (WIMP). WIMPs emerge as new particles
in several Beyond the Standard Model (BSM) theories, for instance, in
Supersymmetry. A signal of cosmic messengers generated by WIMPs in the GeV-TeV
energy range can be detected by different operating experiments such as e.g.
AMS-02, Fermi-LAT, Imaging and Water Cherenkov detectors, SuperK and IceCube
(see e.g. [3]).

The theoretical calculation of the flux of messenger particles produced by DM
annihilation or decay is based on the energy spectra at source111The energy
spectra at source are the one produced directly from DM particles annihilating
or decaying at the astrophysical target of interest. Therefore, they are not
the spectra of particles that reach Earth, which instead need to account for
their propagation from the source, thus involving additional processes, like
oscillation in the case of neutrinos, absorption in the case of $\gamma$ rays
or energy redistribution in the case of charged particles. (hereafter simply
called spectra) of these particles, which are typically calculated using codes
for the generation of high-energy physics collision events also known as
multi-purpose Monte Carlo event generators. One of the most widely used
reference for the source DM spectra is Ref. [4] (hereafter PPPC). The latest
version of PPPC has been produced by using the Pythia event generator version
8.135 to calculate DM spectra for different annihilation channels and masses
from 5 GeV to 100 TeV. PPPC employed a process where DM creates a resonance,
with a center-of-mass energy equal to twice the DM mass, which subsequently
decays into a couple of SM particles. These results can be applied to a wide
range of models and have established themselves as a standard tool used by the
community for indirect DM searches.

For the most widely considered DM candidates, like the WIMP, messenger
particles are mainly produced by three processes: hadronization, leading order
electroweak (EW) processes and EW corrections. The hadronization is initiated
with the production of quarks and/or gluons, either produced directly or
through the decays of heavy resonances, which subsequently generate gluons and
other quarks. After a timescale of the order of a Fermi ($10^{-15}$ m), these
particles hadronize and produce mesons and baryons before decaying into stable
particles at cosmological scales (except, $\bar{p}$ which is stable).
Hadronization constitutes one of the main processes for producing $\bar{p}$,
$\gamma$ rays, through the production and subsequent decay of $\pi^{0}$, and
$\nu$ and $e^{\pm}$ coming from the decay of $\pi^{\pm}$, kaons and hyperons.
Leading order EW processes take into account all EW decays, such as the one of
$\tau$ and $\mu$ leptons or EW gauge bosons. These processes are the main
production mechanisms for electrons and positrons.

EW corrections include initial state radiation (ISR), final state radiations
(FSR) and internal bremsstrahlung (IB). When DM is a singlet under the SM
gauge symmetries, it does not directly emit radiation and thus ISR is not
considered. However, when DM particles couple to the SM gauge bosons, EW ISR
can provide a sizeable contribution to the final spectra, see e.g. [5, 6, 7,
8, 9, 10, 11]. In general, ISR exhibits a degree of model dependence. Since we
focus on a model-independent approach, we will not discuss ISR further in this
work. Instead, we will take into account the radiative emissions from FSR
processes, that can produce final state particles in terms of photons, gluons
and EW bosons. The Pythia standard shower algorithm by default includes FSR of
photons from fermions and gluons from quarks. It includes also the emission of
$W^{\pm}$ and $Z$ bosons, named as EW showers. This latter contribution was
not included in the Pythia version employed in PPPC. Therefore, the authors
have included EW showers semi-analytically to first order, following the
method of Ref. [12].

Figure 1: Feynmann diagram of an event where two DM particles ($\chi$)
annihilate into two $Z$ bosons. Then, the $Z$ particles emit EW radiation,
producing other bosons, which subsequently decay into fermions.

The approach taken in PPPC has represented the state of the art for many
years. However, several new features have been added in Pythia the meanwhile.
Indeed, in the latest Pythia versions, the emission of $W^{\pm}$ and $Z^{0}$
gauge bosons off fermions is an integrated part of the ISR and FSR, and is
fully interleaved with the QCD and QED emissions (see Ref. [13] for details).
More recently, an implementation of a new shower algorithm has been provided
in Vincia [14]. It uses a $p_{T}$-ordered model for QCD+QED/EW helicity-
dependent showers based on the Antenna formalism. Originally developed as an
independent plugin, it has been incorporated into the Pythia 8 source code
since Pythia version 8.3.

The main innovations of Vincia in relation to the standard Pythia shower
process are the following. First, Vincia takes into account the helicity along
all the shower, i.e. Vincia decomposes the shower into distinct terms for each
set of contributing helicities. Secondly, it includes the trilinear gauge
boson interactions $Z^{0}W^{+}W^{-}$, $HW^{+}W^{-}$, $HZZ$ and $\gamma
W^{+}W^{-}$, which are neglected in the standard Pythia shower and that can
give an important contribution at high DM masses (see Appendix A for more
details about the differences between the two shower options). As an example,
in Fig. 1, we show a representative Feynman diagram of an event included in
Vincia but not in standard Pythia. This event involves the production of two
$Z$-bosons from DM annihilation. The two $Z$-bosons produce $W$-bosons. Then
one of the four bosons subsequently radiates off a Higgs boson while the two
$W$-bosons in the lower branch decay into leptons and quarks, respectively.
Exemplarily, we also show gluon radiation off one of these quarks. In addition
to the previous innovations, the authors of Pythia have included new tuning of
the model parameters in the code, which makes the code more compatible than in
the past with collider data.

The main differences between PPPC and Vincia are the following. First, PPPC
includes EW corrections via semi-analytical calculations using the results of
Ref. [12]. These corrections were considered at leading order without
resummations. Second, these corrections are then matched to Pythia to add
further QCD emissions, to handle the decays of heavy resonances and to
hadronise colored final-state particles. Instead, in Vincia, both the EW
corrections and the other QED+QCD emissions are handled in a coherent manner.
Furthermore, the decay of heavy resonances such as the $W$-boson or the top
quark is fully interleaved with the rest of the shower machinery, a fact that
has non-trivial impact on the kinematical distributions of the decay products.
Finally, as we will demonstrate later, in the version of Pythia used in PPPC,
the photon yield for leptonic channels tends to be underestimated at low
energies.

Recently, the authors of Ref. [15] (hereafter HDMS 222The abbreviation derives
from their public repository called HDMSpectra.) have provided decay spectra
for DM with masses significantly above the scale of EW symmetry breaking. In
this approach, they evolve from the scale of the DM mass down to just above
the weak scale using the DGLAP equations by adopting an implementation that
considers all interactions in the unbroken SM phase, as well as a partial
treatment of soft-coherence effects. The handling of the decays of the $W/Z$
bosons and the top quark is done in the unbroken phase while the decay of the
SM Higgs boson is performed with Pythia. They perform a matching by evolving
across a parametrically small region through the weak scale, removing all
particles with EW scale masses. Finally, these results are matched at the EW
scale to Pythia, which is then used to calculate the subsequent showering,
hadronization, and light particle decays. The results of HDMS can lead to
theoretical problems in the matching between the physics generated by the
DGLAP formalism in the unbroken phase and the physics described by Pythia. Our
results are not affected by any of the issues described above because all the
processes happening after the DM annihilation or decay, including
hadronization, EW corrections and particle decays, are calculated internally
and consistently by Vincia. Additionally, our approach allows us to release
reliable results down to DM masses much lower than HDMS, which includes
spectra for values above 500 GeV only. Moreover, we use an improved version of
the Pythia code, tuning the most important parameters to LEP data. In Refs.
[16, 17, 18] (hereafter QCDUnc), some of the authors of this paper have
derived a new set of hadronization parameters (tuning) using Pythia 8.2
through a fit to LEP and SLD data at the $Z$-boson resonance. They have also
calculated a conservative set of uncertainties on the shower and hadronization
model parameters. They have estimated the impact on the energy spectra of
$\bar{p}$, $\gamma$ rays, $\nu$ and $e^{+}$, which is at most of the order of
$10\%$–$20\%$ in the peak region.

In this work, we generate state-of-the-art spectra for the following cosmic
particle messengers produced during DM annihilation and decay: $\bar{p}$,
$\gamma$ rays, $\nu$ (in their three flavor states) and $e^{+}$.333We verified
that the spectra of $e^{+}$ and $e^{-}$ are the same and therefore can be used
for those indirect searches that do not distinguish between electrons and
positrons. We improve the existing results on this topic adopting the
following strategy:

  * •

We use the latest version of the Vincia algorithm, implemented in the Pythia
version 8.309, which we interface with MadDM [19, 20, 21]. We use MadDM for
two reasons. First, it is convenient for having an event file containing all
(in particular the helicity) information to be passed to Vincia. Second, we
need a proper event generator for the 4-body processes as detailed in Sec. 3.1
since MadDM is the only event generator dedicated to DM annihilation. Thereby,
we include contributions from triple gauge boson interaction in the showering
and carry the helicity information throughout the showering process. We also
include the effect of running quark masses. Vincia represents the state-of-
the-art code for including EW corrections for masses up to hundreds of TeV.

  * •

We include annihilation channels with off-shell gauge bosons and new channels
such as $ZH$ and $Z\gamma$.

  * •

We carry out an improved tuning of the hadronization parameters in Pythia
(adopting the Vincia shower as our default) to fit the available data for the
production of particles from the $Z$ resonance measured at LEP. These
parameters are used to make predictions for DM annihilation thanks to the jet
universality.

Our results are publicly available in a github
repository.444https://github.com/ajueid/CosmiXs.git While we perform the
simulations for annihilating DM, upon simple rescaling, the results can be
used for decaying DM as well. Specifically, we cover the case of decaying
scalar, pseudoscalar555In the case of pseudoscalar DM the decay into $HH$ and
$W^{+}W^{-}$ and $ZZ$ are forbidden due to CP violation., vector and axial-
vector DM. A dedicated study of decaying fermionic DM is left for future work.

The remainder of the paper is organized as follows. In Sec. 2, we describe the
hadronization and EW processes that are responsible for the production of
particles relevant to DM studies. The main novelties of our analysis are
detailed in Sec. 3. In Sec. 4, we report the results of the tuning of the
Vincia shower algorithm and the consequent uncertainties in the final particle
spectra. Finally, in Sec. 5 we show the results for the particle spectra and
the comparison with other reference results, and we draw our conclusions in
Sec. 6.

##  2 Hadronization and Electroweak Model

In this section, we briefly discuss the physical description of the particle
production from DM. We begin with a discussion of the general features of
stable particle production from DM annihilation, which encompasses a complex
sequence of phenomena including resonance decay, QED and QCD bremsstrahlung,
EW showers, hadronization and hadron decays. We conclude this section with a
discussion of the composition of particle spectra, with some examples of final
state particles and annihilation channels.

###  2.1 Emission mechanism of particles from dark matter

Let us consider a generic annihilation process of DM particles $\chi$ into a
set of final-state particles:

$\displaystyle\chi\chi\to\underbrace{\bigg{[}X_{1}X_{2}\ldots
X_{N}\bigg{]}}_{\rm Intermediate~{}states}\to\overbrace{\bigg{(}Y_{11}\ldots
Y_{1a_{1}}\bigg{)}\ldots\bigg{(}Y_{N1}\ldots Y_{Na_{N}}\bigg{)}}^{\rm
Stable~{}particles}.$ (2.1)

In the case where the narrow-width approximation holds, we can factorise the
whole process into a production part and a decay part. In general, the first
part consists in the production of $N$ particles ($X_{1},\cdots X_{N}$), which
may be quarks, gluons, leptons or heavy resonances such as the $W/Z/H$ bosons
or the top quark. These particles undergo a series of complex processes that
give rise to various particles that are stable over astrophysical/cosmological
scales, such as photons, positrons, neutrinos or antiprotons. The narrow width
approximation holds true for particles with small decay widths, such as the
Higgs boson. However, off-shell effects provide important corrections for the
pair production of massive gauge bosons, especially below their production
threshold. We stress that the narrow width approximation is not valid in this
case since heavy particles produced either in the annihilation process or the
showering process do not decay until their virtuality reaches a scale that is
close to their offshellness scale. The physics modeling of stable particle
production depends on the nature of the intermediate particles $X_{i}$ and/or
their decay products ($Y_{ij}$).

QED bremsstrahlung occurs when the $X$ or the $Y$ particles are either
electrically charged or include photons. In this case, additional photons
and/or electrically charged particles are produced through $X_{i}\to
X_{i}\gamma$ (Fig. 2a) and $\gamma\to f\bar{f}$ (Fig. 2d). The photon emission
is enhanced for both soft and quasi-collinear regions.666For more details, see
Section 4 of Ref. [22]. Note that the collinear photons can have very high
energies ($E_{\gamma}\to M_{\chi}$) provided that the angle between the parent
particle and the photon is extremely small. On the other hand, fermion pair
production through $\gamma\to f\bar{f}$ can occur with subleading
probabilities but enhanced at low values of photon virtualities:
$(p_{f}+p_{\bar{f}})^{2}/M_{\chi}^{2}\to 0$. If the phase space is permitting,
fermions in general can emit $W/Z/H$ bosons (Fig. 2a–c). The inclusion of the
$W/Z$-boson emissions was included in Pythia since version 8.176. On the other
hand, due to the gauge structure of the SM, massive gauge bosons and the SM
Higgs boson can undergo further weak emissions (Fig. 2e–f). Note that not all
of these branchings are included in Pythia but Vincia includes all of them
through the Antenna formalism. In fact, more than 1000 Antenna functions for
weak showers are implemented in Vincia. In the following, the EW radiation
from bosons is called EWBR.

In the case where the $X$ or the $Y$ particles are colored particles, further
colored particles are produced through QCD bremsstrahlung (Fig. 2g-h). QCD
showers are treated in similar fashion as the case of QED showers, which
reflects the enhancement of probabilities for both soft and collinear
emissions ($X\to Xg$ and $g\to q\bar{q}$) in addition to the $g\to gg$
branching. The probabilities of the QCD shower branching are controlled by the
value of the strong coupling constant.

The decay of short-lived particles is an important source for stable particle
spectra, especially $e^{+}$ and $\bar{\nu}$ (Fig. 3). The contribution of
these sources to the particle spectra is dominant for the high-energy region.
We note that two annihilation/decay channels are very special. First the
$t\bar{t}$ channel is the only channel where we can have resonance decays,
QED+QCD showers, EW showers and hadronization. The $W^{+}W^{-}$ channel is
also very special as both electroweak/QED showers and resonant decays can
occur, producing quarks and leptons which then undergo further QED+QCD showers
before producing hadrons. The default treatment in Vincia is that resonance
decays are interleaved with QED/QCD and EW showers. More details about the EW
showers and the treatment of heavy-resonance decays in Vincia are shown in
Appendix A.

Figure 2: Example of Feynman diagrams for the different parton-shower
branchings. Here we show QED+EW emissions of gauge bosons off fermion lines
(a)–(b), Higgs emission off fermion lines (c), photon splitting into
$f\bar{f}$ (d), gauge boson and Higgs boson emissions off bosonic lines
(e)–(f) (labeled as EWBR), gluon emission off quark lines (g) and $g\to gg$
(h). Figure 3: Examples of Feynman diagrams for weak coupling driven decays of
fermions and gauge bosons. We show the leptonic decays of charged leptons (a),
the semi-leptonic decays of the tau lepton (b), the fermionic decays of the
massive gauge bosons (c) and the three-body decay of the top quark (d).

Any colored particle must be confined inside colour-neutral hadrons above
distance scales of order $10^{-15}$ m. This process which is called
hadronization cannot be modeled using first-principles QCD but only using
phenomenological models. There are two main models: string models [23, 24, 25,
26] and cluster models [27, 28] which are implemented in multipurpose Monte
Carlo event generators. Pythia 8 is based on the Lund string model where the
hadronization is modeled by a left-right symmetric fragmentation function
$f(z)$ given by

$f(z,m_{\perp h})\propto
N\frac{(1-z)^{a_{L}}}{z}\exp\left(\frac{-b_{L}m_{\perp h}^{2}}{z}\right),$
(2.2)

which gives the probability that a hadron $h$ gets a fraction $z\in[0,1]$ of
the remaining energy at each step of the hadronization process. In Eq.(2.2),
$N$ is a normalisation constant, $m_{\perp
h}^{2}\equiv\sqrt{m_{h}^{2}+p_{\perp h}^{2}}$ is the square of the transverse
mass of the hadron $h$, $a_{L}$ and $b_{L}$ are tunable parameters. This basic
picture of the hadronization of the $q\bar{q}$ system is unaffected in the
presence of gluons, since a gluon having a colour and an anti-colour structure
can be seen as a kink in the string. Eq. (2.2) can also be generalized to
include flavour effects: in particular, strange quarks and massive quarks at
the endpoints of the string. The inclusion of heavy mass effects is achieved
by the Bowler modification [29, 23]. Baryons composed of three quarks or three
antiquarks are produced similarly to mesons. The production of baryons can be
achieved by breaking the strings by the production of diquarks-antidiquarks
between the quark and the antiquark at the string endpoints. However, this
basic picture leads to a strong correlation between the produced baryon and
the antibaryon in both the flavour and the angular distributions. These
correlations have been falsified by experimental measurements of
$\Lambda^{0}$–$\bar{\Lambda}^{0}$ angular correlations by the Opal
collaboration [30]. To reduce the degree of the correlation between baryons,
the pop-corn mechanism was introduced [31, 32]. In this mechanism, one or more
$q\bar{q}$ pairs are produced in between the diquark-antidiquark pairs, which
enables the production of one or more mesons between the two baryons and
therefore decreases their correlation. With all these modifications, the
string fragmentation function is given by

$\displaystyle f(z)\propto N\frac{1}{z^{1+r_{Q}b_{L}m_{Q}}}z^{a_{\rm
eff}}\bigg{(}\frac{1-z}{z}\bigg{)}^{a_{\rm eff}}\exp\left(\frac{-b_{L}m_{\perp
h}^{2}}{z}\right),$ (2.3)

with $Q=c,b$ and $a_{\rm eff}\equiv a_{L}+a_{\rm QQ}$. The parameters of the
Lund string fragmentation function are given in Tab. 1. We must stress that
for the purpose of this study, we do not need to tune the Bowler parameters
$r_{c}$ and $r_{b}$ and we use their default values. The produced hadrons
within QCD jets decay into $\gamma$ (Fig. 4a), $\mu^{-}\to
e^{-}\bar{\nu}_{e}\nu_{\mu}$ (Fig. 4b), and $\bar{p}$ (Fig. 4c–d). The
contribution of these decays dominates in the peak and the bulk regions of the
spectra. For the $\bar{p}$ spectra, another important source comes from
hadronization, in which case these $\bar{p}$ are called primary (see Refs.
[17, 18] for more details). As mentioned above, the hadronization mechanism is
only solved by phenomenological models with many free parameters. Therefore,
the uncertainties can be estimated based on the parameters of the
hadronization model. The estimation of these uncertainties was done for the
first time in Refs. [16, 17, 18, 33]. These uncertainties were found to be of
the order of $10$–$30\%$ depending on the annihilation channel, the DM mass
and the final state particles. Furthermore, it was found that the impact of
these uncertainties on the best-fit point of the DM mass and the annihilation
cross-section can be dramatic, especially for heavy DM.

Figure 4: Examples of Feynman diagrams for the main hadron decays within QCD
jets. From top left to bottom right, we show the decay of
$\pi^{0}\to\gamma\gamma$ (a), $\pi^{-}\to\mu^{-}\bar{\nu}_{\mu}$ (b), $n\to
pe^{-}\bar{\nu}_{e}$ (c) and $\Lambda^{0}\to p\pi^{-}$ (d).
---
parameter |  Pythia 8 setting | Variation range | Vincia
$\sigma_{\perp}$ (GeV) | `StringPT:Sigma` | 0.0 – 1.0 | 0.305
$a_{L}$ | `StringZ:aLund` | 0.0 – 2.0 | 0.45
$b_{L}$ | `StringZ:bLund` | 0.2 – 2.0 | 0.80
$a_{QQ}$ | `StringZ:aExtraDiquark` | 0.0 – 2.0 | 0.90
$r_{c}$ | `StringZ:rFactC` | 0.0 – 2.0 | 0.85
$r_{b}$ | `StringZ:rFactB` | 0.0 – 2.0 | 1.15
Table 1:  The main parameters of the Lund fragmentation function in Pythia 8
along with their range and their default values when using the Vincia shower
algorithm.

###  2.2 Composition of particle spectra

In this section, we show how the different emission mechanisms explained in
the previous section compose into the DM spectra. We generate the spectra of
cosmic messengers for 64 masses between 5 GeV and 100 TeV777We plan to extend
the range towards larger masses in future analyses.. We simulate 5 million
annihilation events for each mass and produce the result in terms of
$dN/d\log_{10}(x)$, where $x\equiv E/M_{\chi}$. $\log_{10}(x)$ is provided for
100 logarithmic values between -8 and 0. The spectra are produced for the
following channels: $e^{+}e^{-}$, $\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$,
$\nu_{e}\nu_{e}$, $\nu_{\mu}\nu_{\mu}$, $\nu_{\tau}\nu_{\tau}$, $u\bar{u}$,
$d\bar{d}$, $c\bar{c}$, $s\bar{s}$, $t\bar{t}$, $b\bar{b}$, $\gamma\gamma$,
$gg$, $W^{+}W^{-}$, $ZZ$, $HH$, $ZH$, $\gamma Z$. In addition, for the
channels $e^{+}e^{-}$, $\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$ we also calculate
separately the left-handed and right-handed spectra ($e^{-}_{L}e^{+}_{L}$ and
$e^{-}_{R}e^{+}_{R}$) and for the gauge bosons $W^{+}W^{-}$, $ZZ$ the
longitudinal and transverse polarization spectra ($Z_{L}Z_{L}$ and
$Z_{T}Z_{T}$). We make our results publicly available as tables in a github
repository.4 In particular, Fig. 5 illustrates the spectrum of $\gamma$ rays,
positrons and antiprotons, for different annihilation channels. We focus in
the plot on high DM masses where EW corrections are relevant.

When the annihilation channel involves quarks, the main production of
particles comes from the hadronization of these quarks, which emit several
gluons for QCD bremsstrahlung and produce several particles in the final
states. We show in Fig. 5 the case of $M_{\chi}=100$ TeV annihilating into
$b\bar{b}$. In particular, the production of $\gamma$ rays ($e^{\pm}$ and
$\nu$) is mostly due to the hadronization process, which produces $\pi^{0}$
($\pi^{\pm}$) mesons that subsequently decay into two photons (in $e^{\pm}$
and $\nu$ by muon decays). The vast majority of $\bar{p}$ are produced by the
hadronization of $\bar{u}$ and $\bar{d}$ quarks. We can call this prompt or
primary production. However, there is also a relevant production of secondary
antiprotons, which are produced by the decay of resonances ($\Delta$-baryons),
from neutrons or from hyperons, such as $\Lambda$ and $\Sigma$. The
hadronization process accounts for at least $90\%$ of the total production of
particles. This process is labeled as “Hadronization” in Fig. 5. Other two
processes are responsible for the production of particles at a subleading
order. These are EWBR, for which we have the radiation of bosons from other
bosons (see Fig. 2 d-e), and QED FSR with radiation of a photon from a fermion
(see Fig. 2 a). We show in Fig. 6 the evolution of the contribution of EWBR
and FSR to the $\gamma$-ray total yield as a function of the DM mass. In
particular, the EWBR makes at most about $10\%$ of $\gamma$ rays at 100 TeV
and with a decreasing contribution going at lower DM masses. A similar result
is also obtained for $\bar{p}$, $\nu$ and $e^{+}$. QED FSR, instead,
contributes less than $1\%$ at all energies as exemplified in Fig. 5 for
$b\bar{b}$. Similar results are obtained for the other quarks and for the EW
and Higgs bosons.

Figure 5:  Spectra of $\gamma$ rays, positrons ($e^{+}$) and antiprotons
($\bar{p}$) for different annihilation channels, namely $e^{\pm}$,
$\mu^{\pm}$, $\tau^{\pm}$ and $b\bar{b}$. For each plot, we report the
different production mechanisms and the total spectrum in units of
$\log_{10}(x)$, where $x=E/M_{\chi}$. We display the product of hadronization
from quarks generated by DM annihilation (labeled as hadronization). The EWBR
is due to the boson radiation from other bosons. This process creates quarks
that hadronize into antiprotons and pions, which subsequently decay into
$\gamma$ rays or $e^{\pm}$. We display also the production of photons due to
QED FSR and the result of the hadronization products from quarks generated
through pair production from an energetic photon (labeled as
$\gamma\rightarrow q\bar{q}$). In case of $e^{\pm}$ channel we also show the
direct production of positrons from the initial particles and for the $\tau$
channel the production of its decay into $\gamma$ rays through its decay into
neutral pions.

In the case of leptonic channels, the contribution of EWBR and QED FSR becomes
much more relevant. This is particularly true for the production of $\gamma$
rays and antiprotons. Fig. 5 shows the spectra of $\gamma$ rays for the
$e^{\pm}$, $\mu^{\pm}$ and $\tau^{\pm}$ channels for a DM mass of 1 TeV. We
also show the production of positrons for the $e^{\pm}$ channel. This result
is representative of all the energies above a few hundreds of GeV where EWBR
becomes relevant. In particular, for $e^{\pm}$ and $\mu^{\pm}$ cases the
production of photons is dominated by FSR at all energies except for
$\log_{10}(x)\sim[-5,-1]$ where the EWBR gives the most important
contribution. In the case of the production of positrons, electrons and
neutrinos, EWBR makes the largest contribution at $\log_{10}(x)<-2$. Instead,
at higher energies, these particles are produced directly from the initial
leptons (referred to as “Prompt” in the figure). There is another important
process that produces some of these particles, namely the production of some
quarks from very energetic photons (see Fig.2 c). These quarks hadronize and
produce mesons and hadrons, which subsequently decay into $e^{\pm}$ and $\nu$.
This process, called $\gamma\rightarrow q\bar{q}$, accounts for at most $10\%$
of the total yield. Finally, antiprotons are mainly produced by the processes
EWBR and $\gamma\rightarrow q\bar{q}$. In the case of the $\tau^{+}\tau^{-}$
channel, there is another important process that produces $\gamma$ rays,
positrons and neutrinos, namely the decay of the $\tau$ lepton into charged
and neutral pions, which has a branching ratio of $65\%$. This process
contributes mainly at $\log_{10}(x)>-2$.

We show the variation of the contribution of EWBR to the $\gamma$-ray spectrum
with respect to the total spectrum as a function of DM mass in Fig. 6. For the
leptonic channels, FSR dominates the spectrum for the $\mu$ and $e$ channels,
while its contribution decreases to $20-30\%$ at 100 TeV. In the case of the
$\tau$ channel, on the other hand, the FSR yield remains roughly constant with
the DM mass, with a contribution between $15\%$ and $25\%$. In contrast, the
EWBR contribution becomes more important with increasing DM mass, reaching
between $60\%$ and $80\%$ for $M_{\chi}=100$ TeV.

Figure 6:  Fraction between the number of photons produced through the EWBR
(dashed curves) and QED FSR (dotted curves) and the total yield. We show the
results for the leptonic and $b\bar{b}$ channels with different colors.

##  3 Main novelties of our analysis

In this section, we elaborate on the novelties of our model and assumptions
with respect to the state-of-the-art literature.

###  3.1 Polarization and off-shell contribution

The method used in PPPC and QCDUnc to generate the DM spectra is based on
Pythia and the _resonance approach_. This method assumes that a spinless
resonance with a center-of-mass energy of twice the DM mass is produced and
subsequently decays into the specified final state (i.e. annihilation)
channel.888An example of this method is present in the example main07.cc of
the Pythia code. For example, if the annihilation channel is $b\bar{b}$,
Pythia simulates an $e^{+}e^{-}$ collision of energy $E_{\rm beam}=M_{\chi}$,
which generates a resonance ${\cal R}$ with energy of $2M_{\chi}$. Then,
$\mathcal{R}$ decays into $b\bar{b}$ with an isotropic distribution emission
with respect to the polar angle, i.e. the angle of the incoming $e^{\pm}$.
After being produced, the pair of the SM particles, i.e. $b\bar{b}$, undergo
the QED and QCD shower that produce the final particles.

The resonance approach implemented in Pythia and used in PPPC and QCDUnc does
not take into account some important aspects of the particle physics theory of
DM annihilation and particle production:

  * •

Off shell contribution: The contribution for the annihilation into a pair of
massive gauge bosons, $VV$ ($V={W^{\pm},Z}$), is zero if the DM mass is
smaller than the boson’s pole mass, as the resonance is restricted to decay
into a pair of on-shell particles only. However, gauge bosons can be produced
off-shell (see, e.g., [34]) i.e. with an effective mass different from its
pole mass, since their decay is modeled with a Breit-Wigner function with a
width $\Gamma$ that is 2.09 (2.50) GeV for W${}^{\pm}(Z)$ (see [35] for a
recent review).

  * •

Polarization of the gauge bosons: The resonance approach does not take into
account the spin of the massive gauge bosons produced by the DM annihilation.
Instead, the method works as if the annihilation channel is through the Higgs
boson, which is a scalar. Massive gauge bosons are spin-1 particles, and they
carry three possible polarizations (two transverse states ($\pm 1$) and one
longitudinal state). The boson polarization affects the kinematics (angular
and momentum distribution) of the fermions produced after their decay [36].
This effect cannot be taken into account with the generic resonance from which
all the particles produced are unpolarized.

The above two problems can be solved by using an external code that calculates
the matrix elements in a given BSM model. In this way, the spin of the
particles is included in the hard process and the kinematics of the particles
are generated according to the polarization of the bosons and their off-shell
contribution. We use MadGraph_aMC@NLO [37] and its wrapper MadDM [19, 20, 21].
Specifically, we use the LTS version of MadGraph_aMC@NLO (currently 2.9.16)
while for MadDM we use a custom version that incorporates Vincia and our new
tuning. This version of MadDM is not yet public, but will be made available to
the reader upon request. In Appendix B, we describe the MadDM commands we use
to generate the spectra.

For annihilation into off-shell vector bosons, we generate the full four-
fermion process with MadDM. Here, we consider annihilation via a $s$-channel
scalar mediator, specifically the SM Higgs, $\chi\chi\rightarrow H\rightarrow
VV$. To this end, we employ the Singlet Scalar model with a Higgs portal (see,
e.g., [38]). We generate the off-shell process by including the $V$ decays in
the hard process:

$\chi\chi\rightarrow H\rightarrow VV\rightarrow 4f,$ (3.1)

where $4f$ indicates the production of 4 fermions from the $V$ decay. We will
call the diagrams with four fermions in the final state of the hard process
the four-body diagram. In this process, the gauge bosons are virtual
particles. Therefore, they can have a mass below the pole, meaning that the
channel is kinematically open also for $M_{\chi}<M_{V}$. We consider the four-
body process up to 100 GeV. For higher masses, off-shell effects are fully
negligible. In fact, for a mass of 100 GeV, we have checked that the spectra
of the four-body processes agree with the ones of the on-shell production
processes, $\chi\chi\rightarrow VV$, within Monte Carlo uncertainties. The
method of using the four-body diagrams calculated by MadDM together with the
Vincia shower algorithm provides a consistent framework to take into account
the off-shell contribution and the helicity information of the massive bosons.

We now turn to the discussion of polarization of massive gauge bosons. As we
use MadDM for the generation of the hard process, the spin information of all
final state particles are stored in the Les Houches Event File (LHEF) [39].
Then the LHEF is passed onto Pythia which produces the showering taking into
account the polarization information of the massive bosons. In Fig. 7, we show
the effect of the polarization for a DM mass of 1 TeV. We consider the
$W^{+}W^{-}$ annihilation channel and calculate the spectrum for the
production of $\gamma$ rays, $\bar{p}$ and $e^{+}$. The figure reveals that
the largest differences in the spectra occur in the regions below and above
the peak of the distribution. This is true for all particles messengers. For
$\gamma$ rays and $e^{+}$ spectra the difference can be up to $20$–$30\%$,
while for $\bar{p}$ the deviation can also reach $50\%$ but at the tails of
the distribution. At the peak of the spectra, the results obtained with and
without taking boson polarisation into account are very similar in the range
of a few $\%$. The effect of the polarization of the gauge bosons becomes more
and more relevant the larger the DM mass is.

The difference in the spectra between the case with and without polarization
information is due to the different kinematics of the quarks and leptons
produced in the decay of the $W^{\pm}$ in the two cases. To check this, we
produced the spectra of the quarks for the four-body diagram with the Higgs
bosons channel $\chi\chi\rightarrow HH\rightarrow 4f$ (named as 4-body, $HH$),
the two-body case for the $W$ channel $\chi\chi\rightarrow W^{+}W^{-}$ and
with the default spin-0 resonance method of Pythia (spin0, res). Again, we
exemplarily consider a DM mass of 1 TeV. All three methods lead to the same
result for the quark spectra: the energy distribution is flat, i.e. the
probability of producing a quark with any energy between 0 and the DM mass is
the same. In all three cases, the polarization of the generated boson is not
taken into account. Instead, in the four-body process with the $W$ channel,
there is a peak at about half the DM mass and the distribution is not uniform.
This effect is due to the $V-A$ structure of the charge and neutral current
interactions in the SM, which generates asymmetries, called the forward-
backward and left-right asymmetries. The forward-backward asymmetry has been
studied in $e^{+}e^{-}$ collisions at the $Z$ resonance, which subsequently
decays into a pair of fermions (see, e.g., [40]). The fermions are not
produced symmetrically with respect to the polar angle, which is the angle
relative to the electron beam. This asymmetry is typically 0 for
$\sqrt{s}=m_{Z}$. The left-right asymmetry instead requires polarized beams
and is related to the asymmetry in the cross-section for the production of
fermions with the two chiralities.999In some studies related to the top quark,
it was found that the energy spectra of the charged leptons was strongly
correlated to their angular distributions. The latter is a direct probe of the
spin of the top quark (see for example Refs. [41, 42, 43] for more details).

Figure 7:  Comparison between the spectra generated with the Pythia resonance
method (“Pythia, res", dotted curves) and with MadDM using the Singlet scalar
model (dashed curves) for the $W^{+}W^{-}$ channel. We show the spectra of
$\gamma$ rays (green), $\bar{p}$ (blue) and positrons (red).

Finally, we want to mention that for large DM masses, higher order corrections
in the hard process can become significant. When including these
contributions, a matching scheme between the radiation from the fixed-order
calculation and showering algorithm has to be employed. Note that radiation in
the hard process is model-dependent as, in addition to FSR, it can involve IB
and ISR diagrams. In general, in a fixed-order calculation the inclusion of
all diagrams of a given order in perturbation theory is vital to ensure gauge
invariance and, hence, physically meaningful results. In this work, we refrain
from including higher-order corrections in the hard process, to support
applicability of the generated spectra to a wide range of models. Note that
for intermediate mass scales, the effect on the resulting spectra is small.
Considering, for instance, the Singlet Scalar model and a DM mass of 1 TeV,
the spectra from the four-body processes $\chi\chi\rightarrow VV\rightarrow
4f,\;\chi\chi\rightarrow VV\rightarrow 2f\,2V,\;\chi\chi\rightarrow
VV\rightarrow 4V$ and the two-body process $\chi\chi\rightarrow VV$ (with
showering performed by Vincia) agree within less than $5\%$ in the relevant
energy range.

###  3.2 Effect due to the choice of the BSM model

The resonance method implemented in Pythia is a widely used method in the
literature that accounts for the production of a pair of SM particles from the
decay of a spin-0 resonance. Other physical assumptions leading to the
annihilation channel and associated with a specific BSM model could produce
important differences in the final DM spectra. In particular, the choice of a
DM model plays an important role in the annihilation channels with very
relevant EW corrections. For example, DM simplified models [44] have couplings
between new mediators and fermions that are similar to the $V-A$ structure of
the EW processes in the SM. The exact values of the $V-A$ coupling parameters
$g^{V}$ and $g^{A}$ can produce a preference for the production of certain
helicity states that have different particle spectra [4].

As an example, we consider a DM simplified model with a mediator $Y_{1}$ of
spin-1 and a Dirac fermion $\chi$ as DM particle. The additional Lagrangian
added to the SM one is given by:

$\bar{\chi}\gamma_{\mu}(g_{\chi}^{V}+\gamma_{5}g_{\chi}^{A})\chi
Y_{1}^{\mu}+\sum^{N_{f}}_{i}\bar{F}_{i}\gamma_{\mu}(g^{V}_{ij}+\gamma_{5}g^{A}_{ij})F_{j}Y^{\mu}_{1},$
(3.2)

where $F_{i}$ is the fermionic multiplet and $g^{V}_{ij}$ and $g^{A}_{ij}$ are
the vector and axial couplings between $F_{i}$ and $Y^{\mu}_{1}$, which
resemble the $V-A$ structure of the neutral current in the SM. Therefore,
there are coupling parameters $g^{V}_{ij}$ and $g^{A}_{ij}$ for each of the
fermionic states written above. Typically, these models are taken with only
the vector structure, i.e. by assuming that $g^{A}=0$ [45] (see for caveats
e.g. here [46, 47, 48]). This implies that the annihilation of DM particles
produce both left-handed and right-handed fermions, i.e.
$\sum^{N_{F}}_{i}(g^{V}_{L,ij}\bar{F}_{L,i}\gamma_{\mu}F_{L,j}+g^{V}_{R,ij}\bar{F}_{R,i}\gamma_{\mu}F_{R,j})Y^{\mu}_{1}$,
where $F_{L,i}$ and $F_{R,i}$ are the left-handed and right-handed components
of the fermions.

However, this can change drastically if different assumptions are made for the
values of $g^{A}$ and $g^{V}$. Lets for example assume that $g^{A}=-g^{V}$. In
this case the Lagrangian associated to the production of fermions is given by:

$\sum^{N_{f}}_{i}g^{V}_{ij}\bar{F}_{i}\gamma_{\mu}(1-\gamma_{5})F_{j}Y^{\mu}_{1}=\sum^{N_{f}}_{i}g^{V}_{ij}\bar{F}_{L,i}\gamma_{\mu}F_{L,j}Y^{\mu}_{1}.$
(3.3)

This implies that, as in SM charged current interactions between the $W$ boson
and the fermions, the mediator $Y_{1}$ couples only to the left-handed
fermions and there is no coupling with the right-handed fermions. The model
thus produces completely polarized fermions which have left-handed chirality
and negative helicity. The opposite applies for antifermions. Instead, if
$g^{A}=g^{V}$ the mediator $Y_{1}$ couples only to right-handed fermions and
there is no coupling with left-handed fermions. The model thus produces
completely polarized fermions, which have right-handed chirality and positive
helicity. The opposite applies for antifermions.

Instead, in the case of scalar and pseudoscalar couplings between the mediator
and SM fermions, the messenger particle spectra do not change when different
values of the couplings are considered. This can be demonstrated considering
the DM simplified model with a mediator $Y_{0}$ of spin-0 and a Dirac fermion
$\chi$ as DM particle. This model is described by the following interaction
term in the Lagrangian:

$\bar{\chi}(g_{\chi}^{S}+i\gamma_{5}g_{\chi}^{P})\chi
Y_{0}+\sum^{N_{f}}_{i}\bar{F}_{i}(g^{S}_{ij}+i\gamma_{5}g^{P}_{ij})F_{j}Y_{0},$
(3.4)

where $g^{S}_{ij}$ and $g^{P}_{ij}$ are the scalar and pseudoscalar couplings,
respectively, between $F_{i}$ and $Y_{0}$. Regardless the values of the scalar
and pseudoscalar couplings ($g^{S}_{ij}=0$ or $g^{P}_{ij}=0$ or
$g^{P}_{ij}=g^{S}_{ij}$ or $g^{P}_{ij}=-g^{S}_{ij}$), this model does not
select any specific chirality of the fermions.

Note that the same conclusions hold for decaying DM considering its
interaction with SM fermions to be described by the second terms in Eqs. (3.2)
and (3.4), respectively, for spin-1 and spin-0 DM. In this case the mediator
is replaced by the respective DM particle.

As explained in Sec. 2.2, the EW corrections with the emission of a $W^{\pm}$
gauge boson can be a very important production mechanisms in DM spectra. When
we consider models with mediators of spin 0 or mediators spin 1 with $g^{V}$
or $g^{A}=0$, the fermions produced from the mediator contains both left and
right helicities with the same probability. Therefore, the production of
$W^{\pm}$ and $Z$ is turned on. When, $g^{A}=g^{V}$ the mediator produces only
right-handed fermions which do not produce $W^{\pm}$ for EWBR, instead $Z$ for
EWBR are still produced with the same rate as in the previous case. Finally,
when $g^{A}=-g^{V}$ the mediator produces only left-handed fermions which
produce both $W^{\pm}$ and $Z$ for EWBR. Since in this case fermions have only
the left-handed helicity the production of $W^{\pm}$ is enhanced by a factor
of two with respect to the case with spin 0, $g^{V}$ or $g^{A}=0$, which
contain both helicity states.

In Fig. 8, we show the difference between the DM spectra obtained for
different choices of coupling parameters for the $e^{+}e^{-}$,
$\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$ and $b\bar{b}$ annihilation channels. We
tested the following cases

  * •

SHP model, spectra for both helicity states (labeled as Spin 0).

  * •

DMsimp model spin 1 ($g^{V}=0$ or $g^{A}=0$), spectra for both helicity states
(Spin 1, $F_{L}+F_{R}$).

  * •

DMsimp model spin 1 ($g^{V}=-g^{A}$), spectra for left-handed fermions (Spin
1, $F_{L}$).

  * •

DMsimp model spin 1 ($g^{V}=g^{A}$), spectra for right-handed fermions (Spin
1, $F_{R}$).

As expected, the results obtained with the spin 0 case are the same as the
spin 1 case, $F_{L}+F_{R}$, , within the statistical errors. Therefore, in
Fig. 8 we decide not to show the Spin 0 spectra. In case of the $b\bar{b}$
channel, and for all hadronic channels, the differences between the tested
cases are at most at the level of $5-10\%$. This is due to the fact that the
main production process of particles is hadronization generated from pairs of
quarks or gauge bosons produced from DM annihilation (see, Sec. 2.2). Instead,
the EWBR for hadronic channels contribute, through a secondary hadronization
process, at most with $10\%$ of the total yield for $\gamma$ rays (see, Fig.
6).

For leptonic channels, the EWBR can provide a very important contribution and
thus the choice of the coupling parameters of the BSM model can have a much
larger effect with respect to what found for the hadronic channels. In case of
leptonic channels, for $g^{A}=g^{V}$, which selects only right-handed
fermions, is the one with the lowest spectra at the peak. This is due to the
fact that the mediator produces only right-handed fermions, which cannot
couple with $W^{\pm}$. Therefore, the EWBR with $W^{\pm}$ is not present and
does not contribute to the spectrum, giving thus a reduced yield of final
particles. Instead, the case with Spin 0 or Spin 1 with either $g^{A}=0$ or
$g^{V}=0$ contains fermions both right and left-handed. Therefore, when the
left-handed fermions are produced, which happens on average for half of the
events, the coupling with $W^{\pm}$ is present and the EWBR with these bosons
contribute to the spectrum. Finally, when $g^{A}=-g^{V}$ only left-handed
fermions are produced, so the production of $W^{\pm}$ is twice as large as the
Spin 0 case and the spectrum is much larger at the peak. Since the variation
in the spectra between the tested cases is due to the contribution of EWBR the
differences in the spectra are present at values of $\log_{10}(x)=[-5,-1]$
where this process gives the largest contribution (see, Sec. 2.2).

PPPC provides the spectra for the individual helicity states. The authors have
published this result only for leptons since for quarks, the difference is
minimal. We compare our results for left and right-handed helicity states
spectra with the PPPC spectra in Fig. 8. We obtain very similar results for
the $\mu^{+}\mu^{-}$ and $\tau^{+}\tau^{-}$ at energies $\log_{10}(x)>-5$.
Instead, at smaller energies, our results are systematically larger. We will
discuss more extensively this in Sec. 5. Instead, for the $e^{\pm}$ channel
the differences are not only at low energies but also at the peak of the
distribution where the contribution of the EWBR is the largest.

Figure 8:  Effect of the mediator spin and coupling parameters when
considering the DM simplified model in the spectra of $\gamma$ rays. We show
the results obtained for the $e^{-}e^{-}$, $\mu^{+}\mu^{-}$,
$\tau^{+}\tau^{-}$ and $b\bar{b}$ annihilation channels (from top left to
bottom right). In each figure we report the results obtained when using
$g^{A}=0$ which correspond to the case for which we have both left and right-
handed helicity states (Spin 1, $F_{L}+F_{R}$), $g^{A}=g^{V}$
($g^{A}=-g^{V}$), which corresponds to the case where only right-handed (left-
handed) fermions are produced, Spin 1, $F_{R}$ (Spin 1, $F_{L}$ ). Results
obtained with PPPC are shown for comparison.

###  3.3 Improved shower algorithm

The final very important improvement we introduce in our paper is the fact
that we use the Vincia shower algorithm and a specific tuning of the model
parameters. We already discussed this in Sec. 2 so here we will report the key
points. The first advantage of using Vincia is that it takes into account the
helicity of the particles during all the shower algorithm. The standard Pythia
shower algorithm does not take into account the helicity, and thus it misses
important effects as the ones we showed in the previous section, which are due
to the $V-A$ structure of the SM EW interactions. The second improvement is
related to the fact that Vincia takes into account the trilinear diagrams of
the massive bosons, which starts to be important at DM masses above 1 TeV.
These diagrams are not taken into account in the standard Pythia shower
algorithm. The third important new aspect is that we perform a tuning of the
main Vincia model parameter in order to match the LEP data for the production
of particles at the $Z$ boson resonance, which we expect to be similar with
the DM annihilation process. This is a very significant improvement with
respect to, PPPC which was generated with Pythia version 8.135. This version
of Pythia, using the default parameters, did not have a specific tuning to
collider data. Using the latest Pythia version has remarkable effects, as we
will discuss later in the paper. The most important one is for the photons
produced for QED FSR, which was highly underestimated in PPPC at low energies.

##  4 Tuning of the hadronization model

In this section, we discuss the tuning of the Lund hadronization model
parameters to a set of measurements performed by the experimental
collaborations of LEP at the $Z$–boson pole. We first discuss the technical
setup we have adopted in the fits and then present our results.

---
Parameter | Monash |  Vincia (default) |  Pythia [17, 18] | This work
$a_{L}$ | $0.68$ | $0.45$ | $0.601$ | $0.337\pm 0.015$
$b_{L}$ | $0.98$ | $0.80$ | $0.897$ | $0.784\pm 0.020$
$\sigma_{\perp}~{}({\rm GeV})$ | $0.335$ | $0.305$ | $0.307$ | $0.296\pm 0.003$
$a_{QQ}$ | $0.97$ | $0.90$ | $1.671$ | $1.246\pm 0.082$
$\chi^{2}/N_{\rm df}$ | $1034.52/852$ | $786.11/852$ | $676.69/852$ | $660.21/852$
Table 2:  Tuning results of the parameters of the hadronization model in
Pythia 8 using the Vincia shower plugin. For comparison, we show the results
in the baseline Monash tuning in both Pythia and Vincia as well as the results
of Refs. [17, 18]. Figure 9: Tuning results projected on the different
parameters and using the same measurements as in Refs. [17, 18]. $68\%$,
$95\%$ and $99.5\%$ confidence level contours are shown in turquoise, blue,
and dark blue respectively.

###  4.1 Technical setup

We use Pythia version 8.309 [49] to generate samples at the particle level and
with Vincia shower plugin as the default. The Monash tuning is used as our
baseline for further tuning [50]. We use Rivet version 3.1.7 for the
implementation of the different measurements at the particle level [51].
Professor version 2.3.3 is used to perform the tuning of the parameters [52].
In this work, we only tune four parameters which are $a_{L}$, $b_{L}$, $a_{\rm
QQ}$ and $\sigma_{\perp}$. The default values of these parameters in Vincia
along with their allowed ranges are given in Tab. 1. The Professor toolkit
uses analytical expressions that model the physical dependence of the
observables on the different parameters. This dependence is derived by fitting
the Monte Carlo predictions to a set of points in the four-dimensional
parameter space. Then the best-fit points are derived by a standard $\chi^{2}$
minimisation using Minuit [53] and which is implemented in Professor as the
default option. To assess the quality of the tuning we estimate the goodness-
of-fit (GoF) defined as

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

where ${\rm Data}_{(b)}$ is the experimental value of the observable ${\cal
O}$ at a bin $b$, $f_{(b)}(\\{p_{k}\\})$ is the value of the analytical
expression of the response function that models the theory prediction and is
cast as a polynomial of the parameters (defined in equation 4.4)
$\\{p_{k}\\}\equiv\\{a_{L},b_{L},a_{\rm QQ},\sigma_{\perp}\\}$ and
$\Delta_{b}$ is the total error on the observable ${\cal O}$ at a bin $b$.
There are three types of errors on the observable ${\cal O}$ at bin $b$:
experimental errors on the measurements, MC errors due to the limitation of
the size of the MC samples and the theory errors. Given that we have simulated
2 million events for each point in the parameter space, we find that the MC
errors are the smallest contribution to the error budget. In our analysis, we
do not assume any correlation between the different errors, as this
information was not provided by the experimental collaborations. On the other
hand, following Ref. [50] we add a flat $5\%$ uncertainty for each bin so that
we can avoid overfitting effects and as a sanity limit for the accuracy in
theory predictions for both the perturbative and the non-perturbative effects.
The total error is thus given by

$\displaystyle\Delta_{b}=\sqrt{\sigma_{b,\rm exp}^{2}+\sigma_{b,{\rm
MC}}^{2}+[0.05\times f_{(b)}(\\{p_{i}\\})]^{2}},$ (4.2)

The number of degrees-of-freedom ($N_{\rm df}$) is defined as the total number
of bins minus the number of independent parameters

$\displaystyle N_{\rm df}=\sum_{\mathcal{O}}|b\in\mathcal{O}|-N_{\rm params}.$
(4.3)

A good fit implies that $\chi^{2}/N_{\rm df}$ per number of degrees-of-freedom
to be $\approx 1$. The polynomial dependence of the true MC response is cast
as a fourth-order polynomial:

$\displaystyle
f_{(b)}(\\{p_{i}\\})=\alpha_{0}^{(b)}+\sum_{i=1}^{4}\beta_{i}^{(b)}p_{i}+\sum_{i,j=1}^{4}\gamma_{ij}^{(b)}p_{i}p_{j}+\sum_{i,j,k=1}^{4}\delta_{ijk}^{(b)}p_{i}p_{j}p_{k}+\sum_{i,j,k,\ell=1}^{4}\epsilon_{ijk\ell}^{(b)}p_{i}p_{j}p_{k}p_{\ell},$
(4.4)

with $\alpha,\beta,\gamma,~{}{\rm and}~{}\delta$ being the polynomial
coefficients determined in the fit and $\\{p_{i}\\}$ are the parameters of the
Lund hadronization model. To compute the response function, we have randomly
generated 500 MC samples that correspond to the four-dimensional parameter
space. The order of the polynomial function plays a crucial role in both the
quality of the fits and the consistency of the interpolated results with the
true MC response at the minimum of the model parameters. Performance of the
interpolation procedure is determined through the estimate of the residuals,
which are defined as

$\displaystyle{\rm Residuals}\equiv\sum_{i}\frac{f_{i}(\\{p_{k}\\})-{\rm
MC}_{i}}{{\rm MC}_{i}},$ (4.5)

where the sum runs over all the bins and for all the measurements and ${\rm
MC}_{i}$ represents the true MC response. We have checked that a four-order
polynomial response function is good enough to model the true response, since
that $95\%$ of the residual distribution is within $0.02$.

---
Measurement | Experiment |  $\chi^{2}/N_{\rm bins}$ | Measurement | Experiment | $\chi^{2}/N_{\rm bins}$
$1-T$ |  Aleph [54] | 0.13  |  $C$–parameter |  Aleph [54] | 0.39
$\log(1/x_{p})$ |  Aleph [54] | 0.19  | $\langle N_{\rm ch}\rangle$ |  Aleph [54] | 0.028
$\langle N_{\rm ch}\rangle~{}(|Y|<0.5)$ |  Aleph [54] | 0.012  | $\langle N_{\rm ch}\rangle~{}(|Y|<1.0)$ |  Aleph [54] | 0.028
$\langle N_{\rm ch}\rangle~{}(|Y|<1.5)$ |  Aleph [54] | 0.030  | $\langle N_{\rm ch}\rangle~{}(|Y|<2.0)$ |  Aleph [54] | 0.040
$\pi^{\pm}$ spectrum |  Aleph [54] | 0.67  |  $\pi^{0}$ spectrum |  Aleph [54] | 0.24
$\Lambda^{0}$ spectrum |  Aleph [55] | 1.24  |  $\Lambda^{0}$ spectrum ($2$–jet events) |  Aleph [55] | 1.31
Thrust |  Aleph [56] | 0.097  |  $C$–parameter |  Aleph [56] | 0.35
$N_{\rm ch}$ ($y_{\rm cut}=0.01$) |  Delphi [57] | 5.99  |  $N_{\rm ch}$ ($y_{\rm cut}=0.02$) |  Delphi [57] | 4.88
$\Lambda^{0}$ spectrum |  Delphi [58] | 1.34  | $\langle N_{\Lambda^{0}}\rangle$ |  Delphi [58] | 0.53
$\pi^{0}$ momentum |  Delphi [59] | 0.41  | $\log(1/x_{p})$ |  Delphi [59] | 0.33
$1-T$ |  Delphi [59] | 0.18  |  $C$–parameter |  Delphi [59] | 0.34
$\langle N_{\rm ch}\rangle$ |  Delphi [59] | 0.031  | $\langle N_{\pi^{\pm}}\rangle$ |  Delphi [59] | 0.063
$\langle N_{\pi^{0}}\rangle$ |  Delphi [59] | 0.39  | $\langle N_{\rho}\rangle$ |  Delphi [59] | 3.40
$\langle N_{p}\rangle$ |  Delphi [59] | 2.30  | $\langle N_{\Lambda^{0}}\rangle$ |  Delphi [59] | 1.54
$\langle N_{\rm ch}\rangle$ |  Delphi [60] | 0.005  | $\langle N_{\rm\pi^{\pm}}\rangle$ |  Delphi [60] | 0.10
$\langle N_{p}\rangle$ |  Delphi [60] | 0.05  | $N_{p/\bar{p}}/N_{\rm ch}$ |  Delphi [60] | 0.27
$\pi^{\pm}$ momentum |  Delphi [60] | 0.46  |  $p/\bar{p}$ momentum |  Delphi [60] | 0.43
Thrust (udsc events) |  L3 [61] | 0.34  |  $C$–parameter (udsc events) |  L3 [61] | 0.22
Charged multiplicity |  L3 [61] | 3.39  | $\log(1/x_{p})$ |  L3 [61] | 0.96
$x_{p}$ (udsc events) |  L3 [61] | 0.78  |  |  |
$\langle N_{\rm ch}\rangle$ |  Opal [62] | 0.37  |  $\pi^{\pm}$ spectrum |  Opal [63] | 0.25
$\Lambda^{0}$ scaled energy |  Opal [64] | 1.49  |  $\pi^{0}$ scaled momentum |  Opal [65] | 0.12
All events $\log(1/x_{p})$ |  Opal [66] | 0.38  | $\langle N_{\rm ch}\rangle$ |  Opal [66] | 0.16
$1-T$ |  Opal [67] | 0.10  |  $C$–parameter |  Opal [67] | 0.35
Table 3: Contributions to the $\chi^{2}/N_{\rm df}$ per each measurement that
was included in the tuning.

###  4.2 Results of the tuning

In this section, we discuss the results of the fits of the hadronization-
function parameters. In order to have a good model of hadronization, we not
only include the spectra of photons, neutral and charged pions, and baryons
but also the measurements for the event shapes (in particular the Thrust and
the $C$-parameter), mean identified particle multiplicities, charged
multiplicities and the charged momentum distributions. We have included
measurements performed by Aleph [54, 55, 56], Delphi [57, 58, 59, 60], L3 [61]
and Opal [62, 63, 64, 65, 66, 67] Collaborations (more details can be found in
Appendix A of Ref. [18]). In total, we have included 47 measurements
containing 856 bins. To ensure that this fit has a good convergence behavior,
we have used a total number of 100 scans for the minimization. The results of
the tuning are shown in Tab. 2, where we also show the results of the Monash
tuning with the Pythia 8 shower plugin [50], with the Vincia shower plugin
(Antenna shower) and with the tuning presented in Refs. [17, 18]. We also see
that the tuning of this work leads to a very good, $\chi^{2}/N_{\rm df}$ which
is slightly better than previous Monash tuning.

We display the $68\%$, $95\%$ and $99\%$ CL contours projected on the full
parameter space in Fig. 9. The figure shows some degree of correlations
between the parameters. For instance, one can see that $a_{L}$, $b_{L}$ and
$a_{QQ}$ are highly correlated as expected, while $\sigma_{\perp}$ has a
positive correlation with $a_{L}$ and negative correlation with $b_{L}$ and
$a_{QQ}$. To assess the quality of our tuning we calculate the mean
contribution to the total GoF defined as $\chi^{2}/N_{\rm bins}$ for all the
distributions in Tab. 3. We can see that for most of the measurements –
especially those directly connected to the stable particle spectra – the model
at the best-fit point yields a very good, $\chi^{2}$ which is of the order 1
or less.

Tuning $a_{L}$ $b_{L}$ $\sigma_{\perp}~{}({\rm GeV})$ $a_{QQ}$ Central $0.337$
$0.784$ $0.296$ $1.246$ $1\sigma$ eigentunes Variation $1^{+}$ $0.345$ $0.803$
$0.295$ $1.345$ Variation $1^{-}$ $0.329$ $0.766$ $0.297$ $1.149$ Variation
$2^{+}$ $0.968$ $1.444$ $0.339$ $1.074$ Variation $2^{-}$ $0.047$ $0.482$
$0.277$ $1.326$ Variation $3^{+}$ $0.207$ $0.907$ $0.255$ $1.233$ Variation
$3^{-}$ $0.476$ $0.652$ $0.340$ $1.262$ Variation $4^{+}$ $0.327$ $0.791$
$0.346$ $1.247$ Variation $4^{-}$ $0.346$ $0.778$ $0.250$ $1.246$

Table 4: Hessian variations corresponding to $1\sigma$, defined as
$\Delta\chi^{2}/N_{\rm df}=1$.

Figure 10: Spectra of stable particles in DM annihilation into $q\bar{q}$
where $q$ refers to all the quarks except the top quark and having universal
couplings for $m_{\chi}=1000$ GeV. We show the spectra of $\gamma$, $\bar{p}$,
$e^{+}$ (left panel) and of $\bar{\nu}$ (right panel). For each particle
species, we also estimate the $1\sigma$ uncertainty bands from QCD
hadronization. To have a better visibility, we also show these uncertainties
in the ratio-to-nominal subplots.

We close this section with a brief discussion of the uncertainties that are
related to the hadronization model. While in principle a comparison between
different MC event generators such as Herwig 7 [68] and Sherpa 2 [69] with the
results of Pythia 8 can yield a gross estimate of the theory systematics, we
have found that this is not the case. In Ref. [16], we have shown that the
envelope spanned by the three MC event generators, for distributions that are
relevant for e.g. $\gamma$-rays, cannot define a systematic and conservative
estimate of the uncertainties that are allowed by the data. This is a result
that has already been discussed some time ago in Ref. [70] using slightly
older MC event generators. The situation is completely different for the case
of antiproton spectra. In Refs. [17, 18] it was shown that the relative
differences between the different multipurpose MC event generators can be
extremely large and can reach up to $50$–$60\%$, especially at the extremes of
the baryon spectra. Those differences cannot be defined as uncertainties on
the hadronization model. Therefore, estimating uncertainties within the same
model seems to be the correct method for these analyses [16, 17, 18]. Here, we
only briefly discuss the size of these uncertainties assuming the Vincia
shower algorithm as our default option. We do not provide the uncertainties on
the tabulated spectra, but the interested reader can find them in the github
repository.4 The Professor tuning allows for an estimate of the uncertainties
on the parameters using the Hessian method which is widely used in the PDF
community, see e.g. Ref. [71] for more details. The method consists of
diagonalising the $\chi^{2}$ covariance matrix near the minimum of the
parameter space. The variation around the minimum of the $\Delta\chi^{2}$ can
be expanded as:

$\displaystyle\Delta\chi^{2}=\sum_{i}K_{i}(x_{i})(x_{i}-x_{i}^{0})+\sum_{i}\sum_{j}H_{ij}(x_{i},x_{j})(x_{i}-x_{i}^{0})(x_{j}-x_{j}^{0}),$
(4.6)

where $K_{i}(x_{i})$ is the first derivative of the $\chi^{2}$ which vanishes
near the minimum and $H_{ij}=\partial^{2}\chi^{2}/\partial x_{i}\partial
x_{j}$ is the Hessian matrix. The sum is over all the parameters of the model.
The diagonalisation of $H_{ij}$ leads to the so-called prinicipal directions
(eigenvectors) and the corresponding eigenvalues. We get a set of $2\times
N_{\rm params}$ variations, which in our case correspond to $8$ variations.
Imposing a constraint on the maximum variation with maximum radius equal to
the confidence level of the departure from the minimum, one obtains a one-
sigma variation if $\Delta\chi^{2}/N_{\rm df}=1$, two-sigma variation if
$\Delta\chi^{2}/N_{\rm df}=4$ and so on. The results of the Hessian method for
the $1\sigma$ eigentunes101010The eigentunes or the Hessian variations
correspond to the principal directions obtained from the diagonalization of
the matrix $H_{ij}$ near the minimum. In other words, they are obtained as
corresponding to a fixed change in the goodness-of-fit measure which is found
by imposing a constraint on the maximum variation which is defined on an
hypersphere (called the tolerance $T$). A one-sigma eigentunes corresponds to
$\Delta\chi^{2}/N_{\rm df}=1$, a two-sigma eigentunes corresponds to
$\Delta\chi^{2}/N_{\rm df}=4$ and so on. are shown in Tab. 4. The impact of
these uncertainties on the particle spectra is shown in Fig. 10 where we give
an example for a $1000$ GeV DM annihilating into $q\bar{q}$. We can see that
the uncertainties range from $10\%$–$30\%$ depending on the energy region.
These results are in a good agreement with the findings of Refs. [17, 18].

##  5 Results for the particle spectra

###  5.1 Overview of the particle spectra

Here we show a selection of results for the particle spectra we obtain. In
particular, in Figs. 11 and 12, we show the spectra obtained for all channels
considered and focusing on DM masses between 100 GeV and 100 TeV. We show the
results obtained for the annihilation of DM particles. However, for a scalar,
pseudoscalar, vector and axial vector DM particle 111111Due to the
conservation of the spin quantum number, fermionic decaying DM particles
produce an odd number of SM fermions in the final state. Such channels are not
included in our current work and are deferred to future work. our results can
also be used for decaying DM upon simple rescaling. Specifically, the
annihilation spectra for $m_{\rm DM}^{\rm ann}$ correspond to the decay
spectra of $2m_{\rm DM}^{\rm dec}$ while rescaling the $x$ variable
accordingly. We have explicitly checked that the spectra of DM annihilation
via a mediator with a given spin and coupling structure is equal to the ones
of decaying DM particle with the same properties as the mediator. Therefore,
our tables can be used for DM decays that select a specific lepton chirality
or gauge boson polarization.

In Fig. 13, we show the multiplicity for the production of the cosmic
messenger particles from DM annihilation obtained for the different channels
as a function of the DM mass. For the case of $\gamma$-rays, the channels
involving quarks and gluons produce most of the photons together with the
massive bosons. Instead, leptons and neutrinos produce much less $\gamma$
rays. The trend with the mass is similar for the production of positrons,
neutrinos and antiprotons. The difference is about a normalization factor. In
particular, the multiplicity in the production of $e^{+}$, $\nu_{e}$, and
$\nu_{\mu}$ is around 0.5, 0.5, and 1.0 times the one for $\gamma$ rays,
respectively, whereas the multiplicity for $\bar{p}$ and $\nu_{\tau}$ yields a
fraction of 0.3 and 0.003, respectively, of the multiplicity of $\gamma$ rays
only.

Figure 11: Energy spectra of cosmic messengers produced from DM annihilation.
We show from top to bottom the results obtained for DM mass of 100 GeV, 1 TeV,
10 TeV and 100 TeV. We display from left to right the spectra for $\bar{p}$,
$\gamma$ rays and $e^{+}$.

Figure 12: Same as Fig. 11 for the spectra of the three flavors of $\nu$.

Figure 13: Mean multiplicity per annihilation for the production of cosmic
messenger particles as a function of the DM mass for different annihilation
channels. From top left to bottom right we show the mean yields for
$\gamma$-rays, $\bar{p}$, $e^{+}$, $\nu_{e}$, $\nu_{\mu}$, and $\nu_{\tau}$.

###  5.2 Comparison with PPPC and HDMS

Figure 14: Comparison among the spectra obtained with our analysis using the
Pythia and Vincia shower algorithms with the results of PPPC, QCDUnc and HDMS.

In this section we compare the results of our spectra with the ones obtained
in PPPC, HDMS and QCDUnc.

We start by discussing the case of the particle production for the
$e^{+}e^{-}$ and $\mu^{+}\mu^{-}$ channels and DM masses below 1 TeV. We show
in the top left panel of Fig. 14 the comparison of our results with respect to
PPPC for the $e^{+}e^{-}$ and $m_{\chi}=100$ GeV. We do not show HDMS results
because no spectra are provided from this reference for $m_{\chi}<500$ GeV. We
see that PPPC is systematically lower than our spectra for $\log_{10}(x)<-1$.
The most important difference is that the PPPC results show a prominent cutoff
for $\log_{10}(x)<-5$. The PPPC spectra are systematically smaller than ours
at small $\log_{10}(x)$ also for the $\gamma$-ray, $e^{+}$ and neutrino
production with the $\mu^{+}\mu^{-}$ channel. At a DM mass of 100 GeV the main
$\gamma$-ray production mechanism is FSR while the EWBR process produces a
negligible contribution (see Sec. 2.2). Therefore, the significant difference
between our results and PPPC for the $e^{+}e^{-}$ and $\mu^{+}\mu^{-}$
channels should reside on the computation of the FSR. As an additional check,
we generate the spectra using another BSM model, which is the DM simplified
model with vector mediator and fermionic DM particle. We also use the
main_07.cc provided in the examples of, Pythia which produces the spectra with
a generic spin 0 resonance. Also, with these two additional models, the
results are very similar to the ones we obtain. In order to investigate this
discrepancy further, we run the simulation with the same Pythia version used
in PPPC, which is 8.135. We first use the default setup. The results we obtain
with this Pythia version are very similar to the one reported in PPPC (see top
left panel of Fig. 14). Therefore, the deviation between PPPC and our results
are due to the difference in the FSR calculation between version 8.1 and 8.3
Pythia. In order to investigate this further, we change the value of the
parameter TimeShower:pTminChgL from the default value $10^{-4}$ to $10^{-6}$,
which is the one used in the newer version of the code. We therefore traced
the discrepancy to reside in the lower threshold for this parameter adopted in
the MC: by lowering its value, the discrepancy with our results vanishes,
except for a deviation of the order of $20\%$ around $x\sim 10^{-3}$ which is
likely due to the improved tuning. The parameter TimeShower:pTminChgL is very
important for FSR because it represents a cutoff on the tranverse momentum
($p_{T}$) for the QED emissions of photons off charged leptons, as already
mentioned in Refs. [21, 15]. The results obtained by QCDUnc are overall
similar to the ones we find except that at around $\log_{10}(x)\approx-3$ the
QCDUnc spectrum has a deficit of about 20%, which is due to the fact that Ref.
[17] does not include EWBR, which contributes exactly at those energies.

The disagreement in the low-energy part of the spectra for the $e^{+}e^{-}$
and $\mu^{+}\mu^{-}$ channel remains also for DM masses larger than 100 GeV.
This is visible in the right top panel of Fig. 14 where we show the case for
10 TeV. In particular, we can see that our results as well as the ones from
HDMS, at small $x$ do not show a cutoff in the spectrum as in the case of
PPPC. This issue related to FSR in PPPC has a very relevant impact for the
indirect DM searches in $\gamma$ rays with the data of Fermi-LAT or IACTs. In
particular, for analysis of astrophysical objects for which the secondary
emission from inverse Compton scattering is not important, and thus the prompt
emission is the only one considered, e.g. with Milky Way dwarf spheroidal
galaxies. This is particularly relevant for DM masses above 1 TeV for which
the low energy tail is the important part of the spectrum for comparison
between the theoretical predictions and experimental data for the flux.

In the top right panel of Fig. 14 we show the $\gamma$-ray production for a DM
mass of 10 TeV. At these masses the EWBR process starts to be important, and
the spectra obtained with the three methods are different at the peak, which
is located at $\log_{10}(x)$ between $-4$ and $-2$. This is visible comparing
the result we find with default Pythia and Vincia shower algorithm. The
difference obtained between Vincia, PPPC and HDMS is at the level of
$20$–$30\%$ and it is present for the spectra of all particles. These
discrepancies are due to the different way the EW shower of gauge bosons is
implemented in the three analyses. In particular, PPPC is missing subsequent
radiation of gauge bosons and does not include the trilinear diagrams with
three gauge bosons. More importantly, PPPC does not include the helicity of
the particles when producing the shower. As we have seen in Sec. 3.1, this can
change the spectra at the level of $20$–$40\%$ depending on the energy scale
of the process. Instead, HDMS includes all the EW corrections considered also
in Vincia and takes into account the helicity of the particles. However, their
matching procedure may lead to issues for the particle spectra in the low
energy region. There is a last ingredient that can make a difference, which is
related to the fragmentation function parameters of Pythia version used in
PPPC and HDMS which are different from the latest version of Pythia with our
tuning.

For the channel involving quarks, the difference between PPPC or QCDUnc and
our model are at the level of $20$–$30\%$ in the most relevant energy ranges.
This is true for the production of all particles and for all masses. These
deviations are mainly due to the different versions of the Pythia code. We
remind the reader that our results are expected to be more robust because we
use a newer version of Pythia with a more refined shower algorithm that is
tuned to data relevant to DM.121212In our analysis, we employ data from
resonant $Z$-boson production only. As hadronization occurs at long distances
compared to the hard-scattering process, this approach allows us to constrain
the parameters of the hadronization model in a controlled manner. Similar
challenges occur in fits of Parton Distribution Functions. Constraining the
hadronization model using data at higher energies such as the LHC comes with
additional complications and uncertainties. The comparison with the spectra of
HDMS shows much bigger differences. In fact, the spectra of all the particles
show a similar trend in the right-hand part of the peak while in the left-hand
part there is a cutoff in the HDMS results. The reason for this discrepancy in
the low energy range is originated from the way the mass effects are
generated. In the HDMS, the mass effects are accounted for through matching at
the EW scale. This leads to smaller yields for massive particle spectra that
starts at scales of order the particle mass with respect to the EW scale, i.e.
$10^{-2}$ for the antiprotons. This issue has been discussed in details in
Appendix D.2 and shown in Figure 9 of Ref. [15]. This different trend could be
due to differences between the treatment of the HDMS and Vincia ’s EW
showering. These differences include the treatment of spin interference, a
different treatment of soft interference and the matching procedure at the EW
scale, which is not required in Vincia as it performs all evolution in the
broken description of the SM. For example, Vincia includes full soft coherence
effects in QCD processes while HDMS consider it in an approximate way.
Moreover, HDMS uses the massless parton approximation in which all the
particles are massless. Therefore, mass effects in the analysis of HDMS needs
to be taken into account. We refer the interested reader to Ref. [15] where
this point is discussed in detail. Note that our model describes the photon
data for $e^{+}e^{-}\to q\bar{q}$ produced at the $Z$-pole very well. This
data corresponds to the region for which the disagreement between our results
and those of the HDMS are the most dramatic, i.e. photon energies between $1$
and $10$ GeV.

There are other very important differences between our model, PPPC and HDMS
also when considering the spectra for the channel involving $W^{\pm}$, $Z$ and
$H$ bosons. This is visible in the bottom right panel of Fig. 14 for the
spectrum of $e^{+}$ for DM mass of 10 TeV. The high-energy part is the same
for all the three modelings while below $\log_{10}(x)<-2$ they start to differ
significantly. In particular, PPPC is the one that provides the largest yield.
We think that the large differences at the low $\log_{10}(x)$ values are due
to the treatment of resummations showers, EW showers which is absent in PPPC
and to the matching between the analytical calculations and Pythia used in the
HDMS. We stress that more work is needed in order to understand the impact of
matching and resummation on the spectra for very heavy DM.

##  6 Conclusions

We presented an improved prediction of DM annihilation spectra for cosmic
messengers, specifically, $\gamma$ rays, positrons, antiprotons, and the three
flavors of neutrinos, which are relevant for indirect detection of DM. We
employed the Vincia shower algorithm based on the helicity-dependent Antenna
functions implemented in the Pythia code. Vincia includes processes that have
not been included in Pythia before, such as the trilinear boson interaction
and the full soft-coherence in multiple FSR emissions. Furthermore, it takes
into account the helicity of the particles during the entire showering
process. For the most relevant parameters of the string hadronisation model,
we have performed an improved tuning to the LEP data at the $Z$ boson
resonance using measurements of pions, photons, hyperons and event shapes. The
resulting spectra provide a new state of the art, as we demonstrate by a
thorough comparison of our results to previous literature as well as a careful
assessment of the underlying uncertainties. The precision of our prediction
reaches a 10% level in the energy ranges most relevant for DM searches.

We have generated the spectra from 5 GeV to 100 TeV for all annihilation
channels into a pair of SM fermions and bosons including $\gamma Z$, $HZ$ and
off-shell vector bosons previously not considered in publicly available
spectra. To allow for applicability to a wide range of DM models, we have
taken into account fixed helicity states for fermions and polarization for the
EW gauge bosons. However, this approach does not take into account effects
such as ISR and IB, which requires a reevaluation of the hard process in the
specific DM model under consideration. Upon simple rescaling, our results can
also be applied to decaying DM covering the case of (pseudo)scalar and
(axial)vector DM. Our results provide important input for the interpretation
of upcoming searches for DM in the multi-TeV range performed with data from
LHAASO, HAWC and the future Cherenkov Telescope Array Observatory. The
tabulated spectra are publicly available on github.

## Acknowledgments

We thank Christian Bierlich, Peter Skands and all Pythia authors team for the
help provided with the code. We also thank Marco Cirelli, Daniele Massaro,
Olivier Mattelaer, Nicholas Rodd for insightful discussions and comments. N.F.
and M.D.M. acknowledge support from the Research grant TAsP (Theoretical
Astroparticle Physics) funded by Istituto Nazionale di Fisica Nucleare (INFN).
The work of A.J. is supported by the Institute for Basic Science (IBS) under
the project code, IBS-R018-D1. A.J. would like to thank the CERN Theory
Department for its hospitality where this work has been finalized. C.A.
acknowledges support by the F.R.S.-FNRS under the “Excellence of Science” EOS
be.h project no. 30820817. J.H. acknowledges support by the Alexander von
Humboldt foundation via the Feodor Lynen Research Fellowship for Experienced
Researchers and Feodor Lynen Return Fellowship. R.R. acknowledges support from
the Ministerio de Ciencia e Innovación (PID2020-113644GB-I00) and the GVA
Research Project Sabor y Origen de la Materia (SOM) (PROMETEO/2022/069).

##  Appendix A Parton-shower algorithms: Pythia vs Vincia

In Pythia 8, the parton showers are based on the dipole type
$p_{\perp}$-ordered evolution. This algorithm has been available since Pythia
6.3 and used for both Initial State Radiation (ISR) and Final State Radiation
(FSR) [72, 73]. On the other hand, Pythia 8 includes an implementation of the
$\gamma\to f\bar{f}$ splittings and EW massive gauge boson emissions like
$q\to q^{\prime}W$ and $q\to qZ$ [13]. The EW showers are, however, switched
off by default.

Since Pythia 8.219, photon emissions from heavy resonances such as the
$W$–bosons are included in the parton-shower machinery. The default treatment
in Pythia is based on a combination of DGLAP splitting kernels for QED+QCD
radiation with dipole ($2\to 3$) kinematics [72]. Mass effects in the parton
showers are available through matrix-element corrections (MECs) which are
available for ISR [72] and FSR [74]. The strength of the QCD showers in Pythia
8 is controlled by the effective value of the strong coupling constant
$\alpha_{S}$. Note that in Pythia 8 the value of $\alpha_{S}$ at the $Z$-pole
is not equal to $\alpha_{S}(M_{Z})^{\rm\overline{MS}}=0.118$. There are two
reasons for such a choice. In the soft-gluon emission limit, the dominant
terms for the splitting functions can be universally absorbed into the leading
order splitting kernel by a translation to the Catani-Marchesini-Webber (CMW)
scheme (also called the MC scheme) [75]. This results in an increase of
$\alpha_{S}(M_{Z})$ by about $10\%$. Furthermore, in the recent tuning of
Pythia 8 an agreement with experimental data for the measurements of the
$e^{+}e^{-}\to 3~{}{\rm jets}$ is reached if $\alpha_{S}$ is increased by
another $10\%$, see [76, 50]. Note that in Pythia 8 there is the possibility
to choose different values of $\alpha_{S}$ for ISR, FSR, MPI or the hard-
scattering process. The default value of $\alpha_{S}(M_{Z})$ in Pythia 8 for
FSR is $\alpha_{S}(M_{Z})=0.1365$ and the default choice of the RGE running is
the one-loop order. We do not change these options in this study.

It is accepted that a scale proportional to the shower evolution scale
($p_{\perp}$) to be used for the evaluation of $\alpha_{S}$ at each branching
(called thereafter the renormalisation scale). Uncertainty estimates can be
performed by variation of the renormalisation scale by a factor of two in the
positive and the negative directions:
$\\{\mu_{+},\mu_{-}\\}=\\{1/2,2\\}~{}\mu_{R}$. However, this may destroy some
of the universal corrections obtained in the CMW scheme. To solve this issue,
a framework for the automated scale variations was recently developed in Ref.
[77] and was implemented since Pythia 8.215. The formalism allows for
compensation terms to reduce the effects of large variations while having an
agreement with the CMW scheme. On the other hand, this formalism allows for
variations of the non-singular terms of the splitting functions. The firstt
application of this formalism to DM indirect detection has been done in Refs.
[16, 17, 18] where uncertainties were found to be of order $10$–$20\%$
depending on the annihilation final state and the energy region.

Vincia is a $p_{\perp}$-ordered parton-shower model for QED+QCD+EW emissions
based on the Antenna formalism. This formalism was first introduced in Ariadne
event generator [78] which is initially based on the colour dipole model [79,
80] and notably used for LEP studies. The treatment of QCD showers in FSR is
similar to the one used in Ariadne. In the case of ISR, Vincia extends the
concept of backward evolution to the formalism [81] through coherent Initial-
Initial (II), Final-Final (FF) [82] and Resonance-Final (RF) Antennae [83].
Due to the fact that all these components are coherent and interleaved in a
single sequence of decreasing $p_{\perp}$, Vincia possesses the unique
property of soft coherence for all the physical situations. For the QED
showers, the default Antenna functions include fully coherent multipole soft
interference effects, which are added to the collinear DGLAP structure [84,
85]. Such QED multipole treatment is fully interleaved with the QCD evolution
as well. This feature is very unique to Vincia.

In addition, Vincia contains an implementation of EW showers which includes
all the Higgs boson couplings and all the gauge-boson self-interactions [86,
87] allowing for all the $Z/W/H$ branchings. This EW shower module enables for
inclusion of weak corrections in EW Sudakov form factors and resummation of
multiple massive gauge-boson emissions and branchings. However, when
activating this module only the collinear limits are implemented and not the
full soft interference effects. To use the EW shower module, the helicity
information on the produced partons needs to be provided, since this module is
based on helicity-dependent shower [88, 89]. This can be achieved by either
providing the Les Houches Files (LHEF) as input to add parton showers which
explicity contain the helicity information or internally via a Vincia option
for hard-scattering matrix element calculations.

Vincia also includes interleaved resonance decays which means that short-lived
heavy resonances such as the $W/Z/H$ bosons, the top quark or any beyond-the-
SM (BSM) resonance that are produced either in the hard-scattering process or
emitted in the EW evolution are treated to be stable particles until the
evolution scale reaches a $p_{\perp}$ of order of the off-shellness scale. The
system composed of a shower plus decay is then merged into the upstream system
and the QED+EW+QCD evolution of the system continues starting from the
offshellness scale. This picture is physically intuitive, as heavy resonances
can not emit radiations at frequencies that are lower than the inverse of
their lifetime. This unique feature of Vincia leads to dramatic impacts on the
distributions of the reconstructed invariant mass, as compared to the case
where resonance decays are not interleaved with the shower evolution (like in
Pythia simple shower). The mass effects are properly taken into account for
FSR, and the corresponding Antenna functions have the appropriate limits in
the quasi-collinear regions. On the other hand, Vincia supports the sector of
Antenna showers where the phase space is divided into non-overlapping regions
and in which case every sector receives only contributions from one Antenna
branching function. This feature enables for straightforward inclusion of
higher order corrections and also for multi-jet merging, which is also called
sector merging (for more details, see Refs. [90, 91]). The choice of
$\alpha_{S}$ in Vincia for ISR and FSR is similar to that in the simple shower
in Pythia 8\. On the other hand, Vincia allows for different choices of the
scale factors for ISR and FSR emissions and splittings and different cut-off
scales for the II, IF and FF emissions.

##  Appendix B Dark Matter energy spectra calculated with MadDM

In this section, we report the relevant commands that we use for generating
the spectra with the MadDM code. The version of MadDM that we use is a custom
one which includes the Vincia shower algorithm and the tuning we have derived
in this paper. This version will be released in the comings months.

We remind the reader that we use three specific BSM models, which are:

  * •

the Singlet scalar model with a Higgs portal (SHP) (see, e.g., [38]).

  * •

DM simplified model with a fermionic DM particle and a vector boson mediator
with fermionic DM (DMSimp1) (see, e.g., [45]).

  * •

DMS imp model with a CP-odd scalar mediator model with fermionic DM (DMSimp0)
(see, e.g., [45]).

In particular, we use the SHP model for all the channels except for the
neutrinos, which are assumed to be massless and thus have coupling with the SM
Higgs boson. Moreover, for the cases for which we want to get specific
chirality or polarization states for fermions and bosons we use two cases we
use the DMSimp model. The syntax considered in MadDM to generate the spectra
for the SHP model is the following:

import model ScalarHiggsPortal_NLO_UFO
define darkmatter n1
generate indirect_detection b b
output folder_name
launch folder_name
set indirect = flux_source
set vave_indirect_cont 1e-3
set save_output spectra
set precise
set sigmav_method madevent
set nevents 5e6
set msdm 1000

In the first row the model is imported, then the DM particle is defined, and
we request the calculation of the cross-section and spectra for the bottom
channel. The mode is written in output in a folder and the analysis is
launched in the same folder. With set indirect = flux_source we request the
calculation of the flux at the source, and we set the DM velocity to
$10^{-3}c$. We require that the spectra files are saved in the folder, and we
use the precise method to calculation the cross-section using the madevent
model which fixes the DM relative velocity. Finally, we set the number of
events and the DM mass value in units of GeV.

For the spectra with the neutrino channels, we use the DMSimp1 model by
changing the couplings between the mediator and the neutrino to a value
different from 0. We use for this scope the following commands:

import model DMsimp_s_spin1_MD
define darkmatter xd
generate indirect_detection ve ve
output folder_name
launch folder_name
set indirect = flux_source
set vave_indirect_cont 1e-3
set save_output spectra
set precise
set sigmav_method madevent
set gnu11 0.25
set gnu22 0.25
set gnu33 0.25
set nevents 5e6
set msdm 1000

In order to produce the spectra for the $ZH$ channel, we use the DMSimp0
model. The syntax used for this case is:

import model DMsimp_s_spin0_MD
define darkmatter xd
generate indirect_detection h z
….
where the part after the generate command is the same as in previous cases,
shown above.

We use the same model also in the case we want to generate only left-handed
electrons:

import model DMsimp_s_spin1_MD
define darkmatter xd
generate indirect_detection e+ e-
output folder_name
launch folder_name
set indirect = flux_source
set vave_indirect_cont 1e-3
set save_output spectra
set precise
set sigmav_method madevent
set gVl11 0.25
set gVl22 0.25
set gVl33 0.25
set gAl11 -0.25
set gAl22 -0.25
set gAl33 -0.25
set nevents 5e6
set msdm 1000
where we have specified to take values of the $g_{V}$ and $g_{A}$ that are
opposite. In case we desire to produce the opposite chirality, we have to
reverse the sign of the axial couplings.

As discussed in the main text of the paper, for DM masses below the threshold
of the massive boson, the off-shell contribution of the $W$, $Z$ and $H$ can
be very relevant. In order to take into account this effect, we generate the
diagrams with four final fermions as follows:
import model ScalarHiggsPortal_NLO_UFO
define darkmatter n1
define ferm = u u d d c c s s t t b b e- e+ mu- mu+ ta- ta+ ve ve vm vm vt vt
generate indirect_detection w+ w- > ferm ferm ferm ferm / g+ g- g0 ferm
…..

In the last row, the part written as / g+ g- g0 ferm is added to remove the
goldstone bosons and the fermions from the internal legs in the diagrams. We
proceed with the same syntax also for the $HZ$ channel but with the model
DMSimp0.

Finally, for the case of polarized EW bosons, i.e.
$W_{L}W_{L},W_{T}W_{T},Z_{L}Z_{L},Z_{T}Z_{T}$, the user needs to use the
following commands:

import model ScalarHiggsPortal_NLO_UFO
define darkmatter n1
generate indirect_detection w+{0} w-{0}
...
and for transverse gauge bosons, w+{0} w-{0} needs to be replaced by w+{T}
w+{T}. Similar commands can be used for the case of the $Z$-boson, i.e.
generate indirect_detection z{0} z{0} for $Z_{L}Z_{L}$.

The output of MadDM is then passed to Pythia to add parton showers and
hadronization. Techincally, MadDM generated a LHEF with the kinematics of the
final particles in the diagram and Pythia produce the showering and
hadronization processes starting from these final particles. Below we list the
main commands to activate the Vincia Antenna shower module with EW
corrections. First, the following two commands need to be added

PartonShowers:model = 2
Vincia:ewMode = 3

The first command switch to the Vincia Antenna shower module, while the second
command activates the fully-fledged EW corrections. At the run time, Pythia
will display all the changes with respect to the default configuration that is
based on the Pythia simple shower and the Monash tuning. To use the parameters
of the hadronization model that are obtained in our tuning, the following
needs to be added as well

StringZ:aLund = 0.337
StringZ:bLund = 0.784
StringPT:sigma = 0.296
StringZ:aExtraDiquark = 1.246

There are five particles that are considered to be unstable at astrophysical
and cosmological timescales while they are stable at collider experiments,
i.e. $\pi^{\pm}$, $\mu^{\pm}$, $K^{\pm}$, $K_{L}$ and $n$. To make these
particles unstable, there are two possibilities: either increase the limit on
the default value of the proper $c\tau$ for the particle to decaying:

ParticleDecays:limitTau0 = on
ParticleDecays:tau0Max = 10.

The last command is the default option for which particles with
$c\tau>10~{}{\rm mm}$ are considered long-lived. For DM studies, we can change
the value of $c\tau$ to very high values like $10^{6}$ or something like that.
The other option is to explicitly ask Pythia to decay these particles, i.e.

13:mayDecay = true ! muon
211:mayDecay = true ! pi+-
321:mayDecay = true ! K+-
130:mayDecay = true ! Klong
2112:mayDecay = true ! neutron

## References

  * [1] K. Freese, Review of Observational Evidence for Dark Matter in the Universe and in upcoming searches for Dark Stars, EAS Publ. Ser. 36 (2009) 113–126, [arXiv:0812.4005].
  * [2] J. M. Gaskins, A review of indirect searches for particle dark matter, Contemp. Phys. 57 (2016), no. 4 496–525, [arXiv:1604.00014].
  * [3] R. K. Leane, Indirect Detection of Dark Matter in the Galaxy, in 3rd World Summit on Exploring the Dark Side of the Universe, pp. 203–228, 2020\.  arXiv:2006.00513.
  * [4] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection, JCAP 03 (2011) 051, [arXiv:1012.4515]. [Erratum: JCAP 10, E01 (2012)].
  * [5] M. Cirelli, N. Fornengo, and A. Strumia, Minimal dark matter, Nucl. Phys. B 753 (2006) 178–194, [hep-ph/0512090].
  * [6] P. Ciafaloni, M. Cirelli, D. Comelli, A. De Simone, A. Riotto, and A. Urbano, Initial State Radiation in Majorana Dark Matter Annihilations, JCAP 10 (2011) 034, [arXiv:1107.4453].
  * [7] T. Bringmann and F. Calore, Significant Enhancement of Neutralino Dark Matter Annihilation from Electroweak Bremsstrahlung, Phys. Rev. Lett. 112 (2014) 071301, [arXiv:1308.1089].
  * [8] M. Baumgart, T. Cohen, E. Moulin, I. Moult, L. Rinchiuso, N. L. Rodd, T. R. Slatyer, I. W. Stewart, and V. Vaidya, Precision Photon Spectra for Wino Annihilation, JHEP 01 (2019) 036, [arXiv:1808.08956].
  * [9] M. Beneke, A. Broggio, C. Hasner, K. Urban, and M. Vollmann, Resummed photon spectrum from dark matter annihilation for intermediate and narrow energy resolution, JHEP 08 (2019) 103, [arXiv:1903.08702]. [Erratum: JHEP 07, 145 (2020)].
  * [10] M. Beneke, S. Lederer, and C. Peset, Electroweak resummation of neutralino dark-matter annihilation into high-energy photons, JHEP 01 (2023) 171, [arXiv:2211.14341].
  * [11] M. Baumgart, N. L. Rodd, T. R. Slatyer, and V. Vaidya, The Quintuplet Annihilation Spectrum, arXiv:2309.11562.
  * [12] P. Ciafaloni, D. Comelli, A. Riotto, F. Sala, A. Strumia, and A. Urbano, Weak Corrections are Relevant for Dark Matter Indirect Detection, JCAP 03 (2011) 019, [arXiv:1009.0224].
  * [13] J. R. Christiansen and T. Sjöstrand, Weak Gauge Boson Radiation in Parton Showers, JHEP 04 (2014) 115, [arXiv:1401.5238].
  * [14] N. Fischer, S. Prestel, M. Ritzmann, and P. Skands, Vincia for Hadron Colliders, Eur. Phys. J. C 76 (2016), no. 11 589, [arXiv:1605.06142].
  * [15] C. W. Bauer, N. L. Rodd, and B. R. Webber, Dark matter spectra from the electroweak to the Planck scale, JHEP 06 (2021) 121, [arXiv:2007.15001].
  * [16] S. Amoroso, S. Caron, A. Jueid, R. Ruiz de Austri, and P. Skands, Estimating QCD uncertainties in Monte Carlo event generators for gamma-ray dark matter searches, JCAP 05 (2019) 007, [arXiv:1812.07424].
  * [17] A. Jueid, J. Kip, R. R. de Austri, and P. Skands, Impact of QCD uncertainties on antiproton spectra from dark-matter annihilation, JCAP 04 (2023) 068, [arXiv:2202.11546].
  * [18] A. Jueid, J. Kip, R. R. de Austri, and P. Skands, The Strong Force meets the Dark Sector: a robust estimate of QCD uncertainties for anti-matter dark matter searches, arXiv:2303.11363.
  * [19] M. Backovic, K. Kong, and M. McCaskey, MadDM v.1.0: Computation of Dark Matter Relic Abundance Using MadGraph5, Physics of the Dark Universe 5-6 (2014) 18–28, [arXiv:1308.4955].
  * [20] M. Backović, A. Martini, O. Mattelaer, K. Kong, and G. Mohlabeng, Direct Detection of Dark Matter with MadDM v.2.0, Phys. Dark Univ. 9-10 (5, 2015) 37–50, [arXiv:1505.04190].
  * [21] F. Ambrogi, C. Arina, M. Backovic, J. Heisig, F. Maltoni, L. Mantani, O. Mattelaer, and G. Mohlabeng, MadDM v.3.0: a Comprehensive Tool for Dark Matter Studies, Phys. Dark Univ. 24 (2019) 100249, [arXiv:1804.00044].
  * [22] A. Buckley et al., General-purpose event generators for LHC physics, Phys. Rept. 504 (2011) 145–233, [arXiv:1101.2599].
  * [23] X. Artru and G. Mennessier, String model and multiproduction, Nucl. Phys. B 70 (1974) 93–115.
  * [24] T. Sjostrand, The Lund Monte Carlo for Jet Fragmentation, Comput. Phys. Commun. 27 (1982) 243.
  * [25] B. Andersson, G. Gustafson, G. Ingelman, and T. Sjostrand, Parton Fragmentation and String Dynamics, Phys. Rept. 97 (1983) 31–145.
  * [26] T. Sjostrand, Jet Fragmentation of Nearby Partons, Nucl. Phys. B 248 (1984) 469–502.
  * [27] B. R. Webber, A QCD Model for Jet Fragmentation Including Soft Gluon Interference, Nucl. Phys. B 238 (1984) 492–528.
  * [28] J.-C. Winter, F. Krauss, and G. Soff, A Modified cluster hadronization model, Eur. Phys. J. C 36 (2004) 381–395, [hep-ph/0311085].
  * [29] M. G. Bowler, e+ e- Production of Heavy Quarks in the String Model, Z. Phys. C 11 (1981) 169.
  * [30] OPAL Collaboration, G. Abbiendi et al., A Study of parton fragmentation in hadronic Z0 decays using Lambda anti-Lambda correlations, Eur. Phys. J. C 13 (2000) 185–195, [hep-ex/9808031].
  * [31] B. Andersson, G. Gustafson, and T. Sjostrand, Baryon Production in Jet Fragmentation and $\Upsilon$ Decay, Phys. Scripta 32 (1985) 574\.
  * [32] P. Eden and G. Gustafson, Baryon production in the string fragmentation picture, Z. Phys. C 75 (1997) 41–49, [hep-ph/9606454].
  * [33] C. Bierlich, P. Ilten, T. Menzo, S. Mrenna, M. Szewc, M. K. Wilkinson, A. Youssef, and J. Zupan, Reweighting Monte Carlo Predictions and Automated Fragmentation Variations in Pythia 8, arXiv:2308.13459.
  * [34] A. Cuoco, B. Eiteneuer, J. Heisig, and M. Krämer, A global fit of the $\gamma$-ray galactic center excess within the scalar singlet Higgs portal model, JCAP 06 (2016) 050, [arXiv:1603.08228].
  * [35] J. Erler and M. Schott, Electroweak Precision Tests of the Standard Model after the Discovery of the Higgs Boson, Prog. Part. Nucl. Phys. 106 (2019) 68–119, [arXiv:1902.05142].
  * [36] W. J. Stirling and E. Vryonidou, Electroweak gauge boson polarisation at the LHC, JHEP 07 (2012) 124, [arXiv:1204.6427].
  * [37] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079, [arXiv:1405.0301].
  * [38] M. Di Mauro, C. Arina, N. Fornengo, J. Heisig, and D. Massaro, Dark matter in the Higgs resonance region, Phys. Rev. D 108 (2023), no. 9 095008, [arXiv:2305.11937].
  * [39] J. Alwall, E. Boos, L. Dudko, M. Gigg, M. Herquet, A. Pukhov, P. Richardson, A. Sherstnev, and P. Z. Skands, A Les Houches Interface for BSM Generators, 12, 2007.  arXiv:0712.3311.
  * [40] DELPHI Collaboration, P. Abreu et al., Cross-sections and leptonic forward backward asymmetries from the Z0 running of LEP, Eur. Phys. J. C 16 (2000) 371–405.
  * [41] A. Prasath V, R. M. Godbole, and S. D. Rindani, Longitudinal top polarisation measurement and anomalous $Wtb$ coupling, Eur. Phys. J. C 75 (2015), no. 9 402, [arXiv:1405.1264].
  * [42] A. Jueid, Probing anomalous $Wtb$ couplings at the LHC in single $t$-channel top quark production, Phys. Rev. D 98 (2018), no. 5 053006, [arXiv:1805.07763].
  * [43] A. Arhrib, A. Jueid, and S. Moretti, Top quark polarization as a probe of charged Higgs bosons, Phys. Rev. D 98 (2018), no. 11 115006, [arXiv:1807.11306].
  * [44] J. Abdallah et al., Simplified Models for Dark Matter Searches at the LHC, Phys. Dark Univ. 9-10 (2015) 8–23, [arXiv:1506.03116].
  * [45] C. Arina, Impact of cosmological and astrophysical constraints on dark matter simplified models, Front. Astron. Space Sci. 5 (2018) 30, [arXiv:1805.04290].
  * [46] F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz, and S. Vogl, Implications of unitarity and gauge invariance for simplified dark matter models, JHEP 02 (2016) 016, [arXiv:1510.02110].
  * [47] C. Englert, M. McCullough, and M. Spannowsky, S-Channel Dark Matter Simplified Models and Unitarity, Phys. Dark Univ. 14 (2016) 48–56, [arXiv:1604.07975].
  * [48] J. Ellis, M. Fairbairn, and P. Tunney, Anomaly-Free Dark Matter Models are not so Simple, JHEP 08 (2017) 053, [arXiv:1704.03850].
  * [49] C. Bierlich et al., A comprehensive guide to the physics and usage of PYTHIA 8.3, arXiv:2203.11601.
  * [50] P. Skands, S. Carrazza, and J. Rojo, Tuning PYTHIA 8.1: the Monash 2013 Tune, Eur. Phys. J. C 74 (2014), no. 8 3024, [arXiv:1404.5630].
  * [51] C. Bierlich et al., Robust Independent Validation of Experiment and Theory: Rivet version 3, SciPost Phys. 8 (2020) 026, [arXiv:1912.05451].
  * [52] A. Buckley, H. Hoeth, H. Lacker, H. Schulz, and J. E. von Seggern, Systematic event generator tuning for the LHC, Eur. Phys. J. C 65 (2010) 331–357, [arXiv:0907.2973].
  * [53] F. James and M. Roos, Minuit: A System for Function Minimization and Analysis of the Parameter Errors and Correlations, Comput. Phys. Commun. 10 (1975) 343–367.
  * [54] ALEPH Collaboration, R. Barate et al., Studies of quantum chromodynamics with the ALEPH detector, Phys. Rept. 294 (1998) 1–165.
  * [55] ALEPH Collaboration, R. Barate et al., Inclusive production of pi0, eta, eta-prime (958), K0(S) and lambda in two jet and three jet events from hadronic Z decays, Eur. Phys. J. C 16 (2000) 613.
  * [56] ALEPH Collaboration, A. Heister et al., Studies of QCD at e+ e- centre-of-mass energies between 91-GeV and 209-GeV, Eur. Phys. J. C 35 (2004) 457–486.
  * [57] DELPHI Collaboration, P. Abreu et al., Charged particle multiplicity distributions for fixed number of jets in Z0 hadronic decays, Z. Phys. C 56 (1992) 63–76.
  * [58] DELPHI Collaboration, P. Abreu et al., Production of Lambda and Lambda anti-Lambda correlations in the hadronic decays of the Z0, Phys. Lett. B 318 (1993) 249–262.
  * [59] DELPHI Collaboration, P. Abreu et al., Inclusive measurements of the K+- and p / anti-p production in hadronic Z0 decays, Nucl. Phys. B 444 (1995) 3–26.
  * [60] DELPHI Collaboration, P. Abreu et al., pi+-, K+-, p and anti-p production in Z0 —$>$ q anti-q, Z0 —$>$ b anti-b, Z0 —$>$ u anti-u, d anti-d, s anti-s, Eur. Phys. J. C 5 (1998) 585–620.
  * [61] L3 Collaboration, P. Achard et al., Studies of hadronic event structure in $e^{+}e^{-}$ annihilation from 30-GeV to 209-GeV with the L3 detector, Phys. Rept. 399 (2004) 71–174, [hep-ex/0406049].
  * [62] OPAL Collaboration, P. D. Acton et al., A Study of charged particle multiplicities in hadronic decays of the Z0, Z. Phys. C 53 (1992) 539–554.
  * [63] OPAL Collaboration, R. Akers et al., Measurement of the production rates of charged hadrons in e+ e- annihilation at the Z0, Z. Phys. C 63 (1994) 181–196.
  * [64] OPAL Collaboration, G. Alexander et al., Strange baryon production in hadronic Z0 decays, Z. Phys. C 73 (1997) 569–586.
  * [65] OPAL Collaboration, K. Ackerstaff et al., Photon and light meson production in hadronic Z0 decays, Eur. Phys. J. C 5 (1998) 411–437, [hep-ex/9805011].
  * [66] OPAL Collaboration, K. Ackerstaff et al., Measurements of flavor dependent fragmentation functions in Z0 –$>$ q anti-q events, Eur. Phys. J. C 7 (1999) 369–381, [hep-ex/9807004].
  * [67] OPAL Collaboration, G. Abbiendi et al., Measurement of event shape distributions and moments in e+ e- —$>$ hadrons at 91-GeV - 209-GeV and a determination of alpha(s), Eur. Phys. J. C 40 (2005) 287–316, [hep-ex/0503051].
  * [68] J. Bellm et al., Herwig 7.0/Herwig++ 3.0 release note, Eur. Phys. J. C 76 (2016), no. 4 196, [arXiv:1512.01178].
  * [69] T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert, and J. Winter, Event generation with SHERPA 1.1, JHEP 02 (2009) 007, [arXiv:0811.4622].
  * [70] J. A. R. Cembranos, A. de la Cruz-Dombriz, V. Gammaldi, R. A. Lineros, and A. L. Maroto, Reliability of Monte Carlo event generators for gamma ray dark matter searches, JHEP 09 (2013) 077, [arXiv:1305.2124].
  * [71] J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai, and W. K. Tung, Uncertainties of predictions from parton distribution functions. 2. The Hessian method, Phys. Rev. D 65 (2001) 014013, [hep-ph/0101032].
  * [72] T. Sjostrand and P. Z. Skands, Transverse-momentum-ordered showers and interleaved multiple interactions, Eur. Phys. J. C 39 (2005) 129–154, [hep-ph/0408302].
  * [73] T. Sjostrand, S. Mrenna, and P. Z. Skands, PYTHIA 6.4 Physics and Manual, JHEP 05 (2006) 026, [hep-ph/0603175].
  * [74] E. Norrbin and T. Sjostrand, QCD radiation off heavy particles, Nucl. Phys. B 603 (2001) 297–342, [hep-ph/0010012].
  * [75] S. Catani, B. R. Webber, and G. Marchesini, QCD coherent branching and semiinclusive processes at large x, Nucl. Phys. B 349 (1991) 635–654.
  * [76] P. Z. Skands, Tuning Monte Carlo Generators: The Perugia Tunes, Phys. Rev. D 82 (2010) 074018, [arXiv:1005.3457].
  * [77] S. Mrenna and P. Skands, Automated Parton-Shower Variations in Pythia 8, Phys. Rev. D 94 (2016), no. 7 074005, [arXiv:1605.08352].
  * [78] L. Lonnblad, ARIADNE version 4: A Program for simulation of QCD cascades implementing the color dipole model, Comput. Phys. Commun. 71 (1992) 15–31.
  * [79] G. Gustafson, Dual Description of a Confined Color Field, Phys. Lett. B 175 (1986) 453.
  * [80] G. Gustafson and U. Pettersson, Dipole Formulation of QCD Cascades, Nucl. Phys. B 306 (1988) 746–758.
  * [81] M. Ritzmann, D. A. Kosower, and P. Skands, Antenna Showers with Hadronic Initial States, Phys. Lett. B 718 (2013) 1345–1350, [arXiv:1210.6345].
  * [82] H. Brooks, C. T. Preuss, and P. Skands, Sector Showers for Hadron Collisions, JHEP 07 (2020) 032, [arXiv:2003.00702].
  * [83] H. Brooks and P. Skands, Coherent showers in decays of colored resonances, Phys. Rev. D 100 (2019), no. 7 076006, [arXiv:1907.08980].
  * [84] R. Kleiss and R. Verheyen, Final-state QED Multipole Radiation in Antenna Parton Showers, JHEP 11 (2017) 182, [arXiv:1709.04485].
  * [85] P. Skands and R. Verheyen, Multipole photon radiation in the Vincia parton shower, Phys. Lett. B 811 (2020) 135878, [arXiv:2002.04939].
  * [86] R. Kleiss and R. Verheyen, Collinear electroweak radiation in antenna parton showers, Eur. Phys. J. C 80 (2020), no. 10 980, [arXiv:2002.09248].
  * [87] H. Brooks, P. Skands, and R. Verheyen, Interleaved resonance decays and electroweak radiation in the Vincia parton shower, SciPost Phys. 12 (2022), no. 3 101, [arXiv:2108.10786].
  * [88] A. J. Larkoski, J. J. Lopez-Villarejo, and P. Skands, Helicity-Dependent Showers and Matching with VINCIA, Phys. Rev. D 87 (2013), no. 5 054033, [arXiv:1301.0933].
  * [89] N. Fischer, A. Lifson, and P. Skands, Helicity Antenna Showers for Hadron Colliders, Eur. Phys. J. C 77 (2017), no. 10 719, [arXiv:1708.01736].
  * [90] S. Höche, S. Mrenna, S. Payne, C. T. Preuss, and P. Skands, A Study of QCD Radiation in VBF Higgs Production with Vincia and Pythia, SciPost Phys. 12 (2022), no. 1 010, [arXiv:2106.10987].
  * [91] J. M. Campbell, S. Höche, H. T. Li, C. T. Preuss, and P. Skands, Towards NNLO+PS matching with sector showers, Phys. Lett. B 836 (2023) 137614, [arXiv:2108.07133].