Proton–proton total cross-sections at VHE from accelerator data

Up-to-date estimates of proton–proton total cross-sections, σtotpp, at very high energies in the literature were obtained from cosmic rays (>1017 eV) by approximations using the measured proton–air cross-section at these energies. As σtotpp are measured with present day high energy colliders up to nearly 2 TeV in the centre of mass (∼1015 eV in the laboratory), several proven theoretical, empirical and semi-empirical parametrizations for interpolation at accelerator energies were used to extrapolate these measured values to get reasonable estimates of cross-sections at higher cosmic ray energies (∼1017 eV). The cross-section estimates from these two methods disagree by a discrepancy beyond statistical error. Here we use a phenomenological model based on the ‘multiple diffraction’ approach to successfully describe data at accelerator energies. Using this model, we then estimate σtotpp at cosmic ray energies. The model free-parameters used in the fit depend on only two physical observables: the differential cross-section and the parameter ρ. The model estimates of σtotpp are then compared with total cross-section data. Using regression analysis, we determine confidence error bands, analysing the sensitivity of our predictions to the data used in the extrapolations. This work reduces the width of the confidence band around ‘multiple diffraction’ model fits of accelerator data. With the data at 546 GeV and 1.8 TeV, our extrapolations are compatible with only the Akeno cosmic ray data, predicting a slower rise with energy than do other cosmic ray results and other extrapolation methods. We discuss our results within the context of constraints expected from future accelerator and cosmic ray experimental results.


Introduction
Recently a number of difficulties in uniting accelerator and cosmic ray values of the hadronic total cross-sections for proton-proton (p-p), σ pp tot , and antiproton-proton, σp p tot , in the light of the best recent data have been summarized [1]. A united picture would be of the highest importance for the interpretation of results of new cosmic ray experiments such as the HiRes [2] and in designing proposals that are currently in progress as the Auger Observatory [3], as well as in designing detectors for future accelerators, such as the CERN pp Large Hadron Collider (LHC) [4]. Although most accelerator measurements of σ pp tot and σp p tot at centre of mass energies √ s 1.8 TeV are fairly consistent (see figure 1 where the 1.8 values just embrace the 85% bands), unfortunately above √ s > 6 TeV cosmic ray experiments agree among themselves only because of their large uncertainty bands (section 3); 7.5-23% for three experiments at 30 TeV and (3-19)% for 2 experiments at 40 TeV. It is more disturbing that they are consistent with different predictions from the extrapolation of accelerator data up to cosmic ray energies, again only because of their own uncertainty and that of the extrapolations. Some extrapolations predict smaller values of σ pp tot than those of cosmic ray experiments (e.g. [5,6]); others agree at some specific energies with cosmic ray results (e.g. [7]). Dispersion of cosmic ray results is associated mainly with the strong model-dependence of the relationship of the basic hadron-hadron crosssection and the hadronic cross-section in air. The latter determines the attenuation length of hadrons in the atmosphere, which is usually measured in different ways, and depends strongly on the rate (k) of energy dissipation of the primary proton into the electromagnetic shower observed by the experiment. Such a cascade is simulated by different Monte Carlo techniques 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT implying additional discrepancies between different experiments. Furthermore, σ pp tot in cosmic ray experiments is determined from σ inel p-air using a nucleon-nucleon scattering amplitude which is frequently in disagreement with most accelerator data [1].
Although in principle QCD theory gives the exact description of strong interactions, its practical application in the study of hadronic interactions is still limited. In the absence of a pure QCD description, phenomenological models are used to compare experimental data with theoretical schemas. These models are based on general principles, such as 'unitarity', 'analycity' and 'cross-symmetry' and have proven to be successful in the comprehension and prediction of the hadronic amplitude behaviour at high energies [8]. Consequently, a wide range of possibilities is open for the development of models of this kind. Several classifications of phenomenological models have been made according to the general principles on which they are based (e.g. [9]), geometrical scaling models [10], diffraction-dominance models [11], factorized eikonal type models [12] and Reggeon-field-theory models [13]. At present, the more frequently used models (e.g. [14]) are (a) 'analytic amplitude' models based on solutions of the derivative dispersion relations; (b) 'eikonal' models. In the context of (a) outstanding work has been done by several authors (e.g. [6], [15]- [21]), however, the most complete and systematic work, to our knowledge, has been developed in association with the COMPETE Collaboration, in accord with the Regge-Pomeranchuk-Heisenberg type parametrizations and the 'general principles' previously mentioned [22]- [28]. Such systematization has been carried out on the selection and the cleaning of data as well as the application of models and fitting techniques, refining the conventional way (χ 2 ) of measuring the goodness-of-fits, by means of seven statistical indicators related to the data and to the process used to obtain them [25]. The results of this series of work have been of great relevance in the context of hadronic interactions. The present work falls within the framework of models of the type (b) in the impact parameter representation, using 'multiple diffraction' theory [29] to first-order and parametrization of the hadronic form factors and the elemental dispersion amplitude contained in the eikonal. We use an alternative statistical technique to evaluate error bars in the extrapolation process, called the 'forecasting' method. This work produces an improved accelerator extrapolation which both lowers the predicted p-p total cross-section curve and narrows its confidence bands, making it less consistent with most cosmic rays showers [9] and most of the [30] results as feared.
Thus, we dispose of many parametrization models (purely theoretical, empirical or semiempirical based) that fit the accelerator data fairly well. Most of them agree that at the energy of the future LHC (14 TeV in the centre of mass) or higher, the rise with energy of σ pp tot will continue, although the predicted values differ from model to model. Both the cosmic ray and accelerator approaches should complement each other in order to draw the best description of the p-p hadronic cross-section behaviour at ultra high energies. However, the present status is that since the interpolation of accelerator data is nicely obtained with most parametrization models, it has been hoped that their extrapolation to higher energies might yield high confidence values. The accelerator parametrizations are usually based on a small number of fundamental parameters, in contrast to the difficulties found in deriving σ pp tot from cosmic ray results [1]. With the aim of elucidating the problem, we first briefly analyse the way estimates are made for σ pp tot from accelerators in section 2 and from cosmic rays in section 3. We find serious discrepancies in both estimation methodologies. In section 4 we describe the 'multiple diffraction' model and the method used to evaluate the model parameters on the basis of only two data-based physical observables: the differential cross-section and the parameter ρ. For the goal of the present study, we neglect crossing symmetry. In section 5, we present a suitable parametrization to high energies 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT of the free energy-dependent parameters of the model, and we discuss its physical significance. In section 6 we describe the procedure for the determination of error bands as per appendix A. In section 7, on the basis of the 'multiple diffraction' model applied to accelerator data, we predict σ pp tot values with high confidence levels (CL). In section 8, we discuss our results in terms of the hypothesis σp p tot = σ pp tot , and finally in section 9, we conclude with a discussion of the implications of extrapolations within the framework of present cosmic ray estimates.

