The energy spectrum of cosmic rays above 6 PeV as measured at the Pierre Auger Observatory

. Since 2004, the Pierre Auger Collaboration has measured, with unprecedented precision, one of the most important features of ultra-high-energy cosmic rays, the energy spectrum. Located in the Southern hemisphere, the Observatory comprises an array of 1660 water-Cherenkov detectors covering 3000 km² overlooked by 27 fluorescence telescopes. Five sets of measurements have been used to reconstruct the spectrum from 6 PeV to beyond 100 EeV. The highest-energy events, recorded by surface detectors on a 1500 m triangular grid, are reconstructed di ff erently depending on their inclination (vertical events below 60° and inclined ones above 60°). A cross-check of these measurements is made using ‘hybrid events’, in which there are simultaneous detections at both the fluorescence detectors and at least one surface detector. Events with energies below 3 EeV have been studied using a nested array with a spacing of 750 m and using the Cherenkov light recorded with three high-elevation telescopes. In this contribution, updated methods for the reconstruction of all events are discussed, together with their associated uncertainties. A combination of all data sets is reported to construct a spectrum above 6 PeV to the highest energies. A detailed discussion of the spectral features will be presented.


Introduction
More than a century ago, cosmic rays have been observed for the first time and decades ago, the first cosmic ray with an energy larger than 10 EeV has been detected. However, at the highest energies ever observed, the origin of these extraordinary particles is still a source of debate and unknowns among the scientific community. One of the keys to unlock the path to the origin of ultra-high-energy cosmic rays (UHECRs) is to reconstruct the energy spectrum. Thus, experiments have been built or are still operating to measure as precisely as possible the energy spectrum of cosmic rays. When it comes to the highest energies ever observed on Earth, a precise measurement of the spectrum has been reported by the Pierre Auger Observatory [1].
The Pierre Auger Observatory is, at the time when this proceeding is written, the largest ground-based observatory ever built dedicated to the indirect detection of cosmic rays through their interaction in the atmosphere of the Earth and the creation of Extensive Air-Showers (EAS). It is located in the Southern hemisphere, in the province of Mendoza, Argentina, at an altitude around 1400 m above sea level. It consists in two main parts: a Surface Detector (SD) with 1600 water-Cherenkov detector stations, deployed on a triangular grid with a 1500 m spacing, covering an area larger than 3000 km 2 . The SD is overlooked by 24 telescopes located at the edges of the SD, thus forming the Fluorescence Detector (FD) and giving the main characteristic of the Pierre Auger Observatory: its hybrid * e-mail: spokespersons@auger.org design. Both detectors, SD and FD, have been designed to work independently from each other, thus offering the possibility to simultaneously observe the lateral and longitudinal profiles of a shower initiated by a cosmic ray. While FD is working only during moonless and clear nights, reducing its duty cycle ∼15%, SD is operating 24 hours per day in almost all weather conditions. Thanks to this, the Pierre Auger Observatory is able to collect a large statistic of events with the SD, for which the calibration in energy is ensured by a use of a subset of events that can be reconstructed by both detectors at the same time, as it will be described later.
Initially designed to observe UHECRs with an energy above 3 EeV to get a precise measurements of the ankle at ∼5 EeV and the flux at the highest energies with ≳50 EeV, the Pierre Auger Observatory experienced several upgrades. In particular, an extension towards lower energies with the aim to observe the region of the second knee at ∼0.1EeV. For the SD, this extension consists in building inside the existing array, a denser region with 61 water-Cherenkov detectors with a spacing of 750 m and covering an area of 24 km 2 . This infilled array is then overlooked by 3 High-Elevation Auger Telescopes (HEAT), observing the higher altitudes and detecting events dominated by the Cherenkov emission. Several detection and reconstruction techniques have thus been used to reconstruct the energy spectrum of cosmic rays from 6 PeV to the highest energies. Alongside with the data sets used, these techniques will be described below.

