## A new evaluation of the antiproton production cross section for cosmic ray studies

Mattia di Mauro 1 2 , , Fiorenza Donato , Andreas Goudelis , Pasquale Dario Serpico 1 2 2 1 Physics Department, Torino University, and Istituto Nazionale di Fisica Nucleare, Sezione di Torino, via Giuria 1, 10125 Torino, Italy and 2 LAPTh, Univ. de Savoie, CNRS, B.P.110, Annecy-le-Vieux F-74941, France

Theoretical predictions for the cosmic antiproton spectrum currently fall short of the corresponding experimental level of accuracy. Among the main sources of uncertainty are the antiproton production cross sections in cosmic ray inelastic interactions. We analyse existing data on antiproton production in pp scattering, including for the first time the measurements performed by the NA49 Collaboration. We compute the antiproton spectrum finding that in the energy range where data are available (antiproton energies of about 4-550 GeV) different approaches lead to almost equivalent results, with an uncertainty of 10-20%. Extrapolations outside this region lead to different estimates, with the uncertainties reaching the 50% level around 1 TeV, degrading the diagnostic power of the antiproton channel at those energies. We also comment on the uncertainties in the antiproton production source term coming from nuclei heavier than protons and from neutrons produced in pp scatterings, and point out the need for dedicated experimental campaigns for all processes involving antiproton production in collisions of light nuclei.

PACS numbers: 13.85.-t,13.85.Ni98.70.Sa

## I. INTRODUCTION

Cosmic ray (CR) antiprotons are a remarkable diagnostic tool for astroparticle physics. The bulk of the measured flux is certainly consistent with a purely secondary origin in CR collisions onto interstellar medium gas, but additional primary components are not excluded, either of astrophysical origin (see for instance [1]) or of exotic nature, such as dark matter annihilation or decay [2]. At very least, antiprotons provide a consistency check for the current understanding of galactic CR modeling and can narrow down propagation parameters (see e.g. [3-5]).

This tool is however only as sharp as the uncertainties entering the background (i.e. the secondary component) and signal (i.e. the primary component) computations are robust. Statistical and systematic errors reported by the PAMELA collaboration [6] are already at the 10% level up to the 10 GeV scale, below the theoretical error. In a short time, AMS-02 [7] is expected to provide significantly higher precision, calling for a reassessment of the theoretical predictions.

The contribution of different processes to the ¯ secp ondary yield has been studied in the past, see e.g. [810]. In [9], for instance, the uncertainties on the production cross sections were estimated to be ∼ 25%, and already identified as the limiting factor in theoretical predictions (see also [10] for similar considerations). In practice, nuclei heavier than protons and helium only contribute at a few percent level (see e.g. [8]), thus playing a very marginal role, either as projectiles or targets, in the antiproton production. Reactions involving helium (p-He, He-p, He-He) represent a sizable fraction of the total yield, easily reaching ∼ 50% at low energies [9].

While for processes involving helium nuclei no data is available, the situation is different for the proton-proton case, where there are several experimental studies. The latest re-evaluation of the antiproton production yield in pp collsions was reported in [11], while the Tan &amp; Ng parameterization [12] is still largely used, despite being more than 30 years old. The reason is that, until recently, the available dataset was limited to data only collected in the sixties and seventies. In the last decade, however, two more experimental datasets have become available: the BRAHMS data [13] and-more important for the energies of interest for AMS-02 applications-the NA49 results collected at the CERN Super Proton Synchrotron (SPS) [14].

Given the importance of these nuclear data for new measurements in astroparticle physics, it seems thus timely to re-evaluate the antiproton production cross section in pp collisions in light of this newly available information. In this paper we engage ourselves in this task, in order to provide the community with a parametrization for the inclusive antiproton production cross section as well as with a reliable assessment of the corresponding uncertainties that should be taken into account in CR studies.

The outline of the paper goes as follows: in section II we set up the relevant formalism, present the experimental data that will be used subsequently and describe the analysis methods. In section III we present our results. We begin by validating our analysis framework by reproducing existing results in the literature and then move on to evaluate the inclusive antiproton production cross section, first relying solely on the novel NA49 data and then on the full set of available measurements. In section IV we briefly comment on the impact of other contributions entering the secondary antiproton source term, namely of antineutrons and helium nuclei. Finally, in section V we discuss our key results and present our conclusions. Two appendices follow, where we present some-standard but useful-kinematics used in our analysis as well as a new evaluation of the total, elastic and inelastic pp scattering cross sections that we performed for the energy range of

## II. FRAMEWORK, DATA AND METHODS

## A. Theoretical framework

CR protons interact with the interstellar medium (ISM) and may produce secondary antiprotons. Different channels are involved, with the dominant one being the CR proton flux collisions with the target hydrogen gas ( pp ). The corresponding source term is the convolution of the antiproton production cross section dσ p p → ¯ p dE ¯ p ( E , E p ¯ p ) and the interstellar CR proton energy spectrum

$$q _ { \bar { p } } ^ { p p } ( E _ { \bar { p } } ) & = \int _ { E _ { \text{th} } } ^ { + \infty } \frac { d \sigma _ { p \, p \rightarrow \bar { p } } } { d E _ { \bar { p } } } ( E _ { p }, E _ { \bar { p } } ) n _ { H } ( 4 \pi \Phi _ { p } ( E _ { p } ) ) d E _ { p }, \\ \text{when} \ \infty \ \sim \ \text{is} \ \mathbb { m } \ \mathbb { m } \ \text{when} \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m } \ \mathbb { m }$$

where n H is the ISM hydrogen density, Φ p is the CR proton flux, E p and E ¯ p are the CR proton and antiproton energies, and E th the production threshold energy equal to 7 m p . The overall ISM composition is H:He:C=1:0.1:5 × 10 -4 cm -3 [15]. Whenever needed for illustrative purposes, we will fix Φ p to the fit to the preliminary AMS-02 data [16] reported in [17].

The differential cross section dσ pp → ¯ p /dE ¯ p is in turn the integral over the angle ϑ between the incoming proton and the final state antiproton momenta

$$\frac { d \sigma _ { p p \to \bar { p } } } { d E _ { \bar { p } } } ( E _ { p }, E _ { \bar { p } } ) = 2 \pi p _ { \bar { p } } \int _ { \vartheta _ { \min } } ^ { \vartheta _ { \max } } E _ { \bar { p } } \frac { d ^ { 3 } \sigma } { d p _ { \bar { p } } ^ { 3 } } d ( - \cos \vartheta ), \ \ ( 2 ) \quad \text{(with} \quad \\ \text{TeV}, \quad \text{$\dots$}$$

where p ¯ p = √ E 2 ¯ p -m 2 p , ϑ min = 0 ◦ and the expression for ϑ max is given by Eq.(A15) in Appendix A. The integral in Eq.(2) is computed in the galactic frame at fixed antiproton energy E ¯ p . Its integrand represents the 'Lorentz-invariant distribution function' 1 for the process p + p → ¯+ p X , i.e. the inclusive antiproton production. Inclusive cross sections for processes of the form a + b → c + X can be described in terms of the Lorentzinvariant distribution function

$$f ( a + b \to c + X ) = E _ { c } \frac { d ^ { 3 } \sigma } { d p _ { c } ^ { 3 } } = \frac { E _ { c } } { \pi } \frac { d ^ { 2 } \sigma } { d p _ { L } d p _ { T } ^ { 2 } } = \frac { d ^ { 2 } \sigma } { \pi d y d p _ { T } ^ { 2 } } \,, \quad \begin{array} { c } \text{tion c} \\ 2 3 \text{], r} \\ \text{energy} \\ ( 3 ) \text{cover} \end{array}.$$

where p L , p T and y are respectively the longitudinal and transverse momentum and the rapidity of particle c . Traditionally, the independent variables most frequently used to parameterise this quantity are

- · the centre-of-mass (CM) energy √ s = √ 2 m E p ( p + m p ), which is uniquely fixed by the total incident proton energy in the lab frame, E p ;

1 For simplicity, in the rest of this work we will interchange this term with 'cross section'.

- · p T , the antiproton transverse momentum;
- · the so-called 'radial scaling' variable x R , defined as

$$x _ { R } = \frac { E _ { \bar { p } } ^ { * } } { E _ { \bar { p }, \max } ^ { * } } \,, \quad \quad \quad ( 4 )$$

where E ∗ ¯ p is the antiproton energy and E ∗ ¯ max p, is the maximal energy it can acquire, both defined in the CM frame. The maximal antiproton energy is (see Appendix A for details):

$$E _ { \bar { p }, \max } ^ { * } = \frac { s - 8 m _ { p } ^ { 2 } } { 2 \sqrt { s } } \,,$$

which-from the condition E ∗ ¯ max p, ≥ m p -also implies the threshold energy for the incident proton in the lab frame, E p ≥ 7 m p .