Hadronic σ pp tot from accelerators
Ever since the first results of the Intersecting Storage Rings (ISR) at CERN arrived in the 1970s, it is well established that σ pp tot rises with energy [31,32]. The CERN SppS collider found this rise valid for σp p tot as well [33]. Later, the Tevatron at Fermilab confirmed that for σp p tot the rise continues at 1.8 TeV, even if different experiments disagree as to the exact value [34]- [36]. A full discussion of these problems may be found in [37,38]. The amount of the rise of the total cross-section at those energies is still to be determined. Let us resume the standard technique used by accelerator experimentalists [6].
Using a semi-empirical parametrization based on Regge theory and asymptotic theorems, experimentalists have successively described their data from the ISR to the SppS energies. It takes into account all the available data for σ pp tot , ρ pp , σp p tot and ρp p , where ρ pp,pp , is the ratio of the real to the imaginary part of the (pp,pp) forward elastic amplitude at time zero. The fits are performed using the once-subtracted dispersion relations: where C s is the subtraction constant. The expression for σ pp,pp tot where −, + stands for pp (pp) scattering. Cross-sections are measured in mb and energy in GeV, E being the energy measured in the laboratory frame. The scale factor s 0 has been arbitrarily chosen to equal 1 GeV 2 . The most interesting term is the one controlling the high-energy behaviour, given by a ln 2 (s) factor, being compatible, asymptotically, with the Froissart-Martin bound [39]. The parametrization assumes σ pp tot and σp p tot to be the same asymptotically. This is justified by the very precise measurement of the ρp p parameter at 546 GeV at the SppS collider, ρp p = 0.135 ± 0.015 [40], which implies that at present there is no sizeable contribution of the odd under crossing part of the forward amplitude, the so-called 'Odderon hypothesis'. This hypothesis predicts a value of ρp p > 0.17-0.20 [41,42], that is, greater than the SppS value. The eight free parameters are determined by a fit which minimizes the X 2 function: The fit has proved its validity predicting from the ISR pp andpp data (ranging from 23 to 63 GeV in the centre of mass), the σp . σ tot exp from accelerator and from cosmic rays: the solid line indicates the best χ 2 fit obtained from a semi-empirical parametrization [6]. The two dashed lines delimit the region of uncertainty. Cosmic ray estimations from [43]. figure 1. We have also plotted the cosmic ray experimental data from AKENO (now AGASSA) [43] and the Fly's Eye experiment [44,45]. The curve is the result of the fit described in [6]. The increase in σ pp tot as the energy increases is clearly seen. Numerical predictions from this analysis are given in  8 mb). Also, recent analysis using both cosmic ray and accelerator data calculated relatively wide extrapolated error bands, for example 5-7 mb at 14 TeV [46]. It follows that ways to reduce the uncertainties and hence improve extrapolations are needed.