Five data sets from two types of detectors
To reconstruct the energy spectrum of cosmic rays starting at 6 PeV, several instruments and methods have been applied. Thus, not less than five different data sets are used by the Pierre Auger Collaboration: two data sets for the FD, the hybrid and Cherenkov, and three data sets for the SD, SD1500-vertical for events with a zenith angle θ below 60 • , SD1500-inclined for events with 60 • < θ < 80 • and observed using the stations distant from each other by 1500 m, and SD750 for events observed with the denser array. In the following, main characteristics of the data sets and their reliability will be presented.

Fluorescence Detector
The fluorescence detector of the Pierre Auger Observatory is designed to observe the emission of the fluorescence in the atmosphere due to the passage of the particles, giving access to the longitudinal profile of the shower. Four sites with six telescopes each are thus placed at the four corners of the SD and are observing the sky at elevations from 0 • to 30 • . At one of these sites, three other telescopes are constituting the extension of the observations at lower energies, HEAT and are overlooking the infilled region of the SD. They are observing the sky at an higher elevation, from 30 • to 60 • , detecting showers developing at higher altitudes. From these two sets of telescopes, observation of events dominated by fluorescence light or Cherenkov emission lead to the construction of two data sets called hybrid and Cherenkov, respectively.

Hybrid data set
As its name is indicating, the first data set constructed from the observations with the FD is using the hybrid design of the observatory. In an event dominated by the fluorescence light, the shower is seen from the side and a line of triggered pixels is observed on the camera of at least one telescope. The details about the telescopes and its triggering system are presented in Ref. [2]. The known pointing direction in the sky of the triggered pixels is used to determine the shower-detector plane containing the axis of the shower and the location of the considered telescope. Then the timing information of each pixel is used to fit the inclination of the shower axis inside the shower-detector plane. At the energies considered by the hybrid data set, the observed FD events is combined with the detection of the shower by at least one station of the SD on the ground. The timing information and position of the station are then used to improve the constraints on the reconstructed geometry of the shower (location of the core on ground, zenith, and azimuth angles) and thus greatly improving the accuracy of the geometrical reconstruction. With the geometry determined, the next step consists in the reconstruction of the longitudinal profile of the shower to estimate the energy of the cosmic rays. This step requires to take into account the detector efficiency, the fluorescence and Cherenkov yields, and attenuation and scattering in the atmosphere of the emitted photons. The last point is possible thanks to the extensive monitoring of the atmosphere for each night of the measurements. Once this profile is constructed from the data of the telescopes, it is fitted with a Gaisser-Hillas function, which integral gives an estimate of the calorimetric energy deposited in the atmosphere. However, for the total deposited energy, we need to correct the energy (∼15% of the total at 10 18 eV and decreasing to ∼12% at 10 20 eV) for the contribution carried away by muons and neutrinos. The details about the correction of the invisible energy are presented in Ref. [3] but it is worth mentioning that this correction is based on a datadriven analysis.
Selection criteria are applied to guarantee the quality of the spectrum resulting from the use of this data set. Thus only events with θ < 60 • , with good reconstruction of the geometry and the longitudinal profiles are selected. The depth of the maximum development of the shower also has to be observed in the field of view of the detector. Cuts on the atmospheric conditions for which the information on the aerosol content is known, on the energy resolution, etc., complete the selection process [4]. Moreover, given the trigger probability of a single station of the SD and the systematic uncertainty on the exposure, the measurements of the energy spectrum is performed for energies above 10 18 eV. For the reconstruction of the spectrum from the hybrid data, 14675 events are selected in the period from 01/2007 to 12/2017. The resolution of the FD energy for these hybrid events is driven by three contributions: the detector (calibration of the telescopes and reconstruction of the profile), the knowledge of the atmosphere at the time of the observation (molecules and aerosols concentration), and the invisible energy. The results of the study performed in Ref. [5] are reported in Fig. 1. The total resolution of the FD energy is stable with the energy of the cosmic rays at a level of around 7 to 8%.

