Hadronic Footprint of GeV-Mass Dark Matter

GeV-scale dark matter is an increasingly attractive target for direct detection, indirect detection, and collider searches. Its annihilation into hadronic final states produces a challenging zoo of light hadronic resonances. We update Herwig7 to study the photon and positron spectra from annihilation through a vector mediator. It covers dark matter masses between 250 MeV and 5 GeV and includes an error estimate.


Introduction
The fundamental nature of dark matter is the biggest particle physics question of our time. It follows directly from the success of quantum field theory in describing the properties of elementary particles, as well as the standard cosmology after the Big Bang. A problem is that the term particle dark matter is very loosely defined, covering very light particles essentially giving a background wave function all the way to primordial black holes. Theory motivations are widely used to support certain mass ranges, but given the generally modest success of models for physics beyond the Standard Model they should be taken with a truck load of salt [1].
The defining feature of dark matter particles -closely related to the one actual measurement of the relic density [2] -is the dark matter production mechanism. It needs to explain the observed relic density in agreement with the largely known thermal history of the universe. For wide classes of dark matter models this implies an interaction with SM particles beyond an obviously existing gravitational interaction. Dark matter masses around the weak scale and down to the GeV scale can be produced thermally through freeze-out or freeze-in, or through a relation with the baryon asymmetry in the universe. These mechanisms require more or less strong couplings to the SM matter particles, i.e. to leptons or quarks. In this mass range there exists a wealth of relatively model-independent direct detection constraints [3], and their extension to lighter dark matter particles at and below the GeV-scale is one of the most interesting experimental directions [4][5][6]. For such GeV-scale dark matter especially the couplings to quarks are not well constrained.
Indirect searches for dark matter are another way to directly probe the properties of the dark matter in the universe. A leading signature are photons produced in the annihilation of dark matter in dense regions of the sky [7][8][9][10]. These photons can constrain dark matter interactions with the Standard Model on the lepton and the hadron side. The reason is that they can be produced directly and through radiation from any charged annihilation product. Just as for positrons, we know the photon spectrum from many particle physics experiments over the recent decades. They are available for instance through the Pppc4dm tool [11] based on Pythia [12]. Standard tools like micrOMEGAs [13,14], MadDM [15,16], or DarkSusy [17,18] include similar spectra based on multi-purpose Monte Carlo generators. A major technical problem with dark matter annihilation into hadrons is that its description is not available through Pythia once the dark matter masses drop below around 5 GeV. The only exception is the recent Hazma [19] tool for dark matter masses below 250 MeV [20]. This leaves dark matter annihilation to the leading hadronic final states for masses between 250 MeV and 5 GeV essentially uncovered.
Technically, GeV-scale dark matter annihilation through light scalar or light vector mediators is very different. If we assume that a new scalar couples to SM particles with Yukawa couplings roughly reflecting the SM mass hierarchy, increasingly weak GeV-scale dark matter will annihilate into charm quarks and tau leptons, followed by muons and eventually pions and electrons. From a hadronic physics point of view the more interesting scenario are vector mediators, where the SM interactions are generation-universal. In that case we will observe a wealth of hadronic annihilation channels below the bb threshold. These annihilation channels will have distinct photon and lepton spectra, which we will focus on in this study.
Finally, the proper description of dark matter annihilation to hadronic final state is plagued by large uncertainties, as for instance pointed out for Pppc4dm [11] in relation to Pythia [12] and Herwig [21]. Consequently, dedicated comparisons between Pythia and Herwig have been published for dark matter annihilation to tau leptons, bottom and top quarks, and weak bosons [22]. More recently, this comparison has been updated [23] to the most recent versions of Pythia8 [12,24] and Herwig7 [25,26]. A detailed analysis of the the Pythia predictions can be found in Ref. [27]. All of these studies target relatively heavy dark matter annihilation, in line with the common weakly interacting massive particle (WIMP) hypothesis.
In this paper we provide the first proper description of photon and lepton spectra from GeV dark matter annihilating into hadronic final states based on Herwig with an updated fit to electron-positron data, including several new final states. They become relevant when we reduce the dark matter scattering energy below the Pythia limit. We update the fit to electron-positron data as the input to the Herwig description and add the necessary new hadronic final states with up to four hadrons. Especially for the photon spectrum we observe a complete change in the shape of the spectrum when we reduce the dark matter mass, starting from typical hadron decay chains to continuum multi-pion production. In addition we provide a first estimate of the impact of the input-data fit uncertainties on the output spectra.
The paper is structured the following way: after introducing our toy model in Sec. 2 we review the established implementations in Sec. 3. We show how their reliability starts to fade once we go below dark matter masses of 5 GeV and the tools start to extrapolate beyond their common Pythia input. In Sec. 4 we show the results from our new Herwigbased implementation. We focus on shape changes in the photon and lepton spectra when we reduce the dark matter mass towards into the continuum-pion regime. We also show the error bands on the photon and positron spectra from the fit uncertainties to the electronpositron input data. In the Appendix we provide all details about our new fit, the underlying parametrizations, the best-fit points, and the error bands.