Hadronic σ pp tot from cosmic rays
Cosmic ray experiments give us σ pp tot indirectly from extensive cosmic ray air shower (EAS) data. As summarized in [1] and widely discussed in the literature, the determination of σ pp tot is a rather complicated process with at least two well-differentiated steps. In the first place, the primary interaction involved in and determined through EAS is proton-air, yielding the pinelastic cross-section, σ p-air inel , using some measure of the attenuation of the rate of showers deep in the atmosphere, m : The k factor parametrizes the rate at which the energy of the primary proton is dissipated into electromagnetic energy. A simulation with a full representation of the hadronic interactions in the cascade is needed to calculate it. This is done by means of Monte Carlo simulations [47]- [49]. Secondly, the connection between σ p-air inel and σ pp tot is model dependent. A theory for nuclei interactions must be used, which is usually Glauber's theory [29,50]. The whole procedure makes it difficult to get a general accepted value for σ pp tot . Depending on the particular assumptions made, the values may range by large amounts, from as low as 122 ± 11 at √ s = 30 TeV quoted by the Fly's Eye group [44,45] to 133 ± 10 mb by the 'Akeno Collaboration' [43], also at √ s = 30, to 162 ± 38 mb at nearly √ s = 30 (around 160-170 mb √ s = 40 TeV) [30] and even as high as 175 +40 −27 at √ s = 40 TeV [9]. It should be realized that the three values at 30 TeV do overlap within their errors as do the two values at 40 TeV. In the 40 TeV cases, even taking into account the large quoted errors, the values for σ pp tot are hardly compatible with the values obtained from the extrapolations of current accelerator data.
These results do not offer cosmic ray estimations of σ pp tot much help in constraining extrapolations from accelerator energies. Conversely we could ask if a reliable extrapolation based on accelerator data could be used to constrain cosmic ray interpretations.

A 'multiple diffraction' approach for evaluation of σ pp tot
Let us tackle the mismatch of accelerator and cosmic ray estimates using the multiple diffraction model applied to hadron-hadron scattering [51]. The elastic hadronic scattering amplitude for the collision of two hadrons A and B, neglecting spin, is described as where ξ(b, s) is the eikonal, b the impact parameter, J o the zero-order Bessel function and q 2 = −t the four-momentum transfer squared. In the first 'multiple diffraction' theory, the eikonal in the transferred momentum space is proportional to the product of the hadronic form factors G A 7 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT and G B (geometry) and the averaged elementary scattering amplitude among the constituent partons f (dynamics), and can be expressed to first order as ξ(b, s) = C A,B G A G B f , where the proportionality factor C A,B is the free parameter known as the 'absorption factor', and the brackets denote the symmetrical two-dimensional Fourier transform. A connection between theory and experiment may be obtained by means of hadronic factors and elementary partonparton amplitudes, which are not physical observables. However with the help of the optical theorem σ pp tot may be evaluated in terms of the elastic amplitude F(q, s): σ pp tot is a physical observable. The other two physical observables are the differential elastic cross-section and ρ expressed respectively as and 'Multiple diffraction' models differ from one another by the particular choice of parametrizations made for G A and G B and the elementary amplitude f . Hereafter we will adopt the model proposed in [7], which has the advantage of using a small set of five free parameters which are in principle energy dependent. Two of them (α 2 , β 2 ) are associated with the form factor The other three (C, a, λ) are associated with the elementary complex amplitude f where and so that with the opacity (b, s) given as  56 indicating that σ tot tends towards zero up to centre-of-mass energies ∼60 GeV. Data points are from [55]. whose explicit expression is where k 0 , k 1 , k ei , and k er are the modified Bessel functions, and E 1 -E 6 are functions of the five free parameters. The p-p total cross-section is directly determined by the expression which was numerically solved in [52,53]. It should be noted that, according to the principle of 'analyticity', the scattering forward amplitude for particle-particle and particle-antiparticle both come from the same analytical function [8]. The total cross-sections of both reactions are assumed to behave in one of the following ways: up to the ISR energies the differences are attributed to Regge contributions, which are expected to disappear at higher energies [5,54], or the differences are interpreted in terms of the 'maximal Odderon hypothesis', which predicts that they increase as the energy passes the highest energy of the ISR. At this point, the following understanding must be emphasized as essential in the rest of this study. Firstly, we quote the success in the prediction of σp p tot at the SppS collider [33] from the ISR data (mainly pp) using expressions (1) and (2), where σ pp tot and σp p tot were taken to be asymptotically equal [15]. Secondly, according to [55] the difference σ = σp p tot − σ pp tot , tends toward zero as s −0.56 up to energies 2000 GeV in the laboratory (∼60 GeV in the centre of mass; figure 2). And thirdly, we are aware that it can be argued that σ pp tot and σp p tot are different for higher energies, but as indicated in section 2 current evidence indicates the contrary.
In the face of this triple line of evidence, in our 'multiple diffraction' analysis, the same behaviour for both σ pp tot and σp p tot at high energy is assumed. So, here after σ pp tot = σp p tot = σ tot . It is noteworthy that some parametrization models, as RRP [56], predict the same value for both cross-sections at √ s > 70 GeV and the same value for the corresponding ρ at √ s > 110 GeV.

