Cosmic-Ray Positrons Strongly Constrain Leptophilic Dark Matter

Cosmic-ray positrons have long been considered a powerful probe of dark matter annihilation. In particular, myriad studies of the unexpected rise in the positron fraction have debated its dark matter or pulsar origins. In this paper, we instead examine the potential for extremely precise positron measurements by AMS-02 to probe hard leptophilic dark matter candidates that do not have spectral features similar to the bulk of the observed positron excess. Utilizing a detailed cosmic-ray propagation model that includes a primary positron flux generated by Galactic pulsars in addition to a secondary component constrained by Helium and proton measurements, we produce a robust fit to the local positron flux and spectrum. We find no evidence for a spectral bump correlated with leptophilic dark matter, and set strong constraints on the dark matter annihilation cross-section that fall below the thermal annihilation cross-section for dark matter masses below 60 GeV and 380 GeV for annihilation into $\tau^+\tau^-$ and $e^+e^-$, respectively, in our default model.


Introduction
The AMS-02 experiment, onboard the International Space Station, has provided extremely precise observations of the flux, composition and spectrum of local cosmic-rays [1]. These data yield new insights into the energetics of the Galaxy's most powerful sources, and the properties of Galactic transport. Additionally, searches for unexpected spectral features can potentially provide evidence for new physics, such as anisotropic diffusion or even dark matter annihilation.
In this paper, we consider the constraints that can be placed on DM in scenarios that attribute the bulk of the rising positron fraction above 10 GeV to Galactic pulsars. We make use of the high-statistical precision of the AMS-02 data to calculate strong constraints on DM annihilation into highly-peaked leptonic final states, which have spectra that are inconsistent with the high-energy positron spectrum. A schematic example of the scenario we examine is shown in Figure 1 Figure 1. A schematic diagram of the scenario we examine. We model a positron flux that includes a primary contribution from pulsars, secondary contributions from the hadronic interactions of primary cosmic-rays, and a highly peaked contribution from DM annihilation into leptonic final states, which we attempt to constrain. DM signals detectable by AMS-02 generally range from a few GeV to a few hundred GeV.
produced by the interaction of protons and He nuclei with the interstellar medium, while at higher energies, most positrons are produced in electron-positron pairs accelerated by pulsars. Against this background, our analysis probes a subdominant contribution from annihilating DM. DM annihilation into hard-spectrum leptonic final states are well-motivated in the context of leptophilic dark matter models [12,47]. Originally discussed in the context of fitting cosmic-ray excesses from PAMELA, AMS-02, ATIC and DAMA/LIBRA [48,49], leptophilic models have also been theoretically motivated, for example, in the context of Z' models [47], and have been found to produce unique signatures in particle colliders [50][51][52]. Moreover, models with significant leptonic couplings are generically expected in well-motivated scenarios such as dark photon models [53,54]. We also note that signals from hard-spectrum leptonic final states may also be dominant even in models with loop-level hadronic interactions [55].
Our result builds on a similar analysis by Bergström et al. (2013) [56]. Their work tests two discrete models for the astrophysical background, as well as contributions from DM particles between 5-300 GeV annihilating into leptonic final states (τ + τ − , µ + µ − , e + e − γ and e + e − ). Their strongest limits are obtained for the direct e + e − channel, where they set upper limits on the DM velocity-averaged annihilation cross section of about 10 −28 cm 3 s −1 at 5 GeV and about 10 −25 cm 3 s −1 at 250 GeV.
In this paper, we update and improve on the Bergström et al. results in three ways: (1) Our analysis utilizes a larger dataset provided by the latest AMS-02 measurements of positron, proton, and He fluxes [5,[57][58][59]. (2) We use the cosmic-ray propagation code Galprop [60] to create a full cosmic-ray propagation model for the positron, proton and He data. Our fit allows a large variety of parameters, such as the diffusion and particle injection spectra in the fit to float, while Bergström et al. build their background from a handful of pre-defined Galprop models. (3) We implement a new solar modulation model [61] to describe the effect of solar activity on low-energy cosmic-ray particles, rather than approximating solar modulation by a force-field potential.
The paper is outlined as follows: In Section 2, we describe the methodology of this analysis and explain the data selection (Sec. 2.1), the computational model (Sec. 2.2), the DM model (Sec. 2.3) and the statistical analysis (Sec. 2.4). The results are presented in Section 3, and our conclusions are discussed in Section 4.

Methodology
In this section, we describe our selection of the AMS-02 data, and present our methodology for constructing astrophysical models that fit this data. We also describe our statistical analysis to constrain the DM contribution.
Our analysis generally proceeds as follows. To compute the DM limits, we perform a fit of the positron, proton and He data from AMS-02 to constrain the parameters of our astrophysical model. While the high-energy positron spectrum is dominated by primary positrons that are created by pulsars as electron-positron pairs, low-energy positrons are produced as secondaries from proton interactions with a subdominant contribution from He (about 10%) [59]. After obtaining an astrophysical model for the positron flux in the absence of DM, we add a DM contribution and compute the 95% confidence upper limit on the DM annihilation rate. Notably, we allow critical astrophysical parameters to be re-fit for each DM mass and cross-section, producing a more realistic constraint on the DM contribution. In the following sections, the procedure is described in more detail.

Data Selection
Our analysis takes advantage of the recently published AMS-02 data for positrons [5], protons [57,58], and He [59]. A summary of the data sets is given in Table 1. We note that our constraints on leptophilic DM only directly depend on our analysis of the positron data, but we include protons and He in the astrophysical model to obtain better constraints on the model parameters. We also note that protons and He have the same electric charge as positrons, which decreases the relative uncertainties from solar modulation.
For positrons, we consider energies in the range from 2 to 1000 GeV in our fit. We do not include energies below 2 GeV due to their large uncertainties, especially coming from solar modulation [62], and because the DM contribution is small at low energies compared to the very bright contribution from secondary production (see e.g. Figure 1 for the different contributions to the positron flux). To achieve the best statistical precision, we take a large dataset that includes the AMS-02 data from May 2011 to November 2017 [5]. For protons, we consider energies in the range from 2 to 1800 GeV, for which we combine two data sets: From 2 to 60 GeV, AMS-02 provides time-dependent data [57], and we select the data from February 2015 to May 2017 (Bartels rotations 2477-2506), after the polarity flip of the solar magnetic field was completed, to avoid additional uncertainties from the effects of the solar modulation. At energies from 60 to 1800 GeV, we use the time-independent data in [58], recorded from May 2011 to November 2013. At these energies, the effect of solar modulation is negligible. For He, we consider the energy range from 2 GeV to 1000 GeV and use the timeindependent data from May 2011 to November 2013, which includes both 3 He and 4 He [59]. As the conversion from rigidity to energy depends on the number of nucleons, we estimate the ratio of 3 He to 4 He in this flux by linearly extrapolating the ratio given in [63] to fluxes at higher energies.

Computational Model
To simulate galactic cosmic-ray propagation, we use Galprop v.56 [60,[64][65][66][67][68][69][70] 1 , a state-ofthe-art tool to compute cosmic-ray propagation, that takes into account numerous processes, such as diffusion, spallation, energy losses and reacceleration. For our analysis, it is sufficient to solve the propagation equation in 2D, and to include only cosmic-ray contributions up to 4 He, while ignoring heavier nuclei. Galprop simulations begin by injecting a significant flux of primary cosmic-rays (most importantly for our analysis, protons and 4 He) following the standard supernova spatial distribution [71]. The subsequent propagation and hadronic interactions of these primary particles produce a secondary flux, which includes more protons as well as new components (most importantly, positrons and 3 He). We supplement these contributions from supernovae with an additional component of primary electrons and positrons powered by pulsar activity, as described below.
We tune the following Galprop parameters in our model (for a list, see Table 2). We model the diffusion spectrum with a diffusion coefficient D 0 set at a reference rigidity of 4 GV and a diffusion break D break , where δ 1 and δ 2 are the spectral indices below and above the break, respectively. The proton injection spectrum is modeled by three free parameters, including two spectral indices and a critical energy describing their transition. This break primarily affects the lowest energy protons and positrons. Similarly, for He we introduce a spectral break and spectral indices γ He 1 and γ He 2 below and above the break, respectively. At higher energies, most of the positron contribution comes from pulsars, where the free parameters are the spectral index γ psr , the cutoff energy E psr cut of the spectrum and the pulsar formation rateṄ 100 (discussed below). Reacceleration is taken into account by the Alfvén velocity which is most important for positrons and has a smaller effect on protons. The change in convection velocity perpendicular to the Galactic plane dv/dz (with convection velocity v = 0 at z = 0) and the solar modulation parameters (φ 0 and φ 1 , discussed below) are significant at lower energies. Lastly, there is an overall normalization of the spectra computed by Galprop which is related to uncertainties in the Milky Way supernova rate [72]. Since secondary positrons are produced mostly by protons, the normalization is the same for both fluxes.
A particularly important parameter in the model is the halo height, z. The DM contribution is proportional to the halo height, and therefore we examine this parameter individually. We compute the astrophysical model for two fixed values of the halo height: setting z = 5.6 kpc as a standard height (default model), and z = 3 kpc as the lower limit on the halo height (small-height model) [73][74][75].
In our model, we introduce pulsars as a new source of primary positrons by modifying the Galprop code. For this, we adopt the standard pulsar distribution profile from [76,77], to obtain the 3-dimensional pulsar distribution in the Galaxy, where z is the halo height of the Milky Way, and R 0 = 1.528 kpc and z 0 = 0.330 kpc. We model the e + e − spectrum from pulsars as a power-law with an exponential cutoff [6], where γ psr is the spectral index and E psr cut the cutoff energy. For the normalization, we calculate C such that the average total pulsar spindown luminosity is 10 49 erg [6], and the average efficiency for the conversion of pulsar spindown power into e + e − pairs is 1%, a value which is consistent with previous empirical fits to the positron data [78] and forms an intermediate value between the ∼10% injection efficiency of mature TeV halos such as Geminga [29] and the <1% efficiency of young pulsars such as Vela [25]. In our model, we let γ psr , E psr cut anḋ N 100 free, fitting the spectrum and normalization of the pulsar contribution.
Cosmic-ray measurements are affected by the heliosphere of the Sun, causing particle losses at low energies ( 10 GeV). This effect strongly depends on the orientation of the solar magnetic field and solar activity. Often (e.g., as in [56]), the solar modulation is approximated by a force-field potential [62]. Here, we implement the solar modulation model from [61] that takes into account the time-, charge-and rigidity dependence of the solar modulation. In the solar modulation model, we choose time periods that match the data analysis periods from AMS-02, and keep the parameters φ 0 and φ 1 , that describe the heliospheric potential, as free parameters.
Additionally, there are several important parameters that we fix to their standard Galprop values. Most notably, we utilize a Galactic magnetic field model based on Model C of Ref. [79], utilizing the re-fit value for the radial scale height of the regular magnetic field of 13 kpc. This model yields a magnetic field at the solar position of 5.7 µG. All parameters for our astrophysical models can be found in Table 2 for the default model and small-height model.

Dark Matter Model
For the DM model, we assume a local DM density of 0.4 GeV/cm 3 and a core radius of the DM halo of 20 kpc (c.f. [80]), see Table 3, and describe the DM distribution in the Galaxy using a generalized NFW profile [81]. We utilize a non-standard NFW slope of 0.5, but note that this choice has almost no effect on the local DM annihilation rate -producing limits that are more conservative by a few percent. We consider four different leptophilic final states for the annihilation of Majorana DM: χχ → τ + τ − , χχ → µ + µ − , χχ → e + e − , and χχ → φφ → e + e − e + e − , where φ is a light mediator. For the τ + τ − , and µ + µ − final states, we use DM spectra obtained with DarkSUSY v.5.1.1 [82]. For the e + e − e + e − and e + e − states, we utilize analytic calculations of the DM spectrum that we directly implement in Galprop. Specifically, in the e + e − e + e − case, the contribution to the differential positron flux is constant and integrates to the DM mass over the full energy range, following a similar  Table 3. List of parameters that are fixed in the model. calculation in [83]. The e + e − final state is represented by a delta function and we smear out the contribution to the positron flux over the neighboring energy bins in Galprop. Here, we neglect electroweak corrections as these are in the order of 1%.

Statistical Analysis
To calculate the best-fit values for our model parameters, we use two minimization codes: pyMultinest v. We create our astrophysical model in two steps: First we fit positrons and protons over a wide range of parameters. For this we use Multinest as it can handle the large number of free parameters well and is reliable at finding the global minimum. Then we add He to the fit and re-fit the previous parameters together with the He injection spectrum. For this, we use iMinuit as it finds the exact values of the minimum more precisely than Multinest. Because their contribution to the positron flux is negligible, we do not include nuclei heavier than He in Galprop.
Throughout our analysis, we compare our model with the AMS-02 particle fluxes in multiple energy bins. We note that the AMS-02 data includes uncertainties that stem from both statistical and systematic errors. The systematic errors are likely to be correlated between energy bins, but to date, the AMS-02 collaboration has not released a covariance matrix that accounts for the correlation between systematic errors in nearby energy bins. Thus, in this paper, we treat the systematic and statistical errors in each bin independently, adding them in quadrature. This leads to an overestimate of the AMS-02 error budget, and allows our models to achieve goodness of fit values less than 1. Notably, however, this preserves the same minimum value for the best-fit, and still allows us to compare models in terms of their ∆χ 2 .
For the first step, we use Multinest to fit the astrophysical parameters corresponding to proton and positron propagation discussed in Section 2.2. We show the intermediate results of this procedure in Appendix A.
For the second step, we take the best-fit values and perform an additional fit by adding the He data to the model, using iMinuit. We set the 4 He injection spectrum to follow a twice-broken power-law, with a fixed spectral index below 5 GV of 1.8, and allow the two higher spectral indices and the location of the higher spectral break to float. We allow the normalization of the primary He flux to float. We note that there is also a small primary 3 He flux, with a spectrum set to the primary proton spectrum. The relative normalization of these terms, 7.199 · 10 4 to 9.033, makes the effect of primary 3 He injection negligible. Note that in Galprop, these abundances are defined relative to the (arbitrarily chosen) proton normalization of 1.06 · 10 6 . We perform the fit to the He data using iMinuit, also allowing the parameters of the previous step to float. As we want to focus on the positron fit of our model, we introduce an additional systematic uncertainty of 2% into the He data, accepting that our model does not include several uncertainties that are specific to He. More importantly, the systematic uncertainties in the He flux are correlated, but here we treat them as uncorrelated.
For the third step, we investigate the DM contribution. While we include protons and He when fitting our astrophysical background model to the data, to obtain more realistic constraints on the positron flux, we consider the combined (DM and astrophysical contribution) to only the positron flux when we set DM limits. This always produces more conservative upper limits because DM does not directly contribute to the proton or He fluxes (leaving the secondary positron contribution approximately the same), and thus the combined (pro-ton+He+positron) χ 2 fit typically provides stronger limits if DM induced adjustments to the diffusion parameters perturb these datasets.
To determine the DM limits, we study DM masses between 5 and 2000 GeV for different annihilation cross sections σv . For each DM mass and final state, we set up a grid of annihilation cross sections using logarithmic steps in the mass range of interest and an iterative process to obtain values for the annihilation cross section, and fit the model using five parameters that are particularly relevant for the DM contribution to the positron flux (D 0 , δ 2 ) or are degenerate with the DM contribution (E psr cut , γ psr andṄ 100 ), and let them float as well for each grid point. Then we compute the χ 2 profile of σv and use this to calculate the 95% upper confidence limit compared to a null DM contribution (i.e., we take χ 2 DM = χ 2 DM=0 + 3.84 which corresponds to the 95% upper limit on a χ 2 distribution with one degree of freedom).

Results
In Table 2 we show the best-fitting parameters of our astrophysical background model to the combined AMS-02 positron, proton and He data (Step 2). In general, our best fit parameters for the diffusion coefficients, proton injection spectral indices, the z-dependence of the convection velocity, Alfvénic velocity and solar modulation parameters are consistent with previous analyses [73,[90][91][92]. We note that our pulsar formation rate is about an order of magnitude smaller as e.g. in Ref. [6]. However, this value depends strongly on the assumptions made about the pulsars (e.g., total spin-down luminosity and initial period), as well as the spectrum, which controls the fraction of the pulsar spin-down power that is injected into the critical 10-1000 GeV range where the pulsar contribution is significant. Several of these parameters (e.g., the total spin-down luminosity) have order of magnitude uncertainties. Another notable result in our fit is that the break in the proton injection spectrum is positioned at a few GeV, which is very close to the lower limit of our model. This means that a break in the proton injection spectrum is not heavily favored in our model.
In Table 4, we show the resulting fits for both of our models to the combined AMS-02 dataset, as well as the individual fits for each cosmic-ray species. For the default model with a standard halo height of 5.6 kpc, the total χ 2 is 88.60 over 141 degrees of freedom. The χ 2 per degree of freedom for positrons, protons and He are 0.88, 0.43, and 0.57, respectively. The χ 2 /d.o.f. is less than 1 for each model and species, indicating that our model fits the data very well. We note that we treat the systematic uncertainties of nearby energy bins as uncorrelated, even though they are likely to be correlated, which results in an overestimation of the uncertainties and thus a χ 2 /d.o.f. below 1. The total χ 2 is better in the default model than the small-height model by ∼ 3σ, which means that 3 kpc approximately represents the minimum halo height that is consistent with AMS-02 data.
In Figure 2 we show our best-fit model for positrons, protons and He, comparing our results to the AMS-02 data for the best-fit model parameters given in Table 2. In the bottom of each plot in Figure 2, we give the residuals ((data-model)/model) of each fit, showing that our model fits the data to within a few percent, providing a very good fit to the cosmic-ray data that are relevant to the positron flux. We observe some systematic uncertainties for low-energy positrons (below 6 GeV) and protons and He at higher energies (above 600 GeV). This may correspond to additional secondary proton reacceleration in supernova remnants that is not taken into account in Galprop [92]. Note that the gray shaded regions below 2 GeV and above about 800 GeV, 1050 GeV and 1000 GeV for positrons, protons and He, respectively, are not included in the fit. This means that our model has larger uncertainties at energies near these limits.
We point out that our analysis does not take into account the correlation of the systematic uncertainties between nearby energy bins in the AMS-02 data. As can be observed in the residual plots in Figure 2, especially for protons and He, the best-fit residual values are highly correlated between nearby energy bins, pointing towards a strong correlation in the systematic uncertainties. We are not able to take these correlations into account because the covariance matrices are not provided by AMS-02, which leads to a miss-estimation of the uncertainties. We additionally note that our calculation of the χ 2 /d.o.f. for He is artificially decreased by the additional 2% systematic error added into the He uncertainties. However,  we stress that the best-fit value of our parameters, as well as the value of ∆χ 2 between our models, is preserved in light of these systematic uncertainties.
In Appendix B, we show an example of the positron flux including a DM contribution for annihilations into e + e − at 5 GeV where the excess is largest.
In Figure 3, the 95% confidence upper limits on the DM annihilation cross-sections are presented for the various final decay states. The left panel shows the results of the default model and the right panel the small-height model. For reference, we also include the thermal cross section from [93] (for a more recent analysis, see also [94]). We derive the strongest constraints for DM annihilation directly into e + e − pairs with an annihilation cross section of 6 · 10 −29 cm 3 /s for 12.6 GeV in the default model, and an annihilation cross section of 2.6 · 10 −29 cm 3 /s for 10.0 GeV in the small-height model, both of which lie below the thermal annihilation cross section by about 3 orders of magnitude.
Our limits show a small bump at low energies for all final states. Specifically, in the default model for the τ + τ − final state, the bump peaks at 10 GeV and corresponds to a preference for DM at a local significance of 2.8σ, for µ + µ − the peak is at 8 GeV with a local significance of 2.9σ, for e + e − e + e − the peak is at 6.3 GeV with a local significance of 2.8σ, and for e + e − the peak is at 5 GeV with a local significance of 3.1σ. In the small-height model, for τ + τ − the excess peaks at 13 GeV and 2.7σ, for µ + µ − at 8 GeV and 3.9σ, for e + e − e + e − at 6 GeV and 3.6σ and for e + e − at 5 GeV and 3.6σ. These all correspond to the same spectral feature in the low energy positron fit of our model, which can be seen in Figure 2 as a positive discrepancy (about 2%) between our model and the AMS-02 data below about 5 GeV. We note that this region of our model is prone to larger uncertainties, as these energies are close to the lower limit of our model (2 GeV). Additionally, lower energies are highly affected by uncertainties in solar modulation. Note that the significance of these excesses are local, and we have not implemented a proper trials correction here (as e.g. in [95]) which would reduce the local significance of these excesses. Following the method in [95] requires the covariance matrices of the experimental data to simulate pseudo-experiments, but these covariance matrices are not supplied by the AMS-02 experiment. Nonetheless we estimate the correction by counting the number of upcrossing in the residuals of the positron flux, following the method in [95]. We count 17 upcrossing which reduces the significances from about 3.0σ to about 2.0σ, and from 3.64σ to 2.83σ. However, we note that the significance of this result is likely much lower, due to the fact that this excess appears very close to the boundary of our modeling, where our models of solar modulation and cosmic-ray propagation    [5], protons [57,58] and He [59], displayed in red. The left panels show the default model (5.6 kpc halo height) and the right panels the small-height model (   are likely beginning to break down.
We point out that recent studies [73,96] of the cosmic-ray data suggest an additional break in the diffusion spectrum or the particle injection spectra at a few hundred GV that we do not consider here but could potentially strengthen our limits.

Discussion and Conclusions
In this paper, we calculated constraints on the DM annihilation cross section for DM masses between 5 and 2000 GeV for leptonic final states using the local positron flux. We exclude the thermal cross-section for annihilation into τ + τ − below 60 GeV, for µ + µ − below 160 GeV, for e + e − e + e − below 240 GeV and for e + e − below 380 GeV for our default model with a halo height of 5.6 kpc.
We repeated our analysis with a conservative halo height of 3 kpc in the small-height model, for which we obtain similar constraints: Our limits are below the thermal cross section for annihilation into τ + τ − below 50 GeV, for µ + µ − below 170 GeV, for e + e − e + e − below 430 GeV and for e + e − below 500 GeV. We note that our constraints at higher energies are stronger in the small-height model than the default model. This is because our two models have a different diffusion coefficient at high energies, meaning that the diffusion distance is slightly longer in the small-height model which results in the stronger constraints at high masses ( 300 GeV). We conclude from our results that a change in halo height does not significantly affect the local positron flux, contrary to limits on the local antiproton flux that depends on the halo height [73]. This is due to the fact that the average positron travels a shorter distance than antiprotons at similar energy (because their energy losses are faster)making the edges of the diffusion zone less important for leptonic DM limits.
Our analysis improves on the constraints by Bergström et al. by about a factor of 2 for both models. This is despite the fact that our model gives significantly more freedom to the astrophysical background fit -and thus includes more conservative and robust limits. Note that Bergström et al. use a halo height of 4 kpc. We are also able to extend the limits up to DM masses of 2000 GeV as more cosmic-ray data is available at higher energies.
Our limits also set strong constraints on DM models with only subdominant annihilations into the four leptonic final states considered in this work. At about 30 GeV, our limits fall below the thermal cross section by about two orders of magnitude for annihilation into e + e − . This means that our limits even rule out subdominant contributions from DM annihilations (i.e., ∼1% annihilation into e + e − pairs at 30 GeV), and similarly ∼10% for annihilations into µ + µ − . Notably, DM models with significant annihilation fractions into µ + µ − final states are among the most difficult thermal DM models to constrain with current experiments [97]. Compared to the results of Ref. [97], which only ruled out thermal annihilation to µ + µ − up to a DM mass of 30 GeV, we rule out annihilation to µ + µ − up to masses of ∼150 GeV. We note that a significant fraction of this difference pertains to less conservative model assumptions utilized in this work, compared to that of Ref. [97]. For example, they use a local DM density of 0.25 GeV/cm 3 and local magnetic field strength of 8.9 µG, both of which serve to significantly decrease the local positron flux from DM at a given annihilation cross-section. On the other hand, our limits benefit from the use of a physical model for the astrophysical contributions from pulsars and secondary production, rather than producing an empirical fit to the data. Our limits also profit from the better data by AMS-02, which dominates the improvement of our limits at higher energies. Thus our results can also improve generic constraints on DM annihilation to mixed final states.
Our limits are highly complementary with γ-ray constraints on DM annihilation. In particular, direct annihilations to e + e − , µ + µ − and light mediator final states are difficult to probe with γ-ray instruments, because their decays produce few γ-rays at tree-level. On the other hand, γ-ray studies are highly sensitive to DM annihilation into hadronic final states, which are not tested here. Our limits, in conjunction with current limits on DM annihilation to hadronic states (e.g., in dwarf spheroidal galaxies) place strong limits on DM annihilation regardless of final state.

A Astrophysical Model with Multinest
In this section, we present the results from the first step of creating our astrophysical background models, in which we fit our Galprop model to the AMS-02 data for positrons and protons (see Table 1 for the specific AMS-02 data sets) with Multinest [84], i.e., we do not fit the He spectrum in this step. We choose to include this step because Multinest is very reliable at finding the global minimum of the fit, whereas iMinuit is more likely to fall into local minima in large multi-dimensional parameter spaces. The free parameters can be found in Tables 5 and 6 for the default and small-height model, respectively, and include the following parameters: the diffusion spectrum, the proton injection spectrum (with a spectral break and a spectral index below and above the break), the change in convection velocity, the Alfvén velocity, the Galactic pulsar formation rate and spectrum (with a spectral index and energy cutoff), the solar modulation parameters and an overall normalization of the primary particle abundance of 1.06 · 10 6 for protons in Galprop.
We present our results in Tables 5 and 6 for the default and small-height model, respectively, where we give the range in which we allow the parameters to float, the type of prior (i.e., linear or logarithmic), the best-fit value, the posterior value and the posterior uncertainty. In particular, we scan the priors for D 0 , D break , the proton injection spectral break and E psr cut in logarithmic space, while we use a flat distribution for the other parameters. We obtain results that are relatively consistent between the default and small-height model. Figure 4 shows our fit for positrons and protons for our default model and the smallheight model, together with the AMS-02 data for the local fluxes. The fractional residuals between the data and our model are given in the bottom of each panel, showing that our models fit the AMS-02 data to within a few percent for the data that is relevant to the positron flux.
A summary of the χ 2 -values is given in Table 7  The best-fit values (rather than the posterior means) obtained in this step serve as the initial values for the iMinuit fit in the second step, allowing us to start the fitting procedure close to the global minimum. Comparing our Multinest fit with the second iMinuit fit (see Table 2), we remark that the positron and proton fit only worsen slightly (∆χ 2   The left panels show the default model (z = 5.6 kpc), and the right panels show the small-height model (z = 3 kpc). The AMS-02 data is given in red for positrons [5] and protons [57,58]. The total flux of our model is presented in blue (solid), primary components in orange (dashed) and secondary components in green (dashed-dotted). The fractional residuals between the data and our model are given in the bottom of each panel, showing that our models for positrons and protons fit the data to a few percent accuracy for both halo heights.  Table 5. List of free parameters of our Multinest fit for the default model (z = 5.6 kpc). We give the range in which we let the parameters float, the type of prior, the best-fit value, the posterior and the posterior uncertainty. This model provides a good fit to the data with a total χ 2 of 61.54 over 101 degrees of freedom.   Table 7. Summary of the χ 2 -values of the Multinest fit for the default model (halo height 5.6 kpc, top rows) and the small-height model (halo height 3 kpc, bottom rows) for the total flux, as well as the individual positron and proton flux. All χ 2 per degree of freedom are below 1, indicating very good agreement of our astrophysical model with the AMS-02 data.

B Positron Flux with Dark Matter Contribution
In Figure 5 we show an example of the local positron flux as obtained from our simulations with a contribution from 5-GeV DM annihilating into the e + e − for our default model.