Features of the energy spectrum of cosmic rays above $2.5{\times} 10^{18}$ eV using the Pierre Auger Observatory

We report a measurement of the energy spectrum of cosmic rays above $2.5{\times} 10^{18}$ eV based on $215,030$ events. New results are presented: at about $1.3{\times} 10^{19}$ eV, the spectral index changes from $2.51 \pm 0.03 \textrm{ (stat.)} \pm 0.05 \textrm{ (sys.)}$ to $3.05 \pm 0.05 \textrm{ (stat.)}\pm 0.10\textrm{ (sys.)}$, evolving to $5.1\pm0.3\textrm{ (stat.)} \pm 0.1\textrm{ (sys.)}$ beyond $5{\times} 10^{19}$ eV, while no significant dependence of spectral features on the declination is seen in the accessible range. These features of the spectrum can be reproduced in models with energy-dependent mass composition. The energy density in cosmic rays above $5{\times} 10^{18}$ eV is $(5.66 \pm 0.03 \textrm{ (stat.)} \pm 1.40 \textrm{ (sys.)} ) {\times} 10^{53}~$erg Mpc$^{-3}$.

Although cosmic rays having energies above 10 19 eV were first detected nearly 60 years ago [1,2] and are being investigated by the two largest-ever built detectors, the Pierre Auger Observatory [3] and the Telescope Array [4], the question of their origin remains unanswered. Only recently has the belief that such particles are of extragalactic origin been demonstrated experimentally with the discovery of significant directional anisotropies above 8×10 18 eV [5]. These data are well-described by a dipole pattern, the amplitude of which increases from 6% to 10% as the energy rises to 4×10 19 eV [6]. An important observable for an understanding of ultrahigh energy cosmic rays (UHECRs) is the energy spectrum. We report a measurement above 2.5×10 18 eV based on 215,030 events, over 10 times that used in [7]. Over 16,000 events have energies beyond 10 19 eV. This spectral determination is unique in making no assumptions about the mass composition or the hadronic physics. Full details are reported in [8].
UHECRs can only be studied through the detection of the showers of particles (extensive air-showers) they create in the atmosphere. A calorimetric estimate of the energy carried by the primary particle is possible using telescopes to collect the fluorescence light emitted by atmospheric nitrogen excited by the shower. The on-time of this technique is limited to nights with low-background light while, by contrast, an array of particle detectors deployed on the ground can be operated with a duty-cycle close to 100%. The traditional method of assessing the energy of the primary cosmic-ray from observations made with the particle detectors requires assumptions about its mass and the hadronic processes that control the cascade development. This is clearly unsatisfactory as the mass is unknown and the centre-of-mass energy reached at the LHC corresponds only to that of a proton of 10 17 eV colliding with a nitrogen nucleus. Also, details of pion interactions, which play a key role in shower development, are lacking. The presence of unknown processes could also lead to hidden systematic uncertainties.
To circumvent these limitations, the energies are obtained by making use of a subset of events detected simultaneously by the fluorescence detector (FD) and the particle surface detectors (SD). This "hybrid" approach allows a calorimetric estimate of the energy for events recorded during periods when the FD cannot be operated. A spectrum can thus be derived that is free from assumptions about primary mass or hadronic physics.
The Pierre Auger Observatory is such a hybrid system [3]. It is sited near the city of Malargüe, Argentina, at latitude 35.2 • S with a mean atmospheric overburden of 875 g cm −2 [3]. The SD comprises 1600 water-Cherenkov detectors deployed on a 1500 m triangular grid, covering about 3000 km 2 . The array is overlooked from four stations, each containing six telescopes used to detect the emitted fluorescence light. Comprehensive atmospheric monitoring, particularly of the aerosol content and the cloud cover, is undertaken [3,9].
The SD samples the shower particles that reach the ground. Signals in the individual detectors are quantified in terms of their response to a muon travelling vertically and centrally through it (a vertical equivalent muon or VEM). The signals are used to determine the impact point of the shower axis, the arrival direction and the shower size. For the latter, the signal at 1000 m from the shower axis, S(1000), is used. For the grid spacing of 1500 m, this is the distance that minimizes the uncertainty arising in S(1000) from the imperfect knowledge of the functional form describing the fall-off of signal with distance from the shower axis in individual events [10].
Showers detected by the SD arrive from a range of zenith angles, and they are attenuated according to how much atmosphere is traversed. Accordingly, for each event, S(1000) is adjusted to a reference value, S 38 , the magnitude that it would have had, had the cosmic ray arrived at the median zenith angle of 38 • . The long-established procedure for making this correction, the Constant Intensity Method [11], relies on the quasiisotropy of cosmic rays in zenith angle given the small anisotropy contrasts in celestial coordinates [8]. The large number of events has made it possible to refine the original approach and quantify the change in shower absorption as a function of energy. Such an evolution is anticipated as, at a given zenith angle, the ratio of the muon to electromagnetic components falls as the energy increases, even for an energy-independent composition.
For showers detected by the FD, it is possible to measure the deposition of energy lost to ionisation of the atmosphere using a fit to a modified Gaisser-Hillas profile [12]. The integration of the profile yields a calorimetric measure of this loss. The energy of the primary particle, E FD , is then obtained by the addition of an energy-dependent correction of less than 14%, driven by data [13], to allow for the "invisible energy", carried into the ground by muons and neutrinos. The resolution in E FD is well-described by σ FD (E)/E 7.4% over the whole energy range [14].
Hybrid events are thus used to develop a calibration curve such that every estimate of S 38 can be assigned a valuation of E FD . Here 3,338 hybrid events surviving rigorous quality cuts [8] are used to obtain a relationship between S 38 and E FD of the form E FD = A S 38 B , where A = (1.86 ± 0.03)×10 17 eV and B = 1.031 ± 0.004. No zenithal dependence of A or B has been found, further validating the use of the Constant Intensity Method [8]. Such a simple dependence is sufficient to describe the data in full detail. The energies of the hybrid events range from 2.5×10 18 eV to 8×10 19 eV. The most energetic event, detected at all fluorescence stations, has an energy E FD = (8.5 ± 0.4)×10 19 eV, derived from a weighted average of the four independent estimates of the calorimetric energy. For this event S 38 = 354 VEM so that the energy deduced from the calibration curve is E SD ≡ AS 38 B = (7.9 ± 0.6)×10 19 eV. The systematic uncertainty in the energy assignment is about 14% over the whole energy range [15]. This benefits from the high-precision AIRFLY measurement of the fluorescence yield [16] and from an accurate data-driven estimation of the invisible energy [13]. Other contributions to the uncertainty are related to the estimation of the A and B parameters, the characterization of the atmosphere, the reconstruction of the longitudinal profile and the FD calibration, which provides the largest contribution.
To derive the energy spectrum, we use events recorded by the SD with the largest-signal station not located on the boundary of the array, with zenith angle θ < 60 • and energy ≥ 2.5×10 18 eV. These selection criteria not only ensure adequate sampling of the shower but also allow the evaluation of the aperture of the SD in a purely geometrical manner in the regime where the array trigger is fully efficient and independent of the mass or energy of the primary particle [17]. The resulting SD data set consists of 215,030 events recorded between 1 January 2004 and 31 August 2018, from an exposure, E, of (60,400 ± 1,810) km 2 sr yr. The determination of E, dependent only on the acceptance angle, the surface area and the live-time of the array, is discussed in detail in [17].
The procedure for extracting the spectrum from the observations, fully discussed in [8], is summarised here.
The energy spectrum, typically a power law (∝ E −γ ) with spectral index γ in a given energy interval, is estimated as J i = c i N i / (E∆E i ), with N i the number of observed events in differential bins of width ∆ log 10 E i = 0.1 and c i the correction factors required to eliminate the biases caused by the finite energy resolution. The size of the bins is such that it corresponds approximately to the energy resolution in the lowest energy bin, which starts at 2.5×10 18 eV.
The correction factors are needed because, as the spectrum is steep, the finite resolution causes migration between bins, particularly from lower to higher energies, artificially enhancing the flux. At the lowest energies, the correction depends also on the behaviour of the detection efficiency in the energy region where the array is not fully efficient as well as on the bias in the energy due to trigger-selection effects.
A forward-folding approach is used to determine the correction factors. It consists of finding the model of the energy spectrum folded for detector effects that best describes the data, and then using this model to calculate the values of c i . The SD efficiency can be estimated from the fraction of hybrid events that also satisfy the SD trigger conditions, because above 10 18 eV, the hybrid trigger efficiency is 100% independent of primary mass [18]. The energy resolution of E SD , and the bias in its estimate, are found from a study of the distributions of E SD /E FD . The resolution improves from ≈ 20% at 2×10 18 eV to ≈ 7% at 2×10 19 eV and is constant thereafter. The bias is zero above 2.5×10 18 eV and increases smoothly going to lower energies and larger zenith angles: at 10 18 eV it is ≈ 10% at 0 • and ≈ 30% at 60 • .
Thanks to the hybrid measurements, the correction factors are estimated avoiding any reliance on model and primary mass assumptions. The factors are maximal at the lowest energies, ≈ 8%, and less than 5% at the highest energies available. Further details are given in [8].
The model of the energy spectrum that we used for over a decade is a series of two power laws followed by a slow suppression. With the current exposure, this model turns out to describe the data poorly, as the reduced deviance is found to be 35.6/15 [8]. Consequently, we adopt a more complex function with a sequence of four power laws with smooth transitions [19], with j = i + 1 and ω ij = 0.05. The ω ij factors control the widths of the energy intervals over which the slope transitions occur [8]. This model describes the data with a reduced deviance 17.0/12, which allows us to disfavor the previous parameterization with 3.9σ confidence [8].
The resulting differential energy spectrum and the fitted function are shown in Fig. 1. The normalization is J 0 = (1.315±0.004±0.400)×10 −18 km −2 sr −1 yr −1 eV −1 . The ankle is described by a rollover at E 12 = (5.0 ± 0.1 ± 0.8)×10 18 eV, marking a hardening of the spectrum from γ 1 = 3.29±0.02±0.10 to γ 2 = 2.51±0.03±0.05. At E 23 = (13±1±2)×10 18 eV, the spectrum softens from γ 2 to γ 3 = 3.05 ± 0.05 ± 0.10. Finally, the spectrum softens further above a suppression energy of E 34 = (46 ± 3 ± 6)×10 18 eV with γ 4 = 5.1 ± 0.3 ± 0.1, confirming with higher precision previous reports of the strong attenuation of the flux at the highest energies [7,20,21]. The feature at E 23 , calling for a 2-step suppression, is a new observation. For all parameters and observables presented in the text, the first error is statistical and the second systematic. From the measured energy spectrum one can infer the differential energy density per dex 1 , obtained as ln (10) (4π/c) E 2 J(E). It provides a measurement of the energy density of the local Universe attributable to cosmic rays. Above the ankle, a range in which UHECRs are of extragalactic origin [5], the integration over energy results in (5.66 ± 0.03 ± 1.40)×10 53 erg Mpc −3 . This translates into constraints on the luminosity of the sources, as discussed below.
A detailed examination of the systematic uncertainties of the energy spectrum is reported in [8]. The uncertainty in the flux amounts to 30 -40% near 2.5×10 18 eV, 25% at 10 19 eV and 60% at the highest energies. The uncertainties include contributions from the absolute energy scale (the largest), the exposure, the unfolding procedure and the S(1000) reconstruction. No indication of further systematic uncertainties has been found from a comparison of the spectra calculated over different time periods, seasons and ranges of zenith angle. The wide declination range covered, from δ = −90 • to δ = +24.8 • , allows a search for dependencies of energy spectra on declination. For this, we have divided the sky into three declination bands of equal exposure. In each band, the estimation of the spectrum is made as for the whole field of view, but using unfolding-correction factors relevant to the band in question. We report in 1 dex indicates decade in log 10 E, following the convention of [22]. Table I the parameters characterizing the spectral features for each declination range. They are seen to be in statistical agreement. There is thus no obvious dependence with declination over the energy range covered. A trend for the intensity to be slightly higher in the Southern Hemisphere is observed [8], consistent with the anisotropy observations [6]. We therefore claim a second new result, namely that the energy spectrum does not vary as a function of declination in the range accessible at the Auger Observatory other than in the mild excess from the Southern Hemisphere expected in line with the known energy-dependent anisotropies above 8×10 18 eV. A comparison of the spectrum with that of Telescope Array measured in the Northern hemisphere is discussed in [8] and references therein.
Astrophysical implications of the features of the energy spectrum. We now examine the validity of models proposed to explain features of UHECRs using the new information given here and the data on mass composition and arrival directions recently reported [5,6,[23][24][25][26][27][28]. If UHECRs are produced throughout the Universe, to reach Earth they must cross the background photon fields permeating the extragalactic space. In particular, the cosmic microwave background photons induce pionproduction with protons colliding at around 5×10 19 eV and photo-disintegration of heavier nuclei at a roughly similar threshold, leading to the expectation of a spectral steepening (the Greisen-Zatsepin-Kuzmin (GZK) effect [29]). Depending on the energy and chemical composition of the UHECRs, higher-energy background photons, such as infrared light, may also be responsible of interactions producing the flux steepening.
A popular framework has been that what is observed comes from universal sources, uniformly distributed, that accelerate only protons. As a consequence the ankle is then explicable by energy losses of protons through pair production across greater distances [30][31][32] so that the ankle region would be proton-dominated. However, recent results [26] strongly contradict this expectation: in the ankle region, (3 − 5)×10 18 eV, a pure proton composition, or one of only protons and helium, is excluded at the 6.4 σ level. A second consequence [33] concerns the energy, E 1/2 , at which the integral intensity falls by a factor two with respect to a power-law extrapolation from lower energies. The prediction in this framework is that E 1/2 = 5.3×10 19 eV though this number may change with fluctuations of source luminosities and densities that shape the GZK feature [31,34], and with the maximum achievable energy in the sources. The value found here, (2.2 ± 0.1 ± 0.3)×10 19 eV, is at variance with the prediction because of the new feature of the spectrum at ≈ 10 19 eV, which is absent in the popular paradigm that is thus disfavored.
Relaxing the universality of the source spectra, the steepening at ≈ 10 19 eV could stem from the distinctive spectrum of a local source that emits protons and contributes significantly to the total intensity. At these energies, diffusive propagation of protons from a nearby source is excluded by limits set on extragalactic magnetic fields from rotation measures [35]. Approximately, protons would thus arrive to the Galaxy as a uniform, parallel beam that may subsequently be focused or defocused while propagating in the Galactic magnetic field. As seen from the Earth, the image of the source is expected to be shifted and broadened, with the effect growing with decreasing energy. Also, multiple broad images may be produced if uncorrelated regions of the magnetic field are experienced by the particles [36][37][38]. Such a scenario would thus imply the observation of an anisotropy at intermediate angular scales, the size of which depends on the model of turbulence for the magnetic field [39]. Spectral differences would also consequently appear in some parts of the sky. The softening at ≈ 10 19 eV, in particular, would not be expected in every declination range. The absence of such dependence accordingly disfavors the interpretation that the steepening is due to a source in the local Universe emitting protons. Furthermore, the interplay between the luminosity of a given source and its flux attenuation with distance requires fine-tuning to make viable a scenario in which several sources would emit protons with a distinctive spectrum while at the same time no directional effect would be imprinted upon the observed intensity. Energy density [erg Mpc By contrast, our results fit a scenario in which several nuclear components contribute to the total intensity and in which the electromagnetic fields permeate source environments where nuclei are accelerated to a maximum energy proportional to their charge (Z). This scenario, e.g., [40][41][42][43], provides a natural framework to explain the tendency towards heavier masses with increasing energy as inferred from recent works [23][24][25]. To illustrate the main physics aspects without distraction by the many details a full model scenario would require, we consider here, as in [43], several nuclear components injected at the sources with a power-law spectrum and with the maximal energy of the sources modeled with an exponential cut-off. For simplicity, the sources are assumed to be stationary and uniform in a co-moving volume. We show in Fig. 2 the best reproduction of the data by simultaneously fitting the energy spectrum above 5×10 18 eV and the distribution of the depths of the shower maximum (X max ), which is mass-sensitive (using EPOS LHC [44] as model of hadronic interactions in their interpretation). The abundance of nuclear elements at the sources is dominated by intermediate-mass nuclei accelerated to ≈ 5 Z×10 18 eV and escaping from the source environments with a very hard spectrum. In this scenario, the steepening observed above ≈ 5×10 19 eV results from the combination of the maximum energy of acceleration of the heaviest nuclei at the sources and the GZK effect. The steepening at ≈ 10 19 eV reflects the interplay between the flux contributions of the helium and carbonnitrogen-oxygen components injected at the source with their distinct cut-off energies, shaped by photodisintegration during the propagation. We note that the ratio E 34 /E 23 is 3.4 ± 0.3, matching the mass/charge ratio of CNO to He, as expected from the benchmark scenario shown in Fig. 2.
Some cautionary comments on the illustrative model considered here are in order. The presence of a subdominant light component at the highest energies is not excluded by our data, see e.g. [45]. Also, viable source scenarios can be found without resorting to a mixed composition with a rigidity-dependent maximum energy if, for instance, predominately heavy (Si to Fe) nuclei are accelerated and photo-disintegrate in the source environment [46] or en route to Earth [47,48]. Scenarios with a predominantly light composition [32,49] can fit our X max data as well as those of Telescope Array [50] in the ankle energy range, but these scenarios are at odds with measurements of the correlation of particle densities at ground and X max [26]. At ultra-high energies, a significant re-adjustment of current hadronic interaction models would be required [51] to fit our data with a p/He-dominated model while the data of Telescope Array, because of limited statistical power above 10 19 eV, cannot yet be used to draw reliable conclusions about composition in this energy range [52].
Interactions of the accelerated nuclei in the environment of the sources may give rise to copious fluxes of nucleons below the ankle energy, produced through photodisintegration. Neutrons escaping from the magnetic confinement regions may then explain the observed flux of protons deduced from X max measurements [24,25] in this energy range, due to neutron decay during propagation [46,[53][54][55][56][57]. To make up the all-particle spectrum and to fit the composition data below ≈ 5×10 18 eV, an additional component is further required (see e.g. [58][59][60] for discussions). This could be the high-energy tail from the sources emitting the bulk of Galactic cosmic rays of lower energies or, as in the "B-component scenario" [61], further explored in [62,63], this additional high-energy component is produced by different sources in the Galaxy.
Finally, within this scenario, our data constrain the luminosity density that continuously emitting sources must inject into extragalactic space in UHECRs to supply the observed energy density. This amounts to ≈ 6×10 44 erg Mpc −3 yr −1 above 5×10 18 eV at a redshift of zero, in line with the value of 2×10 44 erg Mpc −3 yr −1 that can be inferred dividing the measured energy density by the typical cosmic-ray energy loss time O(1 Gpc/c) (3.3 Gyr) [64]. Classes of extragalactic sources that match such rates in the gamma-ray band include active galactic nuclei and starburst galaxies [65]. The flux pattern from these objects also provides an indication of anisotropy in UHECR arrival directions [27,28].