The inclusive antiproton production cross section cannot be computed from first principles. Our primary goal in this work is to provide reliable estimates for the magnitude and the uncertainties of the invariant distribution (3). Our results will be presented mostly in the form of suitable fitting functions. However, we also want to test how reasonable the Ansatz of the chosen functional form(s) is. To that purpose, we will also compare the fitted functions to an 'agnostic' spline interpolation of the data, which only requires a smooth, piecewise functional dependence. We will mainly focus on antiprotons with energies ranging from a few GeV and up to O (1) TeV, with the upper limit of this interval corresponding roughly to the highest energy that can be probed by AMS-02 and the lower one to the point where astrophysical uncertainties become so large that they constitute the dominant limiting factor in CR studies, a point which we will also briefly comment upon in section IV.

## B. Experimental data

In order to estimate the inclusive antiproton production cross section, we consider the datasets [13, 14, 1823], reported in Table I in terms of the centre-of-mass energy √ s and ( p T , x R ) regions (often not rectangular!) covered by each experiment. Note that not all experiments provide data in terms of these kinematic variables; in those cases, the data were first converted in terms of { s, p T , x R } . We report in Appendix A the straightforward but somewhat lengthy derivation of the transformation equations. The data are graphically illustrated in Fig.1. In the left panel, the cross section is shown as a function of E LAB ¯ p , for different combinations of p T and x R . In the right panel, the same data are seen in the p T -x R plane. The NA49 data cover wide ranges in both p T and x R , and describe lab antiproton energies from about 8 GeV up to 70 GeV.

Compared to the previous works [11, 12], the analyses of the NA49 [14] and BRAHMS [13] datasets are new to

TABLE I. Datasets used in our analysis along with their corresponding √ s values and ( p T , x R ) regions.

| Experiment                        | √ s (GeV)                              | p T (GeV)         | x R                 |
|-----------------------------------|----------------------------------------|-------------------|---------------------|
| Dekkers et al , CERN 1965 [18]    | 6 . 1, 6 . 7                           | (0 ., 0 . 79)     | (0 . 34 , 0 . 65)   |
| Allaby et al , CERN 1970 [19]     | 6 . 15                                 | (0 . 05 , 0 . 90) | (0 . 40 , 0 . 94)   |
| Capiluppi et al , CERN 1974 [20]  | 23 . 3, 30 . 6, 44 . 6, 53 . 0, 62 . 7 | (0 . 18 , 1 . 29) | (0 . 06 , 0 . 43)   |
| Guettler et al , CERN 1976 [21]   | 23 . 0, 31 . 0, 45 . 0, 53 . 0, 63 . 0 | (0 . 12 , 0 . 47) | (0 . 036 , 0 . 092) |
| Johnson et al , FNAL 1978 [22]    | 13 . 8, 19 . 4, 27 . 4                 | (0 . 25 , 0 . 75) | (0 . 31 , 0 . 55)   |
| Antreasyan et al , FNAL 1979 [23] | 19 . 4, 23 . 8, 27 . 4                 | (0 . 77 , 6 . 15) | (0 . 08 , 0 . 58)   |
| BRAHMS, BNL 2008 [13]             | 200                                    | (0 . 82 , 3 . 97) | (0 . 11 , 0 . 39)   |
| NA49, CERN 2010 [14]              | 17 . 3                                 | (0 . 10 , 1 . 50) | (0 . 11 , 0 . 44)   |

this paper. Note that the BRAHMS centre-of-mass energy corresponds to an incident proton energy of roughly 21 TeV in the lab frame, which lies somewhat beyond the energy region of interest for our work. Given the absence of data for incident proton energies above ∼ 200 GeV, however, we have included this dataset since it can help in guiding the high-energy extrapolation of the fit to physical values. It is worth stressing that in the more interesting tens of GeV region for the antiproton laboratory energy, the major impact will be provided by far by the NA49 data.

Another important conceptual issue concerns the possibility to combine data-whose quality and robustness of error assessment is very diverse-in a global fit. There is no simple answer to this question: on one hand there are some systematic effects that are certainly present in the old data and hard to estimate and correct for. A known example is provided by the feed-down effect. A significant fraction of antiproton production (easily of O (20%)) comes from strange hyperon (Λ, Σ) decays, whose decay lengths are comparable to or larger than length scales of current micro-vertex detections or precision tracking. This effect was taken into account in the NA49 data analysis, where the contribution from hyperons has been subtracted from the measured yields . 2 For older experiments, no such correction was performed: while in some cases-as for the CERN ISR-it may be argued that reasonable estimates make the expected correction negligible, for fixed-target experiments covering an extended range of lab momenta the situation is somewhat more complicated. No a priori correction has been applied in the following for this effect, especially since ex-novo simulations of trajectories through the detectors and the collimators would be needed for robust estimates. For a semi-quantitative discussion, we address the reader to [14]. However, in deriving global fits, we shall allow for experiment-dependent renormalizations, which may account (at least in an averaged way) for such a correction, see below.

2 H. G. Fischer, private communication

On the other hand, relying only on contemporary data, notably NA49, means having the invariant cross section at only one point in √ s , i.e. at one beam energy. Then, in order to obtain the general cross section, one has to supplement the data with some additional theoretical assumption, such as the scaling hypothesis [24], namely that the cross section only depends on p T and x R . While this behaviour is expected to be approximately respected, notably at high √ s , its quantitative accuracy can only be gauged by comparison with experimental data. For this reason we decided to apply both strategies and to use either fits or interpolations, to all datasets or to NA49 only, with or without the scaling hypothesis, to assess the importance of these effects.

## C. Method

Our fits were performed with the MINUIT minimization package. Let us denote with k = 1 , . . . , L the different experimental datasets, with i k the i -th point of the dataset k and let C be the vector of the cross section parameters. The fitting procedure consists of varying the values of the cross section parameters C , comparing the theoretical cross section F s ( i k , x i k , p i k ; C ) with the data f i k ( s i k , x i k , p i k ) and finally finding the minimum of χ 2 ( C ) function defined below. This procedure gives the best-fit configuration C best with the corresponding 1 σ errors σ C . We define the χ 2 ( C ) function to be minimised in the fitting procedure in the following way:

$$\chi ^ { 2 } ( \mathbf C ) = \chi ^ { 2 } _ { \text{stat} } ( \mathbf C ) + \chi ^ { 2 } _ { \text{sys} }$$

where

$$\chi _ { s y s } ^ { 2 } = \sum _ { k = 1 } ^ { L } \frac { ( \omega _ { k } - 1 ) ^ { 2 } } { \epsilon _ { k } ^ { 2 } } \,,$$

and

$$\chi _ { \text{stat} } ^ { 2 } ( \mathbf C ) = \sum _ { k = 1 } ^ { L } \sum _ { i _ { k } } \frac { ( \omega _ { k } f _ { i _ { k } } - F ( s _ { i _ { k } }, x _ { R _ { i _ { k } } }, p _ { T _ { i _ { k } } }, \mathbf C ) ) ^ { 2 } } { \omega _ { k } ^ { 2 } \sigma _ { i _ { k } } ^ { 2 } } \,.$$

5

FIG. 1. The data on d σ 3 p p → ¯ p /dp 3 employed in our analysis are reported as a function of E LAB ¯ p (left panel) and in the p T -x R space (right panel). For the data details, see Table I.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-1dfe0be9-6dd9-445c-8f24-7ab3a1a58d08.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-b1e8af41-d585-48c4-bb2f-b54d0aa51a80.jpg)

In the equations above, glyph[epsilon1] k is a systematic overall scale error of the dataset k (either quoted in the experimental paper, or assumed conservatively to be of the same order of the statistical one if this information is not available, notably for older references); the parameter ω k renormalises the dataset k and is determined consistently by the global fit: of course, large renormalizations with respect to glyph[epsilon1] k are disfavoured by the large penalties to be paid in the global analysis; σ i k is the statistical error on the data point i k , while the factor ( ω f k i k -F /ω ) k is the difference between experimental values f (accounting for a possible renormalization ω k , unique for each dataset) and the fitting function F , which depends on the independent variables described above, and on the vector of fitting parameters C .

This method is the most natural generalization of the unbiased one presented in [25] (see therein equation (3) and section 4) and it has already been successfully used in the past for other astroparticle physics analyses involving combinations of different datasets, as for instance nuclear reaction rates in primordial nucleosynthesis [26].

Passing on to the data interpolations against which we will be comparing our fitting procedure results, one difficulty lies with the fact that 3-dimensional interpolation of scattered data is a non-trivial problem in contemporary numerical analysis, with very few (if any) relevant publicly available tools. In order to tackle this issue, in our interpolations we will be making the assumption that the invariant distribution (3) scales with √ s only through an overall multiplicative dependence on the inelastic pp scattering cross section σ in . Under this assumption, by dividing the experimental data with σ in we obtain a √ s -invariant set of points for which only a 2-dimensional interpolation is needed. Besides, as a by-product of our analysis we have re-evaluated the inelastic cross section as described in Appendix B.

