Combined analyses of the antiproton production from cosmic-ray interactions and its possible dark matter origin

Recent cosmic-ray (CR) studies have claimed the possibility of an excess on the antiproton flux over the predicted models at around $10$ GeV, which can be the signature of dark matter annihilating into hadronic final states that subsequently form antiprotons. However, this excess is subject to many uncertainties related to the evaluation of the antiproton spectrum produced from spallation interactions of CRs. In this work, we implement a combined Markov-Chain Monte Carlo analysis of the secondary ratios of B, Be and Li and the antiproton-to-proton ratio ($\bar{p}/p$), while also including nuisance parameters to consider the uncertainties related to the spallation cross sections. This study allows us to constrain the Galactic halo height and the rest of propagation parameters, evaluate the impact of cross sections uncertainties in the determination of the antiproton spectrum and test the origin of the excess of antiprotons. In this way, we provide a set of propagation parameters and scale factors for renormalizing the cross sections parametrizations that allow us to reproduce all the ratios of B, Be, Li and $\bar{p}$ simultaneously. We show that the energy dependence of the $\bar{p}/p$ ratio is compatible with a pure secondary origin. We find that the energy dependence of the evaluated $\bar{p}/p$ spectrum matches the AMS-02 data at energies above $\sim3$GeV, although there is still a nearly constant $\sim10\%$ excess of $\bar{p}$ over our prediction. We discuss that this discrepancy is more likely explained from a $\sim10\%$ scaling in the cross sections of antiproton production, rather than a component of dark matter leading to antiprotons. In particular, we find that the best-fit WIMP mass ($\sim 300$ GeV) needed to explain the discrepancy lies above the constraints from most indirect searches of dark matter and the resultant fit is poorer than with a cross sections scaling.