Toy model
The standard interpretation framework for weak-scale dark matter is thermal freeze-out production or the WIMP paradigm. Embedding GeV-scale dark matter searches in a global analysis [28] provides an excellent illustration of the many cosmological constraints and their model dependence. Above masses around 10 GeV, FERMI constrains these scenarios using photons in dwarf spheroidal galaxies [29,30], while AMS covers leptonic final states [31,32]. In addition, precision measurements of the Cosmic Microwave Background (CMB) [2] are sensitive to the total ionizing energy either directly (electrons and muons) or indirectly. Finally, Big-Bang Nucleosynthesis (BBN) does not allow WIMPs below around 10 MeV [33,34]. The main difference between these different analyses is that some rely on assumptions on the thermal history of dark matter.
To define a toy model for our hadronization study we assume that the observed dark matter density is somehow produced through thermal freeze-out, but with a light mediator. We assume the dark matter candidate to be a Majorana fermion χ, but our results apply the same way to asymmetric dark matter where the dark matter fermion has to be different from its anti-particle. A simple mediator choice starts from an additional U (1) gauge symmetry, where we gauge one of the accidental global symmetries related to baryon and lepton number [54][55][56][57][58]. For our purpose of testing dark matter annihilation into light-flavor jets with a limited number of photons from leptonic channels the most attractive combination is B − 3L µ [59]. This gives us the annihilation channel To avoid strong biases from an underlying model we also show results for Z couplings similar to the Standard Model case for low energies. For consistent field theory models the annihilation to SM quarks will always occur at the loop level, even if they are suppressed at tree level [60]. As our benchmark model we therefore assume an approximately on-shell annihilation The coupling strength of the DM to the mediator can be chosen arbitrarily, because we are only interested in the form of the energy spectra from the hadronic final states. For light dark matter masses the relevant quarks are u, d, s, c, possibly the bottom quark. The charm quark plays a special role, because threshold region is poorly understood. All we can do is rely on the spectra included in Pythia or Herwig, with little improvement on the modelling side.
For the three lightest quarks there exists a wealth of measurements which we can use to constrain dark matter annihilation into hadrons. We decompose a quark DM current J µ DM = q=u,d,s a qq γ µ q into isospin components and a separate ss contribution, where a q are the couplings of the light vector mediator to the light quarks, q = u, d, s. The mediator couplings to quarks are fixed to a q = 1/3 for any anomaly-free B − L model. Depending on the mediator coupling structure to quarks, one or the other isospin current might vanish. As a consequence, some resonance contributions to the channels might vanish, or even more drastically, pure isospin I = 0 channels, for example χχ → πππ, ηπ, . . . .
are absent. We also choose to include the isospin breaking contribution from ω → π + π − in the I = 1 current for simplicity. The general matrix element for DM annihilation can be written in the form with the DM-mediator coupling a DM and the vector mediator propagator d DM νµ . In our toy model we always assume m Z = 2m χ , but given the non-relativistic nature of DM annihilation the mediator mass should only have negligible impact on our spectra. Since the mass of the mediator determines the width of the mediator, we calculate the width in the hadronic resonance region within Herwig through its decays to all kinematically allowed hadronic final states listed in Tab. 2 of the Appendix.