Evaluation of the model parameters
For the evaluation of the free parameters it was assumed, following [7] that two of them are constants: a 2 = 8.2 GeV 2 and β 2 = 1.8 GeV 2 . The other three parameters are: C(s) and α −2 (s) which determine the imaginary part of the hadronic amplitude (equation (11)) and consequently the total cross-section (equation (6)), as well as λ(s) which controls the real part of the amplitude (equation (12)). Setting λ(s) = 0 makes the amplitude purely imaginary, so that a zero (a minimum) is produced in the dip region, where only the real part of the amplitude becomes important [57]. We fit experimental data (figure 3) of the elastic differential crosssection, dσ/dq 2 , with equation (7), choosing the set (C, α −2 ) for which the theoretical and the experimental central values are equated at the precise t (GeV/c) 2 value where the data show its first minimum (the 'dip'position), which coincides with the minimum of the imaginary part of the elastic amplitude. Throughout we use units where c = 1. Data of pp differential cross-section at 13.8 and 19.4 GeV were taken from [58], for 23.5-62.5 GeV from [59], for 546 GeV (pp) from [60], with (−t < 0.5 GeV 2 ) and from [61] with (0.5 −t 1.55 GeV 2 ), and for 1800 GeV from [62] with (0.034 −t 0.65 GeV 2 ). The procedure is carried out at each energy where there are available accelerator data for dσ/dq 2 in the interval 13.8 √ s 1800 GeV, as illustrated in figures 4 and 5. Because data error bars at the dip position are small the fitting procedure is based on the central values. It must be emphasized that, it is precisely because the minimum of the imaginary amplitude is produced in the dip region, that data at 1800 eV are quite suitable for our procedure, since the predicted minimum falls at t = 0.585 GeV 2 where data are available   Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Table 2. Values of the parameters C, α −2 and λ at each energy. They are obtained by equating the accelerator data and the model prediction for the elastic differential cross-sections and for the parameter ρ in the interval 13 √ s 62.5 and 546 √ s 1800 GeV. (dσ/dt = 0.0101). The next datum at t = 0.627 GeV 2 shows a slight increase. Furthermore, the shift of the dip region towards lower values of t as energy increases is qualitatively consistent with the expectation relative to the shoulder at 546 GeV (t = 0.9 GeV 2 ). Once the values of C(s) and α −2 (s) are determined for each energy they are introduced along with the constant parameters a 2 and β 2 into equation (8).
We determine values for λ(s) which produce the experimental value of ρ(s), for each energy, √ s, where there are available accelerator data in the interval 13.8 √ s 1800 GeV. The central values of the three energy-dependent free parameters, C(s), α −2 (s) and λ(s) obtained are listed in table 2. Therefore, the method employed to evaluate the model parameters only requires dσ/dt and ρ(s) data.

