Dark matter or correlated errors? Systematics of the AMS-02 antiproton excess

Several studies have pointed out an excess in the AMS-02 antiproton spectrum at rigidities of 10-20 GV. Its spectral properties were found to be consistent with a dark-matter particle of mass 50-100 GeV which annihilates hadronically at roughly the thermal rate. In this work, we reinvestigate the antiproton excess including all relevant sources of systematic errors. Most importantly, we perform a fully realistic derivation of the correlations in the AMS-02 systematic error which could potentially"fake"a dark-matter signal. The dominant systematics in the relevant rigidity range originate from uncertainties in the cross sections for absorption of cosmic rays within the detector material. For the first time ever, we calculate their correlations within the full Glauber-Gribov theory of inelastic scattering. The AMS-02 correlations enter our spectral search for dark matter in the form of covariance matrices which we make publicly available for the cosmic-ray community. We find that the global significance of the antiproton excess is reduced to below 1 $\sigma$ once all systematics, including the derived AMS-02 error correlations, are taken into account. No significant preference for a dark-matter signal in the AMS-02 antiproton data is found in the mass range 10-10000 GeV.


Introduction
Since their discovery about 40 years ago [1,2] cosmic-ray antiprotons have been used as a sensitive probe of exotic cosmic-ray sources in our Galaxy, such as dark-matter annihilation. As a matter of fact, already their first measurement exhibited an excess over the expected astrophysical background and has stimulated speculations about a dark-matter contribution [3,4]. While significant theoretical and experimental progress has been made since, today the situation appears similar, although on an entirely different level of precision. By now several ten thousand antiproton events have been reported by the AMS-02 experiment on-board the International Space Station rendering statistical uncertainties subdominant over a large range of rigidities [5]. In this range systematic errors are at the level of a few percent and constitute the limiting factor in data analyses, which nonetheless allows us to search for a dark-matter contribution potentially as low as ∼ 10% of the total antiproton flux.
Recently, several groups have reported an excess over the expected antiproton background in the rigidity range 10-20 GV in the AMS-02 data which is compatible with a dark-matter annihilation signal [6][7][8][9][10][11][12][13]. While the significance of the excess is highly controversial (ranging from 1 − 5 σ in the aforementioned studies), a common picture of the preferred dark-matter properties has emerged. It hints at a particle of mass m χ = 50−100 GeV which annihilates into hadronic final states with roughly a thermal cross section, σv ∼ 10 −26 cm 2 s −1 . Intriguingly, dark matter with similar properties has been considered in the context of the galactic center gamma ray excess [14].
The key ingredient to test the dark-matter interpretation of the antiproton excess is a careful modeling of those systematic effects which could, alternatively, have caused the observed spectral feature. To this end, strong efforts have been made to improve the prediction and to quantify the uncertainties of antiproton production by cosmic-ray scattering [15][16][17][18][19]. The updated cross-section modeling entering the antiproton background was, indeed, found to somewhat reduce the antiproton excess [9][10][11]. However, one major piece is missing in all previous studies: the correlations of systematic errors in the AMS-02 data which have so far not been reported by the collaboration. Not surprisingly, these are of paramount importance given that correlated systematics can induce unwanted features in the data. In [11], it has been shown that correlations, in this case modeled by simple covariance functions, potentially have dramatic effects on the significance of the antiproton excess (see also [20]).
In this study we carefully derive estimates for the correlations in the AMS-02 data and investigate their implications for the tentative dark-matter signal. For this purpose, we collect all publicly available information to split the systematic error into its components which we then address individually. In the rigidity range 10-20 GV the dominant systematics in the antiproton flux andp/p ratio arise from uncertainties in the cross sections for (anti)proton absorption in the AMS-02 detector for which the measured fluxes are corrected.
As a first step we undertake a detailed re-evaluation of the uncertainties of the involved nucleon-carbon absorption cross sections (the AMS detector is dominantly composed of carbon). To this end, for the first time in the literature, we perform a global fit of the absorption cross sections within the full Glauber-Gribov theory of inelastic scattering [21][22][23][24][25]. It links the nuclear absorption cross section to the nucleon-nucleon scattering cross sections and nuclear density functions which are as well subject to experimental measurements. The fit allows us to reliably extract the correlations in the cross-section uncertainties which we can then map to the systematic error in the antiproton flux. The second largest contribution to the correlated error stems from the effective acceptance, which the AMS-02 collaboration obtains from a comparison of their detector response between data and Monte Carlo simulation. We estimate the corresponding correlations from the shape (i.e. the "wiggliness") of the correction function employed in an AMS-02 analysis. The full covariance matrix of errors in the AMS-02 antiproton andp/p data, which we derive in this work, is made available in the ancillary files on arXiv.
In the final step, we perform a spectral search for dark matter in the AMS-02 antiproton data, where we fully include the systematic error correlation. The cosmic ray fits are performed independently in two complementary cosmic ray propagation setups (following [9] and [11], respectively). This allows us to draw solid conclusions eliminating further controversies in the assessment of the significance of the antiproton excess.
The remainder of this paper is organized as follows. In Sec. 2 we investigate the nucleoncarbon absorption cross sections and derive the corresponding correlation matrix. In Sec. 3 we collect the various sources of systematic errors to build the overall covariance matrix for the AMS-02 antiproton flux andp/p flux ratio. Finally, in Sec. 4 we derive the implications for the antiproton excess following the two setups mentioned above before drawing our conclusions in Sec. 5. Appendices A and B, respectively, provide additional details on the input cross sections of the Glauber-Gribov model and the derivation of error correlations for further cosmic-ray species used in the analyses. Appendix C summarizes the best-fit values of all involved cosmic-ray propagation parameters in the two setups considered.