Established tools
Different public tools generate energy spectra for different DM annihilation channels to SM particles. They are limited in DM masses by their approach and by their back-end. We summarize • Pppc4dm [11] provides tabulated energy spectra for indirect detection. The e ± ,p,d, γ, and ν e,µ,τ fluxes are generated with Pythia8.135 [12] down to m χ = 5 GeV. We use the provided interpolation routine to extrapolate the results to m χ = 2 GeV.
• micrOMEGAs [13,14] uses tabulated Pythia spectra for γ, e + ,p, ν e,µ,τ and extrapolates down to m χ = 2 GeV. In the manual of version micrOMEGAs2.0 it is mentioned that the strategy for calculating spectra is analogous to that of DarkSusy and that spectra extrapolated to masses below 2 GeV should be taken with care.
• MadDM [15,16] provides two ways of calculating the energy spectra both based on Pythia [24]. The 'fast' calculation is based on the numerical tables provided by Pppc4dm.
In the 'precise' mode, events are generated with MadGraph and then passed to Pythia for showering and hadronization. In this mode it is possible to calculate the fluxes of any final states based on the UFO model implementation.
• DarkSusy [17,18] provides tables down to 3 GeV for energy spectra of two-particle SM final states based on Pythia6.426 [61]. The tool can interpolate and extrapolate the γ, e + ,p,d, π 0 , ν e,µ,τ , µ fluxes for all quark final states. In addition it includes annihilation to µµ, τ τ , gluons, and weak bosons. Dark matter annihilation into e + e − pairs appears to not be included.
• Hazma [19] is a Python toolkit to produce energy spectra in the sub-GeV range. It is based on leading order chiral perturbation theory and is valid in the non-resonance  From this list it is clear that for dark matter masses in the GeV range all public tools are based on Pythia, one way or another. Multi-purpose Monte Carlo tools, such as, Pythia or Herwig can calculate the energy spectra for many hard scattering processes, followed by hadronization or fragmentation and hadron decays. In the range we are interested in these spectra are usually extracted from data, as discussed in the Appendix. The advantage of the Monte Carlo tools is that we can extract the cosmologically most relevant photon, lepton, and anti-proton spectra for each hard dark matter annihilation process. We assume an annihilation process of the kind given in Eq.(2), but allow for any kinematically allowed SM final state. For the numerical results we rely on spectra from the processes e + e − → SM pairs, at a given energy m ee = 2m χ . In Fig. 1 we compare the corresponding Pythia-like spectra from the standard tools discussed above. We show the photon and positron spectra from DM annihilation into muon, tau, and light-quark (u, d, s) pairs and compare them to the standard Herwig output for an alternative description. Starting with the left panels of Fig. 1 we see a flat photon spectrum from soft-enhanced radiation and a triangular positron spectrum from the µ + -decay with a three-particle final state. For taus the hadronic decays produce neutral and charged pions, where for instance the decay π 0 → γγ dominates down to x ≈ 10 −3 . Below that we again find the flat photon spectrum from soft emission. The dominant contribution to the position spectrum is the hadronic decay chain τ + → π + → µ + → e + , with a sub-dominant contribution from the leptonic β-decay τ + → e + . Next, light-flavor quarks u, d, s form a range of hadrons which then decay to π 0 → 2γ. The positron spectrum from these light quarks includes a soft neutron β-decay, which gives rise to the secondary maximum around x ≈ 10 −4 . The neutron decay is not included in our default version of micrOMEGAs, but can be easily added. Finally, moving to DM annihilation into charms we see that the photon and positron spectra are the same as for the light quarks.
In Fig. 2 we show the same spectra, but for a slightly lower dark matter mass of 2 GeV. This value is slightly beyond where Pythia output can be used in a straightforward manner. Essentially all radiation and decay patterns remain the same as for 5 GeV, but the different curves start moving apart. This is an effect of individual extrapolations from the Pythia output. The only interesting feature appears in the annihilation χχ → cc. Here the extrapolated results from Pppc4dm and DarkSusy still include a secondary peak corresponding to the neutron decay in the light quark channel. However, the lightest charm baryon is Λ c has a mass of 2.29 GeV, so at m χ = 2 GeV it cannot be produced on-shell. What we see is likely an over-estimate of off-shell effects or an extrapolation error from the 5 GeV case, which illustrates the danger of ignoring the explicit warning not to use for instance Pppc4dm or DarkSusy below their recommended mass ranges. For micrOMEGAs the spectrum is significantly softer than from the dedicated MadDM call to Pythia and from Herwig.
Altogether we find that for m χ = 5 GeV there is a completely consistent picture, where the Pythia-based results are in excellent agreement with Herwig. Going to m χ = 2 GeV leads to an increased variation between the different tools and illustrates why we might not want to use the standard tools outside their recommended mass ranges.