Parametrization for extrapolation to high energies
Despite other successes of the theory, the microscopic basis, in terms of QCD, of the physics behind the parametrization and fitting procedures have not been completely developed because the hadron diffractive phenomenon in question is essentially a non-perturbative problem with which confinement is still unsolved. However, the closeness of our fits to the collider data (ISR, SppS and Tevatron) suggest that the basic approach of the limited model used here is correct, as described for instance in [51]. Indeed, even with the approximation made in ignoring crossing symmetry, the parametrization process followed here could not be arbitrary. That said, for interpolations and extrapolations, we made parametrizations of the three energy-dependent free parameters. Using the values obtained for those parameters, as described in table 2 and following the procedure to be described in section 6, a second-order fit of the values of C(s) and α −2 (s) and a exponential fit of λ(s) have been obtained from the following analytical expressions: α −2 (s) = 1.8956 − 0.03937 ln s + 0.01301 ln s 2 ,  Results are displayed in table 2 and illustrated in figures 6-8 as the central-solid lines. However, the reliability of the functional parametrizations for extrapolations must have physical support, since it is clear that several parametrizations that may correctly describe the data in the range 13.8 √ s 1800 GeV, may not remain consistent but differ substantially when extrapolated to high energies. Thus parametrization selections should be restricted according to the physical information available. As mentioned before, the fits of C(s) and α −2 (s) in the limited experimental range were based on experimental data of the differential cross-section and ρ(s), yielding values that increase with energy (table 2) with positive curvature (figures 6 and 7). Experimentally, total cross-sections increase with energy as ln s or ln 2s in the concerned energy range, and soft processes are expected to have a ln s behaviour. From the optical theorem, the interdependence of the free parameters and the physical observable (equation (6)) may be connected with the unitary condition, for which lowest order cross-sections within the frame of gauge field theories have ln s terms [63]. The fits, extrapolations, and constraints led naturally to ln s terms appearing in the two-energy dependent parameters C(s) and α −2 (s), and therefore the hypothesis of polynomial functions of ln s seems quite reasonable.
As for the parameter λ(s) a basic property of the 'Glauber multiple diffraction model' is to associate elastic scattering cross-sections of nucleons with the scattering amplitude of their composite partons [51]. Within this framework, the ratio of the real and imaginary parts of the parton-parton amplitudes (equation (12)), λ(s) = Re f(q, s)/Im f(q, s) behaves on the partonic scale as ρ(s) behaves on the hadronic scale. The influence of λ(s) on the hadronic amplitude has been empirically analysed in [7], showing that if λ(s) increases (or decreases), ρ(s) also increases (or decreases) (figure 9), and λ(s) = 0 at the same energy value where ρ(s) = 0. However, due to the lack of ρ(s) data above the experimental energy range used in this work, the parametrization of λ(s) at high energies is based on the conventional assumption that beyond √ s ∼ 100 GeV, ρ(s) has a maximum and then goes asymptotically to zero [37], the rate of convergence depending on the particularities of the model. Considering this and the empirical behaviour of λ(s), shown in table 2, we propose the parametrization given in equation (19), where s 0 = 400 GeV 2 is the value at which λ(s) converges to zero, and the numerical coefficients control the maximum and asymptotic behaviour.
Recall that blackening and expansion are very well-known properties of elastic scattering. Within the context of the present empirical analysis, blackening and expansion are related to the elementary parton-parton amplitude and the hadronic form factors through the energydependent parameters C(s) and α(s) respectively. Since in the straight forward direction the scattering amplitude is basically of diffractive nature and the eikonal becomes purely imaginary, for hadron-hadron ξ(b, s) = C(s) G 2 Im f(q, s) = Im (b, s), so that in terms of equations (9) and (11) the opacity satisfies (b, s) 1, and the free parameter C(s) behaves as an absorption factor (optical theorem). On the other hand, the free parameter α(s) 2 may be connected to the hadronic radius [7] as R 2 (s) = −6[dG(q, s)/dt] t=0 . Then, from equation (9) R 2 (s) = 0.2332 1 α 2 (s) Therefore, from equation (18) and the adopted value β 2 = 1.8 GeV 2 the radius is an increasing function of energy and such a behaviour expresses the expansion effect. The result is that hadrons become blacker and larger as energy increases, consistent with the so-called 'Bell effect' [64]. Since the hadronic scattering amplitude is purely imaginary, the free parameters may be associated with the physical observable by means of equations (6) and (7).