Introduction
Cosmic rays (CRs) constitute a valuable tool for studying the galactic environment (properties of interstellar plasma [1], magnetic field configurations [2,3], astrophysical sources and their surroundings [4,5], etc). After being injected from astrophysical sources (mostly supernova remnants [6][7][8]) CRs propagate in the Milky Way and produce other nuclei (called secondary CRs) when they hit gas in the interstellar medium (ISM). While secondary CRs are used to characterize the diffusion process that CRs experience along their way toward the Earth (due to collision-less interactions with plasma waves), other secondary particles like gamma-rays and neutrinos provide useful information about the sources that accelerate CRs, since the trajectories of these particles are not affected by magnetic fields and plasma waves.
On top of this, CR antinuclei and positrons represent an useful tool for revealing the presence of the evasive dark matter (DM), since it is theoretically expected that they annihilate or decay into standard model particle-antiparticle pairs. Particularly interesting candidates are WIMPs [9], which have a mass sufficient to be able to produce nuclei-antinuclei pairs. In fact, many previous studies have tried to explain unexpected discrepancies between the predicted and measured spectra of positrons and antiprotons by including the expected production of these antiparticles from WIMPs. Nevertheless, there are yet large uncertainties related to the production of these antiparticles that could explain these discrepancies without any need of adding a dark matter component. The best example of this fact is the high-energy positron "excess", which was once considered evidence of dark matter annihilation, but is now more commonly considered to be a product of primary e + e − acceleration from pulsars [10][11][12].
Several studies have recently claimed an excess of antiprotons over the expected flux at 10 − 20 GeV, assuming they are only originated as secondary particles, and have explained it by adding annihilation or decay from a WIMP [13][14][15][16][17]. This excess seems to be compatible with a 50 − 100 GeV WIMP with a thermally-averaged cross section near the thermal cross section ( σv ∼ 3 × 10 −26 cm 3 /s) annihilating to hadronic final states (qq pairs, and mainly to bb). Interestingly, many authors claimed that this WIMP seems to be compatible with the properties of the Galactic Center Excess (GEC) of gamma-rays observed by Fermi-LAT [18,19] and has no tension with other direct or indirect DM searches [20,21]. Nevertheless, a recent study shows that this signals can not be compatible with the most recent measurements [19].
The total source term for secondary production of antiprotons from spallation reactions of CRs with the interstellar gas has the following form [22]: where the key ingredients are the CR fluxes (essentially H and He) and thep production cross sections. In particular, the determination of the cross sections ofp production is subject to important uncertainties, and results in 12 − 20% uncertainties in the evaluation of the local p spectrum [23][24][25]. Further uncertainties in the evaluation of the antiproton flux come from the determination of the propagation parameters (essentially, the diffusion coefficient used to characterize the diffusive movement of CRs and the effective Alfvén velocity). Besides, the uncertainties related to the cross sections for secondary CR production are high (> 20%), which leads to large uncertainties in the determination of the diffusion coefficient [26][27][28], since this is usually determined from the flux ratios of secondary-to-primary CRs. Also, correlated errors in the AMS-02 data seem to lead to looser constraints in the determination of the propagation parameters, reducing the significance of the antiproton excess even below 1σ level and implying that current models are consistent with the AMS-02 antiproton data [29,30]. The problem is that the correlation matrix is not yet public and the only way to proceed is to guess its terms. Additionally, the uncertainties associated to solar modulation can affect the significance of this signal [31]. Then, once antiprotons are produced, they also interact with the interstellar gas as they propagate through the Galaxy, resulting either in annihilation or loss of a fraction of their energy. Non-annihilation inelastic interactions of antiprotons with interstellar protons yield an extra source of lower energy antiprotons, which we refer to as a "tertiary" source of antiprotons, but this represents a very subdominant contribution to the total spectrum of antiprotons [32] so that the uncertainties related to this tertiary antiprotons is negligible for the evaluation of the antiproton spectrum.
In addition to the secondary antiproton production, dark matter can also produce a significant antiproton flux. the source term of antiprotons from DM annihilation (or decay) mainly depends on the DM density profile in the Galaxy ρ( x), on the velocity-averaged DM annihilation cross section (or on the decay rate) into pairs of Standard Model particles (qq pairs, W + W − , etc.), and on the DM particle mass, as shown in the following equation [33]: In this paper, we investigate the antiproton spectrum produced from CR collisions with the interstellar gas by employing an analysis that combines thep/p spectrum with the flux ratios of secondary CRs B, Be and Li to constrain the propagation parameters and accounting for systematic uncertainties in the cross sections parametrizations of production of these secondary nuclei. This combined analysis allows us to reduce the grammage needed to reproduce the spectra of B, Be and Li, therefore enhancing the production of antiprotons from CR collisions with gas. In this way, we examine thep/p ratio in the context of a combined analysis that also matches other ratios of secondary CRs. We find that our prediction greatly reproduces the energy dependence of thep spectrum reported by AMS-02 above ∼ 3 GeV albeit with a constant 10% offset compared to the data. We find that the use of the new AMS-02p/p data is very relevant here and conclude that thep/p spectrum is compatible with a pure secondary origin of antiprotons. On top of this, we discuss that the origin of this 10% discrepancy can be more plausibly explained by a scaling of thep cross sections, which offers a better fit to data with respect to a WIMP signal producing antiprotons.