Cherenkov data set
In contrast to the previous data set, in which the light collected by the telescopes is dominated by the fluorescence emission, the second data set constructed from the FD measurements considers events in which the Cherenkov emission is dominating and pointing directly towards the telescopes. While this allows for an extension of the observations to lower energies, it comes at the cost of having no counterparts from the SD because of the low distance between the FD detector and the location of the core of the shower. Thus the previous hybrid reconstruction procedure cannot be applied anymore for the data recorded at these energies at the HEAT telescopes. The reconstruction of the geometry of the shower is then performed using the Profile-Constrained Geometry Fit (PCGF). All the geometries compatible with the timing information of each triggered pixel are tested and for each one of them the longitudinal profile is reconstructed assuming a Gaisser-Hillas function. The geometry that gives the best expectation of the profile of energy deposit in the atmosphere is then selected.
To decrease the energy threshold of events recorded by the FD and reconstructed with a PCGF method, two data sets have to be distinguished. The first one consists of events for which the full trigger procedure of the telescopes is satisfied. However, this trigger procedure has been designed to reduce the number of events dominated by Cherenkov light [6]. To monitor the evolution of the HEAT telescopes, the Pierre Auger Observatory is randomly storing 10% of events that do not pass the last step of the trigger chain. These events, called minimum-bias events, lead to increase of the statistics of events at the lowest energies considered and thus decrease the energy threshold to 10 15.8 eV and a selection of 123159 events in the period from 06/2012 to 12/2017. As for the hybrid data sets, a selection is applied to the events to minimize the introduction of biases. Cuts on the quality of the reconstruction using the PCGF method, requiring that the maximum of the shower is in the field of view of the telescopes, requiring that the minimum fraction of the Cherenkov light is above 50% etc., are applied. This criteria enable us to reduce the energy bias below 5% and ensure a resolution evolving from 12% at 10 15.8 eV to 6% at 10 18 eV.

Surface Detector
The 1660 detectors of the surface array, whether they are deployed on the triangular grid with a spacing of 1500 m or 750 m, are water-Cherenkov detectors sensitive to both the muonic and electromagnetic components of the showers. With the increase of the inclination of the shower, the electromagnetic component gets more and more absorbed in the atmosphere and asymmetry effects due to the geomagnetic field have to be taken into account to perform an unbiased reconstruction of the characteristics of the cosmic ray initiating the shower. To reconstruct the events observed using such detectors, two techniques have been developed, thus splitting the events in two categories: vertical for which θ < 60 • and inclined for 60 • < θ < 80 • . However this distinction is relevant only for the 1500 m array. The size of the infilled region, 24 km 2 , and the spacing between the detectors, allows a fully triggered detection of EAS only with zenith angles up to 40 • .