The interpolations were performed by means of the Python routine SmoothBivariateSpline contained in the scipy library, choosing piecewise cubic polynomials as interpolating functions. Note that the routine does not actually perform an exact interpolation, rather finds a compromise between the smoothness of the interpolating function and the closeness to the experimental data.

Estimating a statistically meaningful uncertainty band in this approach is fairly tricky. What we did was to consider each experimental best determination and error as the mean and standard deviation of a gaussian probability distribution of the cross section at that point. We then sampled these distributions accordingly, thus creating a large number of pseudo-experimental points. Each set of points is then interpolated (and, depending on the variable one is interested in, eventually integrated over cos ϑ and the proton incident energy/proton flux) to obtain an 'envelope band' for the quantity of interest. The average between maximum and minimum of the envelope at each point then defines an 'average interpolation curve'.

## III. RESULTS

## A. Validation of fitting method

As a preliminary exercise, and in order to validate our kinematical data conversion and fitting routines, we checked if the fit of Eq. (6) in [11] is reproduced, of course restricting ourselves to the datasets available at the time [18-23]. The parameterization of the invariant cross sec-

tion is in this case:

$$E \frac { d ^ { 3 } \sigma } { d p ^ { 3 } } & = \sigma _ { \text{in} } ( s ) ( 1 - x _ { R } ) ^ { D _ { 1 } } e ^ { - D _ { 2 } x _ { R } } \quad \text{ or } \\ & \quad \left [ D _ { 3 } ( \sqrt { s } ) ^ { D _ { 4 } } e ^ { - D _ { 5 } p _ { T } } + D _ { 6 } e ^ { - D _ { 7 } p _ { T } ^ { 2 } } \right ] \,, \quad ( 9 )$$

where σ in is the total inelastic cross section for pp collisions for which here, and only here, we used the parameterization adopted in [11]

$$\sigma _ { \text{in} } ( s ) = \sigma _ { 0 } \left [ 1 - 0. 6 2 \, e ^ { - \frac { E _ { \text{inc} } ( s ) } { 2 0 0 } } \sin \left ( \frac { 1 0. 9 } { E _ { \text{inc} } ^ { 0. 2 8 } ( s ) } \right ) \right ] \,, \ \ ( 1 0 ) \quad \text{less} ),$$

where E inc ( s ) is the incident kinetic energy in MeV defined as E inc ( s ) = s/ (2 m p ) -2 m p and σ 0 = 44 40 mbarn. . We show in Tab. II our best fit values and 1 σ errors for the cross section pararemeters D i in Eq.(9). A modest disagreement with the results of [11] was initially found only for D 1 and D 2 , which eventually we could attribute to a typo in the fitting parameter values reported in their Table V. If we invert D 1 with D 2 , not only do we obtain a very good agreement with our results, but also a reduced chi-squared χ 2 ν glyph[similarequal] 3 6, the same value the authors quote . in their paper. By insisting in interpreting literally the values of their table V, we would find χ 2 ν glyph[similarequal] 9 9, . clearly inconsistent with the value of 3.6. In Fig. 2 we display the comparison of the best fit and 3 σ uncertainty band of the source term derived with our best fit values of the parameters in Tab. II, and the best fit source term with the results reported in [11] with D 1 and D 2 inverted. The two results are essentially indistinguishable.

FIG. 2. Comparison of the source term for antiproton production in pp collisions as obtained in [11] (see however text for a correction in their table) and in this work, by re-fitting the same datasets with the same functional form, with our nominal 3 σ statistical error band. Vertical dot-dashed lines show the domain of energy actually covered by the experiments analysed.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-33615ecf-16b5-4cdc-92e2-44b0475464d0.jpg)

## B. Analysis of NA49 data

Once our routines validated, we proceed first with fits to the NA49 dataset only. We use the functional form

$$E \frac { d ^ { 3 } \sigma } { d p ^ { 3 } } & = \ \sigma _ { \text{in} } ( s ) ( 1 - x _ { R } ) ^ { C _ { 1 } } e ^ { - C _ { 2 } x _ { R } } \\ & \quad \left [ C _ { 3 } \, e ^ { - C _ { 4 } p _ { T } } + C _ { 5 } \, e ^ { - C _ { 6 } p _ { T } ^ { 2 } } \right ] \,, \quad \ \ ( 1 1 )$$

where σ in ( s ) is defined in Appendix B, Eq.(B2). This functional form is a simplified version of the standard parametrization proposed in [27] (it has four parameters less), which we found to provide an accurate and more compact description of the data. Note that we implicitly assume some form of scaling, in that the only dependence on s is given by the overall multiplication with the inelastic cross section. The best-fit values and the 1 σ errors are reported in Table III, with the corresponding fit having a reduced chi-square χ 2 ν = 1 3 for 137 degrees of freedom. . The comparison between data and fitted function (along with the corresponding 3 σ bands) is presented in Fig. 3. We see that the data are well represented by the fitting function, Eq.(11), for all the p T and x R values.

Next, we checked that the chosen fitting formula does not impose too strong a theoretical bias. To that purpose, as described in paragraph II C, we performed an 'educated' interpolation of the data by dividing the datapoints by σ in ( s NA49 ) and assuming that the resulting function is independent of s . The final function which is obtained by re-multiplying by σ in ( s ) thus still has a dependence on s , albeit a trivial one, via the overall factor σ in ( s ). The comparison between our fitting and interpolating procedures is shown in yellow in Fig. 4. The vertical lines correspond to the equivalent antiproton energy spanned by the NA49 experiment, where an interpolation is meaningful. In order to obtain a reasonable interpolation outside this interval, we supplemented the datasets with 'fake' points at the boundaries of the interpolation grid, with very large errors not to artificially influence the curve, yet sufficient to prevent the numerical routine from being driven to extreme functional form interpolations (for example, negative cross sections). No reasonable error can be however assigned outside the region covered by the data, apart for a lower limit that should be at least as large as the maximum relative width of the yellow band. The fact that the average interpolation curve is always within ∼ 3 σ of the best-fit previously obtained suggests that this 3 σ band is roughly representing the maximum uncertainty (at least where data exist), accounting not only for statistical errors, but also for possible theoretical biases, acting as additional systematics, related to the choice of the fitting function.

## C. Global analysis

Finally, we proceed to the global analysis of all available data on pp → ¯+ p X listed in Table I. In this case,

TABLE II. Best-fit values and 1 σ errors for the parameters D i in Eq.(9) resulting from a fit to the [11] dataset.

| D 1 (error)    | D 2 (error)          | D 3 (error)    | D 4 (error)      | D 5 (error)      | D 6 (error)        | D 7 (error)      |
|----------------|----------------------|----------------|------------------|------------------|--------------------|------------------|
| 4 . 22(0 . 66) | 3 . 435(0 . 016) 0 . | 0067(0 . 0014) | 0 . 510(0 . 050) | 3 . 609(0 . 021) | 0 . 0209(0 . 0010) | 3 . 086(0 . 083) |

TABLE III. Best-fit parameters and 1 σ errors to the NA49 data [14] with Eq.(11) .

| C 1 (error)    | C 2 (error)        | C 3 (error)      | C 4 (error)    | C 5 (error)        | C 6 (error)      |
|----------------|--------------------|------------------|----------------|--------------------|------------------|
| 7 . 56(1 . 15) | 0 . 245(0 . 148) 0 | . 0164(0 . 0025) | 2 . 37(0 . 13) | 0 . 0352(0 . 0020) | 2 . 902(0 . 059) |

given that we wish to describe data referring to quite different √ s values and covering different ( p T , x R ) regions, it is expected that we will have to introduce some complication with respect to the previous paragraph. In this spirit, we tried numerous different functional forms, essentially variations of the standard parametrization proposed in [27]. We present here results on our two most successful attempts, which also provide interesting insights on the extrapolation to regions where data are either scarce or altogether unavailable, a point that we shall discuss in more detail in section V A.

As a first step, we used an improved version of Eq.(11) introducing an explicit dependence on s , namely

$$E \frac { d ^ { 3 } \sigma } { d p ^ { 3 } } & = \sigma _ { \text{in} } ( s ) ( 1 - x _ { R } ) ^ { C _ { 1 } } e ^ { - C _ { 2 } x _ { R } } & ( 1 2 ) \quad \text{$non o$} \\ & \quad \left [ C _ { 3 } ( \sqrt { s } ) ^ { C _ { 4 } } e ^ { - C _ { 5 } p _ { T } } + C _ { 6 } ( \sqrt { s } ) ^ { C _ { 7 } } e ^ { - C _ { 8 } p _ { T } ^ { 2 } } \right ]. \quad \text{$omit$} \quad \text{$in $th$}$$

This parametrization of the cross section is similar to the one proposed in [27] except for the absence of a √ s exponent in the (1 -x R ) C 1 term. The fit gives a reduced chi-square of χ 2 ν = 4 16, . with a number of degrees of freedom of 385. The best-fit values and uncertainties are reported in Table IV. We have also checked that considering the exact form as in [27] we obtain an even worse fit to the data ( χ 2 ν = 5 6 with 385 degrees of freedom). .