Simulations setup
The propagation parameters have been usually inferred from the B/C spectrum [34], nevertheless, nowadays AMS-02 data also provides unprecedented precision for the measurements of the flux ratios of the secondary CRs Be and Li [35][36][37][38][39][40][41][42]. The difficulty in obtaining a precise evaluation of these parameters, instead, resides in the large uncertainties related to the cross sections for the production of secondary CRs. These cross sections are incorporated in CR propagation codes via parametrizations of the many reaction channels involved in the CR nuclear network [43]. Many recent works have studied the impact of these uncertainties and stressed their importance for correctly evaluating the spectra of secondary CRs [26,[44][45][46][47]. A combination of different CR observables would allow us to reduce the effect of uncertainties in the determination of the propagation parameters [48].
In this work, we implement a cylindrically symmetric (2-dimensional) model of the Galaxy. The general formulae describing CR propagation in the Galaxy are given by: which takes into account CR diffusion in space (characterized by the diffusion coefficient, D) and in momentum (also known as reacceleration, characterized by D pp ), convection or advection, energy losses, injection from sources, decays and collisions with the ISM. These equations are numerically solved with the recent DRAGON2 code 1 . More information about the code and the propagation equation can be found on refs. [49] and [32]. The diffusion equation employed in this work has a rigidity-dependence modeled as: where the reference rigidity is set to R 0 = 4 GV and the values of the rigidity break (R b ), the change in spectral index (∆δ = δ − δ h , where δ h is the spectral index for rigidities above the break), and the smoothing parameter s are: ∆δ = 0.14 ± 0.03, R b = (312 ± 31) GV and s = 0.040 ± 0.0015, as determined in Ref [50] from a detailed analysis of the proton and helium AMS-02 fluxes. The spatial diffusion coefficient and diffusion coefficient in momentum space are related by the Alfvén velocity, V A [51,52]: .
To simulate the solar modulation effect on CRs reaching the Earth, we make use of the Force-Field approximation [53], characterized by the Fisk potential, φ. Charge-sign modulation is applied here, following the approach developed in Ref [54]. In this case, the solar modulation potential is modified to have the form: choosing F (R/R 0 ) ≡ R 0 R and R 0 = 1 GV, similarly to what was done in Ref [13]. We set φ + 1,AM S−02 = 0 and explore values of φ − 1 around 0.9 GV, close to that derived in [13]. φ 0 is the usual Fisk potential, which we set to be 0.61 GV, since it allows us to reconcile the Voyager-1 [55,56] and AMS-02 data in the period 2011-2016 (period of data collection for the AMS-02 data on secondary-to-primary ratios, He, C and O). This value is also chosen to be consistent with the data from the NEWK neutron monitor experiment 2 (see [57,58]). The AMS-02 proton data was collected in the period of 2011-2013, leading to a difference in the Fisk potential ∼ 0.1 GV, so that to compute the proton flux at Earth we use φ 0 = 0.60 GV.
The new AMS-02 data on thep/p spectrum [59] is obtained from the period 2011-2018, so that we employ a value of φ 0 = 0.58 GV for the evaluation of its spectrum (∆φ is calculated using NEWK data from both periods, following [57,58]).

Analysis setup
In this work, we perform combined analyses of the secondary CRs: B, Be and Li, combined with antiprotons (p), in order to study the predictedp/p spectrum from CR interactions with the ISM gas and to study a possible signal of a WIMP-producedp in the Galaxy. The propagation parameters (H, D 0 , V A , η and δ) are inferred from the flux ratios of these CR species by means of a Markov chain Monte Carlo (MCMC) procedure, that relies on Bayesian inference. This analysis determines the probability distribution functions for the propagation parameters studied in this work to describe the AMS-02 data and their confidence intervals. This algorithm was presented in [27], but it has been upgraded to be able to include the halo height in the procedure. The analysis also incorporates nuisance parameters (scale factors) to further adjust the normalization of the spallation cross sections of B, Be and Li production (S B , S Be , S Li ), which is the main source of uncertainty in the evaluation of secondary-toprimary flux ratios. The injection parameters of primary CRs 1 H, 4 He, 12 C, 14 N, 16 O, 20 Ne, 24 Mg and 28 Si are adjusted to reproduce the AMS-02 data. These parameters are readjusted after a set of propagation parameters is obtained from the MCMC analysis, and we repeat again the analysis. We iterate in this way until the algorithm converges. This is mainly relevant for the evaluation of thep/p spectrum, for which a good adjustment of the proton and He spectra are important. For more details on the MCMC analysis, we refer the reader to Ref [27]. We take into account the production of antiprotons from nuclei heavier than He (namely, C, N, O, Ne, Mg and Si) by scaling the proton cross sections (both, cross sections of interactions with p and He as targets in the ISM) by a A 2/3 factor, where A is the mass number of the CR nuclei colliding with ISM gas. The contribution from these nuclei implies a ∼ 5 − 6% of the total antiproton flux between 10 GeV/n and 100 GeV/n, in agreement with what was found in Ref. [23]. Furthermore, we stress that we employ here the new AMS-02 data for thep/p spectrum [59].
Recently, the authors of Ref [27] showed that the secondary-to-primary flux ratios of B and Be infer roughly identical values of the normalization (D 0 ) and spectral index (δ) of the diffusion coefficient (eq. 2.2) using the DRAGON2 parametrization after rescaling their cross sections according to the AMS-02 Be/B, Li/B and Li/Be flux ratios. Nevertheless, the ratios of Li show significant differences in the determination of the propagation parameters with respect to B and Be ratios. In addition, it was demonstrated that this parametrization yield the best agreement to AMS-02 data in a combined analysis of B, Be and Li flux ratios. Therefore, the DRAGON2 parametrizations are used for the production of heavier secondary CRs [26] and the analyses are focused on the secondary-to-primary ratios of B and Be. The cross sections derived in Ref [60] for thep production cross sections are used for the evaluation of the antiproton spectrum.
Different analyses have been performed: The "Standard analysis", which includes the B/C, B/O, Be/C, Be/O flux ratios (which are able to constrain the values of D 0 /H, V A , η and δ), the Be/B, Li/B, Li/Be ratios (which allow setting the values of the scale parameters) and the 10 Be/Be and 10 Be/ 9 Be flux ratios (allowing us to constrain the halo height, H [61]). Then, we carry out another analysis including also thep/p ratio (taking into account antiprotons only produced from CR interactions with the ISM gas, eq. 1.1). Nevertheless, it must be mentioned that in these analyses we do not have into account correlations in the AMS-02 data, since the public data from the Collaboration do not include the full covariance matrix [29,45].
The combination of all these CR observables allow us to constrain the galactic halo height and the rest of propagation parameters, evaluate the impact of the spallation cross sections uncertainties in the determination of the antiproton spectrum and test the excess of antiprotons at ∼ 10 GeV. Other works have also included the uncertainties related to the spallation cross sections, although they only used the B/C ratio and a refit to the very uncertain cross sections measurements of B production (with 1σ uncertainties > 20%, in average [26]). Here, in turn, we make use of the B and Be CR species and the scale factors calculated by combining the high energy part of the flux ratios among B, Be and Li with the cross sections data. Furthermore, these works fix the value of the halo height and do not include it in their analyses (mainly due to the huge uncertainties related to its determination).

Computation of a WIMP contribution
Finally, we also perform an analysis including the possible annihilation of a dark matter particle (a WIMP) into hadronic states that producep (eq. 1.2), in order to infer the WIMP mass and thermally averaged annihilation cross sections ( σv ) needed to explain the discrepancies found assuming a pure secondary origin of antiprotons. This analysis is focused in the χ → bb annihilation channel. Two different DM density profiles are analysed to assess the influence of this component in the determination of the WIMP mass and thermally averaged annihilation cross sections. Among the most used DM profile densities, derived from N-body simulations and cosmological studies [62,63], the most common one is the Navarro-Frenk-White (NFW) density profile [64], which we use as reference. Other common density profiles are, the Einasto profile [65,66], Isothermal profile [67], Moore profile [68], or the Burkert profile [69,70]. In order to compare with the NFW profile, we use the Burkert profile, since they present the largest differences, as can be seen from figure 1 of ref [71]. These profiles are generally given by: with a characteristic halo radius r h ∼ 20 kpc , and a characteristic halo density ρ h , used to normalize the profiles from observations [72]. In particular, we normalized these profiles to a value ρ = 0.43 GeV at the Solar System position, r = 8.3 kpc. We use DarkSusy [73] to compute the antiproton production from WIMPs at Earth.

Standard analysis
This analysis focuses on the combination of the ratios of secondary CRs mentioned above. The diffusion parameters obtained from this analysis are given in table 1, along with their statistical uncertainties. As we see, the best fit halo height, H, is around 7 kpc, compatible with recent analyses [26,46,74] and the ratio D 0 (10 28 cm 2 /s)/H (kpc) is close to 1. We must mention here that the determination of the H value is very uncertain, since it is determined from the ratios of 10 Be, for which most of experimental data are located around a few hundred MeV/n. This energy region is particularly affected by the uncertainties related to the solar modulation effect and the diffusion coefficient at such energies -the low energy part of the diffusion coefficient is the most uncertain, since this energy region is particularly complex, given the effect of damping and dissipation of plasma waves that seem to become important at sub-GeV energies [75,76] and the fact that there are very few data points in the secondary-toprimary spectra below 1 GeV/n. On top of this, at energies around a few hundreds of MeV/n the production of secondary CRs is dominated by nuclear resonances. The value found for the effective Alfvèn velocity, V A , is ∼ 22 km/s, in agreement with the expected values derived in Ref [77]. In addition, the value of the spectral index of the diffusion coefficient, δ, is ∼ 0.43, well compatible with other recent analyses [47,78]. The determination of this parameter is also affected by an extra systematic uncertainty around ±0.01, since we fixed the value of the spectra index after the break (and break position), as discussed in [27]. Interestingly, the scale factors derived in the analyses, as nuisance parameters, are very small: a 5 − 6% scaling is found for the cross sections of B production and a ∼ 1% scaling for Be and Li. We note that, while the related 1σ statistical uncertainty is of ±0.1, the systematic uncertainty related to these scale factors has been estimated to be around 10%, with a ∼ 5% uncertainty associated to multi-step reactions and a ∼ 3% due to inelastic cross sections [27]. The gas distribution employed in the simulations can also add an additional 4% uncertainty in the determination of these scale factors.
As we can see from panels a, b and c of figure 1, the ratios of the secondary CR nuclei can be fitted simultaneously within the AMS-02 uncertainty (< 5%). Nevertheless, the evaluation of thep/p with these propagation parameters yield a 20 − 30% discrepancy with respect to the AMS-02 data, as shown in 1, d. Experimental data are taken from the ASI Cosmic Ray Data Base [79]

Combined analysis of thep/p spectrum
In this section, we report the results obtained from the combined analysis of the flux ratios included in section 3.1, adding also thep/p ratio above 4 GeV in the analysis, without including additional nuisance parameters for the cross sections of production ofp (i.e. not modifying the original cross sections). We have chosen 4 GeV to reduce the impact of the solar modulation uncertainties in the evaluation of thep/p ratio and performed the analysis also with thep/p spectrum from energies from 3 and 5 GeV, obtaining no significant differences in the evaluated propagation parameters above 1σ. As we see from figure 2, the ratios of secondary CR nuclei are still compatible with AMS-02 data within 2σ statistical uncertainty. The scale factors play a major role in this analysis, allowing us to readjust the grammage and diffusion parameters (mostly determined by the ratios B/C, B/O, Be/C and Be/O) to improve our predictedp/p spectrum while keeping a good agreement to the other secondary ratios. In fact, these parameters change by up to 10% with respect to the values predicted from the Standard analysis (see table 1). We remind the reader that the uncertainties associated to the cross sections measurements are > 20% for the main channels of production of Be and  B. We also observe that the value of δ inferred from this analysis is around 0.49, closer to a Kraichnan spectrum for the turbulence. In addition, we notice that the Alfvén velocity tends to be compatible with no reacceleration (V A = 0 km/s) and that the best-fit halo height is slightly different from that found in the Standard analysis (∆H ∼ 1 kpc), which means that the antiproton spectrum is also affected by this parameter (and not only to the ratio D 0 /H). Figure 3 shows thep/p spectrum evaluated with these propagation parameters, and using the values φ 0 = 0.58 GV and φ n = 0.9 GV for the Fisk potential (eq. 2.4). From this evaluation we obtain an important conclusion: the resultant residuals with respect to AMS-02 data are roughly constant and at the level of ∼ 10%, below the uncertainty associated to the cross sections of production ofp [23]. Therefore, we have also added to this figure a dashed line that represents the same predicted spectrum but scaled by 10% in order to show that it is able to fit the AMS-02p/p experimental data within the reported 1σ errors above ∼ 3 GeV. We have varied the value of φ n from 0.8 to 1 GV, finding that the residuals at 10 GeV remain the same within ±2%. We argue here that, while the fact that the residuals are so low is mainly due to the inclusion of the scaling factors for B, Be and Li production cross sections in the analysis, the new data makes the residuals to be flatter than with the previousp AMS-02 data, as shown in more detail in appendix A. We argue that the uncertainty related to the contribution of CR species heavier than H and He could also affect this discrepancy by ∼ 3%.  On top of this, the same spectrum but scaled by 10% is also shown as a dashed line, allowing us to see that a simple 10% scaling would make us reproduce experimental data within 1σ errors.

Testing possible antiproton production from a WIMP
The fact that these residuals are so flat above 3 GeV could be related to the effect of cross sections uncertainties, instead of an extra source of antiprotons (such as antiprotons produced from dark matter), whose contribution would have a bump-like structure. In fact, it is very reasonable to think of a 10% rescaling of the cross sections of antiproton production, either from the component of prompt production of antiprotons (direct creation, CR + gas →p) or from the component of antineutron and antihyperons decay (CR + gas → {n,Λ,Σ} →p), or a scaling in both (see, e.g. eq. 5 of ref [60]). Nevertheless, we also evaluate the spectrum of antiprotons produced from a WIMP, in order to test if this contribution could explain the excess found in thep/p spectrum. To do so, we test the mass and thermal annihilation cross sections ( σv ) that yield the best fit with AMS-02 data. This analysis consists of evaluating the spectrum of antiprotons produced from a WIMP annihilating in the bb channel and sum this contribution to the diffusep produced from CR collisions (with propagation parameters found in the combined analysis ofp/p spectrum). The computation of the antiproton spectrum generated from the WIMP is performed with the DarkSusy code, for the NFW and Burkert dark matter profiles. Then, we use our MCMC algorithm to find the values of WIMP mass and σv that provide the best fit on thep/p ratio.
In the left panel of figure 4, we display the result of the analysis of the NFW profile (the WIMP spectrum looks the same for the Burkert profile) in comparison to AMS-02 data. Here, we can see that the component of antiprotons produced from DM annihilation does not seem consistent with the energy dependence of the thep/p spectrum, although it can nearly reproduce this spectrum (although antiproton production from two, or more, annihilation channels could fit well with experimental data). We must highlight that the χ 2 value of this fit is ∼ 57 (evaluated over the 44 data points of thep/p spectrum above 4 GeV), a 1): %85.(57 Figure 4. Results of the best fit WIMP candidate able to reproduce the discrepancy with respect AMS-02 data. Left panel: thep/p spectrum produced from the best WIMP candidate annihilating in the bb channel is shown as a dotted line. The uncertainty associated to the determination of mass and σv is shown as an orange band. The diffusep/p spectrum from CR collisions with interstellar gas, found in the combinedp/p analysis, is shown as a dashed line, with its associated statistical uncertainty as a yellow band. Finally, the sum of both contributions is shown as a solid line and the total statistical uncertainty is shown as a green band. Right panel: Triangle plot showing the mass and thermally averaged annihilation cross sections of the WIMP candidate that provides the best fit to AMS-02 data. The 68% and 95% uncertainty regions for the Burkert and NFW dark matter density profiles are shown. The best-fit mass value is found to be 304.35 +8.37 higher value than from scaled spectrum, which is ∼ 31. In the right panel of figure 4, the best fit values of WIMP mass and σv are shown as a triangle plot with the 68% and 95% uncertainty regions. The best fit values are of ∼ 300 GeV for both profiles and a value of the thermal annihilation cross sections of ∼ 1.5 × 10 −25 cm 3 /s for the NFW profile and ∼ 4.7 × 10 −25 cm 3 /s for the Burkert profile. These values are much larger than the values usually quoted in indirect searches of DM from antiprotons (as seen in appendix A this is mainly related to the use of the newp AMS-02 data) and well above with calculations of the thermal relic cross sections [81]. In addition, these values are in tension or above the Fermi-LAT constraints of Dwarf galaxies [82] and Milky Way satellites [83], and are not compatible with the possible dark matter signal from the Galactic Center Excess (GCE) [18]. However, these values are neither discarded from Planck CMB constraints nor from AMS-02 constraints on positrons, as shown in ref [82]. It is also important to mention that, even though other dark matter annihilation channels would change the mass and σv predicted, the other hadronic channels (tt, cc and ss) produce a very similar flux of antiprotons at 10 GeV and the W + W − channel leads to a slightly larger σv to reproduce thep/p excess. The direct inputs and outputs from the DRAGON2 runs with the propagation parameters obtained from both analyses are available at this repository 5 , as well as the best-fit spectra of proton, helium and antiproton. Moreover, we include in the repository the tables of antiproton flux produced from the annihilation of WIMPs for diverse masses and σv values in the bb channel for the two DM density distributions studied in this work and show a comparison between the flux of antiprotons produced from annihilation from all the hadronic channels and the W + W − channel at the best-fit mass and σv found for the NFW distribution.

Discussion and conclusions
In this work, we have studied the diffusep/p spectrum produced in CR interactions with the interstellar gas. We first carry out a combined analysis, without including thep/p ratio, in which we have taken into account different secondary CRs and their ratios, as well as nuisance parameters that allow us to renormalize the cross sections of production of B, Be and Li. These observables allow us to determine the propagation parameters (H, D 0 , V A , η and δ) to provide a simultaneous fit of the ratios of the different CR species. We show that thep/p spectrum predicted from this analysis differ from the AMS-02 data by a factor of 20 − 30%.
Then, we include thep/p spectrum in our combined analysis, making use of the new AMS-02 data for this spectrum. This analysis yield a much closer prediction of thep/p spectrum with respect to recent AMS-02 data, producing discrepancies that are constant above 3 GeV and are of ∼ 10%. Interestingly, the new AMS-02 data lead to flatter residuals than with the previous dataset. We discuss that this discrepancy seems to be plausibly explained by just a 10% scaling of the cross sections ofp production, which makes our model match completely the AMS-02p/p experimental data within 1σ uncertainties. To test alternative explanations of this excess, we evaluate if an extra component of antiprotons from annihilation of a WIMP could reproduce this discrepancy. We find that, although this component can be close to reproduce thep/p spectrum, the values of mass and thermally averaged annihilation cross sections (the best fit values found are: WIMP mass ∼ 300 GeV and σv ∼ 1.5 − 4.7 × 10 −25 cm 3 /s) are in tension or ruled out by several other indirect dark matter searches. Also the χ 2 value obtained from the fit with a 10% scaling is almost a factor of two lower than the fit considering DM annihilating in antiprotons.
Therefore, we conclude that taking into account all the sources of uncertainties in the evaluation of the secondary antiprotons produced from CR interactions with gas allows as to explain thep/p spectrum without any need of an extra source of production ofp. Specifically, it seems that the energy dependence of thep/p spectrum is well reproduced assuming a pure secondary origin of antiprotons and that the excess found is plausibly explained by a rescaling of the cross sections ofp production. On the contrary, it seems difficult that the component of WIMP annihilation alone could be the responsible of this excess, although this source of antiprotons could be still present, below the high level of uncertainties associated to the evaluation of thep spectrum.
In this way, we provide a set of propagation parameters and scale factors for renormalizing the cross sections parametrizations that allow us to reproduce all the ratios of B, Be, Li andp simultaneously.
Finally, we highlight the most important points revealed of this work: • An analysis that takes into account uncertainties in the production cross sections of the main secondary CRs B, Be and Li leads to residuals in thep/p spectrum always lower than 12% at 10 GeV. The total uncertainties are, at least, a factor of two greater.
• The use of the newp/p AMS-02 data has an important impact in the conclusions reached, since the residuals found are much flatter than when the old dataset is employed.
• The shape of the reported excess does not resemble a bump, as would be expected from dark matter annihilation. In addition, the best-fit WIMP candidate obtained with this new dataset seems to be incompatible with most of the indirect searches of dark matter currently done. The fact that we obtain a much larger WIMP mass than the one predicted in previous works is mainly due to the use of the new AMS-02p/p data.

Comparison with other works
Many recent works have been dedicated to the study of the antiprotons produced from CR collisions and the possible signature of antiprotons produced from DM. The main difference with respect to previous analyses is the use of the scale factors for renormalizing the cross sections of production of the main secondary CRs and the use of the new AMS-02 data. Previous works have found an excess of antiprotons with maximum residuals around 10 GeV, with a bump-like structure similar to what we find in the analysis with the old AMS-02 data ( fig. 5). Examples of these studies are refs. [15,17,24]. However, they did not take into account the full systematic uncertainties in the evaluation of thep production due to the uncertainties in spallation cross sections of B (since they use only the B/C ratio). In addition, they predict a WIMP mass close to 100 GeV to explain the antiproton excess, which would have changed significantly (towards larger WIMP mass) using the new AMS-02 dataset.
Some other works claim that this excess is not significant at all under the systematic uncertainties. This is the case of, for example, refs. [29,60,84], where the authors claim that thep spectrum is compatible with a pure secondary origin when taking into account all the systematic sources of uncertainty. The two last references also include possible correlations in the AMS-02 data from a guessed correlation matrix. This is something we do not include in this study, since we employ only the errors reported by the AMS-02 collaboration, and would have the effect of broadening the uncertainties associated to the determination of the propagation parameters.
A Results with the previous AMS-02p/p dataset In this appendix, we show the results of the combined analysis including the old AMS-02 data for thep/p ratio (data collected in the period 2011-2015). We apply the same analysis as in section 3.2, with the difference that the Fisk potential employed here is set to φ 0 = 0.61 GV also for thep/p spectrum (corresponding to the period of data collection of the previousp/p dataset).
Relevantly, the propagation parameters found in the analysis with the oldp/p dataset are roughly identical to those obtained using the new AMS-02 dataset (see table 1), as expected.
In figure 5, we show the predicted spectrum evaluated with the Winkler cross sections and the propagation parameters labeled as "p/p old data" in table 1. In the leftpanel, the experimental data are those from the previous AMS-02 dataset (data collection from 2011-2015), while in panel b they are the newly published data (data collection from 2011-2018). A clear change in the shape of the residuals is obtained, even though the evaluated spectrum is the same in both cases. The new dataset leads to a flattening of the shape of the residuals and the WIMP signal which is able to explain this discrepancy is shifted towards larger masses.
In figure 6, we repeat the exercise shown in figure 4, but for the WIMP signal that achieves the best fit the oldp/p dataset. Panel a shows the evaluated spectrum with the contribution of a WIMP producing antiprotons (WIMP annihilating into the bb channel) which provides the best fit to this experimental data. In the right panel, we show the 1 and 2σ significance regions of the probability distribution function of the inferred mass and σv , for the Burkert and NFW dark matter profiles in the Galaxy. While the inferred σv value (equal to ∼ 1.2 cm 3 /s for the NFW profile and ∼ 3.7 cm 3 /s for the Burkert profile) is slightly lower than the one obtained in the analysis with the new dataset, the inferred mass is significantly lower (from ∼ 300 GeV to ∼ 130 GeV). In fact, the WIMP signal necessary to reproduce thep/p spectrum with the new experimental data is much less credible than with the previous dataset.

B Summary of the MCMC results: propagation parameters
In this appendix, we report the diffusion parameters obtained for the standard and combined p/p analyses in table 1, which contains the median values (maximum posterior probability value) ± the 1σ uncertainty related to their determination and the actual range of values contained in the 95% probability of the distribution. We also show the probability distribution functions (PDFs) of each propagation parameter obtained from the different analyses in figure 7.

Propagation parameters and scale factors
H (kpc) D 0 (10 28 cm 2 /s) v A (km/s) η δ S B S Be S Li Figure 7. Probability distributions of the considered propagation parameters obtained for the main analyses performed in this work. The contour plots highlight the 68% and 95% credible regions.