Vertical data sets
Covering large arrays at small costs implies a sparse deployment of the detectors, thus the lateral profile of the shower observed on the ground is sampled only by a finite number of detectors for which the time and amplitude of the signal they observe is extracted. While for the reconstruction of the geometry of the shower, i.e. the location of the core on the ground and the arrival direction (zenith and azimuth angles), we exploit the timing at which the stations are triggered by the passage of the particles, the first step of the reconstruction of the energy of the vertical data of the SD1500 and SD750 data corresponds to the reconstruction of the Lateral Distribution Function (LDF) [7,8]. Thus, using the amplitude of the signal of each station participating in an event, the signal at a particular distance from the shower axis is extracted. In the case of the SD1500, this optimal distance is 1000 m and it has been shown in Ref. [9] that the fluctuations of the signals at this distance are minimal. This signal, S (1000), is then the more accurate estimator of the shower size and by extension, of the energy of the primary. Recent studies of the impact of this procedure, especially the choice of one optimal distance, are reported in Ref. [10]. For the SD750, the same procedure is applied, extracting S (450), the signal at 450 m. However, these estimators of the shower size are dependent on the inclination of the shower since with an increase of the inclination, the electromagnetic component is more and more absorbed. To overcome this, the Constant Intensity Cut (CIC) method is applied, assuming that the flux of cosmic rays is isotropic above a defined threshold in energy. Based on that assumption, for which the validity has been scrutinized in Ref. [11], attenuation curves are established at different intensities (∼energies) of cosmic rays as shown in Fig. 2 for SD1500. The estimator of the shower size is reduced by a factor 2 when comparing S (1000) at 0 • and 60 • . According to the model parametrised from these curves, S (1000) and S (450) are corrected in such a way that the energy estimators are defined as S 38 and S 35 , i.e. the signals that would be observed if inclination of the shower was 38 • and 35 • for SD1500 and SD750, respectively. Finally, it is important to mention that these energy estimators are then corrected to take into account effects from the atmospheric conditions at the time of the observations [12] and effects from the geomagnetic field [13]. Thanks to these procedures and the detection methods deployed in the field of the observatory, the SD array is fully efficient for vertical events for energies above 10 18.4 eV with a selection of 215030 events in time period from 01/2004 to 08/2018. For the SD750 array, the extension to the lower energies is first assured by the smaller distance between two adjacent detectors but, on top of that, the use of triggers dedicated to the detection of low-energy events ensure a full efficiency above 10 17 eV [14]. In the period from 01/2014 to 08/2018, a total of 545090 events were selected to reconstruct the spectrum of cosmic rays from the SD750.

Inclined data set
For the more inclined showers, θ > 60 • , the signals observed on the ground are dominated by the muonic component. Moreover, at such inclinations, the circular symmetry around the shower axis is broken by attenuation and geometrical effects and by the influence of the geomagnetic field on the propagation of muons. While the geometry can be reconstructed using the timing of the station, the procedure described in the previous section cannot be applied anymore. For such events, instead of a LDF maps of the density of muons on the ground are built from Monte-Carlo simulations of EAS [15]. Under the assumption that the normalisation of this distribution of muons depends only on the energy of the shower, the scaling parameter N 19 is used as an estimator of the energy of the shower. It corresponds to the scaling factor needed to reconstruct the map of the density of muons for a simulated shower initiated by a cosmic ray with an energy of 10 19 eV to the data. The dependency on the zenith and azimuth angles of the showers is thus incorporated in the density map, leaving the estimator of the energy N 19 free of such dependencies. These methods of detection and reconstruction ensure a full efficiency of SD1500 for the inclined events above 10 18.6 eV with the selection of 24209 events over a period of acquisition from 01/2004 to 08/2018.

Calibration in energy
To convert the estimators of the energy, S 38 , S 35 , and N 19 (for the SD1500-vertical, SD750, and SD1500-inclined, respectively), the hybrid design of the Pierre Auger Observatory is playing a central role, making it possible to perform a calibration in a data-driven way, free of any assumptions about the mass composition of the cosmic rays or about the models of hadronic interactions. For each of the three SD data sets, a subset of events has also been observed and reconstructed by the FD. Applying a set of selection criteria on the quality of the reconstruction for both the SD and FD, on the knowledge of the atmospheric conditions at the time of the observation, etc., the energy calibration of the SD estimators is performed comparing the SD estimators with the calorimetric energy obtained from the FD reconstruction. The calibration of the SD1500 data sets (vertical and inclined) is ensured by a comparison with the standard FD telescopes while the one of the SD750 is also counting on the use of the HEAT data. This procedure is shown in Fig. 3. For each data set, the correlation with the FD energy, is fitted using whereŜ is the considered estimator of the energy and the parameters A and B are left as free parameters. The results of the fit is reported in Fig. 3 by the black lines. The evaluation of the energy resolution of the SD data sets is performed using the same subsets of events seen by both detectors and thus building the distributions of ratio E SD /E FD as a function of energy. This ratio is shown in Fig. 4 for the SD1500-vertical data set at different energies. Knowing the FD energy resolution [5], ∼7 to 8%, the SD energy resolution is estimated decreasing from ∼20% to less than 10% in the energy range of each data set. The bias is inferred using the same ratio for the SD1500 vertical and inclined data sets, while for the SD750 data set a use of MC simulations complete the study at the lowest energies considered [8,11,16].
The systematic uncertainties in the energy scale of the Pierre Auger Observatory are derived the same as for the energy resolution, using hybrid events [5]. The main contributions to the systematic uncertainties are coming from the uncertain knowledge of the atmosphere (evolving with  energy from 3.4% to 6.2%), the calibration of the FD telescopes (9.9%), and the method to reconstruct the longitudinal profile (decreasing with energy from 6.5% to 5.6%). Adding the uncertainties arising from the invisible energy, fluorescence yield, and the total uncertainties in the energy scale results in ∼14%.