Herwig4DM spectra
To extend the range of valid simulations of dark matter annihilation to quarks we start with the standard Herwig7 [25,26,62] implementation. We then add a set of additional final states and update some other spectra, as discussed in the Appendix. This allows us to cover DM masses down to twice the pion mass. Below the threshold m Z = 2m π ≈ 250 GeV the annihilation to hadrons will be suppressed and the annihilation to electrons and photons will dominate. In Fig. 3 we show the photon and positron spectra from the annihilation process for decreasing DM masses from m χ = 2 GeV to 250 MeV.

Spectra
Most photons and positrons in hadronic processes come from neutral and charged pion decays, respectively. These pions are either directly produced or are the end of a decay chain of all forms of hadronic states listed in Tab. 2 in the Appendix. In a few cases, photons can also be directly produced in DM annihilation, for instance χχ → ηγ, πγ .
In the left panel of Fig. 3 we see how photon production channels drop out when we reduce the DM mass or center-of-mass energy of the non-relativistic scattering process. Whereas for m χ > 1 GeV all possible hadronic final states contribute to the round shape of the spectrum, for lower energies only photons and positrons from very specific processes give a characteristic energy spectrum. For example for m χ = 500 MeV or equivalently a center-of-mass energy of 1 GeV we expect two kaons from the φ resonance to provide most photons through consecutive decays of kaons to pions to photons. This leads to a triangular shape of the photon spectrum. If we go down to 250 MeV, the only remaining annihilation channels are χχ → π 0 γ, ππ, 3π .
Of those, the photons mainly come from the π 0 γ final state, so one photon is produced directly with an energy around the DM mass. It leads to the sharp peak around x ≈ 1. The additional photons come from the π 0 -decay and are responsible for the distribution to roughly x ≈ 10 −1 .
The same applies for m χ = 375 MeV with an additional bump-like contribution from neutral pions in the 3π and 4π channels as well as additional photons from the dominantly neutral ηγ → (2γ)γ, (3π 0 )γ decay including a direct photon.
The basic shape of the positron spectrum is given by the neutron pair production threshold. Above threshold, we observe an additional peak slightly above x ∼ 10 −4 resulting from positron production in the neutron β-decay. For m χ < 1 GeV, all positrons come from charged pion decays. The peak position depends on how early that charged pion decay occurs for the dominant processes at the respective center-of-mass energy. For example, for m χ < 500 MeV, charged pions are mainly produced directly in ππ, 3π, 4π production and hence the peak of the spectrum is shifted towards x = 1.
As mentioned in Sec. 2, the composition of the DM current changes with the way the mediator couples to quarks. In any (B − 3L)-like model with equal couplings to quarks, the isospin I = 1 contribution vanishes and consequently some resonance contributions as well as all channels listed in Eq.(5) vanish. For m χ = 250 MeV this implies that without the ππ channel, π 0 γ becomes the dominant annihilation mode. The direct photon production lifts the photon spectrum, as seen in the upper panels of Fig. 3. This is accompanied by a drop in the positron spectrum that only receives contributions from the subdominant 3π final state. If we choose a center-of-mass energy below the 3π threshold, positron spectra from quarks would be completely absent. For m χ = 375 MeV with an increasing 3π contribution towards the ω(782) resonance, the position spectra are lifted. For higher energies and the contribution from several channels, the (B − 3L)-like spectra resemble the SM-like case. Although their sources are not identical channel by channel, the way the photons and positions are produced is similar.