The extrapolation procedure
The procedure followed to obtain predictions of σ pp tot at high energies with confidence intervals based on the 'forecasting' statistical method described in appendix is as follows. (ii) Using the least-squares method in matrix formalism, as described in the appendix (equations (A.10)-(A.12)) we obtained the constantsβ 0 ,β 1 ,β 2 for each of the parameters. The autocorrelation constant was determined as described in [65,66].   16), but from the substitution of the extreme values of the error bands of the three energy-dependent parameters into equations (15) and (16), followed by the corresponding fits (figures 10 and 11).

Results
The total cross-sections with their respective errors are summarized in table 4. For √ s 62.5 GeV data were taken from [59] and for 546 GeV, from [33]. For the value at 1800 GeV there exist three different measurements: the value of the 'CDF Collaboration' (80.3 ± 2.24 mb) [35], the value of the 'E710 Collaboration' (72.8 ± 3.1 mb) [34] and the value of the 'E811 Collaboration' (71.7 ± 2.02 mb) [36].
There is controversy over whether the correct value is that of the 'CDF Collaboration' or of the lower values of the 'E710 Collaboration' and the 'E811 Collaboration', mainly related to the estimation of the different backgrounds. For instance, the 'CDF Collaboration' has decided to use only its σ pp tot value to quote luminosity values for all of its physics programmes [70] whereas the 'Tevatron Collaboration', D0, has adopted the average of the three measurements [70]. For this work, one of us (J V) suggested following this second option to use the arithmetic weighted mean of the three values (74.91 ± 1.35 mb). The effect on our results of the different values of the three collaborations is mentioned in the next section. We quote in  Figure 10 represents graphically the first case together with cosmic ray data. For further discussion, we include in table 4 predicted values from two standard accelerator-based data extrapolations: the first one, σ [21] 546 makes a fit to all σ pp tot and σp p tot data in the interval 10 √ s 546 GeV using expression (2) [21], and the second one, σ [6] ρ (figure 1), in the same energy range, makes the simultaneous fit to all σ pp tot , σp p tot , ρ pp tot and ρp p tot data through the method described in section 2 [6]. Table 4. Central, upper and lower values for σ tot obtained with the procedure of section 6 and the method described in the appendix. The σ 1800 tot , σ 546 tot , σ 62.5 tot columns include data up to 1800, 546, 62.5 GeV respectively. The quoted σ [21] 546 and σ [6] ρ values are from [21] and [6] respectively.

Discussion
Let us first examine what happens when our main assumption, the asymptotic equality of σ pp tot and σp p tot , is not used. Analysis of figure 11 shows that if we limit our fitting calculations to the accelerator domain √ s 62.5 GeV, where σ pp tot data exist, the extrapolation at cosmic ray energies produces an error band so large that potentially any cosmic ray result becomes compatible with results at accelerator energies. It can be seen that, in this case, the σ pp tot values obtained when extrapolated to ultra high energies seem to confirm the highest quoted values of the cosmic ray experiments [9,30]. Also it can be noted that such extrapolation to ultra high energies may claim not only agreement with the analysis carried out in [30] and the experimental data of the Fly's Eye [9], but even with the Akeno Collaboration [43], because their experimental errors are so big that they overlap with the errors reported in [9], and of course fall within the predicted error band for that case √ s 62.5 GeV ( figure 11). If true, that would imply that the extrapolations cherished by experimentalists are meaningless. But the prediction shown in figure 11 gives σ pp tot = 69 mb at the CERN SppS Collider (546 GeV), and 91.6 mb at the Fermilab Tevatron (1.8 TeV). Comparing with table 1 we see that the measured σp p tot at 546 GeV is smaller than the predicted σ pp tot by nearly 8 mb, and that at 1.8 TeV by more than 15 mb, which no available model is able to explain [38]. However, when the initial hypothesis, σp p tot = σ pp tot asymptotically is used, then the existing σp p tot data at higher accelerator energies may safely be included. This permits enlargement of the lever arm for the extrapolation by a great amount, and both the predicted values and the error band change considerably. This can be clearly seen in figure 10 [9,30] by several standard deviations, although not so different from the Fly's Eye or Akeno results, nor from the predicted value in [6].
The quoted error bands, obtained as described in section 6 and the appendix illustrate that the inclusion of 'residual correlations' in the 'forecasting' method produce a neat fit, irrespective of the parametrization model. This seems to be an advantage with respect to other statistical techniques. However, it must be emphasized that it is not valid to compare different statistical techniques when using different parametrization models, or when running the same parametrization with different input values. In spite of this, we would like to make a strict qualitative observation, in the sense that other parametrizations [6,21] with similar input quantities (for instance ρ and dσ/dt) lead to central values that are only slightly higher than ours, whereas the quoted errors are larger than ours, some of them by nearly a factor of 3, as can be seen in table 4, or in figures 4 and 10. For instance, the quoted error in the fit σ 1800 tot at 100 TeV in the centre of mass (which corresponds to 10 19 eV in the laboratory), 147 +9.63 −7.68 , is comparable to (or even better than) the error obtained at much more lower energies in other works as can be seen in table 4, at 16 and 14 TeV ( 10 17 eV in the laboratory), corresponding to the energy range of future LHC CERN Collider. Also, it can be appreciated from table 2 in [6], that the half width of the cross-section uncertainty bands, 7.2 and 10%, at 16 and 40 TeV respectively, is higher than the equivalent in this work, 6.4 and 7.9% (the values quoted in column 5 of table 4); in both cases the values of σ pp tot have been estimated using as lever arm the value at 546 GeV. Obviously, these comments are only suggestive but important for the questions they raise: though it is clear from statistical theory that the forecasting technique is a highly precise method (section 6, appendix and references therein), its advantage over other techniques is, however, in principle only a second order effect (as inferred from equations (A.1) and (A.8)) when applied to the same conditions. Does the goodness-of-fit improve for certain parametrizations and input because the parametrization model is better? The answer must depend on a deep quantitative analysis, applying the forecasting method to different parametrization models of the total crosssection, in order to discern which model yields the best confidence estimates, with the goal of drawing more reliable physical inferences. Such an analysis is beyond the scope of the present paper and will be the subject of a forthcoming work.

