AMS-02 antiprotons' consistency with a secondary astrophysical origin

TL;DR

This paper presents a detailed prediction of antiproton fluxes from secondary astrophysical processes, demonstrating their consistency with AMS-02 data and highlighting challenges in achieving high-precision dark matter searches.

Contribution

It introduces an improved method for calibrating cosmic-ray propagation models and provides a robust flux prediction accounting for uncertainties, supporting a secondary origin of antiprotons.

Findings

AMS-02 antiproton data are consistent with secondary production.

Uncertainties and correlations significantly impact flux predictions.

Achieving percent-level precision in predictions remains challenging.

Abstract

The AMS-02 experiment has ushered cosmic-ray physics into precision era. In a companion paper, we designed an improved method to calibrate propagation models on B/C data. Here we provide a robust prediction of the $\bar{p}$ flux, accounting for several sources of uncertainties and their correlations. Combined with a correlation matrix for the $\bar{p}$ data, we show that the latter are consistent with a secondary origin. This paper presents key elements relevant to the dark matter search in this channel, notably by pointing out the inherent difficulties in achieving predictions at the percent-level precision.

Figures

Tables

Table 1. Table 1: Respective p 𝑝 p -values for different sources of errors. We take dof = 57 absent 57 =57 , i.e. the number of p ¯ ¯ 𝑝 \bar{p} data. Total errors on data are defined to be σ tot = σ stat 2 + σ syst 2 subscript 𝜎 tot superscript subscript 𝜎 stat 2 superscript subscript 𝜎 syst 2 \sigma_{\rm tot}=\sqrt{\sigma_{\rm stat}^{2}+\sigma_{\rm syst}^{2}} .

Error considered	${\chi }^2/\mathrm{dof}$	p-value $\left({\chi }^2\right)$	p-value (KS)
${\sigma }_{\mathrm{stat}}$	23	0	0
${\sigma }_{\mathrm{tot}}$	1.69	$8.3\times {10}^{-4}$	0
${𝒞}^{\mathrm{data}}$	0.85	0.79	0.97
${\sigma }_{\mathrm{stat}}$ and ${𝒞}^{\mathrm{model}}$	1.32	0.05	0.99
${\sigma }_{\mathrm{tot}}$ and ${𝒞}^{\mathrm{model}}$	0.37	1.0	0.01
${𝒞}^{\mathrm{data}}$ and ${𝒞}^{\mathrm{model}}$	0.77	0.90	0.86

Table 2. Table 2: Best-fit parameter values and χ 2 / dof superscript 𝜒 2 dof \chi^{2}/{\rm dof} on H, He, C, and O AMS-02 data for the three benchmark transport models BIG , SLIM , and QUAINT   [ 58 ] . The first line recalls the modulation parameter value fixed by the analysis of AMS-02 B/C data [ 58 ] , which is also taken for H, He, C, and O AMS-02 data. See text for discussion.

	BIG	SLIM	QUAINT
${\varphi }_{\mathrm{FF}}^{\star }$ [GV]	0.731	0.734	0.725
${\chi }_{\min }^2/{\mathrm{dof}}^{†}$	0.90	0.71	1.93
Source parameters
${\alpha }_H$	2.37	2.41	2.34
${\alpha }_{\mathrm{He}}$	2.31	2.34	2.30
${\alpha }_{Z>2}$	2.33	2.36	2.33
${\log }_{10}{\left(q_H\right)}^{‡}$	-3.98	-3.92	-4.05
${\log }_{10}{\left(q_{\mathrm{He}}\right)}^{‡}$	-4.44	-4.39	-4.46
${\log }_{10}{\left(q_C\right)}^{‡}$	-5.86	-5.82	-5.86
${\log }_{10}{\left(q_O\right)}^{‡}$	-5.75	-5.71	-5.75
High-rigidity break parameters
${\Delta }_h$	0.17	0.26	0.13
$R_h$ [GV]	358.	303.	436.
$s_h$	0.04	0.15	0.002

Equations

$${\chi }^2=(\mathrm{data}-\mathrm{model}{)}^{\mathrm{T}}({\mathcal{C}}^{\mathrm{model}}+{\mathcal{C}}^{\mathrm{data}}{)}^{-1}(\mathrm{data}-\mathrm{model}).$$

$${\sigma }_{\mathrm{inv}}^{\mathrm{tot}}={\sigma }_{\mathrm{inv}}(2+{\Delta }_{\mathrm{IS}}+2{\Delta }_{\Lambda }),$$

$$\bar{\Lambda }/\bar{p}=c_1+\frac{c_2}{1+(c_3/s{)}^{c_4}},$$

$${\Delta }_{\mathrm{IS}}=c_0(x+c_2{)}^{c_3}\exp (-x/c_1)$$

$${\mathcal{C}}_{ij}^{\alpha }={\sigma }_i^{\alpha }{\sigma }_j^{\alpha }\exp \left[-\frac{1}{2}{\left(\frac{{\log }_{10}(R_i/R_j)}{{\ell }_{\alpha }}\right)}^2\right],$$

$${\mathrm{c}}_{ij}^{\alpha }=\frac{{\mathcal{C}}_{ij}^{\alpha }}{\sqrt{{\mathcal{C}}_{ii}^{\alpha }}\sqrt{{\mathcal{C}}_{jj}^{\alpha }}},$$

$${\mathcal{C}}^{\mathrm{model}}=\sum _{\mathrm{a}}{\mathcal{C}}^{\mathrm{a}}={\mathcal{C}}^{\mathrm{xs}}+{\mathcal{C}}^{\mathrm{transport}}+{\mathcal{C}}^{\mathrm{parents}},$$

$${\mathcal{C}}_{ij}^{\mathrm{a}}=\frac{1}{N}\sum _{n=1}^N\left({\Phi }_{i,n}^{\mathrm{a}}-{\mu }_i^{\mathrm{a}}\right)\left({\Phi }_{j,n}^{\mathrm{a}}-{\mu }_j^{\mathrm{a}}\right).$$

$$\mathcal{C}={\mathcal{C}}^{\mathrm{data}}+{\mathcal{C}}^{\mathrm{model}},$$

$${\chi }^2=\sum _{i,j}{x}_i{\left({\mathcal{C}}^{-1}\right)}_{ij}{x}_j\equiv {\mathbf{x}}^{\mathrm{T}}{\mathcal{C}}^{-1}\mathbf{x},$$

$${\tilde{x}}_i=U_{ij}{x}_j\text{while}\tilde{\mathcal{C}}=U\mathcal{C}{U}^{\mathrm{T}}.$$

$${\chi }^2=\sum _i{\tilde{z}}_i^2.$$

$${\tilde{R}}_i=\sum _j{U}_{ij}^2{R}_j,$$

Full text