Error bands
Uncertainties on the energy spectra are dominated by the uncertainties from the fits to electron-positron data discussed in the Appendix. We define ranges of model parameters to cover bands in the e + e − -annihilation cross sections as a function of the energy and propagate those parameter ranges through the hadronic currents into the energy spectra. This means that the error on a given spectrum corresponds to the uncertainty of the dominant channel at the corresponding energy.
In the upper panels of Fig. 4 we see that the photon spectrum at m χ = 250 MeV inherits large uncertainties from the poorly measured dominant π 0 γ channel in that energy range. For m χ = 375 MeV the more precisely measured 3π channel suppresses the π 0 γ channel, but still leaves us with visible error bands. For even higher energies several channels contribute to the uncertainty of the photon spectrum. We observe the smallest error bands for spectra that benefit from precisely measured dominant processes, for instance peak regions such as the φ resonance at 1 GeV in the KK channel, the ρ resonance in the 2π decay, or generally well-measured channels such as 4π. Positron spectra with their dominant 2π, 3π, 4π channels are always well measured. The only exception is m χ = 1 GeV spectrum, especially the lower peak around ∼ 10 −4 , which comes from the neutron β-decay. As discussed in the Appendix, the nn channel is poorly measured and leaves us with larger uncertainties in that regime.
In (B − 3L)-like models, we will not get any contributions from well-measured 2π and 4π final states. This means the uncertainties on the position spectrum for m χ = 250 MeV are slightly larger than in the SM-like case, see the lower panel of Finally, we want to ensure that our error estimates are conservative. In Fig. 4 we use the uncertainties on the individual channels bin-wise, add all contributions up and normalize by the sum of their corresponding cross-sections. For channels with large cross-sections that are also giving the main contribution to the total amount of photons/positrons in the spectrum, the error bars can completely cancel for the normalized spectra. This way, we only get sizable uncertainty bands for spectra where one channel is dominating the shape of the spectrum, but is playing a sub-dominant role in the total cross-section. An example is the π 0 γ final state for the SM-like photon spectrum at m χ = 250 MeV or the lower bump in the 1 GeV positron spectrum caused by nn. This assumption can be considered somewhat aggressive in a situation where we do not have full control of the full error budget. Instead, we can maximize and minimize all spectra channel by channel and separately normalize them by the smallest and largest total cross-section possible. This way there will be no cancellation for single-channel spectra, and in Fig. 5 we indeed see much increased uncertainties. Obviously, the real error bands are going to be somewhere between the results shown in Fig. 4 and Fig. 5 determined by analysis details beyond the scope of this first analysis.

Outlook
We have studied the positron and photon spectra from non-relativistic dark matter annihilation in a dark matter mass range from 250 MeV to 5 GeV (with the exception of the poorly understood region near the charm threshold). We consider a light vector mediator with general couplings to SM fermions. For the photon spectra we see a smooth interpolation from typical hadron decay chains with their round spectra down to the pion continuum with a triangular shape. For positrons the main feature is the secondary neutron decay above threshold.
Because we are relying on an updated fit to electron-position input data to Herwig we can also propagate the uncertainties from poorly measured channels into the photon and positron spectra. Already for relatively heavy dark matter the positron spectrum shows sizeable error bars. In the case of photons, smaller dark matter masses with fewer and less well measured annihilation channels are also plagued by significant error bars, eventually covering an order of magnitude for m χ = 250 MeV.
Our new implementation closes the gap between standard Pythia-based tools such as Pppc4dm, MicrOMEGAS, MadMD, or DarkSUSY and the comparably simple smallmass continuum regime and should allow for a reliable study of GeV-scale dark matter even if it dominantly interacts with SM quarks. *

Acknowledgments
First, we want to thank Patrick Foldenauer for his early contributions to Figs. 1 and 2. TP is supported by the German Research Foundation DFG under grant no. 396021762-TRR 257 and would like to thank Stefan Gieseke for help starting this project. Peter Reimitz is funded by the Graduiertenkolleg Particle physics beyond the Standard Model (GRK 1940). Peter Richardson is supported by funding from the UK Science and Technology Facilities Council (grant numbers ST/P000800/1, ST/P001246/1), and benefited from the European Union's Horizon 2020 research and innovation programme as part of the Marie Sk lodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104). * The code we have used to produce these results will be available in a future version of Herwig7. If there is sufficient interest we will also think about providing the output as a cool and fast neural network.

A Updated fits with error envelopes
If we limit ourselves to dark matter annihilation through a vector mediator we can relate the dark matter annihilation process to the corresponding and measurable process e + e − → hadrons.
Its matrix element has the form The electromagnetic quark current J µ em = q=u,d,s e qq γ µ q can be decomposed into its isospin components I = 0, 1 and its strange-quark content, with We study all hadronic states which appear in the total cross section σ(e + e − → hadrons) in the MeV to GeV range. An list of all channels, their parametrizations, their data fits, and their threshold values is given in Tab Table 3: Parameters of the nucleon form factor from our fit using the model describing pp production from Ref. [92].
meson dominance [116]. In that case the hadronic current had|J µ em |0 can be described by a momentum-dependence and a form-factor that include all resonances allowed under certain isospin symmetry assumptions. The parametrization and fit values for the form-factors for the πγ, ππ, 3ππ, ωπ, and ηγ final states are taken from Refs. [67,69,72], as implemented in the event generator Phokhara [117,118], and the Born cross section formulae from the SND measurements [63,73,93]. For all other channels, we provide new fits.