20
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Conclusions
It has been shown in this work that high CL extrapolations to high energy σ pp tot values are strongly dependent on the energy range covered and the available number of data values. In particular, if we limit the input σ pp tot accelerator data to the range √ s 62.5 GeV, for the extrapolation to cosmic ray energies, the results are compatible with most of the cosmic ray experiments, and other predictive models [6], because the predicted error band is wide enough to cover their quoted errors ( figure 11). However, as data to be extrapolated reach higher energies, including σp p tot data up to 1.8 TeV, that is, when all experimental available data are taken into account, the curve is lowered by this last datum and the extrapolated error bands are much reduced. Accepting the accelerator extrapolations as indicative of the most probable σ pp tot , we conclude that they might help normalize the interpretation of cosmic ray experiments, perhaps, as an example, suggesting that the parameter k be kept free. The k value thus computed could help to tune the complicated Monte Carlo calculations used to evaluate the development of the showers induced by cosmic rays in the upper atmosphere. Extrapolations from our parametrization model would imply that σ inel p-air should be smaller than usually considered, which would have important consequences for the development of high energy cascades.
Although found at the Tevatron could be made. Additionally, experiments such as the HiRes [2] and the Auger Observatory [3] will bring new light to the extrapolation problem, i.e. whether the real values of σ tot are described by parametrizations consistent with a fast rise with high energies, to be called the 'cosmic result' (e.g. [71]) or with a slower rise with energy, to be called the 'accelerators result', the result presented in this paper. Instead of looking for a new physics based on 'exotic' events to explain the fast rise in energy of σ tot , what we need is highly reliable data at intermediate energies (2 √ s 15 TeV) to be used as confidence lever arms for extrapolations, to more accurate σ tot at cosmic ray energies. In should be kept in mind that extrapolations from accelerator data must be included in interpreting cosmic ray phenomena and determining crosssections. Although there are basic differences of this work and the COMPETE Collaboration [22]- [28] regarding the model employed and the techniques for extrapolation, the statistical method for calculation of the confidence intervals and the value of the Tevatron collider crosssection data used, nevertheless, it must be noted that the conclusions reached are very similar (e.g. [26]). Another parametrization with central value predictions not so different from ours, is that of the 'Ensemble A', discussed in [67,68].
Finally, it is worth mentioning that this work reduces the width of the confidence band around 'multiple diffraction' model fits of accelerator data, strengthening suspicions that this extrapolation disagrees with cross-sections derived from cosmic ray showers. The reasons for the disagreement can be in the theory and method of analysis of the showers, in misunderstandings of the measurements at both ends, in the incompleteness of our understanding of p-p cross-sections, etc. The recognition and quantification of the discrepancy contributed here is a necessary step. As the reduction in uncertainty is largely due to the inclusion of higher energy accelerator data 22 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT where S 2 d is the variance, defined as the square of the standard deviation (usually called the standard error of estimate),ȳ n+m is the corresponding central prediction and t {k−p} {δ/2} denotes the probability density function (pdf) known as the Student's t-distribution for the k values of the independent variables with p degrees of freedom, and a confidence coefficient δ/2. However, the previous generalization does not incorporate effects related to the position of the experimental values of y around the employed central value, that is the correlation among residuals that generates the regression model [74]. In other words, no interaction terms are assumed, implying that each of the independent variables affects the response y independently of the other independent variables and the error associated with any one y value is independent of the error associated with any other y value. Such an omission of the residuals' correlation leads to a notorious modification of the statistical estimators.