The AMS-02 experiment has ushered cosmic-ray physics into precision era. In a companion paper, we designed an improved method to calibrate propagation models on B/C data. Here we provide a robust prediction of the $\bar{p}$ flux, accounting for several sources of uncertainties and their correlations. Combined with a correlation matrix for the $\bar{p}$ data, we show that the latter are consistent with a secondary origin. This paper presents key elements relevant to the dark matter search in this channel, notably by pointing out the inherent difficulties in achieving predictions at the percent-level precision.

††preprint: LAPTH-017/19, ULB-TH/19-04, LUPM:19-056

Introduction —

The spectrum of Galactic cosmic-ray (CR) antiprotons ($\bar{p}$’s) originating from interactions of CR protons on the interstellar (IS) gas was first calculated fifty years ago [1, 2]. CR $\bar{p}$’s were discovered a few years later in balloon-borne experiments [3, 4], but in excess compared to this expected secondary background; this was soon interpreted as the presence of an extra contribution, possibly from primordial black holes [5] or dark matter (DM) annihilation [6, 7]. Over the past thirty years, a lot of progress has been made both theoretically and experimentally.

From the 90’ to the 00’s, balloon flights (MASS, IMAX, CAPRICE, HEAT-pbar, and BESS [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]) and AMS-01 onboard the shuttle [18] collected hundreds of $\bar{p}$’s, up to $\sim 10$ GeV. Improvements in the predictions, accounting for tertiary production [19], more accurate cross sections [20], and reacceleration [21] paved the road to obtain compatibility of the data with modern diffusion models [22, 23], the transport parameters of which were fitted on B/C [24].

In the following decade, the next generation of instruments—BESS-TeV/Polar [25, 26, 27] and the PAMELA satellite [28, 29, 30]—measured several thousands of $\bar{p}$’s up to a few hundreds of GeV. Updates of the model parameters and corresponding predictions, using better measured $p$ and He CR spectra, yielded a secondary flux in agreement with the data [e.g. 31].

The state-of-the-art AMS-02 experiment has been operating on the ISS since 2011. It has already recorded several tens of thousands of $\bar{p}$’s [32], up to TeV energies. For the first time, the measurement is dominated by systematic uncertainties (at a few percent level), implying new challenges for their interpretation. Here we address this issue, going beyond previous analyses of preliminary AMS-02 $\bar{p}$’s data [33, 34, 35].

We underline that CR $\bar{p}$’s are one of the most sensitive astroparticle probes of annihilating/decaying DM in the GeV-TeV range [36, 37, 38, 39, 40], and any constraint on DM candidates depends on how well the astrophysical background is controlled. This is especially important as claims for $\bar{p}$ excesses attributed to DM are being debated [41, 42, 43, 44, 45].

Methodology —

The flux of CR $\bar{p}$’s at Earth depends on (i) the cross sections entering their production, scattering, and annihilation, (ii) the CR propagation model, (iii) the IS spectrum of the (most abundant) CR nuclei, and (iv) modulation of fluxes in the Solar cavity.

Uncertainties on relevant nuclear cross sections are among the dominant ones for the $\bar{p}$ flux calculation [22, 33]. Recently, new data have been taken [46, 47, 48], leading to improvements in the cross section parametrisations [49, 50, 51, 52, 53]. A crucial ingredient needed for the calculation is the Lorentz invariant cross section of prompt $\bar{p}$’s produced in $pp$ interactions. We make use of the parametrisation proposed in Ref. [52] and improved in Ref. [53] (Param II). For nucleon-nucleon interactions, we use the scaling relation B proposed in Ref. [53] which provides the best agreement with LHCb data. We use the covariance matrices of errors on the parameters to propagate the uncertainties to the $\bar{p}$ flux calculation. Antineutrons ($\bar{n}$’s) and antihyperons are produced in hadronic interactions and decay into $\bar{p}$’s; their contribution has to be included in the total cross section. The yield of $\bar{p}$’s produced via $\bar{n}$’s is larger than in the prompt $\bar{p}$ channel, an energy-dependent effect studied and parameterised in Ref. [52], in agreement with the very scarce experimental data. We introduce an exponentially decreasing energy-dependent parametrisation in order to reproduce the 1$\sigma$ CL interval of Fig. 8 in Ref. [52]. Regarding the anti-hyperon contribution, we make use of the parametrisation determined in Ref. [52] and the covariance matrix of the corresponding parameters to propagate the uncertainties on the $\bar{p}$ flux calculation. More details on the nonprompt channels and on the way we propagate the related uncertainties are given in the App. A.1 and App. A.2, respectively. Inelastic annihilating and non-annihilating interactions of $\bar{p}$’s with the IS medium (ISM) are treated following the procedure described in [38].

We use the 1D diffusion model implemented in usine v3.5 [54] which assumes a thick diffusion halo size $L$, and a thin disc containing the sources and the gas. Computationally more expensive scenarios with extended geometries [55, 56, 50] and time dependence [57] have been considered in the literature. Yet, should simple 1D models prove sufficient to account for the data in a secondary production scenario, more complex scenarios would a fortiori work as well. We use the most generic transport model defined in [58] (called BIG) where the transport parameters are fitted on B/C following the methodology described in [59], i.e. a model of the covariance matrix of B/C AMS-02 errors has also been incorporated in the fit. Using the minuit [60] package to fit the B/C, we compute the best-fit values of the free parameters (transport and nuisance parameters) and their covariance matrix. We fully propagate the transport uncertainty to the $\bar{p}$ flux from the covariance matrix. We have compared the distribution of the best fit parameters from the covariance matrix with the distribution obtained by resampling the B/C data to check that the Gaussian treatment of the uncertainties is valid in this context.