pp (update)
The data and the fit function for this channel are given in Tab. 2. We updated the data set used for our fit since from the input to the the previous fit [92] Ref . [119] is superseded by Ref. [84], Ref. [120] by Ref. [89], and Ref. [121] by Ref. [75]. For asymmetric data uncertainties we symmetrize statistical and systematic uncertainties separately and then add both in quadrature. We refrain from a more sophisticated error analysis for instance including correlations between systematic uncertainties, since in most cases detailed information about the systematic uncertainties is either missing or the statistical uncertainty dominates. For the fit, we get χ 2 /n.d.f = 1.069, and the best-fit values are shown in Tab. 3.

ηππ, η ππ (update)
The fit function for the ηππ and η ππ hadronic currents are based on [96]. We re-fit the fit function to more recent data sets [94,95] compared to those used in [96]. The fit values can be found in Tab Table 5: Parameters for the description of KK production from our fit using the model of Ref. [67]. All masses and widths are given in GeV, all other parameters are dimensionless

KK (update)
We parametrize the hadronic current for the K 0K0 and K + K − channels in the same way as done in Ref. [67]. Unlike Ref.
[67], we do not fix all masses and widths of the ρ, ω and φ states to their PDG values but let them float in the fit. Furthermore, we use an updated data set for the fit, as mentioned in Tab. 2 and included the τ − → K 0 S π − ν τ data from Ref. [122] to better constrain the I = 1 component of the current. The fit values are listed in Tab. 5. The last coupling of each resonance is calculated via Eq.(16) in Ref. [67], and we keep η φ = 1.055, γ ω = 0.5 and γ φ = 0.2 fixed such as in Ref. [67]. For the simultaneous fit to K 0 K 0 and K + K − data we obtain χ 2 /n.d.f = 1.621.

4π (update)
For the 4π channel, we use the parametrization of Ref. [72] and fit it to more recent rate measurements for e + e − → 2π 0 π + π − and e + e − → 2π + 2π − from BaBar [70,71]. We obtain a χ 2 /n.d.f = 1.28 and the fit values are listed in Tab. 6.  Table 6: Parameters for the 4π channel for our fit using the model from [72]. All masses and widths are in GeV; couplings β j i , (j = a 1 , f 0 , ω and i = 1, 2, 3) as well as c ρ are dimensionless; c a 1 and c f 0 in GeV −2 and c ω in GeV −1 .  Table 7: Fit values for the ηφ, ηω, and φπ channels.

ηφ, ηω, φπ (new)
Our first new fit is to the processes e + e − → ηφ, ηω, φπ, where the momentum-dependent Born cross sections are where α em (s) is the fine structure constant, P f (s) = q 3 cm,X the final-state phase space, q cm,X the final-state particle momentum and F is the respective form factor. The resonant contributions are simply parametrized by where we take the s-dependent width Γ(s) from Ref.
[100]. All parameters and fit values for ηφ, ηω, and φπ production are listed in Tab. 7.

ωππ (new)
Next, for the ωππ channel, we use for the hadronic current. In our energy range we only need to consider one vector meson mediator ω , namely the ω(1650) meson. For the f i mediator we have where m σ and Γ σ are the mass and width of the σ meson and using the Flatté parametrization [124] BW f 0 (m ππ ) = g ω ωf 0 (980) m f 0 (980) with Γ ππ = g ππ q π (m ππ ) for the f 0 (980) meson, with parameters from Ref. [125]. If not mentioned otherwise, the parameter are set to their PDG values [123]. The σ meson contribution can be viewed as a phase space contribution to the ωππ channel more than resonant contribution. Therefore, the width is chosen to be large, see Tab. 8.