Motivated by the relatively poor quality of the fit, we tried an extended version of Eq.(12), namely

The improved value of χ 2 ν is obtained at the expense of some rescaling of the datasets. With respect to our best result given by Eq.(13), the measurements reported in [13, 14, 18-23] are renormalised respectively by factors ω k of { 0 87 . , 1 04 1 16 0 98 . , . , . , 0 95 . , 1 13 . , 1 02 1 16 . , . } . Therefore the NA49 data [14], which represent the bulk of the fitting procedure, are renormalised by a negligible ∼ 2% while [13, 18, 19, 23] by more than 10%. Interestingly, the largest renormalization value is 16% for the BRAHMS dataset [13] giving a factor of 1.16, still not particularly significant given the statistical errors, yet perhaps indicative of some 'theoretical error' effects which become more prominent when an agreement over a large energy range is demanded. We display in Fig. 5 the comparison of the cross section best fit and 3 σ uncertainty band according to Eq.(13) with the datasets [14, 18-23]. We omit the comparison with the BRAHMS results, because in this case the cross section has only been measured along a line in the ( p T , x R ) space (see Fig. 1). Nevertheless, the difference between our best fit cross section and the data in [13] is at most ∼ 30%. We see that most of the data are well reproduced by the fitting function of Eq.(13) at all p T values. This is true in particular for the NA49 data, except for a slight overestimation at the lowest p T value. We have however checked that a 20% shift in the differential cross section for p T &lt; 0 15 has a . negligible effect on the antiproton source term (less than 5%).

$$E \frac { d ^ { 3 } \sigma } { d p ^ { 3 } } = \sigma _ { \text{in} } ( s ) ( 1 - x _ { R } ) ^ { C _ { 1 } } e ^ { - C _ { 2 } x _ { R } } \left | [ C _ { 3 } ( \sqrt { s } ) ^ { C _ { 4 } } e ^ { - C _ { 5 } p _ { T } } + \begin{array} { c } \text{We} \\ \text{In this} \end{array} \\ C _ { 6 } ( \sqrt { s } ) ^ { C _ { 7 } } e ^ { - C _ { 8 } p _ { T } ^ { 2 } } + C _ { 9 } ( \sqrt { s } ) ^ { C _ { 1 0 } } e ^ { - C _ { 1 1 } p _ { T } ^ { 3 } } ] \right |, \quad ( 1 3 ) \quad \text{there} \\ \text{points} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$}$$

where the absolute value simply prevents the function from becoming negative in some corners of parameter space. Compared to the previous function, this one further contains the exponential of a cubic function of p T and an additional s -dependence. The best-fit parameters for Eq.(13) are reported in Table V. This parametrization yields a somewhat better χ 2 ν = 3 30 for 382 degrees . of freedom. In order to test the validity of the scaling hypothesis, we obtain that the fit to the same data with Eq.13 containing no dependence on s but not in σ in ( s ), gives a reduced χ 2 ν = 4 82. .

We then repeated the interpolation analysis, previously only performed for NA49, for the entire dataset. In this case, the parameter space coverage is such that there is no need to supplement the dataset with 'fake' points, as previously done for the NA49 data alone. The spline method results in this case are, thus, fully datadriven, modulo our implicit assumption concerning the cross section √ s scaling according to an overall factor σ in ( s ).

In Fig. 6 we compare the results obtained for the antiproton source term through our fits according to the equations (13) and (12), as well as the estimate based on our spline interpolation. The energy range where pp data (except for BRAHMS) are available is bracketed by the vertical lines. We see that above 10 GeV, and within the region where experimental data are available, all three

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-272560ee-b52e-4dd2-8d3a-291871d544a3.jpg)

FIG. 3. Comparison between NA49 data with the fitting function of Eq.(11), see Table III, with 3 σ error bands. For clarity, the data and the theoretical curves at each p T value have been multiplied by a factor of 0 9 . n p T , where n p T is the integer counting the p T , from lower to higher ( i.e. : for p T = 0.60 GeV/ c 2 the rescaling is 0 9 . 6 ).

TABLE IV. Best-fit values and 1 σ errors for the parameters C i in Eq.(12) derived with a fit to all datasets.

| C 1 (error)   | C 2 (error)   | C 3 (error)      | C 4 (error)   | C 5 (error)   | C 6 (error)    | C 7 (error)   | C 8 (error)   |
|---------------|---------------|------------------|---------------|---------------|----------------|---------------|---------------|
| 4.499(0.040)  | 3.41(0.11)    | 0.00942(0.00083) | 0.445(0.027)  | 3.502(0.018)  | 0.0622(0.0086) | -0.247(0.049) | 2.576(0.027)  |

TABLE V. Best-fit values and 1 σ errors for the parameters C i in Eq.(13) derived with a fit to all datasets.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-41d278d5-aa8d-4d32-a8b2-8b54124fd782.jpg)

| C 1 (error)   | C 2 (error)   | C 3 (error)               | C 4 (error) C   | 5 (error) C 6 (error)     |    |
|---------------|---------------|---------------------------|-----------------|---------------------------|----|
| 4.448(0.035)  | 3.735(0.094)  | 0.00502(0.00036)          | 0.708(0.019)    | 3.527(0.014) 0.236(0.024) |    |
| C 7 (error)   | C 8 (error)   | C 9 (error)               | C 10 (error)    | C 11 (error)              |    |
| -0.729(0.036) | 2.517(0.027)  | - 1 . 822(0 . 009)10 - 11 | 3.527(0.022)    | 0.384(0.021)              |    |

methods yield compatible results. At lower and higher energies, however, there is a significant departure of the three estimates. We will discuss the implications of these results in much more detail in section V A.

In order to compare the results derived in this section to previous published proton-proton cross section estimates, we show in Fig. 7 the best fit and 3 σ uncertainty band source term calculated with our results and the best fit source term derived with the parametrizations adopted in [11, 12]. In the range of antiproton kinetic energy where data exist, our 3 σ band is marginally compatible with the parameterization in [11, 12], which is overestimated (underestimated) below (above) about 20 GeV.

fects of nuclear projectiles and targets in the collisions, and the yield coming from antineutron production. An exhaustive treatment of both subjects goes beyond our current purposes. For completeness, however, in the following we summarise the re-scalings of the yield from the pp → ¯ process that are usually adopted in the literature p to account for both the processes, and some of the issues involved.

## IV. CONTRIBUTIONS FROM NEUTRONS AND NUCLEI

In order to obtain the total antiproton source term, two more effects should be taken into account: the ef-

Concerning nuclear enhancements (effects of protonnucleus, nucleus-proton, and nucleus-nucleus collisions), unfortunately very little data are present, notably none for the most important channels which are the ones involving helium. One possible strategy to deduce cross sections for reactions involving helium is to constrain those of nuclear species for which some data are available, and extrapolate from heavier species to lighter ones, see e.g. [11]. Given that helium is quite light, however, it has often been considered reliable to deduce the relevant cross sections from rescaling the pp ones either with semiempirical formulae or via hadronic models, see e.g. [9]. The most recent dedicated studies were performed on

FIG. 4. Comparison between the function Eq.(11) fitted to NA49 data, with 3 σ error band (solid curve with cyan/blue shaded band), and interpolated curve (dashed red), with the interpolation envelope band, yellow/light shading. The dashed vertical lines correspond to the equivalent antiproton energy sampled by the NA49 experiment, where an interpolation is actually justified.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-af26c106-8115-4ec0-ae53-7ec4f6c30c1d.jpg)

the basis of the Monte Carlo (MC) model DTUNUC in [28] and in [9]. The models implemented in the software are based on the Dual Parton Model [29] and the Gribov-Glauber approach for a unified treatment of soft and hard scattering processes. The former are parameterised according to Regge phenomenology whereas the latter rely on lowest order perturbative QCD. Eventually, DTUNUC formed the basis of/merged into DPMJET-III (see [30] and refs. therein for further details) 3 .

Afit of the nuclei enhancement yield of antiprotons found in [28] compared to the one in pp collision, is given in [10]:

$$Q _ { t o t } / Q _ { p p } = 0. 1 2 ( T _ { \bar { p } } / G e V ) ^ { - 1. 6 7 } + 1. 7 8 \,, \quad ( 1 4 ) \quad \text{anti} \text{.}$$

with the above expression assuming 10% density ratio of He to H nuclei. Note that this ignores different spectral indices and species-dependent spectral breaks, which have been reported by some experiments but have not been confirmed by preliminary results of AMS-02. To gauge their possible effect at high-energies, we address the reader to the brief discussion in [32].