Important inputs for the $\bar{p}$ calculation are the IS fluxes for all progenitors (mostly H, He, C and O, but also all nuclei up to Fe, see App. B.3 for their detailed contributions; we also account for heavy elements in the ISM, see App. B.4). There are two options to propagate the associated uncertainties: a standard approach in semianalytical models [22, 31, 33, 52] is to demodulate top-of-atmosphere (TOA) data, in order to extract the parameters describing the IS fluxes and their uncertainties, which allows us to yield the $\bar{p}$ flux with a fast procedure and to easily propagate the progenitors uncertainties. The preferred approach in fully numerical propagation codes is to directly use the IS fluxes calculated at the propagation stage, ensuring that the calculated fluxes fit the data [23, 61, 35, 55]. In this global fitting, the transport and source uncertainties are determined simultaneously. Although using a semianalytical model, we follow this second approach but do not rely on a global fitting of the transport and source parameters. Instead, we follow and extend the two-step procedure detailed in [58] (see Sec. III therein): We start from the best-fit parameters obtained from the B/C analysis (BIG) and then perform a simultaneous fit of H [62], He, C, and O [63] AMS-02 data to determine the source parameters, i.e. four normalisations (1H, 4He, 12C, and 16O) and three slopes ($\alpha_{\rm H}$, $\alpha_{\rm He}$, and a universal source slope $\alpha_{\rm Z>2}$ for all other species). Source isotopic fractions are fixed to Solar System values [64]; the abundances of non-fitted elements up to $Z=30$ are fixed such as to match the HEAO-3 data [65] at 10.6 GeV/n. We fit at the same time the high-rigidity diffusion break parameters which are better constrained by elemental fluxes than by the B/C [58]. We include in the fit the covariance matrices of errors on H, He, C, and O data, see App. B.1. The outputs of this fitting procedure are discussed in App. B.2. Note that we fix the Fisk potential 111The force-field approximation [66, 67, 68, 69] has a single free parameter, and it was found to describe well modulated fluxes [70]. $\phi_{\rm FF}$ for H, He, C, and O data to the value yielded by the B/C data fit [58]. While published AMS-02 data are from the same time period for B/C, He, C, and O (05/2011 to 05/2016), they originate from a shorter period for H (05/2011 to 11/2013), so that the associated modulation level should be slightly different. However, as the fit is restricted to data above the $\bar{p}$ production threshold ($E_{k/n}>6\,m_{p}$ [22]), the impact is feeble, and in any case negligible compared to other uncertainties.

State-of-the-art modulation models are 2D or 3D charge-dependent diffusion and drift in the Solar cavity [71]. They go beyond the force-field approximation [72] and at variance with the latter, they predict different modulations for positively- and negatively-charged particles. Charge-sign dependent effects on secondary antiprotons have been estimated by Ref. [73] to be smaller than $\sim 10\%$ below 5 GeV and negligible above. However, these models have a large number of parameters and are still under study. Moreover, the impact of the drift effect is strongly dependent on the data-collection period. The AMS-02 $\bar{p}$ data [32] were analysed from a time period of 4 years (May 2011 to May 2015) longer than that of H and shorter than that of He, C, and O data. As illustrated in [74], drift models can be used to derive effective and different Fisk potentials for nuclei and $\bar{p}$’s, with $\phi_{\rm FF}^{\bar{p}{}}$ which is not clearly smaller or larger than $\phi_{\rm FF}^{\rm nuc}$ on a sufficiently long period. Assuming the same modulation for $\bar{p}$ and H, He, C, and O actually already gives a satisfactory description at low rigidity. For this reason, we did not include further uncertainty on $\phi_{\rm FF}^{\bar{p}{}}$. If anything, this would slightly enlarge the uncertainty of the prediction at low rigidity, improving further the consistency between our calculated $\bar{p}$ flux and the data.

Uncertainties in the parameters discussed above should be propagated to the $\bar{p}$ flux calculation. To proceed, for each source of uncertainty (production cross section, transport in the Galaxy, fit to parents fluxes) we draw $\approx 10000$ realizations of the parameters from their covariance matrix and we compute the corresponding $\bar{p}$ fluxes. We then calculate the 1$\sigma$ confidence-level (CL) envelope for each source of uncertainty as well as the covariance matrix of the model $\mathcal{C}^{\rm model}$ (see App. C). This enables us to soundly assess the compatibility of the model with the data.

In the context of the B/C analysis, we stressed the importance of using a realistic covariance matrix of the data errors to avoid misleading conclusions [59]. We anticipate that the same is true for $\bar{p}$’s. However, since this matrix is not directly provided by the AMS-02 collaboration, we build it from the published systematic errors and associated description of their physics origin, in the same spirit as in [59].

The various contributions to the AMS-02 systematics are broadly described in [32]. For instance, the text “This [selection] uncertainty amounts to 4% at 1GV, 0.5% at 10GV, and rises to 6% at 450 GV” is interpreted as a piecewise power-law behaviour, and is shown as an orange thin line in Fig. 1. Thus we build the seven sources of systematics quoted by the AMS-02 collaboration (coloured thin lines), where the quadratic sum of all contributions leads to the black thin line. In order for the sum to match the total systematic errors provided in [32] (black thick line), we rescale for each rigidity point our separate contributions by the ratio of the thick to the thin black lines. This leads to our model for the AMS-02 $\bar{p}$ systematics (coloured thick lines). The covariance matrix associated with these systematics is then built based on a choice of their correlation length, $\ell$. More details on this procedure are given in App. C and, for the B/C analysis, in [59]. For $\bar{p}$’s, we take as educated guesses for the correlations lengths (in unit of energy decade) $\ell_{\rm Acc.}=0.1$ (acceptance), $\ell_{\rm Cut}=1.0$ (rigidity cut-off), $\ell_{\rm Scale}=4.0$ (rigidity scale), $\ell_{\rm Templ.}=0.5$ (template fitting), $\ell_{\rm XS}=1.0$ (cross sections), $\ell_{\rm Unf.}=1.0$ (unfolding), and $\ell_{\rm Sel.}=0.5$ (selection).

Results —

The top panel of Fig. 2 shows our baseline $\bar{p}$ flux prediction (not a fit) obtained from the best-fit values for the $\bar{p}$ production cross sections, the transport (BIG) and the associated parents fluxes, compared with AMS-02 data with errors taken as the quadratic sum of systematic and statistical errors (black crosses). The ‘standard’ residuals with respect to the baseline model are displayed in the middle panel. Note that the points do not include the model uncertainties, nor correlations in the data uncertainties. We also show on the same plot the 68% total confidence band for the model (grey band). This band could release the tension with the data, even before accounting for the information on the correlations in rigidity bins. The respective contributions of parents, cross sections and transport are also plotted. At tens of GV, the errors from transport and cross sections are almost constant and close to 10%. At larger rigidities, the errors from transport and parents increase because of the increasing experimental uncertainty in the B/C ratio and parent fluxes, respectively. At low rigidity, the error from transport grows for the same reasons and encompasses the uncertainty in the prediction of the low-rigidity behaviour (see App. E). However, we remind the reader that a visual comparison can be deceiving: The presence of non-diagonal values in the covariance matrices is responsible for a better agreement between the model and data than perceived in the residuals (see Table 1).

To test the actual compatibility of our prediction with the $\bar{p}$ data, we present two statistical tests which boil down to probabilistic statements in terms of $p$-value (see App. D). First, we propose a $\chi^{2}$ test, with the help of a covariance matrix of errors on both data and model:

[TABLE]

The covariance matrices of the data ${\cal C}^{\rm data}$ and the model ${\cal C}^{\rm model}$ are given by the sum of the different contributions previously detailed (see SM, Sec. II-A and III). We find $\chi^{2}\approx 44$, and identifying the number of degrees of freedom (dof) with the number of $\bar{p}$ data points (57), we infer a corresponding $p$-value of 0.9 which is reported in the last line of Table. 1.