Nucleon-nucleus absorption cross sections
In this section we describe the computation of the nucleon-nucleus absorption cross section forpC and pC which is the key ingredient in the assessment of the AMS-02 systematic error. We perform the computation in the framework of the Glauber model [21][22][23]. The theory is formulated based on the eikonal approximation and provides a successful theoretical description for the scattering of moderately relativistic particles off nuclei.
In the following we first introduce the nuclear density function used in Sec. 2.1. We then detail the computation within the Glauber model in Sec. 2.2 while shadowing effects due to inelastic screening are discussed in Sec. 2.3. Parametrizations of the elementary nucleon-nucleon cross section serving as the input for the Glauber-model computations are presented in Sec. 2.4. Finally, we perform a global fit of all input parameters (from the nucleon-nucleon cross sections, the nuclear density function and the inelastic screening) to the respective data in Sec. 2.5 and derive a correlation matrix from the fit as detailed in Sec. 2.6.

Parametrization of nuclear densities
In this work we employ the harmonic oscillator shell model density [26][27][28] which provides a good description for the light nuclei, 3 ≤ A ≤ 16 [29]. It reads where r 2 ch p and r 2 ch A are the mean square charge radii of the proton and nucleus, respectively. We take r 2 ch p = 0.7714 fm 2 and r 2 ch 12 C = 6.1 fm 2 (for carbon) as nominal values [29,30]. We estimate the uncertainty on r 2 ch 12 C to be 0.44 fm 2 by comparing the above value to the one obtained when taking into account nucleon-nucleon repulsion with an expulsion radius of d = 0.9 fm (see [29]).

Computations within the Glauber model
The absorption cross section of a nucleon N on a nucleus A is obtained by subtracting the respective elastic and quasielastic (pA → pA * ) part from the total cross section (see e.g. [31]): Within the Glauber model [21][22][23], neglecting Coulomb effects and spin-orbit interactions [32,33], it is described by (see e.g. [34] for a recent account on the subject): where b is the impact parameter and χ N (b) the nuclear phase-shift function. The last expression in Eq. (4) corresponds to the optical approximation valid in the limit of large A (see e.g. the discussion in Ref. [35]) which is, however, not used in the numerical analysis. 1 The phase-shift function reads where σ NN j is the total nucleon-nucleon cross section, α NN j is the ratio of the real-to-imaginary part of the forward scattering amplitude and β NN j is the slope of the differential inelastic cross , with t being the four-momentum transfer squared. In general, these three quantities depend on the momentum of the incoming nucleon, p lab , in the lab-frame. Furthermore, r T is the transverse part of r.
For a spherical symmetric nuclear density function, ρ(r), the integrals in Eq. (5) can be rewritten to where J n denotes the Bessel function of the first kind. Note that the two integrals in Eq. (6) are real such that the term proportional to α NN j does not contribute to Imχ N . For the harmonic oscillator shell model density, introduced in Sec. 2.1 the integrals in Eq. (6) can be solved analytically leading to: which we use to solve Eq. (4) numerically.