Lacking empirical information for the most relevant channels involving helium, it is hard to assess the accuracy of the previous models. The overall uncertainty (on the total source term yield Q tot ) was estimated in [28] to reach 40%, from the dispersion of predictions based on different prescriptions, but this conclusion is overly pessimistic, since not all the models/evaluations have the same reliability (some were based on obsolete prescriptions, for instance). In [9] the error estimate was closer

3 Recently, some new theoretical evaluation appeared in the Appendix of [31].

to 20 -25%, provided that the pp cross section does not depart from the Tan &amp; Ng parameterization [12] by more than 10%, which seems to be only marginally compatible with our results.

Since all these prescriptions do not include subtle nuclear effects, it is also likely that the uncertainty at low energies (where they are expected to play a larger role) is significantly larger than at high energy. If, in addition, one considers the more complicated astrophysical propagation effects at low-energies (reacceleation, convection, solar modulation) 4 and the need to correctly account for (catastrophic and non-catastrophic) energy losses, themselves affected by errors, it is clear that below a few GeV's the lower the energy the less reliable is the theoretical prediction. Most likely, this energy window cannot be used (but very crudely) for astroparticle physics constraints.

Another correction which is needed to infer the total antiproton flux from σ pp → ¯ p consists in accounting for the antiproton flux coming from antineutron production. Traditionally, it has been assumed that σ pp → ¯ n = σ pp → ¯ p , so that the previous results have been simply multiplied by a factor 2. However, the NA49 collaboration itself [34] has reported an isospin-dependence in studies of secondary yields in np and pp collisions: in pp reactions, there is a significant preference of the positively charged pn ¯ combination over ¯ pn (the opposite being true for neutron projectiles). This results in σ pp → ¯ n = κσ pp → ¯ p with κ glyph[similarequal] 1 5 . around x F ∼ 0 (see also Fig. 3 in [34]; x F is defined in Appendix A), although the effect depends on x F to some extent. Given the still rudimentary knowledge of these effects, an energy independent rescaling κ glyph[similarequal] 1 3 . ± 0 2 should encompass the data and be a . better approximation than the usually assumed κ = 1. It is clear that addressing these issues is of paramount importance for further reducing the uncertainties in the antiproton source term.

## V. DISCUSSION AND CONCLUSIONS

## A. Discussion

We now discuss our findings focusing on the global analysis outlined in Sect. III C. As we can see in Fig. 7, for antiproton energies lying roughly within the interval (10 300) GeV, we find that our results on the antipro-, ton source term from proton-proton scattering are consistent with previous estimates. They are moreover stable with respect to reasonable changes in the parametrization choice and in agreement with data-driven methods. These findings can be understood considering that the

4 It has been clearly shown as different propagation setups can be responsible for a ∼ 50% min-max dispersion in the flux expected at low-energy [33 ? ].

FIG. 5. Differental cross section for antiproton production in pp scattering, as a function of x R , for different p T values. The curves refer to the 3 σ uncertainty band around the best fit obtained with a fit of Eq.(13) to the datasets in Tab. I. The data, from top left to bottom right, are from [14, 18-23]. For the sake of clarity, the data from [19] and [14] and the relevant theoretical curves at each p T value have been rescaled by a factor 0 6 . n p T and 0 9 . n p T , respectively, as described in Fig.3.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-53f2890f-003e-4277-b268-b1106792b946.jpg)

FIG. 6. Comparison between fitted function of Eq.(13) with 3 σ band (solid curve with cyan/blue shaded band), of Eq.(12) with 3 σ band (dot-dashed curve with green/light green shaded band) and interpolated curve (dashed red), with the interpolation envelope band, red/orange shading. The dashed vertical lines correspond to the equivalent antiproton energy sampled by the global dataset, where an interpolation is in principle meaningful.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-1ada9c8a-53da-4e66-a70f-2d665523890d.jpg)

FIG. 7. The best fit and 3 σ uncertainty band source term derived with the fit of Eq.(12) and Eq.(13) to all datasets in Tab. I is shown together with the source term obtained using [11, 12] cross section parametrizations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-545960a8-1dff-423e-84e7-3bfc8f3d1b87.jpg)

majority of the data lie in the T ∈ (10 300) GeV range, , where the most reliable estimates of the distribution in Eq.(3) can be obtained and which, even prior to the NA49 and BRAHMS measurements, were already discreetly populated with data. In this sense, given that the NA49 data are not in contradiction with previous experimental results, it is expected (and verified) that the estimates presented in Tan &amp; Ng [12] and Duperray et al [11] are in good agreement with our findings for this energy range. Moreover, as long as a reasonable functional form is adopted for the invariant distribution, it is more or less bound to predict a comparable source term within this energy range. The small discrepancies of our spline interpolation and fitting approaches could be likely attributed to the fact that the interpolation essentially neglects the scaling violation , while the fits do allow for some flexibility (extra dependence on s ) to accommodate it.

On the other hand, at low and high energies, the relatively small amount of available data essentially implies extrapolations of the fits performed principally for T between 8 and 300 GeV. Consequently, moderately different assumptions can yield significantly different results. This is demonstrated by the fact that adopting two slightly different parameterizations while using the same dataset changes the high-energy source term prediction quite dramatically. Moreover, these findings are insensitive to the inclusion or not of the BRAHMS data in the analysis, which means that the results in [13] are not sufficient to constrain the high-energy behavior of the invariant distribution and, hence, the antiproton source function. This is due to the fact that the data of [13], only cover the exponentially suppressed highp T region (similarly to the ones of [23]), see Fig. 1. In this sense, both the lowenergy and high-energy behavior of the invariant distribution remain highly uncertain. Given that both the spline method and the fit with Eq.(13) demonstrate a similar trend at high energies, we believe that making any conclusive statement concerning the high-energy behaviour of the antiproton inclusive cross section would be risky. This is all the more the case since spline interpolations can be notoriously misleading when extrapolated outside data-covered regions.

Whereas in the low-energy regime this point is not very important, given that in any case the secondary antiproton flux is dominated by huge uncertainties coming from astrophysical sources (solar modulation, propagation parameters, antiproton scattering cross sections), it is plausible that in the region of several hundreds of GeV and higher the main uncertainty is still due to the antiproton production cross section.

We summarise in Table VI the pp -induced source term along with the associated percentage uncertainties resulting from our analysis of the NA49 data according to Eq.(11), our global analysis according to Eqs.(12) and (13), our spline interpolation method of the full dataset, and the previous estimates in [12] and [11], for a few representative values of the antiproton energy. This table simply illustrates the results reported in Figs. 6 and 7: with increasing energy, the different approaches turn from marginally compatible (at the lowest energies, few GeV) to fully compatible until, towards the end of the region for which experimental data are available, they yield very different results.

Concerning the error estimates, we also point out that the nominal 1 σ error band obtained from the χ 2 minimization procedure is underestimated, for several reasons. In some case where χ 2 /dof is close to 1, as in the fit to NA49 data only with a simple fitting formula, we showed how the agreement with an interpolation method is only meaningful if roughly a 3 σ band is used as typical estimate of the error. This is the choice we presented in

TABLE VI. Best-fit values and corresponding percentage relative errors for the pp -induced source term (in GeV -1 cm -3 s -1 ), for some representative antiproton energies and different approaches in the data analysis.

|   T (GeV) | Eq.(11) (% error)        | Eq.(12) (% error)        | Eq.(13) (% error)        | spline (% error)             | Tan & Ng         | Duperray et al   |
|-----------|--------------------------|--------------------------|--------------------------|------------------------------|------------------|------------------|
|         5 | 1 . 23 · 10 - 30 (4 . 9) | 1 . 47 · 10 - 30 (6 . 1) | 1 . 67 · 10 - 30 (5 . 4) | 1 . 38 · 10 - 30 (2 . 7)     | 1 . 42 · 10 - 30 | 1 . 40 · 10 - 30 |
|        10 | 4 . 31 · 10 - 31 (4 . 2) | 4 . 87 · 10 - 31 (3 . 0) | 5 . 17 · 10 - 31 (4 . 8) | 4 . 34 · 10 - 31 (2 . 5)     | 4 . 96 · 10 - 31 | 4 . 74 · 10 - 31 |
|       100 | 1 . 70 · 10 - 33 (5 . 9) | 1 . 82 · 10 - 33 (8 . 7) | 1 . 77 · 10 - 33 (6 . 8) | 2 . 03 · 10 - 33 (3 . 2)     | 1 . 82 · 10 - 33 | 2 . 04 · 10 - 33 |
|       500 | 2 . 42 · 10 - 35 (6 . 2) | 2 . 82 · 10 - 35 (9 . 5) | 3 . 39 · 10 - 35 (8 . 8) | 3 . 26 · 10 - 35 (5 . 2)     | 2 . 38 · 10 - 35 | 3 . 27 · 10 - 35 |
|      1000 | 3 . 13 · 10 - 36 (6 . 9) | 4 . 16 · 10 - 36 (11)    | 6 . 83 · 10 - 36 (10)    | 7 . 02 · 10 - 36 (5 . 8) 3 . | 29 · 10 - 36     | 4 . 93 · 10 - 36 |