Such a test does not directly assess a possible overestimate of the errors, and also relies on the notion of number of dof (which may be a shaky concept in some circumstances, see e.g. the discussion in [76]). Thus, we also perform a Kolmogorov-Smirnov (KS) test, which obviates the above limitations. We compute the distribution of the ‘eigen residuals’ ($\tilde{z}$-score) corresponding to the residuals of the eigen vectors (data-model) of the total covariance matrix. In the bottom panel of Fig. 2 we show these ‘eigen residuals’ as a function of ‘rigidity’ (actually, the one rotated in the eigen basis, see App. D), and in the inset, we compare the corresponding histogram with a Gaussian. The KS test leads to a $p$-value of 0.27 which is also very good and is clearly consistent with the hypothesis that $\bar{p}$’s are of secondary origin. For completeness, we also report in Table. 1 the $p$-values when considering different combinations of errors: (i) If there was no uncertainty in our baseline model, the covariance matrix of data errors alone (${\cal C}^{\rm data}$) would already give enough freedom to allow for a very good agreement between the data and the secondary flux prediction; (ii) Considering only the statistical uncertainties in the data and the uncertainties in the model ($\sigma_{\rm stat}$ and ${\cal C}^{\rm model}$), this prediction is marginally consistent with the data at the 2$\sigma$ level, with the KS test leading to an even better $p$-value. Also note the relevance of the KS test (as opposed to the $\chi^{2}$ test) to spot error overestimates, in the case of $\sigma_{\rm tot}$ and ${\cal C}^{\rm model}$; (iii) In the most realistic case considering both ${\cal C}^{\rm data}$ and ${\cal C}^{\rm model}$, $p$-values are very good for both the $\chi^{2}$ and KS test. Thus, not only is a secondary origin for the locally measured $\bar{p}$’s statistically consistent with the data, but, as shown by these considerations, it is also robust with respect to error mismodelling in either model or data errors.

Conclusions —

Percent-level details in the model predictions now matter, as do more subtle aspects of the data error treatment. In this paper we have presented a major upgrade of the $\bar{p}$ flux prediction and analysis by: (i) using the latest constraints on transport parameters from AMS-02 B/C data, (ii) propagating all uncertainties (with their correlations) on the predicted $\bar{p}$ flux, and (iii) accounting for correlated errors in $\bar{p}$ data. The multi-component nature of the systematic error, with different $R$-dependencies and correlation lengths, has a crucial impact on the analysis, and was not captured in more simplified treatments as in Ref. [44]. With these novelties, we unambiguously show that the AMS-02 data are consistent with a pure secondary astrophysical origin. We stress that this conclusion is not based on a fit to the AMS-02 $\bar{p}$ data, but on a prediction of the $\bar{p}$ flux computed from external data. Our results should hold for any steady-stade propagation model of similar complexity, as they all amount to the same “effective grammage” crossed to produce boron nuclei (on which the analysis is calibrated), with roughly the same grammage entering the secondary $\bar{p}$’s. We have checked that this conclusion is robust with respect to a variation by a factor of a few of the correlation lengths of the AMS-02 systematic uncertainties. Also, recent analyses of Fermi-LAT data are suggestive of a spatial dependent diffusion coefficient, notably different in the inner Galaxy [77]. Moving to more complex scenarios containing the 1D framework considered here as limiting case would broaden theory space, but would not alter our conclusions on the viability of secondary production to explain antiproton data. On the technical aspects, more computationally expensive methods could allow one to go beyond the quadratic assumption (i.e. assuming multi-Gaussian error distributions) embedded in the covariance matrix of errors. For more advanced applications, sampling techniques like Markov chain Monte Carlo could be used (e.g., [78]). However, a significant improvement in our perspectives for DM searches in the $\bar{p}$ flux can only be achieved by simultaneously reducing the systematics in the data and the errors of the modelling. On the data side, a covariance matrix of errors directly provided by the AMS-02 collaboration would definitively be an important improvement to fully benefit from the precision achieved by AMS-02. On the modelling side, the next step would be to combine more secondary-to-primary ratios (Li/C, Be/C, and B/C) to further decrease the propagation uncertainties. Of course, better data and modelling on $\bar{p}$ and $\bar{n}$ production cross sections is also required, and the sub-leading error due to primary source parameters could be reduced by combining AMS-02 data with higher energy data from CREAM, TRACER and CALET [79]. In the current state of our analysis, we can anticipate that, from the frequentist point of view, a clear statistical preference for an additional feature in the data is unlikely. However, this conclusion must rely on a quantitative analysis that we postpone for a forthcoming paper.

Acknowledgements — MB is grateful to Michael Korsmeier and Martin Winkler for very useful discussions. We are grateful to all the members of the Cosmic Rays Alpine Collaboration. This work has been supported by the “Investissements d’avenir, Labex ENIGMASS”, by Univ. de Savoie AAP “DISE”, and by the French ANR, Project DMAstro-LHC, ANR-12-BS05-0006. The work of Y.G. is supported by the IISN, the FNRS-FRS and a ULB ARC. We also acknowledge a partial support from the Agence Nationale pour la Recherche (ANR) Project No. ANR-18-CE31-0006, the Origines, Constituants, et EVolution de l’Univers (OCEVU) Labex (No. ANR-11-LABX-0060), the CNRS IN2P3-Theory/INSU-PNHE-PNCG project “Galactic Dark Matter”, and the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreements No. 690575 and No. 674896.

Appendix A Antiproton production cross sections

Most of the $\bar{p}$’s are produced via decays of antihyperons and $\bar{n}$’s. These contributions are not included in the prompt $\bar{p}$ cross section measured by the collider experiments used in Ref. [53]. Ref [80] proposed energy-dependent parametrisations based on the most up-to-date collider data. The total $\bar{p}$ production cross section is

[TABLE]

where $\sigma_{\rm inv}$ is the Lorentz invariant cross section of promptly produced $\bar{p}$’s, $\Delta_{\rm IS}$ and $\Delta_{\Lambda}$ the isospin asymmetry and anti-hyperon corrections, respectively.

A.1 Antihyperons

A sizeable fraction of $\bar{p}$’s is produced via the decay of antihyperons $\bar{\Lambda}$ and $\bar{\Sigma}$. The ratio of hyperon-induced to promptly produced $\bar{p}$’s is $\Delta_{\Lambda}=(0.81\pm 0.04)(\bar{\Lambda}/\bar{p})$, where $\bar{\Lambda}/\bar{p}$ is the ratio of $\bar{p}$’s produced via the decay of $\bar{\Lambda}$ to the total yield. Here, it is assumed that antihyperons decay equally into $\bar{p}$’s and $\bar{n}$’s. The energy dependence of $\bar{\Lambda}/\bar{p}$ is given by

[TABLE]