The combined spectrum and a new feature
For each of the five data sets previously described, the raw spectrum of the cosmic rays is obtained by dividing the rate of events in a particular bin of energy by the width of this bin and its corresponding exposure.

Exposure
In the case of the three data sets constructed from the SD, the computation of the exposure is reduced to a geometrical problem. Indeed, at energies above the full-efficiency threshold, the total exposure is calculated by integrating in time the aperture of all working cells of the surface array, i.e. hexagonal cells in the case of the triangular grid. This calculation is possible using the monitoring database of the SD, and because for each event in the SD1500-vertical,inclined, and SD750 data set, it is required that the six stations around the station recording the highest amplitude of the signal are active at the time of the detection. The exposure, computed in such a way, is then independent of any MC assumptions and of the energy. Its uncertainties are kept below 3%. The results are presented in Fig. 5. Thanks to its large area of detection and its more than 15 years of data taking, the Pierre Auger Observatory has reached exposures of 60400 km 2 sr yr, 17500 km 2 sr yr, and 105 km 2 sr yr for the SD1500-vertical, -inclined, and SD750 data sets, respectively. The computation of the exposure for the FD data sets is less straightforward and has to be performed simulating a large number of MC events that are taking into account the status of the atmosphere and of the telescopes [17].  Due to the constraints inherent in the use of MC simulations, the mass composition chosen is following the one derived from the Global Spline Fit model [18] and simulations assuming pure composition are then used to derive the associated uncertainties related to the lack of knowledge of the true composition reaching Earth at the highest energies. The resulting exposure is thus dependent on the energy of the cosmic rays and is ∼2200 km 2 sr yr at 10 19 eV and ∼10 km 2 sr yr at 10 17.5 eV for the hybrid and Cherenkov data sets, respectively. The uncertainties arising from an exposure computed using MC simulations are described in Ref. [6].