our plots. A similar prescription was found to be more indicative of the real uncertainty, once global fits were performed. In this case, the inadequacy of the nominal 1 σ error band was already hinted to by the relatively large reduced χ 2 , never smaller than χ 2 ν = 3 30. . We attribute these results to a combination of factors: i) underestimated experimental errors, notably in (some of) the older datasets, due to effects that were neglected as the feed-down we mentioned. ii) inadequacy of any simple functional form tested to describe faithfully the data, especially on a large dynamic range; iii) some sort of more or less implicit analytical extrapolation assumption in order to achieve coverage of the 3-dimensional space ( √ s, p T , x R ) starting from a discrete set of points. Note that this also applies to interpolation techniques, which for instance rely on some theoretical assumptions such as scaling. The situation may be certainly improved if high-quality measurements such as the ones provided by NA49 could be extended to a broader dynamic range.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0288v3-figure-c454cf11-9de2-4082-a239-e7bcc269e32b.jpg)

We also stress that outside the regions where data are available, there is no compelling reason for either one of our results according to equations (12) and (13) to be more realistic than the other. Whereas the agreement of all of our computations at intermediate energies hints that the error estimates there is fairly reliable, this is not at all the case at very low and high energies. A more conservative approach is to assume that in this case the error is dominated by the extrapolation uncertainty, for which a proxy is given by the region spanned by the ensemble of our approaches, amounting to about 50% at 1 TeV.

As a practical summary of our analysis, we report in Fig. 8 an estimate for the uncertainties inherent to the production of antiprotons from inelastic pp scatterings. The results are expressed as the ratio of the antiproton source term in Eq.(1) to a reference value. For the blue and the red bands, this reference value has been fixed to the source term obtained by setting the pp production cross section to the best fit to all the data obtained with Eq. (13) (parameters as in Table V). Outside the vertical bands-delimiting the energy range in which data are available-we extrapolate the production cross section by means of the same formula.

The blue band corresponds to considering parametriza-

FIG. 8. Estimate of the uncertainties in the antiproton source term from inelastic pp scattering. The blue band indicates the 3 σ uncertainty band due to the global fit with Eq.(13), while the red band corresponds to the convolution of the uncertainties brought by fits to the data with Eq.(13), Eq.(12) and with the spline interpolation (see Fig.6.). The orange band takes into account the contribution from decays of antineutrons produced in the same reactions. Vertical bands as in Fig.6. See text for details.

tion (13) alone. By simple inspection we can clearly see that the relevant uncertainty is maximally of the order of 10%. The red band is obtained by convoluting the uncertainty bands resulting from fits through Eqs.(13) and (12) and (within the vertical bands) the spline interpolation. This more conservative approach sizes the uncertainties from 20% at the lowest energies to the extrapolated 50% at 1 TeV.

The most conservative estimate is shown by the orange band, where the additional uncertainty on the antineutron production has been taken into account. In this case, the normalization has been fixed to a source term in which the antineutrons produced in pp scatterings contribute with an energy-independent rescaling factor κ = 1 3 (w.r.t. . 1). The relevant uncertainty band has been derived by shifting the (red) previous convolution by an additional factor to account for the antineutron decay, κ glyph[similarequal] 1 3 . ± 0 2, . as discussed in Sect. IV. The orange band indicates that the antiproton source term may

vary by 30% at 1 GeV, and by up to more than 50% at 1 TeV. In the energy range where scattering data are available (between 4 GeV and 550 GeV), the uncertainty is of O (20 -30%) non-symmetrically around the reference source term, with the ambiguity increasing with energy.

Finally, although we did not include an analysis of the uncertainties due to the contribution of nuclei to the yield, it is obvious that the total relative uncertainty cannot be smaller than the one showed above.

## B. Conclusions

In this work we have performed a re-evaluation of the Lorentz-invariant distribution for inclusive antiproton production in pp collisions in light of the recent results from the NA49 and BRAHMS experiments. We have combined these observations with older measurements at different centre-of-mass energies in order to extract a reliable estimate both for the average value of the cross section and for the current status of the corresponding uncertainties. Our main results have been presented in sections III B, III C and IV and summed up as handy fitting functions.

We have paid special attention in quantifying the extent to which the functional form we adopt introduces theoretical bias, by comparing our primary results to those obtained by using different parametrizations, datadriven estimates or different subsets of the data. We are therefore confident that the uncertainties quoted throughout our paper are robust (cf also the discussion in section V A).

Our findings show that, despite the experimental progress in measuring the inclusive antiproton production cross section, uncertainties persist. In most of the well-constrained intermediate antiproton energy range (few GeV - few hundreds GeV), these uncertainties can be as low as 10 -15%. At higher energies, we have shown that our knowledge is much worse, with extrapolations leading to errors larger than ∼ 50% at 1 TeV.

A complete determination of the antiproton yield from pp scattering must include the antineutron decay contribution, usually assumed to be identical to the antiproton one. A recent measurement has reported a significant isospin dependence, which could amount to a slight increase of the antiproton source term and be accompanied by a non negligible uncertainty in the overall flux.

Last but not least, we remind that a significant contribution to the cosmic antiproton flux is due to the reactions involving helium nuclei, both in the incoming radiation and in the ISM. The relevant cross sections have never been measured, and the corresponding antiproton source term can be estimated only by relying on some theoretical models and extrapolations, introducing an additional uncertainty which is hard to quantify.

We expect that for the few years to come, these uncertainties will continue imposing non-negligible limitations on the interpretation of cosmic antiproton data, expected to be measured by AMS-02 with an unprecedented accuracy. As a side-effect, it appears unlikely that any definite conclusion for dark matter indirect detection could be drawn from a relatively featureless 'excess' in the antiproton yield, expected to be be below a few tens of percent of the overall flux, unless perhaps correlated features are found in several different cosmic ray channels. At very low energies, the situation is further aggravated by the poor knowledge of astrophysical parameters, plus additional nuclear and particle physics uncertainties related e.g. to non-catastrophic energy losses of antiprotons.

Despite current limitations, many of these sources of error are not irreducible, but could be addressed with dedicated experimental campaigns. As we stressed, the key advances would come from planning experiments with an extended coverage in the ( p T , x R ) plane, like NA49, with higher energies, with helium target/projectile, and by providing systematic analyses of anti-neutron yields. Another very useful enterprise would be a re-analysis of existing datasets with the aim of assessing the feed-down corrections associated to different experiments. Indeed, the CR antiprotons must include the contributions from hyperon decays, while the current databases are not uniform in that respect.

Cosmic ray antiprotons offer an important tool to address a number of astrophysical and astroparticle questions, as certified by the resources invested in balloon-borne and space-borne detectors in recent years. We believe that a nuclear and particle physics commitment to measure many of these missing ingredients should be identified as a strategic task, to provide astroparticle missions with the crucial laboratory input to fully exploit their results.

Note Added: During the completion of this work, we became aware of a study by R. Kappl and M. W. Winkler, arXiv:1408.0299. On one side, it focuses on NA49 data interpolation rather than performing a global dataset reanalysis with fitting formulae, like ours. On the other hand, it extends to propagated fluxes including other channels as well. As far as a comparison is possible, our results agree with theirs within quoted errors.

## ACKNOWLEDGMENTS

We would like to thank H. G.Fischer for valuable information about NA49 experiment, D. Maurin for stimulating discussions and for careful reading of the manuscript, and R.H. Bernstein for spotting a typo in Eq. (10) and another one in Table (II) in an earlier version of this manuscript. A.G. would like to acknowledge useful discussions with P. Aurenche on issues in hadron physics. M.D.M. would like to thank Andrea Vittino for interesting and useful discussions on antiprotons secondary production uncertainties. At LAPTh, this activity was developed coherently with the research axes supported

by the ANR Labex grant ENIGMASS. F.D. further acknowledges support from AAP-'LeSAdHE' of University of Savoie.

## Appendix A: Useful kinematics

In order to extract the invariant distribution of Eq.(3), we have transformed all data (when needed) in the CM frame, expressed them in the form ( p T , x R , f ( √ s, p T , x R )) and then performed our fit and interpolation as described in the main text. Here we collect for convenience some useful formulae by means of which these conversions can be done. In what follows, all kinematic variables carrying stars are defined in the CM frame whereas those carrying LAB superscripts are defined in the lab frame.

Note that experiments [18, 19, 23] and [14] are fixed-target experiments, whereas [20-22] and [13] are colliding beam ones. With the exception of NA49 [14], fixed-target experiments give their results in the LAB (or, equivalently in this case, target) frame. All other datasets are given in the CM frame.

The experiments [18] and [19] do not provide results for the quantity defined in (3), rather for d σ/d 2 Ω dp which can be recast as [35]