where $\sqrt{s}$ is the centre-of-mass energy, $c_{1}=0.31$, $c_{2}=0.30$, $c_{3}=(146\,\rm GeV)^{2}$, and $c_{4}=0.9$. We determine the uncertainty by randomly drawing 10000 realisations of the parameters from the covariance matrix. This method enables us to reproduce the median and the 1$\sigma$ CL obtained in Ref.[80]. They are shown in the top panel of Fig. 3.

A.2 Antineutrons

Antiprotons are preferentially produced via the decay of $\bar{n}$’s. The energy dependence of the isospin asymmetry proposed in Ref. [80] is motivated by the fact that a pion has to be produced when $\bar{p}$’s are promptly produced. We found that taking the parametrisation and the covariance matrix of parameters from Ref.[80] does not allow to reproduce the median and the 1$\sigma$ CL interval the authors show in Fig. 8. Instead, we introduce the parametrisation

[TABLE]

where $x=\log(\sqrt{s})$, $c_{0}=0.33$, $c_{1}=6.42$, $c_{2}=0.5$ and $c_{3}=1.6$. We randomly draw 10000 realisations of $c_{0}$ using a Gaussian distribution with a standard deviation of $\sigma_{c_{0}}=0.36$ and with the lower bound $c_{0}>0.04$. This way, we recover the median and the 1$\sigma$ CL interval of Fig. 8 in Ref. [80]. The results are shown in the middle panel of Fig. 3. Since we adopt the nucleon-nucleon scaling relation determined in Ref. [53], which depends (though slightly) on a different parametrisation of $\Delta_{\rm IS}$, we rescale the parameter $D_{2}$ entering Eq. (17) of Ref. [53] in order to obtain the same value of $f^{A_{1}A_{2}}$ at the energy of NA49.

The ratio between $\sigma_{\rm inv}^{\rm tot}$ as expressed in Eq. (2) and $\sigma_{\rm inv}$ is shown in the bottom panel of Fig. 3. In the case of vanishing nonprompt contribution and isospin violation, one would expect a constant value equal to 2. Hence, not only the $\bar{p}$ flux receives an upwards correction of $\sim 20\%$ to $50\%$, but also acquires peculiar energy-dependent features, further affected by relatively large uncertainties.

Appendix B Antiproton parents

B.1 AMS-02 data errors

In H, He, C, and O fluxes published by the AMS-02 collaboration [62, 63], uncertainties are broken-down into several contributions: statistical (Stat.), unfolding (Unf.), rigidity scale (Scale), acceptance (Acc.), and even trigger (Trig.) for protons data. As discussed in [59], the information provided in the Supplemental Material of AMS-02 publications allows us to build a first estimate of the covariance matrix of uncertainties $\mathcal{C}^{\alpha}$. For each uncertainty $\alpha$ (Stat., Unf., Scale, Acc., and Trig.), the $i\!j$-th element of the matrix is taken to be

[TABLE]

where $\sigma^{\alpha}_{i}$ is the systematic uncertainty of type $\alpha$ taken at rigidity bin $R_{i}$, and $\ell_{\alpha}$ is the correlation length indicating the range over which rigidity bins are correlated.

As the data analysis of AMS-02 nuclei follows many similar steps, the correlation lengths derived in [59] for B/C data [75] should not be too different from the ones for the proton [62] and He, C, and O data [63]. Consistently with [59], we take:

• •

$\ell_{\rm Stat.}=0$ (statistical uncertainties are uncorrelated);

• •

$\ell_{\rm Scale}=\infty$ (rigidity scale impacts all rigidities);

• •

$\ell_{\rm Unf.}=0.5$ (unfolding acts on intermediate scales) though the exact value is not critical [59];

• •

$\ell_{\rm Trig.}=1$ (intermediate error, close to a normalisation error, for the H data only);

• •

$\ell_{\rm Acc.}$ is ill-defined as it comes from a mixture of uncertainties with different correlation lengths (data/Monte Carlo corrections, cross section normalisation). Following [59], we split it into

• –

$\ell_{\rm Acc.~{}norm.}\sim 1.0$ (normalisation error);

• –

$\ell_{\rm Acc.\,LE}\sim 0.3$ (low-energy rise);

• –

$\ell_{\rm Acc.~{}res.}=0.1$ (residual acceptance).

Note that $\ell_{\rm Acc.~{}norm}$ is not well constrained, while it was found to be critical for the B/C analysis [59]. Its value relies on indirect consistency arguments (analysis of transport parameters), and we refer the reader to [59] for the discussion and caveats about this choice.

The above numbers and various systematics (shown in Fig. 4) are plugged in Eq. (5) to generate the covariance matrices for H, He, C, and O. These matrices are used for the determination of source parameters.

B.2 Primary CR source parameters from the AMS-02 H, He, C, and O data

As discussed in App. C, for each propagation configuration (BIG, SLIM, and QUAINT, see [58] and also Sec. E), we fit the source parameters and the high-rigidity break of the diffusion coefficient on H, He, C, and O AMS-02 data, relying on their covariance matrix of errors introduced in Sec. B.1. As we are only interested in the region where these species produce $\bar{p}$’s, we limit the fit to data above the threshold production $E_{k}/n=6m_{p}$. As the H, He, C, and O nuclei are the dominant $\bar{p}$ progenitors (see Sec. B.3), we check here that good fits are achieved in order to ensure the reliability of our $\bar{p}$ prediction.

Table 2 gathers the outputs of the fits, i.e. the $\chi^{2}/{\rm dof}$ values, best-fit values for the sources parameters (4 normalisations and 3 slopes) and high-rigidity diffusion break parameters (break slope, position, and sharpness). First, the values obtained here for the diffusion break parameters and source slope $\alpha_{Z>2}$ are consistent with those of the analysis carried out in [58] based on the C and O nuclei only (and only using total errors instead of the covariance matrix of errors in the fit). The additional information brought from this analysis is a comparison of the H, He, and $Z>2$ source spectral indices222He, C, and O data were analysed at a later stage by the AMS-02 collaboration, by which time improvements on the alignment and rigidity scale were made. This could slightly change the best-fit break values and source parameter for this quantity.. While this is not the goal of this study, it is worth mentioning that all models prefer a slightly harder spectral index for He than for heavier species, $\sim 0.02$ harder (uncertainties, not shown in the Table, are $\sim 0.007$). We also confirm a different source slope, $\sim 0.06$, between H and He [81, 62], except in the QUAINT propagation model configuration. However, the latter give a bad fit to the primary elemental data333The bad fit of H data in QUAINT is understood as this model has high levels of reacceleration and convection, which affects more strongly $A/Z=1$ species (i.e. H) than $A/Z\approx 2$ species (He, C, O). with $\chi^{2}/{\rm dof}=1.93$.

