Antiproton bounds on dark matter annihilation from a combined analysis using the DRAGON2 code

Early studies of the AMS-02 antiproton ratio identified a possible excess over the expected astrophysical background that could be fit by the annihilation of a weakly interacting massive particle (WIMP). However, recent efforts have shown that uncertainties in cosmic-ray propagation, the antiproton production cross-section, and correlated systematic uncertainties in the AMS-02 data, may combine to decrease or eliminate the significance of this feature. We produce an advanced analysis using the DRAGON2 code which, for the first time, simultaneously fits the antiproton ratio along with multiple secondary cosmic-ray flux measurements to constrain astrophysical and nuclear uncertainties. Compared to previous work, our analysis benefits from a combination of: (1) recently released AMS-02 antiproton data, (2) updated nuclear fragmentation cross-section fits, (3) a rigorous Bayesian parameter space scan that constrains cosmic-ray propagation parameters. We find no statistically significant preference for a dark matter signal and set strong constraints on WIMP annihilation to bb̅, ruling out annihilation at the thermal cross-section for dark matter masses below ∼ 200 GeV. We do find a positive residual that is consistent with previous work, and can be explained by a ∼ 70 GeV WIMP annihilating below the thermal cross-section. However, our default analysis finds this excess to have a local significance of only 2.8σ, which is decreased to 1.8σ when the look-elsewhere effect is taken into account.


Introduction
The Alpha Magnetic Spectrometer (AMS-02) has produced high-precision measurements of the local cosmic-ray antiproton spectrum, which has long been considered a critical channel for indirect dark matter (DM) searches [1,2].Intriguingly, early studies in 2016 [3,4] found statistically significant (∼ 5σ) evidence for an excess over the expected astrophysical background in the cosmic-ray antiproton spectrum at an energy of ∼10 GeV that was best-fit by the annihilation of DM particles with a mass of ∼60 GeV into b b final states.
Compared to the much-discussed cosmic-ray positron excess [5][6][7] first observed by PAMELA [8], and later by AMS-02 [9], the DM interpretation of the antiproton excess has a stronger theoretical motivation due to its consistency with a standard thermal WIMP.On the other hand, the suggested antiproton excess is more challenging to identify because the signal-to-background ratio is much smaller (≲ 10%).
Due to the small relative amplitude of the excess, it is imperative to correctly account for systematic errors.Recent analyses have pointed out a number of effects that individually appear capable of decreasing, or entirely eliminating, the statistical significance of the excess, including: (1) uncertainties in cosmic-ray propagation and secondary production [10][11][12][13][14][15], (2) uncertainties in the production cross sections of astrophysical antiprotons as well as other nuclear production cross-sections that are used as inputs for the cosmic-ray propagation codes [16][17][18][19][20], (3) correlations in the systematic uncertainties in the AMS-02 spectral data, which are undoubtedly present, but have not been publicly quantified by the AMS-02 collaboration [10,21,22].While most studies agree that there is a positive residual in the antiproton flux at an energy of ∼10-20 GeV, the local significance of this excess spans from ∼ 7σ [23] all the way down to a negligible ∼ 1σ [24].
Recently, a new antiproton dataset was released by the AMS-02 collaboration, which included updated antiproton data spanning the period 2011-2018 [25].A few recent studies have analyzed the antiproton spectrum using this dataset, generically finding that the size of the residual has decreased, and has also moved towards a higher antiproton energy.Thus, these analyses have found a statistically insignificant excess in the range of 0.5σ to 2σ [22], with best-fit DM masses in the range of 80 GeV to 120 GeV.
In this paper, we utilize the DRAGON2 code in order to fit the new antiproton measurements from AMS-02, employing an advanced Bayesian analysis that constrains our fit to the antiproton data through simultaneous fits to the secondary spectra of B, Be, and Li.We include the halo height as a free parameter in our fit, along with the other propagation parameters and nuisance parameters that allow us to account for cross sections uncertainties.This comprehensive analysis improves on previous work due to the number of constraints (datasets) employed, the use of refined cross section parametrizations, the treatment of possible correlations in AMS-02 errors and the energy dependence of antiproton cross-section uncertainties.
We remark that there is no previous analysis which fits the different spectra (and ratios to primary particles, such as C and O) of B, Be and Li simultaneously with antiprotons.Recently, efforts in this direction were made in Ref. [22], where the authors map the probability distribution of the halo height -obtained from a previous analysis of Li/C, Be/B and B/C -into their antiproton fit, however keeping the remaining propagation parameters fixed.Furthermore, Ref. [22] is based on a slightly simplified semi-analytic propagation model as opposed to the fully numerical approach implemented in DRAGON2 and we make use of Bayesian inference.
Focusing on DM annihilation to b b final states, we find no statistically significant evidence for a DM annihilation component.The largest local significance for an excess stands at 2.8σ for a DM mass of ∼ 66 GeV.On the one hand, this mass range is consistent with an interesting region of DM models that are also capable of explaining the Galactic Center Excess (GCE) [26].On the other hand, the corresponding global significance the corresponding global significance, which takes into account the look-elsewhere effect, is only 1.8σ.Utilizing these results, we place strong constraints on the DM annihilation cross-section.
In Section 2, we describe our simulation set-up and the methodology we employ in this work.We also compare how our set-up and analysis differs from other recent works in Section 2.3.Then, in Section 3.1 we report our fits to the antiproton spectra for the scenario where we only consider the production of antiprotons from CR interactions.We discuss and show our results for the combined analyses including DM antiproton production in Section 3.2, where we also derive DM bounds.Finally, we summarise and discuss our main findings in Section 4.