$$\frac { d ^ { 2 } \sigma } { d \Omega d p } = p _ { \tilde { p } } ^ { L A B, 2 } \frac { d ^ { 3 } \sigma } { d p ^ { 3 } }. \quad \quad ( \text{A1} ) \quad \text{all}$$

This leads to

$$E _ { \tilde { p } } ^ { L A B } \frac { d ^ { 3 } \sigma } { d p _ { \tilde { p } } ^ { 3 } } = \frac { E _ { \tilde { p } } ^ { L A B } } { p _ { \tilde { p } } ^ { L A B, 2 } } \frac { d ^ { 2 } \sigma } { d \Omega d p } \quad \quad ( A 2 ) \quad \text{Th} \quad \text{res} \quad \text{ }$$

which can be computed straightforwardly from the available data using E LAB ¯ p = √ p LAB, 2 ¯ p + m 2 p .

The results are given as a function of the transverse momentum component p T and ϑ LAB , the ¯ emission anp gle with respect to the incoming proton beam in the lab frame. Since p T remains unchanged in the two reference frames, one needs to compute the x R values which the measurements correspond to. We have that

$$E _ { \tilde { p } } ^ { * } = \gamma ^ { C M } ( E _ { \tilde { p } } ^ { L A B } - v ^ { C M } p _ { \tilde { p } } ^ { L A B } \cos \vartheta ^ { L A B } ) \quad ( \mathbf A 3 ) \quad \text{From}$$

where v CM and γ CM are the velocity of the CM reference frame with respect to the lab one and the corresponding Lorentz boost. These are given by

$$v ^ { C M } & = p _ { p } ^ { L A B } / ( E _ { p } ^ { L A B } + m _ { p } ) \quad \text{(A4)} \quad \text{and} \\ & = p _ { p } ^ { L A B } / ( \sqrt { p _ { p } ^ { L A B ^ { 2 } } + m _ { p } ^ { 2 } } + m _ { p } ) \quad \text{$\text{It is}$} \\ \gamma ^ { C M } & = ( E _ { p } ^ { L A B } + m _ { p } ) / \sqrt { s } \,. \quad \text{$\text{For}$}$$

Then, direct use of the definition (4) allows us to compute x R .

Experiments [20, 21] and [14] provide results for f ( p + p → ¯+ p X ) in the CM frame, as a function of p T and the alternative Feynman scaling variable x F defined as

$$x _ { F } = \frac { 2 p _ { L } ^ { * } } { \sqrt { s } } \simeq \frac { p _ { L } ^ { * } } { p _ { L, m a x } ^ { * } } \,, \quad \quad ( A 5 )$$

where p ∗ L is the antiproton longitudinal momentum and p ∗ L,max is its maximum possible longitudinal momentum. It is thus necessary to find the correspondence between the x R and x F scaling variables. First, from p ∗ 2 L = p ∗ 2 -p ∗ 2 T and p 2 = E 2 -m 2 we get

$$E _ { \bar { p } } ^ { * 2 } = p _ { L } ^ { * 2 } + m _ { p } ^ { 2 } + p _ { T } ^ { * 2 }. \text{ \quad \ \ } ( \text{A6} )$$

whereas E ∗ ¯max p is given by (5). Direct use of (A5) and (4) then gives

$$x _ { R } = \frac { \sqrt { x _ { F } ^ { 2 } ( s / 4 ) + m _ { p } ^ { 2 } + p _ { T } ^ { 2 } } } { E _ { \bar { p }, \max } ^ { * } } \,. \quad \quad ( A 7 )$$

Reference [23] provides results for f ( p + p → ¯ + p X ) in the CM frame as a function of p T and ϑ LAB , also providing the values ϑ ∗ of the angle in the CM frame. The p T values remain unchanged in the CM frame. So all that is left is to calculate x R . Since p ∗ ¯ p = p T / sin ϑ ∗ , by using E = √ glyph[vector] p 2 + m 2 p and equation (A13), we can calculate x R through (4).

The BRAHMS experiment [13] chooses to give its results for the invariant cross section

$$\frac { 1 } { 2 \pi p _ { T } } \frac { d ^ { 2 } \sigma } { d p _ { T } \, d y } = E _ { \tilde { p } } \frac { d ^ { 3 } \sigma } { d p _ { \tilde { p } } ^ { 3 } } \quad \quad ( A 8 )$$

as a function of the antiproton transverse momentum p T and the antiproton rapidity in the CM frame, y ∗ , defined as

$$y ^ { * } = \frac { 1 } { 2 } \ln \left ( \frac { E ^ { * } + p _ { L } ^ { * } } { E ^ { * } - p _ { L } ^ { * } } \right ). \quad \quad \left ( A 9 \right )$$

From this definition and E ∗ 2 ¯ p = p ∗ 2 L + p 2 T + m 2 ¯ p we get

$$E _ { \bar { p } } ^ { * } = \sqrt { m _ { \bar { p } } ^ { 2 } + p _ { T } ^ { 2 } } \cosh y. \text{ \quad \ \ } ( A 1 0 )$$

Then, x R can be calculated by means of this relation and (4).

It is maybe useful to report the computation of the maximal antiproton energy E ∗ ¯max p introduced in 5. For a general inclusive reaction a + b → c + X in the CM frame we have

$$\sqrt { s } = E _ { c } ^ { * } + E _ { X } ^ { * } \. \quad \quad \quad ( A 1 1 )$$

By replacing E ∗ c , E ∗ X through E = √ p 2 + m 2 , squaring the resulting relation, reintroducing the energies in the crossed terms, eliminating E X in favour of √ s and E ∗ c , using p ∗ 2 -E ∗ 2 = -m 2 to eliminate p ∗ 2 and solving for E ∗ c we get

$$E _ { c } ^ { * } = \frac { s + m _ { c } ^ { 2 } - m _ { X } ^ { 2 } } { 2 \sqrt { s } } \,, \quad \quad \ ( \text{A12} ) \quad \ \text{with} \quad \text{where}$$

where we introduced for compactness a slight abuse of notation for the X system by assigning it a mass variable m X . In reality m X simply refers to internal energy of the X system.

Now, for s and m p fixed, we see that E ∗ ¯ max p, is obtained through (A12) once the energy m X of the X system becomes minimal. Conservation of baryon number fixes m X, min = 3 m p (production of 3 additional protons at rest). So we get

$$E _ { \bar { p }, \max } ^ { * } = \frac { s - 8 m _ { p } ^ { 2 } } { 2 \sqrt { s } } \,. \quad \quad ( A 1 3 ) \quad \text{ } T A$$

Finally we derive the maximum antiproton angle with respect to the incoming proton in the laboratory frame ϑ LAB max . This quantity can be derived using the condition x R ≤ 1 and the definition of x R in Eq.(4). Then writing E ∗ ¯ p in the laboratory frame (see Eq.(A3)) the condition x R ≤ 1 can be written as

$$x _ { R } = \frac { \gamma ^ { C M } 2 \sqrt { s } ( E _ { \tilde { p } } ^ { L A B } - v ^ { C M } p _ { \tilde { p } } ^ { L A B } \cos \vartheta ^ { L A B } ) } { s - 8 m _ { p } ^ { 2 } } \leq 1 \,. \quad \begin{array} { c } \text{we are} \\ \text{data} \\ \text{in the} \\ \text{$\text{GeV.}$} \\ \text{avail} \end{array}$$

By using the definition of v CM and γ CM given in Eq.(A4) we get