Figure 5 shows further details with a graphical view of the goodness-of-fit for the H (blue), He (red), C (green), and O (orange) nuclei data in and outside the fit regions (grey shaded area on the left and centre panels). The left panels show $z$-scores, i.e. the residuals measured in unit of $\sigma_{\rm tot}$, (model-data)/$\sigma_{\rm tot}$, for the BIG (top), the SLIM (middle), and the QUAINT (bottom) propagation model configurations: a good fit implies that 68% of the points should be distributed within $1\sigma_{\rm tot}$. However, this representation does not account for correlations in the systematic errors. As discussed in Sec. D, it is more interesting to show the same quantity in a rotated basis in which the covariance matrix of errors is diagonal (called $\tilde{z}$-scores)444The price to pay is that now the matrix of rigidity positions is no longer diagonal, and the value we pick to represent in the middle panels of Fig. 5 is the rotated $z$-score (dubbed $\tilde{z}$-score) defined in Sec. D.. The $\tilde{z}$-scores shown in the second column of Fig. 5 no longer display the spectral features seen in the left panels. This illustrates the importance of properly visualising scores with the covariance accounted for, as we may otherwise be visually biased by these structures. The right panels of Fig. 5 are the histograms of the $\tilde{z}$-score points displayed in the middle panels. They indicate that whenever a good $\chi^{2}/{\rm dof}$ value is obtained, the normalised residuals follow a $1\sigma$-width Gaussian distribution (solid black line). These plots provide further details on the separated contributions from H, He, C, and O data: in all cases, the H data are always less well fitted (broader distribution), especially in the QUAINT model configuration, which explains the bad $\chi^{2}/{\rm dof}$.

In the light of this study, we correctly reproduce the H, He, C, and O AMS-02 data (except for the QUAINT model configuration) which gives good confidence in the self-consistency of the $\bar{p}$ flux calculation.

B.3 Ranking parent contributions

From the BIG model configuration, we show in Fig. 6 the fractional IS $\bar{p}$ contribution per parent element, properly accounting for each CR isotopes (and summed over all ISM components, here H and He). We wish to highlight a few features in this plot: the most significant channels are those involving the most abundant CR fluxes, i.e. H ($\sim 80\%$), He ($\sim 20\%$), then O, C, Fe, Si, and Mg ($\sim 1-3\%$), followed by Ne and N ($\sim 0.5\%$). All the remaining ones give sub-percent contributions, but we stress that all the contributions below 0.2% (grey shaded area) add up to a total of $\sim 1\%$ (solid black line). Also, if we only account the contributions of elements heavier than Si, they add up to a total of $\sim 2\%$ (solid red line).

The contributions of the two most important elements (H and He) are broken down into their isotopic content (blue and green solid lines for H and He isotopes respectively), highlighting the contributions from primary species (1H and 4He) w.r.t. sub-dominant secondary ones (2H and 3He, lighter-coloured lines).

Focusing on the rigidity dependence, the feature at a few GV is related to secondary species: the peaks/dips seen at these rigidities directly reflect the peaks in secondary-to-primary ratios (e.g., B/C, 3He/4He, etc.). Also the spectra of secondary species are one $\delta$ (slope of the diffusion coefficient) steeper than those of primary species, so that their contributions fall off with rigidity. Among the primary species, a difference is observed between H and all heavier elements: H has a steeper source spectrum than the other species (see Table 2), so that its relative contribution w.r.t. these elements decreases with rigidity.

For the sake of illustration, we show in Fig. 7 that these fractional contributions are not very sensitive to the transport configuration considered (BIG, SLIM, and QUAINT setups, see Sec. E). This is not surprising as all three are tuned to match the B/C, H, He, C, and O data. The main difference is observed at high rigidity in the QUAINT configuration, because in this model, contrary to the BIG and SLIM ones, H and other species have a similar source spectral indices (see Table 2 in Sec. B.2).

B.4 Impact of heavy species in the ISM

To be able to perform tens of thousands fits of the H, He, C, and O data together with the associated $\bar{p}$ flux predictions, almost all our calculations were carried out restricting ourselves to CR parents up to 30Si. Moreover, a standard practice in the field is to assume that the IS gas is only made of H and He. However, in the context of high-precision $\bar{p}$ data, it is important to reassess the impact of heavier species in both CRs and in the ISM, to control the accuracy of the $\bar{p}$ flux calculation. For instance, the fraction of CNO in the ISM is estimated to be $\sim 0.5-1\%$ [21, 82, 53], whereas that of heavy species in CRs (e.g. Fe) is estimated to be $\lesssim 1\%$ [53] (see Fig. 6).

Before coming to our prescription to account for these contributions, a few remarks are in order, since ‘heavy’ CRs and ISM enter the calculation at several steps: to extract the transport parameters from the B/C analysis, in the calculation and fit of H, He, C, and O fluxes, and then in the $\bar{p}$ flux calculation. Our reference calculation here is based on the transport configuration BIG (see Sec. 11 of [58] for more details), and the best-fit transport parameters proposed in [58] only consider CR parents up to 30Si with H and He only in the ISM. Hence, for a self-consistent calculation, we must refit the B/C (using the same procedure as in [58]), to extract the transport parameters for the same CR and ISM content as for the $\bar{p}$ propagation. The last missing piece is the detail of the implementation of the ISM composition and cross sections for our calculations. The IS gas composition and species abundances are taken up to Fe from [83], the production cross sections for $\bar{p}$’s were already discussed in the main text and in Sec. A, and for the nuclear inelastic and production cross sections, we rely on a $A_{\rm target}^{0.31}$ scaling relation proposed in [84] for 3He production.

We show in Fig. 8 the ratio of a given TOA calculation to our reference $\bar{p}$ flux. To better understand each effect, we show separately the impact from the ‘propagation stage’ and from the fully consistent calculation: in the former case, we keep the same CRs and ISM as in the reference for the $\bar{p}$ flux calculation, but use the best-fit transport parameters with ‘heavy’ CRs and/or ISM species, whereas for the latter case, ‘heavy’ CRs and ISM species are self-consistently used at all stages of the calculation.

• •

‘Propagation stage only’ calculation (dashed and dotted thin lines): the thin dotted blue line shows that accounting for CRs up to Fe555Heavier species are negligible: Ni contributes less than 0.05% to $\bar{p}$ (see Fig. 6) and the abundances of $Z\geq 30$ elements drop further by several orders of magnitude. impacts on the transport coefficients—and thus $\bar{p}$—at the percent level; adding an extra abundance of heavier ISM species (black dashed line) further increases it to $\lesssim 4\%$ effect, but with a very different spectral shape. Intuitively, these changes go in the right direction, as opening more channels for boron production means less grammage to get the same B/C data, and thus less $\bar{p}$’s. The fine details are difficult to track down as they involve a mixture of increased production (CRs and ISM species) balanced by an increased destruction on the heavier ISM species (more efficient both at lower rigidities and for heavier species).