Five spectra
The raw spectra are suffering from the distortions induced by the finite energy resolution and the shape of the spectrum, causing a bin-to-bin migration of events. Taking into account the detector effects and assuming a spectrum following a multiply-broken power law, a correction factor is established as a function of the energy and applied to the raw spectra. The resulting spectra are shown in Fig. 6 after normalizing (independently of the energy) the spectra by <1%, +5%, −2%, <1%, and +7% for the SD1500- More details can be found in Ref. [11] vertical, -inclined, SD750, hybrid, and Cherenkov spectra, respectively. These normalization values are obtained applying the combination method that will be described later in this section. The discrepancies between all the different reconstructed spectra is well encompassed by the uncertainties in the energy scale (∼15% for the Cherenkov data set, ∼14% for the other data sets). Moreover, it is interesting to note that the uncertainty in the energy scale is the main source of uncertainties in the reconstruction of the spectra observed by the Pierre Auger Observatory. In Fig. 7, the analysis of all possible sources of uncertainties in the spectrum is reported for the SD1500-vertical data set [11]. Systematic uncertainties coming from the reconstruction of the energy estimator, computation of the exposure, or forward-folding procedure are below 3% and thus well below the uncertainties in the energy scale. The domination of the uncertainties in the energy scale in the total uncertainty of the spectrum is found also for all the other data sets used [4,6,8,16]. It is worth noting that at the lowest energies of the spectrum reconstructed from the Cherenkov data set, the uncertainty are dominated by the uncertainties in the exposure due to the assumptions made with the use of simulations.
The combination of all the five spectra is performed using a forward-folding approach, beginning with the evaluation of the uncorrelated uncertainties between the different data sets, a shift in exposure and in the energy calibration is introduced to modify a given spectrum relative to the others. For each spectrum, a bin-to-bin migration of events is taken into account considering the reconstruction biases and resolutions of each data set. The shape of the expected spectrum is assumed as a sequence of six power-laws written as where J 0 is the normalization parameter, γ i are the spectral indices, and E i j the energies at which the transitions occur. The width of the transition between two power-law regions is set by the parameters ω i j fixed as ω 01 = ω 12 = 0.25 and ω 23 = ω 34 = 0.05 (half of the width of one logarithmic bin in energy). Thus a combined likelihood, written as the product of Poissonian terms, is applied to determine for each spectrum the best shift (constrained by a Gaussian penalization included in the likelihood) in exposure and in energy and the best set of parameters is selected to get the combined spectrum shown in Fig. 8. The combined spectrum is fitted using Eq. (2) and the best-fit parameters are reported in Tab. 1 with the transition width fixed as mentioned previously. This spectrum is now covering five orders of magnitude in energy, from 6×10 15 eV to the highest energies ever observed. Positions of the low-energy ankle at 2.8×10 16 eV, firstly reported in Ref. [19], of the second knee at 1.58×10 17 eV, of the ankle at 5×10 18 eV and of the suppression at 4.7×10 19 eV are thus reported. With the increase of the number of events collected by the Pierre Auger Observatory, a new feature has recently been reported at an energy of 1.4×10 19 eV and is named in-step [11]. The interpretation of all these features is now under scrutiny to understand the mechanisms behind the acceleration of the cosmic rays. Such interpretation is presented in Ref. [20] using the data about the mass composition of the cosmic rays. For example, in the scenario explored, the in-step feature appears as a transition between helium-like nuclei and CNO-like nuclei.

Conclusion: a glimpse into the future
The continuous observation of UHECR by the Pierre Auger Observatory in the last 18 years allows the reconstruction of the energy spectrum with an unprecedented precision. However, the story of the energy spectrum is not planned to end soon. Associated to its longevity, the Observatory is constantly improved. To reinforce the measurements at the lowest energies, an array of 19 water-Cherenkov detectors has been deployed in the region of the SD750 with a spacing between two adjacent detectors reduced to 433 m [21]. It is then expected, from its measurements, to build a sixth energy spectrum of cosmic rays starting at 6×10 16 eV for vertical events (θ < 45 • ).
Such scenarios are already tested combining the energy spectrum as presented in this proceeding with the mass composition extracted from the measurement of the shower maximum. A recent update of this analysis performed by the Pierre Auger Collaboration has been reported in Ref. [23]. The transition between Galactic and extra-galactic cosmic rays is thus scrutinized, as well as the end of the spectrum interpreted as a combination of a rigidity-dependent acceleration and propagation effects.
Finally it is not possible to track back the origins of the UHECRs without analyzing the arrival directions of the cosmic rays. Anisotropy studies at different scales have been conducted and led to the discovery of a dipole above 8 EeV and an indication of the Starburst galaxies as a possible sources of UHECRs. The combination of all the characteristics of the cosmic rays: energy, mass, and arrival direction is constituting the next step in solving the question of the origins of such particles and is already under way [24].