$$\cos \vartheta _ { m a x } ^ { L A B } = \frac { E _ { \tilde { p } } ^ { L A B } ( E _ { p } ^ { L A B } + m _ { p } ) - m _ { p } ( E _ { p } ^ { L A B } - 3 m _ { p } ) } { \sqrt { E _ { \tilde { p } } ^ { L A B } - m _ { p } ^ { 2 } } \sqrt { E _ { p } ^ { L A B } - m _ { p } ^ { 2 } } } \,. \quad \text{ the n} \quad \text{to } \text{ca} ] \\ ( A 1 5 ) \quad \text{in } T a l \\ \text{if} \quad \text{$\#\cdot\cdot\cdot\cdot\cdot\cdot\cdot$} }$$

## Appendix B: Inelastic cross section parameterization

The inelastic proton cross section is defined as the difference between the total pp scattering cross section σ pp tot and its elastic counterpart σ pp el

$$\sigma _ { \text{in} } ^ { \text{pp} } = \sigma _ { \text{tot} } ^ { \text{pp} } - \sigma _ { \text{el} } ^ { \text{pp} }. \quad \quad \ ( B 1 ) \quad \text{ag} _ { \text{el} }$$

In order to estimate σ pp in for our energy region of interest, we employ the experimental data provided by the

[1] P. Blasi and P. D. Serpico, Phys.Rev.Lett. 103 , 081103 (2009), arXiv:0904.0871 [astro-ph.HE].

[2] P. Salati, F. Donato, and N. Fornengo, (2010), arXiv:1003.4124 [astro-ph.HE].

Particle Data Group (PDG) on the total and elastic pp cross sections [36]. We fit this data by means of the highest-ranking parametrization of the total proton cross section suggested by the PDG itself, which reads

$$\sigma _ { t o t } ^ { \text{pp} } & = Z ^ { \text{pp} } + B ^ { \text{pp} } \log ^ { 2 } ( s / s _ { M } ) + Y _ { 1 } ^ { \text{pp} } ( s _ { M } / s ) ^ { \eta _ { 1 } } - Y _ { 2 } ^ { \text{pp} } ( s _ { M } / s ) ^ { \eta _ { 2 } } \\ ( \text{A12} ) \quad. \quad \underset { \theta \theta } { \ e \theta } \quad \theta \theta \quad \theta \quad \theta \quad \dots \quad \theta \quad \dots \quad \theta \quad \dots$$

where B pp = π ( glyph[planckover2pi1] c ) 2 /M 2 , s M = (2 m p + M ) 2 , all energies are given in GeV and all cross sections in mb.

TABLE VII. Fit results for the total and elastic proton scattering cross sections according to Eq.(B2).

| Parameter   |   Total | Elastic        |
|-------------|---------|----------------|
| M           |   2.06  | 3.06           |
| Z pp        |  33.44  | 144.98         |
| Y pp 1      |  13.53  | 2.64           |
| Y pp 2      |   6.38  | 137.27         |
| η 1         |   0.324 | 1.57           |
| η 2         |   0.324 | -4.65 × 10 - 3 |

Although this parametrization is given for the total proton cross section, noticing the resemblance in the √ s dependence of both the total and the elastic cross section, we employ the same functional form in order to fit the data on both σ pp tot and σ pp el . In both cases, we include in the fitting procedure only data referring to √ s &gt; 5 GeV. Moreover, given the large amount of relevant data available for low √ s values as opposed to the relative scarcity of measurements for √ s &gt; 80 GeV, we include the highest-energy data contained in the sample in order to capture the cross section behaviour over our entire energy range of interest in a reliable manner.

The fitting parameters for both σ pp tot and σ pp el are shown in Table VII. The resulting functions provide very good fits to the available data, with a χ 2 per degree of freedom of 0 78 for the total cross section and 1 5 for the elastic . . case, where statistical and systematic uncertainties have been added in quadrature. With these two parametrizations at hand, we then compute the inelastic cross section through its definition (B1). We have verified that, modulo small differences due to the specific goals of our analysis described before, our results are in very good agreement not only with the PDG results on σ pp tot but also with the function for σ pp el quoted by all LHC experiments (see e.g. [37, 38]) as well as with the results on both σ pp tot and σ pp el by Block and Halzen [39].

[3] D. Maurin, R. Taillet, F. Donato, P. Salati, A. Barrau, et al., (2002), arXiv:astro-ph/0212111 [astro-ph].

[4] C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, JCAP 0810 , 018 (2008), arXiv:0807.4730 [astro-ph].

- [5] F. Donato, D. Maurin, P. Brun, T. Delahaye, and P. Salati, Phys.Rev.Lett. 102 , 071301 (2009), arXiv:0810.5292 [astro-ph].
- [6] O. Adriani et al. (PAMELA Collaboration), Phys.Rev.Lett. 105 , 121101 (2010), arXiv:1007.0821 [astro-ph.HE].
- [7] http://www.ams02.org/ .
- [8] I. Moskalenko, A. Strong, and O. Reimer, Astron.Astrophys. 338 , L75 (1998), arXiv:astroph/9808084 [astro-ph].
- [9] F. Donato, D. Maurin, P. Salati, A. Barrau, G. Boudoul, et al., Astrophys.J. 563 , 172 (2001), arXiv:astroph/0103150 [astro-ph].
- [10] I. V. Moskalenko, A. W. Strong, J. F. Ormes, and M. S. Potgieter, Astrophys.J. 565 , 280 (2002), arXiv:astroph/0106567 [astro-ph].
- [11] R. Duperray, C.-Y. Huang, K. Protasov, and M. Buenerd, Phys.Rev. D68 , 094017 (2003), arXiv:astroph/0305274 [astro-ph].
- [12] L. Tan and L. Ng, Phys.Rev. D26 , 1179 (1982).
- [13] I. Arsene et al. (BRAHMS Collaboration), Phys.Rev.Lett. 98 , 252001 (2007), arXiv:hep-ex/0701041 [hep-ex].
- [14] T. Anticic et al. (NA49 Collaboration), Eur.Phys.J. C65 , 9 (2010), arXiv:0904.2708 [hep-ex].
- [15] J.-P. Meyer, Apjs 57 , 173 (1985).
- [16] H. Haino and the Ams-02 Collaboration, Talk at the 33rd ICRC Conference (2013).
- [17] M. Di Mauro, F. Donato, N. Fornengo, R. Lineros, and A. Vittino, JCAP 1404 , 006 (2014), arXiv:1402.0321 [astro-ph.HE].
- [18] D. Dekkers, J. Geibel, R. Mermod, G. Weber, T. Willitts, et al., Phys.Rev. 137 , B962 (1965).
- [19] J. Allaby, F. Binon, A. Diddens, P. Duteil, A. Klovning, et al., (1970).
- [20] P. Capiluppi, G. Giacomelli, A. Rossi, G. Vannini, A. Bertin, et al., Nucl.Phys. B79 , 189 (1974).
- [21] K. Guettler et al. (British-Scandinavian-MIT Collaboration), Nucl.Phys. B116 , 77 (1976).
- [22] J. Johnson, R. Kammerud, T. Ohsugi, D. Ritchie, R. Shafer, et al., Phys.Rev.Lett. 39 , 1173 (1977).
- [23] D. Antreasyan, J. Cronin, H. J. Frisch, M. Shochet, L. Kluberg, et al., Phys.Rev. D19 , 764 (1979).
- [24] F. Taylor, D. C. Carey, J. Johnson, R. Kammerud, D. Ritchie, et al., Phys.Rev. D14 , 1217 (1976).
- [25] G. D'Agostini, Nucl.Instrum.Meth. A346 , 306 (1994).
- [26] P. D. Serpico, S. Esposito, F. Iocco, G. Mangano, G. Miele, et al., JCAP 0412 , 010 (2004), arXiv:astroph/0408076 [astro-ph].
- [27] A. N. Kalinovskii, M. V. Mokhov, and Y. Nikitin, American Institute of Physics (1989).
- [28] M. Simon, A. Molnar, and S. Roesler, Astrophys. J. 499 , 250 (1998).
- [29] A. Capella, U. Sukhatme, C.-I. Tan, and J. Tran Thanh Van, Phys.Rept. 236 , 225 (1994).
- [30] S. Roesler, R. Engel, and J. Ranft, , 1033 (2000), arXiv:hep-ph/0012252 [hep-ph].
- [31] M. Kachelriess and S. Ostapchenko, Phys.Rev. D86 , 043004 (2012), arXiv:1206.4705 [astro-ph.HE].
- [32] F. Donato and P. D. Serpico, Phys.Rev. D83 , 023014 (2011), arXiv:1010.5679 [astro-ph.HE].
- [33] M. Cirelli, D. Gaggero, G. Giesen, M. Taoso, and A. Urbano, (2014), arXiv:1407.2173 [hep-ph].
- [34] H. Fischer (NA49 Collaboration), Heavy Ion Phys. 17 , 369 (2003).
- [35] E. Byckling and K. Kajantie, (1971).
- [36] http://pdg.lbl.gov/2013/hadronic-xsections/ hadron.html .
- [37] G. Antchev et al. (TOTEM), Europhys.Lett. 101 , 21004 (2013).
- [38] G. Antchev et al. (TOTEM Collaboration), Phys.Rev.Lett. 111 , 012001 (2013).
- [39] M. M. Block and F. Halzen, Phys.Rev. D86 , 014006 (2012), arXiv:1205.5514 [hep-ph].
Mattia di Mauro Fiorenza Donato Andreas Goudelis Pasquale Dario Serpico
A new evaluation of the antiproton production cross section for cosmic ray studies
Theoretical predictions for the cosmic antiproton spectrum currently fall
short of the corresponding experimental level of accuracy. Among the main
sources of uncertainty are the antiproton production cross sections in cosmic
ray inelastic interactions. We analyse existing data on antiproton production
in $pp$ scattering, including for the first time the measurements performed by
the NA49 Collaboration. We compute the antiproton spectrum finding that in the
energy range where data are available (antiproton energies of about 4-550 GeV)
different approaches lead to almost equivalent results, with an uncertainty of
10-20\%. Extrapolations outside this region lead to different estimates, with
the uncertainties reaching the 50\% level around $1$ TeV, degrading the
diagnostic power of the antiproton channel at those energies. We also comment
on the uncertainties in the antiproton production source term coming from
nuclei heavier than protons and from neutrons produced in $pp$ scatterings, and
point out the need for dedicated experimental campaigns for all processes
involving antiproton production in collisions of light nuclei.