• •

‘Fully consistent’ calculation (thick solid lines): accounting for heavy CRs and heavy ISM species in the transport parameters and in the $\bar{p}$ flux calculation leads to an overall impact of $\sim 1-2\%$, i.e. smaller than in the ‘propagation stage only’ case discussed above. This is not so surprising since as we open more channels for boron production, we also open more channels for $\bar{p}$ production, the two effects mostly cancelling out666The rigidity-dependence of the extra contributions for nuclei and $\bar{p}$’s are not exactly the same, leading to different spectral features. For instance, considering heavier CRs (solid blue line) leads to a $\bar{p}$ flux increasing with rigidity up to $2\%$ at 1 TV, likely to be related to the increased fraction of Fe w.r.t. lighter CRs: lighter elements fragment less than heavy ones at low rigidity.. Finally, accounting for heavier species in the ISM leads to a plateau between 3 and 200 GV ($\sim 2\%$ increase), whose shape is similar to that obtained from the ‘propagation stage only’ calculation (dashed thin black line).

There are three important messages to take away from this analysis.

Reaching a few % accuracy in the $\bar{p}$ flux calculation, even within the framework of a single propagation model, requires to carefully account for heavy CRs up to Fe, and also to account for heavy species in the ISM, notwithstanding the very large uncertainties on the specific cross sections to consider to do so. 2. 2.

These contributions do not amount to a simple normalisation effect, but have subtle spectral features. If heavy species are not consistently taken into account, they may lead to spectral features at the few percent level. 3. 3.

Because we want to have a calculation as accurate as possible, all our $\bar{p}$ fluxes calculated from ‘CRs up to 30Si’ and ‘ISM = H, He’ are rescaled by the solid black curve of Fig. 8; the analysis presented in the main paper thus accounts for the contribution of all relevant heavy species.

Appendix C Correlations of data and model uncertainties

To compute the likelihood of the AMS-02 $\bar{p}$ flux, we first assume the data to follow a Gaussian distribution, with covariance matrix ${\cal C}^{\rm data}$, around the theoretical prediction. The various sources of experimental errors, each contributing to ${\cal C}^{\rm data}$, are discussed in the main text. Here, for each of these, we plot in Fig. 9 the corresponding correlation matrix defined as

[TABLE]

where $i$ and $j$ stand for rigidity and ${\cal C}$ is the covariance matrix. We notice that ${\rm c}_{ij}^{\alpha}$ lies in the range from 0 to 1. The first eight panels show the rigidity correlation of the individual uncertainty contributions. The bottom right panel shows the total correlation matrix for the data, which combines the individual correlation matrices weighted by their associated uncertainties (Fig. 1 of the main text). This latter panel shows that at low- and high-rigidities, the bins are weakly correlated (dominated by statistical uncertainties). At intermediate rigidities, the dominant uncertainties is from the Acc., hence a correlation on short rigidity ranges (bottom centre panel), but contributions from other systematics add significant larger-range correlations (white and light blue shade in bottom right panel).

The predicted fluxes also suffer from uncertainties related to model parameters. The likelihood is built assuming these fluxes also follow a Gaussian distribution, with covariance matrix

[TABLE]

each of the above-mentioned sources contributing independently to the covariance matrix ${\cal C}^{\rm a}$. Strictly speaking, predictions would follow a Gaussian distribution, should they vary linearly with respect to Gaussian-distributed model parameters. As errors are small, this idealized case is not far from reality. But cuts on input parameters, such as the Alfvénic speed or the convective wind, and the presence of insuperable although small non-linearities, have motivated us to go a step further. We still stick to the Gaussian assumption, but we directly explore the $\bar{p}$ flux behaviour at each rigidity bin $i$ while deriving how different rigidities are correlated.

To do so, for each source ${\rm a}$ in Eq. (7), we first draw randomly $N\sim 10^{4}$ realizations of the input parameters, assuming they are Gaussian-distributed with a covariance matrix which we borrow from previous analyses. In the case of CR transport, for instance, we use the best-fit values and covariance matrix of the best-fit low- and intermediate-rigidity parameters as published in [58]. For each set $n$ of input parameters, we calculate the $\bar{p}$ flux $\Phi_{i,n}^{\rm a}$ at AMS-02 $\bar{p}$ rigidity $i$. From the resulting distribution, we compute the $1\sigma$ CL flux envelope. We also infer the average prediction $\mu_{i}^{\rm a}$ and covariance matrix

[TABLE]

The associated correlation matrix ${\rm c}_{ij}^{\rm a}$ is defined as in Eq. (6). It is plotted in Fig. 10 for the three sources of model uncertainties, i.e. $\bar{p}$ production cross section, CR transport and parents. Notice how strongly the different rigidity bins are correlated. This plainly justifies the use of covariance matrices, which turn out to be of paramount importance. The three sources of uncertainties show three different behaviours for their correlation matrices: (i) ‘Parent’ shows a strong correlation between all bins; (ii) ‘Transport’ shows a weaker correlation between low- ($\lesssim 4$ GV) and high-rigidities, probably related to the low-energy break in the transport model which decouples them; (iii) ‘XS’ is in-between, reflecting the various cross section contributions; (iv) ‘Total’, as for the data case, combines the individual correlation matrices weighted by their associated uncertainties (middle panel in Fig. 1 of the main text). It is dominated by ‘Transport’ at low rigidities, it is a mixture of ‘Transport’ and ‘XS’ at intermediate rigidities, and a mixture of ‘Transport’ and ’Parents’ at high rigidities. We stress that the total correlation matrix for the model is very different to that of the data (Fig. 9).

Appendix D The likelihood tests

Our analysis aims at establishing whether or not the $\bar{p}$ flux measured by AMS-02 is compatible with a pure secondary origin. We do not perform any fit of the data. We just gauge how likely they are, once our model for CR transport and secondary $\bar{p}$ production is given. Although our approach is frequentist, priors on the theoretical parameters are incorporated in the definition of the likelihood. The sources of model uncertainties are to be found in $\bar{p}$ production cross section, CR transport and parent fluxes.

The likelihood tests are performed by constructing a $\chi^{2}$ that accounts for both experimental and theoretical uncertainties, with their correlations. The global covariance matrix is defined as the sum

[TABLE]

while the $\chi^{2}$ is taken to be

[TABLE]

where $i$ and $j$ are rigidity indices. The i-th element $x_{i}$ of vector ${\mathbf{x}}$ is the residuals at rigidity $i$, i.e. the difference between the AMS-02 measurement and the baseline theoretical prediction (inferred with the best-fit CR parameters of model BIG of [58] and using full covariance errors for parent fluxes). Once the $\chi^{2}$ is determined, the $p-$value gauges how good is the agreement between data and baseline model.