Methodology
In this section we summarise the main points of the simulation procedure and analysis setup.We note that our procedure is close to that described in Ref. [27], and we refer the reader to that paper for further technical details.

Simulations setup
DRAGON21 is an advanced CR propagation code designed to self-consistently solve the diffusionadvection-loss equation describing CR transport for all species involved in the CR network, including cosmic rays (CRs) of both astrophysical and exotic origin (e.g., from DM annihilations/decays).The transport equation features fully position-dependent and energydependent transport coefficients in either two-dimensional (assuming cylindrical symmetry) or three-dimensional configurations of the Galaxy structure (i.e.considering the spiral arm pattern of the Galaxy).The flexibility of the DRAGON2 code allows for detailed studies of both small-scale and large-scale structures (e.g., the spiral structure of the Galaxy) in steady-state and transient mode, refining the spatial resolution on the regions of interest (e.g., local bubble, GC, or Galactic Plane).
In this work, we solve the coupled system of propagation equations for all the isotopes involved in the CR network [28] up to Z = 14 (Silicon) with a customised version of the DRAGON2 code [29,30] that is optimized for studies of antiprotons and antinuclei2 .For these analyses, we solve these equations assuming cylindrical symmetry in the Galaxy, which makes use of a homogeneous gas distribution implemented according to Ref. [31].We describe the distribution of CR sources (expected to be mainly supernova remnants -SNRs) [32][33][34] using the spatial distribution derived by Ref. [35].We parameterise the source injection as a broken-power law in energy per nucleon.
We set the diffusion coefficient to be spatially independent, with a rigidity dependence described by a broken power-law modeled as: where we fix R 0 = 4 GV, the value of the rigidity break R b = 312 GV, the change in the spectral index ∆δ = δ − δ h = 0.14, and the smoothing parameter s = 0.040, following the values determined in Ref. [36] from a detailed analysis of the proton and helium AMS-02 fluxes.
We consider reacceleration, which is directly related to the spatial diffusion coefficient, and set the Alfvén velocity as a free parameter in our study.However, we neglect convection, since analyses of secondary-to-primary ratios at AMS-02 do not favor any relevant impact of convection on CR propagation [13,[37][38][39].
We calculate inelastic cross-sections using the CROSEC parameterisation [40], but utilize the DRAGON2 inclusive spallation cross-sections for the production of secondary nuclei (Z > 1) [11,41].We note that there have been significant efforts over the last few years dedicated to improving the current cross section parameterisations and to account for their uncertainties in the analysis.Other commonly considered parameterizations include those from GALPROP [42,43] or FLUKA [19] (see also the recent work by [44]).However, the DRAGON2 cross sections have recently been demonstrated to provide a better simultaneous fit to the B, Be and Li flux ratios, compared to other cross sections data-sets [39].For a comparison between different cross section models for some of the main channels involving the production of B, Be and Li, we refer the reader to Appendix A.
We utilize cross sections for antiproton production from CR proton and helium collisions following the calculations in Ref. [45] ("Winkler" option in the DRAGON2 code).In addition, we take into account the tertiary contribution for production of protons and antiprotons.In this work, we also account for the production of p particles from nuclei heavier than He.This is done by scaling the proton spallation cross sections by a factor of A 0.9 [17], where A is the mass number of the CR nucleus colliding with ISM gas.This results in a slight increase (∼3%) in the p production rate with respect to our previous work [27].
Additionally, in our analysis, we do not consider the p/p data below 4 GV, in order to limit the effect of the uncertainties in solar modulation in our results.This cut leads to more conservative limits for low DM masses (below ∼ 100 GeV).In our calculations we implement solar modulation using a modified Force-Field approximation [46] which accounts for the effects of charge-sign dependence [47].In this case, the Fisk potential takes the following form: parameterising its rigidity-dependence as term accounts for the additional energy loss of positively (negatively) charged particles during a negative (positive) solar polarity phase.This additional energy loss occurs because particles with charge opposite to the magnetic polarity access the heliosphere by complicated inward drift along the heliospheric current sheet, while particles whose charge aligns with the polarity can enter on direct trajectories along the poles.Since AMS-02 data were mostly taken during a positive solar polarity phase we set ϕ + 1,AM S−02 = 0 and only consider a non-vanishing ϕ − 1,AM S−02 .Neglecting the short time-period of opposite polaritythe polarity flip occurred near the beginning of AMS-02 operation (around 2012) -should not affect our analysis, in particular as we only include the p/p data above 4 GV.
A reasonable value of the constant Fisk potential, ϕ 0 , was found to be 0.61 GV for the period of 2011-2016 (relevant for the AMS-02 data-sets for the secondary-to-primary ratios, C, O and He) [48], since it allows us to reproduce Voyager-1 [49,50] data and AMS-02 data in that period, while also remaining consistent with the NEWK neutron monitor experiment3 [51,52] (see also [53,54]).In order to adjust the value of ϕ 0 to the data taking period 2011-2018 (relevant for the newly released antiproton data-set [25] as well as for the primary CRs Ne, Mg, Si [55] and the most recent H and He data sets) we consider the time-variation in the neutron intensity between 2011 and 2018 as detected by NEWK [51][52][53][54].In this way we obtain ϕ 0 = 0.58 GV for the 2011-2018-period.Because charge-dependent effects (mainly linked to propagation within the heliospheric current sheet) are necessary to explain the different behaviour of electrons and positrons [56], we allow that negative CR particles can have a non-vanishing ϕ + 1 , finding a good agreement with the data for antiprotons using ϕ − 1 around 0.9 GV, consistent with Ref. [16,21,27].
We stick to the more conventional setup for our analysis, although we note that we neglect other possible astrophysical processes that could significantly affect the propagation of CRs and their production.For example, we neglect the spatial dependence of the diffusion coefficient, while there are both theoretical arguments and γ-ray observations of the diffuse emission (see, e.g.[19,[57][58][59]) supporting a change of the diffusion coefficient in different zones of the Galaxy.It was shown that models that consider possible differences in the transport of CRs in the galactic halo and disk predict only a slightly different antiproton spectrum measured at Earth [60,61].Additionally, the production of secondary CRs within the acceleration region of supernova remnant can change the spectra of secondary CRs [62][63][64][65].These modifications of the standard scenario of Galactic CR transport could affect the astrophysical searches for DM signatures and will be explored in a future work.

Analysis setup
As in our previous work [27], our set-up (similar to the one tested in Ref. [39], where more details are given) relies on a Markov chain Monte Carlo (MCMC) analysis that is based on Bayesian inference [66].We determine the propagation parameters involved in the transport equation (namely H, D 0 , V A , η and δ) from a combined fit to AMS-02 spectral data for the main secondary CR nuclei (B, Be and Li) along with the p spectra.Many previous studies [11-13, 19, 20] have shown the importance of combining different secondary CR observations in order to mitigate the effect of systematic uncertainties (primarily those associated with the production cross-sections of secondary CR) on the propagation parameters of the simulation.Our fitting procedure includes the injection parameters of the primary CRs included in our simulation set-up ( 1 H, 4 He, 12 C, 14 N, 16 O, 20 Ne, 24 Mg and 28 Si).In this work, we fit the B/C, B/O, Be/C, Be/O, p/p flux ratios and the 10 Be/Be and 10 Be/ 9 Be flux ratios (allowing us to constrain the halo height, H [67]), as well as the Li/B flux ratio.We neglect correlations between the nuclear CR data sets (for instance between B/C and B/O) and treat them as statistically independent.This is justified since such correlations only slightly affect the uncertainties in the propagation parameters, and have almost no impact on the DM bounds we compute in this work.
The Fisk potential is adjusted for the period of data collection (2011-2018), as described above.Moreover, in this analysis we include nuisance parameters that allow us to modify the normalization of the spallation cross sections, since spallation cross uncertainties above a few GeV/n mostly affect the overall normalization [11,17,18,68].Extending the analysis of our previous work, a scale factor, that allows us to apply a constant offset in the production rate of antiprotons, has also been added as a nuisance parameter for the p cross sections, so that there are four total nuisance parameters (scale factors) included in this analysis: S B , S Be , S Li and S Ap .These factors allow the model to further adjust the production of these secondary nuclei and are necessary to alleviate the tension between the grammage inferred from analyses of secondary-to-primary ratios of secondary CR nuclei (B, Be and Li) and the grammage predicted by analyses of the antiproton-over-proton flux ratio, which reflects the impact of cross sections uncertainties in the production of these secondary particles, as discussed in Ref. [27].
An important point of this analysis is that we include, in the likelihood definition, a penalty factor associated to each of our nuisance scale factors (see more details in Ref. [39]).In practice, this approach penalizes large variations from the original cross sections, and can be interpreted as an implementation of existing constraints on production cross-sections of B, Be, Li and p.These penalty factors are inversely proportional to the relative variance of the cross section data for the channels that we consider (similar to what is defined in Eq. 4 of Ref. [13]).In this case, for B, Be and Li we consider the variance in the cross section data for the channels of production of each of their isotopes ( 11 B, 10 B, 10 Be, 9 Be, 7 Be, 7 Li, 6 Li) with C and O colliding with p [11].
For the p production cross sections, the amount of experimental data is more limited, and we only include the variance for the data of direct p production from p+p collisions 4 .
Given the difficulty of assessing the systematic uncertainties related to antiproton production we adopt two different analyses that differ in how we treat antiproton cross sections: A simplified analysis, where there is only a scale factor for the antiproton cross sections and its penalty factor has an associated variance of 12% (we call this analysis the "Simplified" analysis) and a more advanced analysis based on Ref. [69] where we include a covariance matrix for the energy dependence of antiproton cross section uncertainties and the penalty associated with the scale factor has a variance of 6.6%.We call this the "Canonical" analysis.In both cases, scaling factors for the cross sections of CR secondary nuclei (namely B, Be and Li) are included.
For both analyses, we study scenarios where we include the contribution of a WIMP annihilating into b b final states that subsequently produces p.We then test whether a DM contribution is preferred and derive constraints on the DM mass and annihilation cross section (⟨σv⟩) along with the rest of transport parameters and nuisance factors.The annihilation yield at production (dN/dE) for every DM mass is computed using the PPPC4DM package [70,71].We use the Navarro-Frenk-White (NFW) DM density profile [72] as a reference with a DM density of ρ ⊙ = 0.4 GeV at the Solar System position, r ⊙ = 8.3 kpc [73].
In our procedure, we incorporate correlations in the AMS-02 systematic errors through covariance matrices computed in Ref. [21], which primarily affects the statistical significance of the reported anomaly, and only slightly changes the diffusion parameters obtained in our fits and the (best-fit) predicted spectra of secondary CRs.We remind the reader that these matrices are not released by the AMS-02 collaboration and only constitute best estimates obtained with the publicly available information.Therefore, our main tests will consist of the Canonical and Simplified analyses considering correlations in AMS-02 systematic errors and, for comparison, we also present the results obtained assuming no correlation in AMS-02 errors.
In order to reduce the impact of uncertainties in solar modulation, we do not consider in the analysis the p ratios below 4 GeV, as explained above.As tested in our previous work, analyses that use alternative energy cuts at 3 GeV or 5 GeV do not significantly change the results.We also test that small variations on the parameter ϕ − 1 do not affect our results.This confirms the robustness of our analysis with respect to solar modulation uncertainties.

Comparison with other recent analyses of the 2018 p data-set
In addition to our previous work, several other groups have analyzed the new AMS-02 p data set over the last two years [22,69,[85][86][87].We briefly review the analysis choices of each of these studies here: In the analyses of Kahlhoefer et al and Balan et al (Refs.[85] and [87]), the authors derive bounds on the DM annihilation rate by a joined fit to the hydrogen and helium spectra and the p/p ratio.The focus of these works is, however, on the implementation of machine learning techniques into CR analysis, rather than on a fully realistic modeling of all relevant systematics (which would require, for instance, the inclusion of additional CR species to constrain the halo height).They show how different diffusion setups (i.e.different ingredients, such as convection or reacceleration, and assumptions on the injection and breaks in the diffusion coefficient) lead to significantly different bounds.
In the analysis of Di Mauro and Winkler [69], they derived bounds on DM annihilation from a joined fit to the AMS-02 antiproton and B/C data.While the modeling of cross sections and systematic uncertainties is similar to our work, Ref. [69] considered a few fixed choices of the halo height rather than deriving it within the fit (which would require the inclusion of radiactive species as in our work).Furthermore, Ref. [69] employs the semi-analytic two-zone diffusion model for CR propagation as opposed to our fully numerical approach with DRAGON2.In Calore et al. [22], the authors derive DM limits by a fit to the AMS-02 antiproton spectrum.Systematic uncertainties in the antiproton cross sections, as well as AMS-02 error correlations are included through covariance matrices similar to Ref. [69] and our work.In order to include the halo height as a fit parameter the authors extract its probability distribution from an earlier CR study on secondary-to-primary ratios.Uncertainties on the remaining propagation parameters are included in a simplified way through a covariance matrix.In contrast to our work, Ref. [22] employs the two-zone diffusion model for CR propagation.This analysis is mainly based on the 2015 p data-set, while we analyze the 2018 dataset, which includes nearly a factor of two improvement in the exposure.
More recently, Ding et al. [86] derived DM bounds from a combined fit to AMS-02 protons and antiprotons (2011-2018 dataset), using the GALPROP code, where the propagation parameters were fixed to the ones that reproduce the B/C ratio.The antiproton cross sections are those from the GALPROP v.54 code and neither cross section uncertainties nor other systematic errors are taken into account.
Our study contains a number of improvements compared to these earlier works.In particular, we perform the first combined fit to the AMS-02 p/p data and a large set of secondary-to-primary and secondary-to-secondary CR ratios.In this way we are able to fully include all relevant propagation uncertainties in our analysis.Specifically, we directly constrain the halo height H, which is a key parameter in dark-matter studies, within our fit.This is in contrast to previous works which, at best, considered some fixed values of H (with the exception of Ref. [22] which incorporated a simplified probability distribution for H).Another new ingredient are the improved DRAGON2 spallation cross sections which resolve some tensions in the fits to nuclear cosmic-ray spectra and allow for a more accurate determination of the CR propagation parameters.Furthermore, we perform a state-of-the art modeling of the relevant cross section uncertainties and of the error correlations in the AMS-02 data -which is missing in some of the previous works.In this light, our analysis represents the most advanced study of the AMS-02 2018 antiproton data set performed so far.
3 Results of the combined antiproton analyses

Secondary antiproton spectrum
Here, we report and compare the results from the combined antiproton analyses without including DM contributions to the antiproton spectrum.In Figure 1, we show the spectra of the main secondary flux ratios (p/p flux ratio, in the left column, and ratios of B/O, Be/O and Li/O, in the right column) derived from the best-fit propagation parameters7 in the scenario of antiproton production solely from CR interactions.In particular, we show p/p spectrum (left panels) and secondary-to-primary flux ratios of B/O, Be/O and Li/O (right panels) evaluated with the transport parameters obtained in our Canonical analyses in a scenario that does not include any DM production of antiprotons.The top row shows the spectra obtained from the analysis including no correlation in the AMS-02 systematic errors while the bottom row shows the spectra for the case when we include correlations.The variation in the flux ratios due to the uncertainty in the determination of the propagation parameters is shown (not including modulation uncertainties) in each case and the error bars on AMS-02 data are the 1σ errors reported by the collaboration.The residuals in the left panels are drawn with respect to AMS-02 data (2018).
the results obtained from the Canonical analysis (where we include the covariance matrix for antiproton cross sections uncertainties) in the case where AMS-02 errors are treated as uncorrelated (top panels) and correlated (bottom panels).The spectra obtained assuming correlated AMS-02 errors from the Canonical analysis and the Simplified one (where we only include a cross sections scaling factor) yield roughly the same spectral fits, although with the Simplified analysis we obtain larger cross sections scale factors.
In the analyses where we assume that the AMS-02 data errors are correlated, the quality of the fit is either χ 2 = 476.1 (in the case where we only include scaling factors for cross sections, Simplified analysis), or χ 2 = 553.4(in the case where the cross sections covariance matrix is included, Canonical analysis).In both cases, there are 329 degrees of freedom.Meanwhile, we obtain a χ 2 ∼ 300 in the uncorrelated cases.The low χ 2 /d.o.f obtained from the uncorrelated case points to the fact that we are overestimating the uncertainties in AMS-02 data if we assume that errors are uncorrelated [10,88].Indeed the AMS-02 measurements are dominated by systematic uncertainties that have been shown to be highly correlated [21].However, given the difficulties on estimating all the factors accounting for these systematic errors in the AMS-02 detector, our modifications to the χ 2 can only be thought of as estimates based on publicly available information.
The probability distribution functions (PDFs) obtained for every parameter are shown in Appendix D. We find very flat residuals in the p/p ratio with respect to the 2018 data-set, such that the p spectra are easily reproduced with a constant scaling of their production cross sections.This result is similar to both our earlier work and that of Ref. [69].For the Simplified analysis, we find a scaling factor of S Ap = 1.11±0.02(a ∼ 11% increase in the cross sections of p production) and S Ap = 1.14 +0.01 −0.02 for the uncorrelated and correlated cases, respectively.We highlight that these scale factors are compatible with the variance associated to the p cross sections data in this analysis (which is ∼12% [17]) and are similar to those obtained in the analyses by Ref. [89].In our Canonical analysis (where we include the covariance matrix for the antiproton cross sections), we find a best-fit scaling of S Ap = 1.07 ± 0.01 for the correlated case, which is still compatible with the 1σ cross-section uncertainties reported by Ref. [69] and a S Ap = 1.08 ± 0.01 for the uncorrelated case.The scaling associated with Be and Li is found to be within 8%, indicating that the cross sections of B production are compatible with very small or no scaling in both kind of analyses.
The PDFs of the inferred propagation parameters in the scenario without DM annihilation are shown in Figure 7, left panel, where we overlap the PDFs obtained in each of the analyses, including correlations in AMS-02 data.The transport parameters found in both analyses are compatible within 1σ, which means that the determination of the propagation parameters is not significantly affected by the cross sections scaling model.
The halo height inferred in these analyses is H ∼ 4 kpc and the ratio D 0 /H is close to 10 28 cm 2 s −1 kpc −1 .The effective Alfvèn velocity is in the range V A ∼ 10 − 17 km/s in the uncorrelated analyses and V A ∼ 15 − 25 km/s in the analysis that includes correlations in the uncertainties of the cosmic-ray data.Finally, the spectral index of the diffusion coefficient ranges from 0.45 to 0.48, compatible with a Kraichnan-like [90,91] turbulence spectrum of the plasma fluctuations below the break (∼ 300 GV) and Kolmogorov-like [92] above.This result is in good agreement other recent studies [12,13,93].
These results confirm our earlier findings that current diffusion models allow us to simultaneously reproduce all the light secondary CRs (see also Ref. [19]) within the context of a model that assumes only secondary production in the interstellar medium and a constant diffusion coefficient throughout the Galaxy.We note that if we decreased the uncertainties on the cross sections of production of secondary CRs, our results could clearly show tensions between the parameters needed to reproduce the different secondary CRs, therefore, future accelerator experiments could provide measurements that would help us to spot tensions and reveal signatures of new physics in astrophysical data.A summary of the propagation parameters inferred is also given, in form of box plots, in Figure 9 and several examples of the best-fit DRAGON2 input files are available at github.com/tospines/Customised-DRAGON-versions/tree/main/Custom_DRAGON2_v2-Antinuclei.
Finally, we remark that the fact that the energy dependence of the different secondary ratios agrees so well indicates that our simple diffusion scenario does not require significant modifications.However in order to further improve the search for new physics signatures, it will be necessary to also explore the implications of inhomogeneous diffusion in the Galaxy (which is suggested by some observations [57,94,95] and supported by theory [58,61]).We leave this analysis for future work.

Search for WIMP features in the antiproton spectrum
In Figure 2 we show the best-fit p/p spectrum obtained in our Canonical analysis in the case where we do not include (left panel), or do include (right panel) correlations in the AMS-02 systematic errors.The best-fit DM contribution to the p/p spectrum is lower in the case where we include error correlations.The best-fit mass and annihilation rate obtained in each analysis, as well as the local significance associated with the DM contribution, are reported in Table 1.
Our results highlight the importance of including the correlations in the AMS-02 systematic errors when assessing the significance of a potential DM signal.If we ignore the correlations, we find a relatively large local excess of 3.5 − 3.8σ, compatible with early analyses [3,4,23].However, the local significance is reduced to 2.8 − 3σ once we account for the systematic error correlations in AMS-02 data.This general trend agrees with the findings of other recent studies [16,21,22,69].In our uncorrelated analyses, the best-fit mass of the DM particle ranges from 51 GeV to 100 GeV when we include the covariance matrix for antiproton cross sections and when we only include a scaling factor, respectively.The best-fit annihilation rate is slightly below the thermal relic value for the Canonical analysis (⟨σv⟩ = 0.89 ± 0.18 × 10 −26 cm 3 /s) and compatible with the thermal relic value in the Simplified analysis (⟨σv⟩ = 2.01 +0.4 −0.35 × 10 −26 cm 3 /s), (see Table 1).We have, furthermore, computed the global significance of the excess by creating mock data samples to determine how often a local 2.8σ excess randomly occurs within the considered mass range (the look-elsewhere effect). 8We obtain a global significance of only 1.8σ, meaning that there is no statistically significant preference for a DM contribution in agreement with [16,21,22,69].In the analyses that include AMS-02 error correlations, we get a best-fit mass of around 66 GeV for an annihilation rate, ⟨σv⟩ ∼ 10 −26 cm 3 /s, somewhat below the thermal relic one.These values are quite similar to those found from the analysis of B/C and the p/p ratio (2015 AMS-02 dataset) using GALPROP in Ref. [21].
Other factors that can impact our results are the cross sections parameterizations that are employed, the evaluation of the correlation matrices for AMS-02 systematic errors, the low energy parameterization of the diffusion coefficient, the set of cosmic-ray observables included in the analysis, and the spatial dependence of the diffusion coefficient.For example, it was found by Ref. [21] that including the antiproton flux instead of the p/p ratio tends to decrease the significance of any antiproton excess.
It is also worth noting that, in scenarios where we include the contribution from WIMP annihilation into b b final states, we still need to scale the cross sections of antiproton production by ∼ 6 − 7% (and as much as 14% in the case that does not include the covariance matrix for antiprotons cross sections).However, the difference in the grammage needed to explain B, Be and Li compared to the one to explain antiprotons is reduced and the cross section scale factors for B, Be and Li are closer to unity (i.e.no scaling).In general, we find that the best-fit cosmic-ray propagation parameters determined in models that include DM annihilation are not very different from those inferred in scenarios that include only p secondary production.This means that the determination of these parameters is mainly driven by the secondary ratios of B, Be and Li.A summary of the propagation parameters, scale factors, WIMP mass and ⟨σv⟩ found in these analyses is given in Figure 10.Future improvements in our understanding of the transport of CRs at low energy and solar modulation will allow us to pin down the astrophysical fluxes at low energy even more precisely and to further improve the sensitivity of our analysis to low-mass DM .

Dark matter bounds
Here, we report 95% confidence upper limits on the dark matter annihilation cross-section, ⟨σv⟩, from our canonical analysis and compare our results to those of other groups.The limits are obtained using the same MCMC analysis described above, with the DM cross-sections computed by logarithmically scanning the WIMP mass in 42 mass bins, from 10 to 1500 GeV.These bounds are shown in Fig. 3, where the left panel displays the limits obtained in our analysis including error correlations in the AMS-02 data, in comparison to the ones obtained by Refs.[22,85].Over most of the considered mass range our limits fall between those of Ref. [22] and of Ref. [87].The differences compared to those previous works arise due to the advancements of our CR analysis, described in Section 2.3.We finally remark that these limits allow us to rule out the thermal relic cross sections for masses below ∼ 200 GeV in the b b-channel.
In the right panel of Figure 3, we show the bounds we obtain for dark matter profiles that follow a contracted NFW profile (favored by the DM interpretation of the galactic center gamma ray excess) with γ = 1.2 and a scale radius of 20 kpc, compared to those obtained from analyses of the antiproton cross sections in Ref. [69].The 2σ contour represents the best-fit DM candidate reproducing the GCE derived by Ref. [96] and the constraints derived from dwarf spheroidal galaxies in Ref. [97].In this plot, we show our upper limits on the annihilation rate for our Canonical and Simplified analyses, including correlations in the AMS-02 errors.
As we can see from this comparison, the DM interpretation of the GCE is in tension with respect to these antiproton analyses for the bb channel (even in the uncrorrelated case) as was shown also in Ref. [69].In fact, the best-fit candidate found by Ref. [96] (which is one of the analyses leading to a lower annihilation rate for the DM candidate fitting the GCE) is around a .Left panel: 95% confidence Upper limits on ⟨σv⟩ derived from our main analyses with a NFW profile, compared to those obtained by Refs.[22,87].Right panel: Comparison of the upper limits obtained for a contracted NFW profile (see details in the text) with other limits derived for the same profile.In particular, we compare our limits with the 2σ contour representing the best-fit DM candidate reproducing the GCE in Ref. [96], the constraints derived from dwarf spheroidal galaxies from Ref. [97] and those from Ref. [69] for a contracted-NFW profile with a halo height of H = 4 kpc (and the scale radius adjusted to 20 kpc).
factor of 3 above our 95% confidence limits.However, we remark that the exploration of other annihilation channels and configurations of the DM profile (mainly contraction index and scale radius), which is rather uncertain close to the Galactic center, is needed to clearly state that the DM explanation for the GCE and the AMS-02 antiproton data are inconsistent (for example, Ref. [69] found that both are compatible in the µ ± channel).Moreover, uncertainties relating to inhomogeneous cosmic-ray propagation, e.g., a strong convective wind near the galactic center, could change the relative intensity of the γ-ray and antiproton signals.Then, we show the limits from Ref. [97], which were shown to be robust to different DM profiles.This also allows us to illustrate the fact that, when using a contracted NFW, antiproton measurements can certainly be the leading observable to constrain the existence of DM.The bounds obtained in the uncorrelated cases, as well as other analyses, are discussed and shown in App. C. In particular, we show the impact of considering uncorrelated AMS-02 systematic uncertainties in the DM bounds by comparing with those including correlations in Fig. 6 (left panel).
In view of these results, we conclude that, even with the simple propagation setup generally assumed, the spectra of the different secondary CRs agree with pure secondary production, and do not require any additional antiproton production from WIMPs.These analyses allow us to limit the contribution from a WIMP particle of mass below ∼ 1 TeV to be, at most, a ∼ 15% of the secondary contribution, for the bb channel.This highlights the great need to further improve cross sections through accelerator measurements in order to (along with the development of a more robust astrophysical modeling) improve our searches for small signals on top of large backgrounds in the antiproton spectrum.For instance, our rescaling factors on nuclear and antiproton cross sections should be revised with future accelerator data.
Until the uncertainties in the relevant CR cross sections and in the CR propagation are significantly reduced it will be difficult to unambiguously identify a DM signal with CR antiprotons -even with more precise AMS-02 data.Therefore, it is of great importance to Correlated AMS-02 errors Uncorrelated AMS-02 errors M W IMP (GeV) ⟨σv⟩ (10 −26 cm 3 /s) local σ M W IMP (GeV) ⟨σv⟩ (10

3.8
Table 1.Best-fit parameters characterizing the WIMP properties, namely mass and annihilation rate, along with the local significance found for the WIMP contribution.In particular, the median of the PDF and the 1σ error are shown.The reported error corresponds to the 1σ uncertainty obtained from the PDF fixing the other parameters to their best-fit values.
investigate DM signals in complementary CR channels with lower astrophysical backgrounds.
In a companion paper, we make use of our antiproton results and update the predictions for antinuclei (antideuteron and antihelium) detection, using newly derived cross sections for their production from CR interactions and WIMP annihilation and discussing the implications of the updated predictions in view of the preliminary detected events reported by the AMS-02 collaboration.

Discussion and conclusions
Over the last decade, many studies have investigated the spectra of CR antiparticles in order to search for hints of new physics.Although there have been indications of significant anomalies in the p spectrum that indicate a possible signal from DM, no clear evidence has yet been uncovered.In this work, we analyzed the spectra of the light secondary CR species B, Be, Li in combination with p in a scenario where antiprotons can also be produced by annihilation of a WIMP into b b final states, completing and extending our previous work.This constitutes the current most complete derivation of DM bounds from AMS-02 antiproton data, characterized by an analysis that combines antiprotons with the rest of light secondary CRs, uses the updated DRAGON2 cross sections for spallation interactions, and accounts for uncertainties in the cross sections parametrizations and possible correlations in AMS-02 data, to search for WIMP signals in the antiproton spectrum.We find that the antiproton spectrum is compatible with pure astrophysical production and obtain a local significance of ≲ 3σ for any additional antiproton component from DM annihilation.This result is obtained within our main analysis that includes correlations between the AMS-02 systematic errors.We have tested that this corresponds to a global significance of < 2σ when accounting for the trials factor from the scan of DM masses, indicating that there is no strong statistical preference for a DM contribution.
Regarding the differences obtained considering and neglecting correlations in the AMS-02 errors, we observe that the significance in our main analyses vary by around 0.7σ, being of ∼ 3.5 − 3.8σ (local) in the uncorrelated case.This shows the importance of including correlations in these analyses in order to prevent from overestimating a potential DM signal.
We have reported DM bounds and showed the effect of treating cross sections uncertainties in different ways, finding that our bounds are compatible with other recent analyses of the 2018 antiproton data-set.In particular, our analyses allow us to consistently rule out the thermal relic cross sections for WIMP masses below ∼ 200 GeV.The mass of the DM particle that produces the largest improvement in the antiproton spectrum is 66 GeV, which agrees well with previous works, and is also interesting because it closely corresponds to the best-fit mass for DM explanations of the galactic center γ-ray excess.However, the bounds on the annihilation rate obtained in our analysis are in tension with the best-fit cross-sections for DM explanations of the GCE in the bb channel.
Let us remind the reader that our analysis takes into account cross section uncertainties (on the secondary production of antiprotons and nuclei) through overall scaling factors and also includes energy-dependent uncertainties on antiproton cross sections through covariance matrices.Given the complexity of our fit, which involves a much larger set of CR species compared to previous analyses, incorporating the energy scaling of cross section uncertainties for each isotope of B, Be and Li goes beyond the scope of this work.We note, however, that the significance of the excess is expected to further (slightly) reduce once the energy-scaling of cross section uncertainties for these secondary CRs is taken into account (see discussion in Refs.[16,21,22,69]).A more detailed investigation of this aspect is left for future work.
Finally, we note that it would be interesting to generalise our analyses to account for possible spatial dependence of the diffusion coefficient, something that could affect slightly both, the spectra of particles produced by CR interactions and the DM signal at Earth.This study is left for a next work.Additionally, in a companion paper, we will use our constraints to re-evaluate the expected fluxes of antideuteron and antihelium at Earth.The version of the code employed in this and our companion work is publicly available at github.com/tospines/Customised-DRAGON-versions/tree/main/Custom_DRAGON2_v2-Antinuclei.

A DRAGON2 cross sections for spallation interactions of B, Be and Li
In this appendix, we show different models of direct cross sections for spallation reactions used in the literature compared to experimental data in some of the most important reaction channels for the production of the light secondary CRs B, Be and Li.Experimental data are taken from various experiments and authors (see Refs. [30,41]): EXFOR (Experimental Nuclear Reaction Data) 9 , the GALPROP cross-section database (isotope_cs.dat)and from various publications and experiments (Bodemann1993, Davids1970, Fontes1977, Korejwo1999, Korejwo2002, Moyle1979, Olson1983, Radin1979, Read1984, Roche1976, W90, W98a and Zeitlin2011).It is important to remind that the FLUKA cross sections are computed by means of the FLUKA code [19] while the other models constitute different parameterizations fitted to data existent at the time of release [98][99][100].
Full details about the different parameterizations can be found in Refs.[11,19,20,42] and references therein.The band of statistical uncertainties is related to the determination of the FLUKA cross sections (see Ref. [19] for more details).

B Other observables included in the analyses
In this appendix we report the predicted fluxes and ratios evaluated from the best-fit parameters obtained in our canonical analysis.In the top panels of Fig. 5 we show a comparison of the different secondary-to-secondary rations of B, Be and Li evaluated from the parameters obtained assuming correlations in AMS-02 errors (top-left panel) and without correlations (top-right panel), compared to the AMS-02 data.The 2σ uncertainty from the determination of the diffusion parameters is also reported as a band for each ratio not including modu-lation uncertainties.As we see, in both kind of analyses, the these models achieve a good reproduction of the experimental data, even at the 1σ level.
In the bottom-left panel of this figure, we show the fluxes of H and He compared to AMS-02 (for the modulated fluxes) and Voyager-1 data (for the unmodulated fluxes) as well as the 2σ uncertainty bands from the determination of the propagation parameters and the Fisk potential (∆ϕ 2σ ∼ 0.12).The obtained flux of these primary CRs is roughly the same in every analysis since the injection parameters for primary CRs are left free in every analysis.Finally, the bottom-right panel shows the 10 Be/ 9 Be flux ratio with the uncertainties related to the determination propagation parameters.In this panel, we have compared our results with data from the ACE [76], IMP [77,78], ISEE [79], ISOMAX [80], Ulysses [81] and Voyager [82] experiments.   .Box-plots describing the distributions for each parameter from the scenario where only secondary production of antiprotons is assumed.In the top panel we show the results obtained from the analysis that no not include correlations in the AMS-02 systematic errors and in the bottom panel those from the analysis including correlations.The median is shown as an orange line inside the boxes, while the mean is represented as a green dashed line.The boxes represent the interquartile region (first and third quartiles of the distribution), while the edges of the whiskers represent the upper and lower extremes of the distribution.As specified above "XS Cov + Scaling" refers to our Canonical analysis and "XS Scaling" to the Simplified analysis.  .Box-plots describing the distributions for each parameter from the scenario where production of antiprotons from WIMP annihilation is taken into account.In the top panel we show the results obtained from the analysis that no not include correlations in the AMS-02 systematic errors and in the bottom panel those from the analysis including correlations.The median is shown as an orange line inside the boxes, while the mean is represented as a green dashed line.The boxes represent the interquartile region (first and third quartiles of the distribution), while the edges of the whiskers represent the upper and lower extremes of the distribution.As specified above "XS Cov + Scaling" refers to our Canonical analysis and "XS Scaling" to the Simplified analysis.
Figure 1.p/p spectrum (left panels) and secondary-to-primary flux ratios of B/O, Be/O and Li/O (right panels) evaluated with the transport parameters obtained in our Canonical analyses in a scenario that does not include any DM production of antiprotons.The top row shows the spectra obtained from the analysis including no correlation in the AMS-02 systematic errors while the bottom row shows the spectra for the case when we include correlations.The variation in the flux ratios due to the uncertainty in the determination of the propagation parameters is shown (not including modulation uncertainties) in each case and the error bars on AMS-02 data are the 1σ errors reported by the collaboration.The residuals in the left panels are drawn with respect to AMS-02 data (2018).

Figure 2 .
Figure 2. p/p spectra evaluated in the scenario where the contribution from WIMP annihilation into b b final states is included, with the best-fit transport parameters obtained from the Canonical analysis assuming uncorrelated AMS-02 systematic errors (left panel) and accounting for the AMS-02 systematic error correlations (right panel).The statistical uncertainty in the determination of the propagation parameters (not including modulation uncertainties) is shown as a yellow band and the uncertainty related to the determination of the best-fit WIMP properties (mass and annihilation rate) is shown as an orange band.

Figure 3
Figure 3. Left panel: 95% confidence Upper limits on ⟨σv⟩ derived from our main analyses with a NFW profile, compared to those obtained by Refs.[22,87].Right panel: Comparison of the upper limits obtained for a contracted NFW profile (see details in the text) with other limits derived for the same profile.In particular, we compare our limits with the 2σ contour representing the best-fit DM candidate reproducing the GCE in Ref.[96], the constraints derived from dwarf spheroidal galaxies from Ref.[97] and those from Ref.[69] for a contracted-NFW profile with a halo height of H = 4 kpc (and the scale radius adjusted to 20 kpc).

Figure 4 .
Figure 4. Most used cross sections models in literature compared to available experimental data for the production of different isotopes of B, Be and Li coming from 12 C (top row) and 16 O (bottom).The band of statistical uncertainties is related to the determination of the FLUKA cross sections (see Ref.[19] for more details).

1 )Figure 5 .
Figure 5. Top row: Comparison of the different secondary-to-secondary rations of B, Be and Li evaluated from the parameters obtained assuming correlations in AMS-02 errors (top-left) and without correlations (top-right), compared to the AMS-02 data.The 2σ uncertainty from the determination of the diffusion parameters is also reported as a band for each ratio not including modulation uncertainties Bottom-left: Modulated (solid lines) and unmodulated (dashed lines) fluxes of H and He compared to AMS-02 and Voyager-1 data.The 2σ uncertainty bands from the determination of the propagation parameters and the Fisk potential are also shown.Bottom-right: Flux ratios of Li/Be, Li/B and Be/B compared to AMS-02 data, including 2σ uncertainty bands from the determination of the propagation parameters.Bottom-right panel: 10 Be/ 9 Be flux ratio with the uncertainties related to propagation parameters determination, compared with data from the ACE, IMP, ISEE, ISOMAX, Ulysses and Voyager experiments.

Figure 9
Figure 9. Box-plots describing the distributions for each parameter from the scenario where only secondary production of antiprotons is assumed.In the top panel we show the results obtained from the analysis that no not include correlations in the AMS-02 systematic errors and in the bottom panel those from the analysis including correlations.The median is shown as an orange line inside the boxes, while the mean is represented as a green dashed line.The boxes represent the interquartile region (first and third quartiles of the distribution), while the edges of the whiskers represent the upper and lower extremes of the distribution.As specified above "XS Cov + Scaling" refers to our Canonical analysis and "XS Scaling" to the Simplified analysis.

Figure 10
Figure 10.Box-plots describing the distributions for each parameter from the scenario where production of antiprotons from WIMP annihilation is taken into account.In the top panel we show the results obtained from the analysis that no not include correlations in the AMS-02 systematic errors and in the bottom panel those from the analysis including correlations.The median is shown as an orange line inside the boxes, while the mean is represented as a green dashed line.The boxes represent the interquartile region (first and third quartiles of the distribution), while the edges of the whiskers represent the upper and lower extremes of the distribution.As specified above "XS Cov + Scaling" refers to our Canonical analysis and "XS Scaling" to the Simplified analysis.
−26 cm 3 /s) local σ Analyses with only secondary production of p