KKπ (new)
Below 2 GeV center-of-mass energy the process e + e − → KKπ is dominated by e + e − → KK * → K(Kπ) where KK * can be either K 0 K * 0 (890) or K ± K * ∓ (890). We can relate the possible final states through their isospin I = 0, 1 and can use the following relations for the corresponding amplitudes A 0,1 [126], In the energy range we are dealing with, we expect the resonances to be φ(1680) and φ(2170) for I = 0 and ρ(1450) and ρ(1700) for I = 1. The lower resonances ρ(770) and φ(1020) are not considered in the energy range of the fit and we set their couplings to zero. Furthermore, we fix the mass and the width of the intermediate K * resonance to m K * = 0.8956 GeV and Γ K * = 0.047 GeV and use a p-wave Breit-Wigner propagator of the form with the s-dependent width where m 1 , m 2 are the decay products of the K * state and their velocity in the rest frame of K * . The K * Kπ coupling is given by Furthermore, we include a small φπ 0 contribution for final states including neutral pions by adding the φπ 0 cross section obtained by the φπ fit and the corresponding branching fractions BR(φ(1020) → K 0 L K 0 S ) = 0.342 and BR(φ(1020) → K + K − ) = 0.489. We perform a simultaneous fit to all possible final states in order to obtain the fit parameters of the amplitudes A 0,1 . The fit values can be found in Tab. 9.
We show all numerical best-fit solutions as blue lines for all final states in Figs.6, 7, and 8. The error bars on the data are dominated by statistical uncertainties. All fits describe the most recent data sets over the entire range shown.

Error bands
In addition to the central values of the relevant parameters describing the e + e − data we also estimate the error bands for the relevant processes. The reason is that some of the channels are rather poorly measured, and it is important to propagate these uncertainties through the analysis. Because most fit parameters are physical parameters appearing in the analytic description of the e + e − cross sections, such as masses or widths or rates, we do not find them suitable for a proper statistical analysis. For instance a total cross section measurement will lead to uncontrolled correlations between widely different phase space regions in the fit, where the different phase space regions are crucial to describe the dark matter spectra for a variable dark matter mass. Examples for the impact of a known form of the energy dependence of the scattering process on poorly measured phase space regions are the ηππ channel in Fig. 6, the ππ channel in Fig. 7, or the 3π channel in Fig. 8.
Instead, we define envelopes by varying a sub-set of fit parameters around their mean value within their uncertainty provided our python IMinuit [127,128] fit or as stated in papers. For poorly resolved peak structures as in the η ππ, φπ, and ηω case or higher resonances as in ηφ and KKπ, we do not vary any widths and only some masses, since are determined from the peak structure and bias the off-peak spectrum through correlations. The contribution of phases to our envelopes is only considered if no other set of parameters is sufficient to describe the measurement uncertainties. For channels with simple parametrizations with fixed masses and widths and floating peak cross sections and phases as in the case of πγ [63] and ηγ [93], we vary all peak cross sections and the phases of the φ and ω resonance, respectively. In these cases, we see that away from the resonance region the error envelopes increase. For precisely measured phase space regions, we consider the full set of parameters describing these regions. These are usually large peak structures such as the φ → KK and ρ → ππ resonances in Fig. 7 or the ω, ρ → 3π peak around 0.78 GeV in Fig 8. Those resolved regions turn out to be well described and are stable against variations of the parameters, so they give only small envelopes.
It can be challenging or nearly impossible to obtain consistent envelopes for some channels, where one parametrization is used for several sub-channels as in the case of KK and pp/nn. As long as the shape of the data is the same as in the case of 4π, KKπ and the φ resonance region in the KK channel, this does not cause any problems. Here we can assume that a parameter and its variation influence the fit curve in the same way. However, for energies above 1.4 GeV in the KK channel, the trend of the data of K + K − and K 0K0 is completely different. Therefore, already the fit to the data is challenging and only possible by allowing for more resonance fit parameters in the parametrization [67]. A variation of the parameter might influence both channels differently and it is not clear that an extremal value in the one case is also extremal in the other. This tension of both data sets causes too small error bands for energies above 1.8 GeV. For the pp/nn channel, we do not have sufficient data for nn to describe this channel properly as already described in Ref. [92].     [100] BaBar, B. Aubert et al., Measurements of e + e − → K + K − η, K + K − π 0 and K 0 s K ± π ∓ cross-sections using initial state radiation events, Phys. Rev. D77 (2008) 092002, arXiv:0710.4451 [hep-ex]. [103] BaBar, J. P. Lees et al., Cross sections for the reactions e + e − → K 0 S K 0 L π 0 , K 0 S K 0 L η, and K 0 S K 0 L π 0 π 0 from events with initial-state radiation, Phys. Rev. D95 (2017) 5, 052001, arXiv:1701.08297 [hep-ex].