A more direct yet qualititive test is provided by the visual inspection of the ${z}$-score which, at rigidity $i$, is defined as $z_{i}={x_{i}}/{\sqrt{{\cal C}_{ii}}}$. However, in the case of non-vanishing correlations between different rigidity bins, the ${z}$-score could be dangerously misleading. Residuals $x_{i}$ need actually to be rotated in a new basis where the covariance matrix becomes diagonal. Such a rotation is performed with the orthogonal matrix $U$ and yields

[TABLE]

By construction, the covariance matrix $\tilde{\cal C}$ is diagonal, with elements $\tilde{\cal C}_{ii}=\tilde{\sigma}_{i}^{2}$. The rotated ${z}$-score, dubbed $\tilde{z}$-score, is defined as $\tilde{z}_{i}={\tilde{x}_{i}}/{\tilde{\sigma}}_{i}$ and is directly related to the $\chi^{2}$ through

[TABLE]

The index $i$ labels now the pseudo (or rotated) rigidity which we define as

[TABLE]

with $R_{j}$ the physical rigidity. In the basis of the eigenvectors of the covariance matrix $\tilde{\cal C}$, the rigidity is no longer diagonal. In our case, the rotation is small, with $U$ close to unity. The pseudo rigidity $\tilde{R}_{i}$ is not very different from the physical value $R_{i}$.

Finally, by virtue of Eq. (12), the $\tilde{z}$-score can be interpreted as the correct distance between the data and the baseline model. With our assumptions (see above), each value $\tilde{z}_{i}$ is expected to be distributed normally, with mean 0 and variance 1. This motivates us to perform a Kolmogorov-Smirnov (KS) test on the AMS-02 data, to supplement the $\chi^{2}$ analysis.

Appendix E Uncertainties from the SLIM/BIG/QUAINT transport configurations

In [58], we introduced three propagation configurations: SLIM (pure diffusion with low- and high-rigidity breaks), QUAINT (‘standard’ convection/reacceleration transport with high-rigidity break), both limiting cases of the more general configuration BIG (low- and high-rigidity breaks, convection, reacceleration). From a statistical point of view, BIG is slightly preferred by the B/C analysis, but the other two also provide satisfactory $\chi^{2}/{\rm dof}$ [58]. Repeating our procedure for the $\bar{p}$ flux calculation (take the best-fit transport parameters, fit high-rigidity break and source parameters on H, He, C, and O data, and calculate $\bar{p}$), we investigate here how sensitive the results are to the propagation configuration.

E.1 $\bar{p}$ from the BIG/SLIM/QUAINT setups

Figure 11 shows the relative difference w.r.t. BIG, for calculations based on the SLIM (dashed) and QUAINT (dash-dotted) best-fit parameters. The top panel shows the relative difference between the best-fit H (blue), He (red), C (green), and O (orange) TOA fluxes, while the bottom panel shows the relative difference for the corresponding IS $\bar{p}$ flux prediction777We recall that the different configurations have different best-fit values for $\phi_{\rm FF}$, see Table 2. The H, He, C, and O fits are performed assuming the same modulation level as obtained in the B/C analysis. See discussion in the main paper.. The grey shaded area in the top panel highlights rigidities below the kinetic threshold of $\bar{p}$ production888The production threshold is $E_{k}=6m_{p}$, which translates into slightly different rigidities for H ($A/Z=1$) and other elements ($A/Z\sim 2$), but for simplicity we use a single region for all elements in the plot.. Both panels share the same rigidity axis, but to interpret these graphs, we recall that $\bar{p}$’s at rigidity $R$ GV are created from parents whose rigidities are ten to thousand times larger. Indeed, the low-rigidity variation of the $\bar{p}$ flux between configurations can be attributed to the configurations themselves, because the few percent variation on the fit parents above the threshold (top panel) cannot provide the $\gtrsim 10\%$ $\bar{p}$ variation below a few GV. At high rigidity, the $5-10\%$ increase in $\bar{p}$ in the QUAINT setup seems to be related to the rise of the H ‘bad’ fit (discussed in Sec. B.2). However, most of the difference observed for the $\bar{p}$ flux in the SLIM configuration may be related to the different high-rigidity break $\Delta_{h}$ (0.26 vs 0.17 in BIG, see Table 2). The latter break in the diffusion coefficient appears only once in the formula of primary fluxes, whereas it appears twice for secondary fluxes: for given fluxes of parents (top panel), the high-rigidity slope of the $\bar{p}$ predictions (bottom panel) will differ by $\Delta_{h}$, asymptotically.

We have two last comments on Fig. 11. First, all the fits to the H, He, C, and O data make use of the covariance matrices of errors introduced in Sec. B.1, except for the blue solid line shown in the bottom panel. The latter was calculated for the same configuration as for the reference case (BIG), but using for the H, He, C, and O data uncertainties the quadratic sum of statistical and systematic errors. This shows that a better description of the data uncertainties for the $\bar{p}$ parents amounts to a $\sim 1-2\%$ impact on the calculated $\bar{p}$ flux. Second, in the range where AMS-02 measured the TOA $\bar{p}$ flux ($1-300$ GV), the calculations differ by $\lesssim 20\%$, meaning that the inter-model variation is already larger or similar to the uncertainties on the $\bar{p}$ flux data. We also stress that the difference is not a pure normalisation one but also a spectral one; this is important to keep in mind as potential percent excesses in $\bar{p}$ data (w.r.t. to the astrophysical prediction) may vanish using another propagation model.

E.2 Intra- vs inter-model comparison and $z$-scores

The top panel of Fig. 12 shows the relative difference of various quantities, (config-ref)/ref, with ref taken to be our baseline configuration with the transport model BIG and where the covariance matrix of errors of AMS-02 H, He C and O have been used. The symbols show the residuals of AMS-02 data, which have spectral features. However, other propagation configurations (SLIM and QUAINT) have different spectral residuals, highlighting once more that bumps and dips in the residuals should not be taken too seriously. The solid blue lines bracket the $68\%$ CL on the transport uncertainties for BIG. It is worth noting that the intra-model transport uncertainty is at the same level as the inter-model uncertainty. This is not surprising since SLIM and QUAINT are actually subcases of the BIG model (see [58]).

The second and third panels of Fig. 12 show the $z$-score in the standard and ‘rotated’ basis respectively (see Sec. D). Indeed, the middle panel shows some spectral features (bumps) that the eyes would be eager to interpret as potentially significant DM induced excesses. However, the middle panel does not account properly for correlations in the data between different rigidity bins, whereas the bottom panel does. There is then no significant correlation left, and the inset shows that the $\tilde{z}$-scores distribution is close to the expected normal distribution (solid black line). This shows that the three transport configurations proposed in [58] are all consistent with the $\bar{p}$ data. For a quantitative analysis (using the BIG transport configuration), see the main paper.