A.2. The forecasting method
The above limitation, can be surmounted by identifying the obtained dataset (x 1i , x 2i , x 3i , . . . , x ki ; y i ) with a time series, and then evaluating the correlation among consecutive residuals (the so-called autocorrelation procedure, of which the simplest one is the autoregressive 1st-order model). The forecasting statistical method is based on autocorrelation models adopted for the evaluation of the correlation among consecutive residuals, in and out of the data set by means of an iterative process [74]. In general, this procedure modifies the fitting constants of the model and the estimated variance of residuals and minimizes the width of the error intervals for interpolation or (extrapolation), in the process increasing the CL of predictions [74]. To quantify the effect that the autocorrelation of residuals has on the regression model and associated estimators, let us represent such a model by a response variable y = E(y) + , where E(y) = β 0 + β 1 x 1 + β 2 x 2 + · · · + β k x k , (A.2) represents the deterministic component of the proposed regression model; x 1 , x 2 , . . . , x k are the independent variables, accepted as assigned and may represent higher order terms and even functions of variables as long as the functions do not contain unknown parameters; β 1 , β 2 , . . . , β k are the unknown coefficients to be determined by the least-squares method, representing the contribution of the independent variables x i , and represents the random error component. For the evaluation of residuals R i = y i − E(y i ) it is assumed they have a normal distribution with mean zero and constant variance. In the autoregression model of first order each residual R i is related with the previous one as where φ is the autocorrelation constant among the residuals (|φ|< 1) [65,66], and r i in this case is a residual called white noise, uncorrelated with any other residual component. The incorporation of the effect of autocorrelation of residuals to the solution of the regression problem through equation (A.2) leads to a modification of the regression constants and the corresponding variance: using the data set, the estimate of the k + 1 regression constants of equation (A.2) and the constant φ is obtained from the autocorrelation model according to the following interpolation equation for the response variableŷ i : y i =β 0 +β 1 x 1,i +β 2 x 2,i + · · · +β k x k,i +φR i−1 , (A.4)

A.3. Matrix approach in regression analysis
This fit method is based on multiple regression theory and consists in creating a prediction equation for a quantity y (dependent variable), which depends on k independent variables (x i ), that is Therefore, the application of a multiple regression model to a given problem leads to a system of n equations with n incognitos, so that its solution is better obtained through a matrix formalism. Denoting by Y the matrix of (n × 1)-dimension of the dependent variables and by X the matrix of [n × (k + 1)]-dimensions of the k independent variables, the row 1, x 11 , x 12 , . . . , x 1k multiplied by the column matrix of the βs determines the value y 1 of the dependent variable, the row 1, x 21 , x 22 , . . . , x 2k multiplied by the column matrix of the βs determines the value y 2 and so on: where (X t X) −1 denotes the inverse matrix of X t X. Essentially, this equation minimizes the quadratic sum of the deviations of points y j with respect to the fitted function (A.9) ( [74] p 783).
With the previous matrices, several statistical estimators are easily determined, such as the sum of square errors (SSE) .13) and the variance required to evaluate the confidence intervals, which is , (A.14) where the denominator defines the number of degrees of freedom for errors, given by the number of β i -parameters. Once the 'central' values are known, we evaluate the confidence interval for a particular value of the response variable, y p , using the matrix of the particular values of the independent variables which determine the estimated value of y p . Such a matrix, namely A, denotes the column-matrix of (k + 1) × 1 dimensions, where elements {1, x 1p , x 2p , . . . , x kp } correspond to the numerical values of the β i appearing in equation (A.9). Therefore, the confidence interval for prediction within the range of data is determined as ( [74], p 795): Interp B =ŷ ± t n−p δ/2 S 2 f A t (X t X) −