Shadowing corrections
The Glauber model as described above accounts for elastic screening, i.e. multiple nucleon scattering where the intermediate particle is the nucleon. While providing a good description to the experimental data at low momenta, p lab 10 GeV [34], for larger beam momenta it is expected to gradually lose its applicability. This is due to the increasing importance of inelastic screening effects, i.e. multiple scatterings with excited intermediate states [24,25,[36][37][38] (see also [39,40] and references therein for more recent accounts on the subject). These effects increase the shadowing correction leading to a reduction of the total cross section.
Here we follow Ref. [41] where inelastic screening has been considered for the case of neutron projectiles using the expression [38]: where d 2 σ/(dM 2 dt)(t = 0) is the differential diffraction cross section for the process N + N j → N j + X evaluated at t = 0, q L = (M 2 − m 2 N )m N /s is the longitudinal momentum transfer in the production of the (exited intermediate) state X with mass M , M 2 0 = (m N + m π ) 2 1.17 GeV 2 , and F is the form factor: For the purpose of solving Eq. (8) we parametrize the differential diffraction cross section by where the first and second term provide an effective description of the resonant and continuum contribution, respectively. Note that the structure of the resonant part is in general more complex [41][42][43]. 2 However, F is a slowly varying function of M 2 for M 2 < 5 GeV 2 , in particular for p lab 10 GeV where inelastic screening effects are relevant. Hence, the exact shape of the resonant part is not relevant for the computation of the integral in Eq. (8). We treat a 1 and a 2 as effective free parameters choosing a 1 = 10.6, a 2 = 4.05 as nominal values and assigning a relative uncertainty of 20%. These values are consistent with the results of [41] obtained from a fit to the Fermilab data [42] on p + p → p + X and p + d → d + X, where the latter was corrected for binding of the deuteron in an approximate way. 3 The final absorption cross section is obtained from

Nucleon-nucleon cross-sections parametrizations
For the nucleon-nucleon cross sections used in the computations within the Glauber model we consider the parametrizations σ pn = σ asymp 1 + c 16 where ∆ s (p lab ) = (s(p lab ) − s(0)) /GeV 2 with and σ asymp taken from Ref. [44]: The expressions in Eq. (12)(13)(14)(15) are valid at p lab 1 GeV. Furthermore, we parametrize the slope of the differential inelastic cross section in the forward direction by: and assume βp n = βp p and β pn = β pp . We determine the best-fit values and covariance matrices for the involved parameters c 1 , . . . , c 20 , d 1 , . . . , d 6 by performing fits to the cross-section data collected in [45] 4 and the data on the slope parameter [46] 5 . The best-fit values and the covariance matrices are given in Appendix A.

Global fit
In order to compute the correlations in the cross-section uncertainties we perform global fits of thepC and pC absorption cross section within the Glauber model (including inelastic screening) varying the input parameters c i , d i as well as r 2 ch 12 C , a 1 , a 2 according to their nominal values and (correlated) uncertainties (using covariance matrices shown in Appendix A). Note that the parameters r 2 ch 12 C , a 1 , a 2 are varied independently forpC and pC, considering the respective parametrization to be an effective description that might account for additional effects not explicitly considered in the calculation.
For the experimental measurements of thepC and pC absorption cross sections we use data from [47][48][49][50][51] 6 and from [49,50,, respectively. Since we observe (mild) inconsistencies in the data, in addition to the reported experimental errors, we included an additional normalization error of 10% for each data point which we took to be fully correlated within the same experiment (but uncorrelated between different experiments).
We sample the 15-(17-)dimensional input parameter space forpC (pC) with the multimodal nested sampling algorithm MultiNest [74]. 7 The corresponding best fits are shown in Fig. 1 as the solid (dark green) curves. The green band denotes the 1σ uncertainty band. For comparison we display the respective results used in the AMS-02 analyses taken from [75] (gray dashed curve and gray shaded error band). The best-fit values for the input parameters in the global fits are displayed in Tab. 1. It is interesting to note that they are very close to their nominal values (not involvingpC and pC data). A similar observation is made by performing the fit without the χ 2 -contribution frompC and pC data (not shown here) which provides very similar results for thepC and pC absorption cross section as in Fig. 1. We, hence, observe that our computation within the Glauber model including inelastic screening effects provides an excellent description of the measuredpC and pC cross section over the whole considered range of momenta. 8 4 We observe that the cross-section data [45] contain underestimated systematic errors. This becomes obvious from the fact that, with errors taken at face value, the cross-section data are poorly fitted by any smooth function. We, therefore, introduce an additional systematic error of 5% for each data point, which we assumed to be fully correlated within one experiment. In the case ofpp-scattering, where inconsistencies are larger, we furthermore had to multiply the errors by a factor of two. With the described procedure we arrived at a realistic goodness of fit of ∼ 1/dof. 5 While the slope parameter is defined at momentum-transfer t = 0, the experimental data [46] were taken at small, but non-vanishing t. This causes small systematic differences in the normalization of data taken at the same laboratory momentum. In order to evade a degradation of the fit, we, therefore, bin the data with a bin size of 0.2 in units of log 10 p lab (p lab in GeV). For each bin we, furthermore, add a systematic uncertainty of 5% to account for the error caused by the non-vanishing t extraction. 6 The data set [48] originally contains a large number of bins in the momentum range 1.6−3.25 GeV. Since systematic errors of 2−5% dominate over statistical errors, their correlations would be important. Due to their unavailability, we combine the data set into three bins, and conservatively take the error to be 5% in each bin (statistical errors are negligible after the combination). By reducing the number of bins, we minimize the impact of the (unknown) error correlations. 7 We use 5000 livepoints, a sampling efficiency of 0.65 and an evidence tolerance of 10 −18 . 8 This holds strictly forpC. However, we observe a significant deviation of a2 in the case of pC towards a reduction of inelastic screening effects. (Such a reduction is not unexpected and can arise from multi-nucleon correlations, see e.g. [76].) While this leads to a slightly larger cross section for p lab 100 GeV, it does not have a significant effect on the correlation matrix derived from the fit.   Table 1: Best-fit values of the input parameters and the respective χ 2 contributions for the global fits of the absorption cross section forpC and pC.

Correlation matrices for AMS-02 cross sections
From the global fit described above we can infer a covariance matrix for σ abs evaluated at the values of p lab that correspond to the rigidity bins chosen by AMS-02. In this way we can incorporate correlations in the AMS-02 flux measurements arising from σ abs uncertainties. Generally, the covariance matrix of a quantity θ can be computed using where w k is the statistical weight (as given by MultiNest) of the kth point in the statistical ensemble (the MultiNest chain) andθ i = k w k θ k i . For the absorption cross section uncertainties θ k i = σ k abs (p lab,i ), where the p lab,i correspond to the bins of AMS-02. The corresponding correlation matrices read  The correlation matrices of the antiproton-carbon and proton-carbon absorption cross-section uncertainties are provided in the ancillary files on arXiv. Note that in the following we make use of these correlation matrices only, conservatively taking the absolute cross section uncertainties as reported by the AMS-02 Collaboration, see Sec. 3.2.

Error correlations in the AMS-02 antiproton data
In order to determine the error correlations in the AMS-02p andp/p data, we follow a two-step procedure: we first split the systematic errors into individual contributions (as described below). Then, in Sec. 3.2, we derive the correlations for each sub-error and build up the full AMS-02 covariance matrices. Our cosmic-ray fits will also require the AMS-02 covariance matrices for the proton and helium fluxes as well as the B/C ratio as an input. Their calculation (which proceeds analogous to the antiproton case) is described in Appendix B.

Systematic errors
In the following, we will denote relative systematic errors in the antiproton flux by ∆p and in thep/p ratio by ∆p /p .

Unfolding error
Detector resolution effects cause the migration of events into neighboring rigidity bins. This must be corrected for through the unfolding procedure. The choice of the migration matrices (characterizing the migration probabilities) is associated with a systematic error. This unfolding error is ∆p unf = 1% at R < 200 GV and 1.5% at R = 450 GV for the antiproton flux. The error partially cancels in thep/p ratio for which it becomes ∆p /p unf = 1% at R = 1 GV and 0.5% at R > 2 GV [5]. Between the stated rigidity intervals we interpolate logarithmically.

Cross-section error
The (rigidity-dependent) AMS-02 acceptance is sensitive to the fraction of cosmic rays which are absorbed in the detector. The survival probability P N of the incoming particle N (N =p, p) with momentum p is estimated as where n A (p) accounts for the amount of detector material with mass number A which has to be traversed by the incoming cosmic ray, while σ N A abs is the corresponding absorption cross section. We note that the material thickness acquires an effective momentum-dependence due to cuts on track length performed in the AMS-02 analysis. For simplicity, we will neglect subdominant material admixtures and assume that the AMS-02 detector is entirely comprised of carbon, as the corresponding cross section error correlations are expected to be very similar.
We can extract the cross-section error ∆p /p xs in thep/p data by (quadratically) subtracting the unfolding error (as derived above) from the acceptance error as given in [77]. Notice that ∆p /p xs is the quadratic sum of the proton and antiproton contribution to the cross-section uncertainties, From Eq. (22), it, furthermore, follows that at linear order. Here, ∆σp C abs and ∆σ pC abs denote the (absolute) uncertainties in the antiproton and proton absorption cross sections on carbon, respectively, which we extract from [75]. By combining Eqs. (23) and (24), we then extract ∆ p xs and ∆p xs . Of course, only ∆p xs contributes to the systematic error in the antiproton flux.

Scale error
The absolute rigidity scale of measured events can be affected by misalignment of the tracker planes and small uncertainties in the magnetic field map of the inner tracker. The AMS-02 collaboration estimated the corresponding systematic uncertainty by comparing electron and positron rigidity measurements in tracker and electromagnetic calorimeter. The scale error ∆p /p scale is directly given in [77]. Since the rigidity scale error affects protons and antiprotons in the opposite way, one finds ∆p scale = 0.5 ∆p /p scale .

Effective acceptance error
A residual systematic error in the effective folded acceptance is estimated by comparing efficiencies in several detector parts as extracted from Monte Carlo simulation with direct measurements. The effective acceptance error amounts to ∆p eff. acc. = 5% at R = 1 GV and 2% at R > 20 GV [5]. Between 1 and 20 GV, we perform a logarithmic interpolation. The effective acceptance error affects antiprotons and protons in the same way, it cancels in thep/p ratio.

Geomagnetic error
In order to reject indirect cosmic rays produced in the earth atmosphere, AMS-02 applies a rigidity cut above the geomagnetic cutoff. The measured cosmic-ray fluxes (at low rigidity) exhibit a small residual dependence on the exact numerical choice of the cutoff factor. The corresponding systematic error is estimated by varying the cutoff factor in the event selection. The geomagnetic error is ∆p geo = 1% at R = 1 GV and vanishes at R > 2 GV [5], in between we interpolate logarithmically. Since the geomagnetic error is significantly smaller for protons [5], we assume ∆p /p geo = ∆p geo . Template shape error and selection error AMS-02 uses templates to separate signal from background events. Systematic uncertainties arise from the choice of the template shape. The template error mostly affects antiprotons. It is dominated by the contribution related to the charge confusion of incoming protons. In addition, the event selection is affected by a systematic error related to the cuts on the track shape which are used to identify a certain cosmic-ray species. The selection error again mostly affects antiprotons. From [77], we can extract the systematic error on the event number which corresponds to the quadratic sum ∆p In order to derive the individual errors, we use the geomagnetic error from above and assume the following relative size of template and selection errors ∆p /p selection /∆p /p template = 0.48 + 6.5/R 0.78 . This function was chosen to reproduce the ratio of the two errors at several rigidities as given in [5]. Since both errors are mainly relevant for antiprotons, we take ∆p selection = ∆p /p selection and ∆p template = ∆p /p template . When (quadratically) summing the individual systematic errors in the antiproton flux, a tiny mismatch with the published overall systematic error is observed. We eliminate the mismatch by rescaling all individual errors with a correction factor which varies at most by a few percent from unity over the full rigidity range. Figure 2 summarizes the resulting systematic errors in the antiproton flux and thep/p ratio as a function of rigidity. 9

Covariance matrices for AMS-02 errors
After splitting the AMS-02 systematic error into its various components, we will now assign correlation matrices ρp a , ρp /p a (a = unf, xs, scale, . . . ) to each of the sub-errors. The leading uncertainty (in the regime where the systematic error dominates over the statistical error) derives from the absorption cross sections. The reported, i.e. absorption-corrected (anti)proton flux scales inversely with the (anti)proton survival probability which was defined in Eq. (22). Therefore, at linear order in the cross-section error, the correlation matrices ρp xs and ρ p xs are identical to the correlation matrices of uncertainties in the absorption cross sections σp C abs and σ pC abs , respectively, (assuming that the same rigidity bins are chosen) which were derived in Sec. 2.6. Note that the correlation matrix ρp /p xs contains the contributions from ρp xs and ρ p xs (weighted by the relative magnitude of the antiproton and proton cross-section uncertainties).
For the remaining uncertainties, we follow the approach of [78,20] and define the following correlation matrix ρp a ij = exp − 1 2 Apart from the cross-section error, the effective acceptance error plays a significant role. It is derived from a data vs. Monte Carlo comparison and may receive contributions from mismodeling of efficiencies in various detector parts or from small errors in the detector composition model. Since it amounts to a collection of different residual errors, it is difficult to gain any intuitive insights into the corresponding error correlations. However, a realistic estimate of the correlation length can be obtained by analyzing the 'wiggliness' of the data/MC correction function employed by AMS-02. In the AMS-02 analysis, the latter is determined from proton data and then assumed to be identical for antiprotons (this is why the effective acceptance correction and the corresponding error cancel in thep/p ratio). We extract the data/MC correction function from the proton flux analysis in the PhD thesis [79] by taking the ratio of effective and geometric efficiency. 10 In Fig. 3 we fit a polynomial of twelfth degree in log 10 R to the error function which well reproduces its overall shape. 11 From this fit, we can directly extract the correlations in the effective acceptance error. By a subsequent fit of the correlations to the form (25), we finally obtain 12 eff. acc. = 0.15 .  As an alternative approach to derive the correlation length of the effective acceptance, we have considered the correction function shown in Fig. 3 to be an estimate for the systematic uncertainty of the effective acceptance itself. To this end, we have defined a likelihood for the parameter eff. acc. from Eq. (25) and two nuisance parameters describing the overall size of the systematic uncertainty and a pure re-normalization uncertainty. A similar strategy was used in [11]. Profiling over the three parameters we obtain eff. acc. ∼ 0.1 which is comparable with the result above.
The remaining uncertainties play a subleading role. They are always subdominant to either the two previously discussed systematic errors or the statistical error. For those errors, we refrain from a detailed analysis and adopt the correlation lengths estimated in [20]: In the next step, we build the covariance matrix for each sub-error by multiplying the entries of the correlation matrix by the AMS-02 errors (as displayed in Fig. 2) in the corresponding bins. The covariance matrices for each sub-error are then added to build the AMS-02 covariance matrix for the full systematic error. Figure 4 illustrates the overall correlations in the AMS-02 antiproton andp/p systematic errors as derived by our method. It can be seen that the systematic error in the antiproton flux is correlated on a shorter length scale compared to the error in thep/p ratio. This is because the effective acceptance error, which has a relatively short correlation length, affects the antiproton flux, but not thep/p ratio. The full AMS-02 covariance also containing the statistical error is provided in the ancillary files.

Implications for the AMS-02 antiproton excess
The AMS-02 error correlations can now be used to gain new insights into cosmic-ray spectra. Of particular interest is the question how the correlations affect the interpretation of the antiproton excess at R = 10−20 GV. The latter has been considered as a possible dark matter signal in a number of studies [6][7][8][9][10][11][12][13]. At the same time, the significance of the excess is rather controversial (ranging from 1 to 5 σ within the mentioned references). The previous studies took the systematic errors in the AMS-02 antiproton flux to be uncorrelated (or modeled correlations in a simplistic way).
In the following, we will reinvestigate the antiproton excess in Sec. 4.3, fully including the derived correlations in the AMS-02 systematic errors. We decided to perform two complementary likelihood analyses on the AMS-02 data. The two analyses differ substantially in the modeling of cosmic-ray propagation and in the considered species. Hence, we can directly verify the robustness of the conclusions we draw on the correlations with respect to the propagation model. Before we describe our analysis methods in Sec. 4.2 we will briefly review the production and propagation of cosmic rays in Sec. 4.1.

Cosmic-ray production and propagation
Cosmic rays are mainly composed of galactic matter which has been energized by supernova shock acceleration. This so-called primary component includes protons, helium and heavier nuclei like carbon and oxygen. Primaries, when they propagate through the Galaxy, induce secondary cosmic rays by scattering processes in the interstellar disc. The source term for a secondary a, which denotes its differential production rate per volume, time and energy, takes the form where i runs over the relevant primary cosmic-ray species (Φ i and T denote their flux and kinetic energy, respectively) and j over the target nuclei in the galactic disc (n j denotes their number density). Secondary production is the main source of certain nuclear cosmic rays like boron, but also plays a major role in the generation of antimatter. For antiprotons, which are the main concern of this work, we will employ the differential production cross sections dσ ij→p /dT (with i, j = p, He) derived in [16,18]. Since production cross sections are only known to a few percent precision, they comprise an important source of systematic error in the modeling of antiproton fluxes. These will be fully included in our cosmic-ray analysis. One of our fits will also include the boron-to-carbon ratio B/C. The boron production cross sections with uncertainties are taken from [9] (see also [81,82]).
In addition to the secondary background, primary antiprotons can be induced by dark matter annihilation. The primary source term reads where ρ χ , m χ and σv stand for the dark matter energy density, mass, and annihilation cross section. The antiproton energy spectrum per annihilation is denoted by dN/dT . In this work we will consider χχ →bb as an exemplary annihilation channel and extract dN/dT from [83]. The scan range for the dark matter mass is taken as m χ =10-10000 GeV. Independent of their origin, cosmic rays follow complicated trajectories on their passage through the galaxy. Scattering on magnetic field inhomogeneities induces a random walk which can be described as spatial diffusion. Convective winds may blow charged particles away from the galactic disc. In addition, cosmic-ray interactions with matter lead to energy losses and annihilation, while Alfvén waves induce reacceleration. The collection of processes is encoded in the diffusion equation. In this work, we will consider two approaches to solve the diffusion equation following our previous studies [9,11]: the semi-analytic solution in the two-zone diffusion model [84][85][86] and a fully numerical solution based on the GALPROP code [87][88][89]. In both schemes, diffusion is assumed to occur homogeneously and isotropically in the galactic halo. Magnetohydrodynamics considerations suggest a power-law form of the diffusion coefficient K (in the GALPROP code, K is denoted as D xx ), where R is the rigidity and β the velocity of the cosmic-ray particle. While η = 1 in the standard implementation (as well as in our previous works [9,11]), here, we will include η as a free fit parameter. This change is partly motivated by recent studies [90][91][92][93] which observed a substantial improvement in the fit to secondary nuclear cosmic rays within the diffusion model with free η. In addition, the freedom of η can be viewed as a conservative choice since it tends to slightly reduce the significance of a potential dark-matter signal (we will return to this point later). From a physical perspective, an increase of the diffusion coefficient (negative η) towards low rigidity can be motivated by wave damping on cosmic rays [94]. Following our previous work [9,11] we also include a break in the power law index δ at R ∼ 300 GV as required to fit nuclear primary and secondary cosmic rays. The high-energy break does, however, virtually not affect our dark-matter analysis. Both, the two-zone diffusion model and the numerical implementation of GALPROP, allow for convective winds perpendicular to the galactic plane. Both include energy losses, annihilation and reacceleration processes, although slight differences in the implementation exist (concerning e.g. the modeling of the interstellar material and the spatial extension of the reacceleration zone). For details we refer to the original references [87][88][89][84][85][86]. Some custom modifications to the default setups have been described in our previous works [9,11].

Methodology
In order to investigate the significance a possible dark-matter signal in the AMS-02 antiproton data, we consider two complementary setups which we describe in the following. Within both setups, we include the derived AMS-02 covariance matrices of systematic errors for all species in the fit.

Setup 1
The first setup implements the approach described in [9]. It is based on the two-zone diffusion model of cosmic-ray propagation. The primary fluxes of p, He, C, N, O, Ne, Mg, Si are taken as an input to predict the secondary fluxes of antiprotons and boron. In addition, a primary antiproton component from dark-matter annihilation is included. Solar modulation is taken into account through the improved force-field approximation [95], which -in order to describe charge-breaking effects -contains one parameter in addition to the Fisk potential.
In setup 1, a simultaneous fit to the AMS-02 B/C ratio, the AMS-02 antiproton flux and the antiproton flux ratio between AMS-02 and PAMELA is performed. The combination of B/C and p is needed to determine the propagation parameters and the significance of a dark-matter excess, while the AMS/PAMELA flux ratio fixes the charge breaking parameter in the solar modulation (see [9] for details). The total goodness of fit is determined by which is evaluated with and without a dark matter component. Uncertainties in the production cross sections of antiprotons and boron are fully included via covariance matrices [9].

Setup 2
The second setup follows the approach of [11] and employs the GALPROP code for cosmic-ray propagation. The primary source terms of proton and helium are modeled as broken power laws. In addition, a primary antiproton source term from dark-matter annihilation is included. From this input, the full network of cosmic-ray propagation, scattering and propagation of secondaries is employed to determine the proton, helium and secondary antiproton flux. A major difference compared to setup 1 is that the propagation parameters are constrained by primary fluxes and thep/p ratio instead of the B/C ratio. Solar modulation is implemented through the standard force-field approximation (with individual Fisk potentials for both charge signs). A more refined  treatment can be evaded as the AMS-02 data is cut at rigidity R = 5 GV. The propagation and solar modulation parameters as well as the significance of a dark-matter excess is determined by a simultaneous fit to the AMS-02 and Voyager data on protons and helium as well as the AMS-02 antiproton-over-proton data. Uncertainties in the production cross sections of antiprotons are again taken into account through the covariance matrix matrix from [9]. 13

Results
The results of our fits are summarized in Tab. 2, which provides the goodness-of-fit, dark-matter parameters and significances for the best-fit points of the setups 1 and 2. The χ 2 values with and without a dark-matter component in the antiproton flux are included in each column, the latter is displayed in parentheses. In order to highlight the impact of the derived AMS-02 correlation matrices, we compare results with and without the error correlations. The corresponding cosmicray spectra and residuals are depicted in Fig. 5 and Fig. 6. The best-fit propagation parameters within the two setups are listed in Appendix C. The stated local significance of the dark-matter signal refers to a ∆χ 2 -test for one degree of freedom. For the global significance we take into account the look-elsewhere effect by evaluating the probability distribution among mock-data sets created under the background-only hypothesis. 14 In order to make the connection to previous studies, let us first turn to the fit results without including error correlations. Up to the extra propagation parameter η (cf. Eq. (31)), this case   corresponds to the default configuration in [9] (for setup 1) 15 and to the configuration "XS uncertainty by covariance matrix" in [11] (for setup 2). Not surprisingly, we qualitatively reproduce the results given in these references. An antiproton excess is observed at R = 10-20 GV. The latter is compatible with a dark-matter particle of mass m χ ∼ 80 GeV and annihilation cross section σv ∼ 10 −26 cm 3 /s into bottom quarks. However, in setup 1, the global significance is only ∼ 1 σ, while it reaches ∼ 2 σ in setup 2. In both setups, the significance is slightly smaller compared to [9,11] which is due to the additional freedom in the propagation. The extra diffusion parameter η allows for a stronger "bending" of cosmic-ray fluxes towards low energy and, -setup 2 - Figure 6: Antiproton-to-proton ratio (top) and proton flux (bottom) of the fit without (left) and with dark matter (right) within setup 2. The solid red and blue curves denote the best-fit spectra at the top of the atmosphere with and without correlations in the AMS-02 errors, respectively. The dashed curves denote the corresponding interstellar fluxes. We display the cosmic-ray measurements of AMS-02 (proton and antiproton-to-proton ratio) and Voyager (proton). The cosmic-ray fit of the AMS-02 data is restricted to rigidities between the dotted black lines. Residuals are shown only for the AMS-02 data points. Error bars denote the statistical and systematic uncertainties (according to the diagonal entries of the total experimental covariance matrix).
hence, absorbs a small fraction of the excess. When we turn to the fits including the correlations in the AMS-02 systematic errors, we observe an overall increase in the χ 2 values. This is to be expected, since -with a realistic modeling of systematic errors -one expects a total χ 2 comparable to the number of degrees of freedom (dof). Setup 1 nicely fulfills this criterion, which gives us confidence that the derived AMS-02 error correlations are indeed realistic. In setup 2 a somewhat higher χ 2 /dof ∼ 1.5 occurs. This could possibly stem from a mild underestimation of systematic errors in the proton and helium data of AMS-02. We wish to point out that the unofficial AMS-02 helium analysis performed in the PhD thesis [98] indeed derived larger uncertainties compared to the published data. Alternatively, it could indicate that we slightly overestimated the correlation length in the proton and helium systematic errors. Even if this is the case, it would not affect our conclusions on the dark-matter excess as we have explicitly verified. 16 Our key result is that the significance of the dark matter excess decreases substantially, once we include error correlations in the AMS-02 data. In setup 1 the preference for a darkmatter signal disappears completely, but even in setup 2 the global significance drops below 1 σ. Correspondingly, the best-fit dark-matter signal is reduced by about a factor of 2 in both setups. It is thus obvious that the systematic errors of AMS-02 provide sufficient freedom to absorb the antiproton excess, once their correlations are properly taken into account.
In both setups, the correlations in the absorption cross section uncertainties, which we derived in Sec. 2, play an important role. However, there is also an interesting distinction: in setup 1, the effective acceptance error in the antiproton flux contributes to reducing the significance of the antiproton excess. In setup 2, thep/p ratio is employed in the fit (instead of the antiproton flux). Since the effective acceptance error cancels inp/p one might think that it is irrelevant in this case. However, we observe, that the fit then takes advantage of the effective acceptance error in the proton flux. This can be seen in Fig. 6, where the excess inp/p partly shifts into the proton channel, once correlations are included.
We finally wish to emphasize that our conclusion, namely that the antiproton excess is explainable by systematic effects, is not build on the systematic error correlations alone. Rather, the full interplay between all uncertainties matters -those in the cosmic-ray propagation, those in the antiproton production cross sections and those in the AMS-02 systematic errors. In order to emphasize this point, let us turn to an example: we have tested the significance of the antiproton excess in setup 2 without including production cross section uncertainties and in a more restrictive propagation setup (without the diffusion parameter η). In this case, the correlations in the AMS-02 systematic errors had a contrary effect: compared to uncorrelated systematic errors, they strongly increased the significance of the antiproton excess to ∼ 5 σ. The same observation has already been made in [11]. We should thus refine our previous statement: the AMS-02 systematic errors alone might not absorb the observed spectral feature in the antiproton flux. However, in combination with the other main uncertainties, the correlated systematic errors can account for the previously seen excess such that there is no longer a preference for a dark-matter signal.
Finally, we remark that our results should not be understood as to exclude the dark-matter interpretation of the antiproton feature at R = 10-20 GV. Any improvement in the description of cosmic-ray propagation and the modeling of antiproton production cross sections will reduce the uncertainties and help to determine whether an antiproton excess exists. Both, the darkmatter interpretation as well as the interpretation as a combination of systematic effects will thus undergo further scrutiny in the next years. Obviously, additional information on systematic errors, provided by the AMS-02 collaboration, would be extremely useful in this regard.

Conclusions
Four years ago, the AMS-02 collaboration has published their first measurement of the cosmic ray antiproton flux. After background subtraction, it seemingly revealed a small residual component in the spectrum at rigidity R = 10−20 GV. Since then, the antiproton excess has been subject to major controversy in the community. On the one hand, its possible explanation in terms of annihilating dark matter has caused a wave of excitement. In particular since the favored mass m χ = 50−100 GeV and (hadronic) annihilation cross section σv = 0.5−3 × 10 −26 cm 3 s −1 could also tentatively be linked to the galactic center gamma ray excess. On the other hand, systematic effects have been a matter of concern.
Since the antiproton excess was first observed, major improvements in the modeling of the antiproton background and the description of cosmic-ray propagation have been made. But one decisive piece has been missing so far: the full covariance matrix of systematic errors in the AMS-02 antiproton data. Correlated systematic errors are so important because they can cause features in the residuals and potentially 'fake' a dark matter excess. Therefore, in this work, we comprehensively derived the correlations in the AMS-02 data with the purpose of further scrutinizing the excess.
In the rigidity range R = 10−20 GV, the systematic error dominantly descends from cosmic ray absorption in the AMS-02 detector. The relevant nuclear cross sections which enter the AMS-02 simulation carry major uncertainties. These directly translate to the measured fluxes of antiprotons, protons and heavier cosmic rays. We carefully computed the mentioned absorption cross sections within the Glauber model and performed a global fit to the data from nucleonnucleus and nucleon-nucleon scattering experiments, the latter of which serve as input to the theoretical model. Including inelastic shadowing effects we obtained an excellent fit to the data. This computation also allowed us to reliably infer the correlations of the absorption cross section error in the AMS-02 systematics. Furthermore, we found a data-driven estimate for correlations in the effective acceptance error and used covariance functions for subdominant contributions related to the AMS trigger and unfolding procedure. From the correlations of individual errors, we constructed the covariance matrix for the total experimental error in the antiproton and proton flux as well as thep/p ratio. Similarly, we also derived the corresponding covariances for He and B/C which also enter our cosmic ray analysis. The final AMS-02 covariance matrices are made available for public use in the ancillary files on arXiv. 17 In the final step, we turned to a reevaluation of the antiproton excess. To this end we performed a spectral search for dark matter in the AMS-02 antiproton data, where we fully took into account the derived correlations in the AMS-02 systematic errors. Our fits also incorporate a detailed modeling of the uncertainties in the antiproton production cross sections which affect the prediction of the astrophysical antiproton background. The results were obtained in two complementary propagation setups with a conservative choice of propagation parameters. In the first setup, we observe that the antiproton excess disappears, once the AMS-02 systematic error correlations are taken into account -without correlations a local (global) excess of 1.9 σ (0.8 σ) had been found. In the second setup, the systematic error correlations reduce the significance of the excess from 2.6 σ (1.8 σ) locally (globally) to 1.8 σ (0.5 σ). We conclude, that with all uncertainties properly taken into account, the AMS-02 data do not prefer a dark matter interpretation at this stage. The fact that we obtain consistent results in both propagation setups gives us further confidence in the robustness of this conclusion.
While our findings emphasize the importance of error correlations in the AMS-02 data, we remark that including the systematic errors in the background model turns out to be equally crucial. We found that ignoring uncertainties in the antiproton production cross sections can even lead to the (wrong) conclusion that AMS-02 systematic error correlations increase the significance of the antiproton excess. Only the interplay between correlated errors in the data, uncertainties in the antiproton production and sufficient freedom in the cosmic ray propagation allows the fit to fully absorb the antiproton excess (without a dark-matter component).
Our result should not be understood as to exclude a dark matter candidate with mass m χ = 50 − 100 GeV and thermal annihilation cross section. Rather, with the full account of systematic errors, our cosmic-ray fit is only weakly sensitive to dark matter of this type. It can simply not distinguish between a dark-matter excess at the 10%-level and a correlated fluctuation in the systematics. Further reducing these systematics must be a main objective in the next years. It is of paramount importance to strengthen the efforts to precisely measure cosmic-ray cross sections at accelerators. At the same time, existing and upcoming nuclear cosmic-ray data will help to further advance the physics of cosmic-ray propagation. In the precision era, which cosmic-ray physics has entered, it is no longer the statistics, but the control of systematic errors which decides on the prospects of indirect dark-matter detection. 2017, funded by MIUR, and 'The Anisotropic Dark Universe' No. CSTO161409, funded by Compagnia di Sanpaolo and University of Turin. Simulations were performed with computing resources granted by RWTH Aachen University under project rwth0085.

A Nucleon-nucleon cross-section fit
In this appendix we display the best-fit values and covariance matrices for the parameterizations of the nucleon-nucleon cross sections and the slopes of the differential inelastic cross sections as described in Sec.
The covariance matrices read: B Error correlations in the AMS-02 proton, helium and B/C data The calculation of the covariance matrices of AMS proton, helium and B/C errors proceeds analogously to those for antiprotons andp/p which was described in Sec. 3. We first split the systematic errors into their components.

Proton flux
In the proton case, the error related to absorption cross sections is 1% at R = 1 GV, 0.6% from R = 10−300 GV and 0.8% at R = 1800 GV [80]. Between the given rigidities we interpolate logarithmically. The effective acceptance error is obtained by (quadratically) subtracting the cross-section error from the acceptance error given in the supplementary material of [80]. The unfolding, scale and trigger uncertainties can directly be taken from the same reference.
Helium flux cross-section errors are taken to be 1% for helium over the full rigidity range [99]. Effective acceptance, unfolding, scale and trigger error are extracted in the same way as for the proton flux from the supplementary material of [99].

B/C
For B/C we employ the effective acceptance error of 1% over the full rigidity range [100]. The cross-section error is obtained by (quadratically) subtracting the effective acceptance error from the total acceptance error given in the supplementary material of [100]. The systematic unfolding, scale, trigger and background subtraction errors are taken from the same reference. The background subtraction error (which was absent in the proton and helium fluxes) is related to the spallation of heavier cosmic-ray species within the AMS detector.
After characterizing the individual error components, we assign a correlation matrix to each of them (following our approach in Sec. 3.2). The cross-section error correlations for the proton flux can directly be taken from our fit (see Sec. 2.6). The cross-section error correlations for helium and B/C are expected to be of similar shape. Therefore, we refrain from performing cross-section fits for the absorption of these species and rather adopt the correlations from the σ pC -fit. However, the natural unit for correlations is the momentum (per nucleon) which (approximately) corresponds to half the rigidity in the case of nuclei. Since the AMS data are provided in terms of the rigidity, we thus have to stretch the correlations from the proton case by a factor of two in order to apply them to helium and B/C. The remaining correlations matrices are again taken to be of the form Eq. (25). We consistently choose the same correlation lengths as in Eqs. (27) and (26). For the trigger error which contains the geomagnetic error, we take trigger = 1 [78]. Similarly, for the background subtraction error which only affects B/C we assume background = 1 (since the contamination should be strongly correlated between a few neighboring bins). Table 3 summarizes the best-fit propagation parameters for the two setups considered. For the definition of all parameters (except for η which was defined in Eq. (31)) and the respective details of the propagation model we refer to [9] (for setup 1) and [11] (for setup 2). In setup 1, a degeneracy between K 0 and L arises and, therefore, L was fixed to the lowest value L = 4 kpc suggested by positron data [9] (the value L = 4 kpc is also motivated from a recent Be/B analysis in the same propagation setup [101], see also [102]). In addition to the parameters given in the table, in both setups a rigidity break in the diffusion coefficient is taken into account in the same way as in [9] and [11], i.e. ∆δ = 0.157, R b = 275 GV, s = 0.074 and δ − δ 2 = 0.12, R 1 = 300 GV, respectively. Note that although some parameters in the two setups are equivalent we choose the very same notation as in the respective references in order to keep the one-to-one correspondence.  Table 3: Propagation, solar-modulation and dark-matter parameters yielding the best fit within the setup 1 (left) and 2 (right) with and without including correlations in the AMS-02 systematic errors. The values given in each column refer to the fits with (without) dark matter.