A measurement of the cosmic-ray energy spectrum above $2.5{\times} 10^{18}$ eV using the Pierre Auger Observatory

We report a measurement of the energy spectrum of cosmic rays for energies above $2.5 {\times} 10^{18}~$eV based on 215,030 events recorded with zenith angles below $60^\circ$. A key feature of the work is that the estimates of the energies are independent of assumptions about the unknown hadronic physics or of the primary mass composition. The measurement is the most precise made hitherto with the accumulated exposure being so large that the measurements of the flux are dominated by systematic uncertainties except at energies above $5 {\times} 10^{19}~$eV. The principal conclusions are: (1) The flattening of the spectrum near $5 {\times} 10^{18}~$eV, the so-called"ankle", is confirmed. (2) The steepening of the spectrum at around $5 {\times} 10^{19}~$eV is confirmed. (3) A new feature has been identified in the spectrum: in the region above the ankle the spectral index $\gamma$ of the particle flux ($\propto E^{-\gamma}$) changes from $2.51 \pm 0.03~{\rm (stat.)} \pm 0.05~{\rm (sys.)}$ to $3.05 \pm 0.05~{\rm (stat.)} \pm 0.10~{\rm (sys.)}$ before changing sharply to $5.1 \pm 0.3~{\rm (stat.)} \pm 0.1~{\rm (sys.)}$ above $5 {\times} 10^{19}~$eV. (4) No evidence for any dependence of the spectrum on declination has been found other than a mild excess from the Southern Hemisphere that is consistent with the anisotropy observed above $8 {\times} 10^{18}~$eV.

We report a measurement of the energy spectrum of cosmic rays for energies above 2.5×10 18 eV based on 215,030 events recorded with zenith angles below 60 • . A key feature of the work is that the estimates of the energies are independent of assumptions about the unknown hadronic physics or of the primary mass composition. The measurement is the most precise made hitherto with the accumulated exposure being so large that the measurements of the flux are dominated by systematic uncertainties except at energies above 5×10 19 eV. The principal conclusions are: 3. A new feature has been identified in the spectrum: in the region above the ankle the spectral index γ of the particle flux (∝ E −γ ) changes from 2.51 ± 0.03 (stat.) ± 0.05 (sys.) to 3.05 ± 0.05 (stat.) ± 0.10 (sys.) before changing sharply to 5.1 ± 0.3 (stat.) ± 0.1 (sys.) above 5×10 19 eV. 4. No evidence for any dependence of the spectrum on declination has been found other than a mild excess from the Southern Hemisphere that is consistent with the anisotropy observed above 8×10 18 eV.

I. INTRODUCTION
Although the first cosmic rays having energies above 10 19 eV were detected nearly 60 years ago [1,2], the question of their origin remains unanswered. In this paper we report a measurement of the energy spectrum of ultra-high energy cosmic rays (UHECRs) of unprecedented precision using data from the Pierre Auger Observatory. Accurate knowledge of the cosmic-ray flux as a function of energy is required to help discriminate between competing models of cosmic-ray origin. As a result of earlier work with the HiRes instrument [3], the Pierre Auger Observatory [4] and the Telescope Array [5], two spectral features were identified beyond reasonable doubt (see, e.g., [6][7][8][9][10][11] for recent reviews). These are a hardening of the spectrum at about 5×10 18 eV (the ankle) and a strong suppression of the flux at an energy about a decade higher. The results reported here are based on 215,030 events with energies above 2.5×10 18 eV. The present measurement, together with recent observations of anisotropies in the arrival directions of cosmic rays on large angular scales above 8×10 18 eV [12] and on intermediate angular scales above 3.9×10 19 eV [13], and inferences on the mass composition [14,15], provide essential data against which to test phenomenological models of cosmic-ray origin. As part of a broad study of directional anisotropies, the large number of events used in the present analysis allows examination of the energy spectrum as a function of declination as reported below.
The determination of the flux of cosmic rays is a nontrivial exercise at any energy. It has long been recognised that 70 to 80% of the energy carried by the primary particle is dissipated in the atmosphere through ionisation loss and thus, with ground detectors alone, one must resort to models of shower development to infer the primary energy. This is difficult as a quantitative knowledge of hadronic processes in the cascade is required. While at about 10 17 eV the centre-of-mass energies encountered in collisions of primary cosmic rays with air nuclei are comparable to those achieved at the Large Hadron Collider, details of the interactions of pions, which are key to the development of the cascade, are lacking, and the presence of unknown processes is also possible. Furthermore one has to make an assumption about the primary mass. Both conjectures lead to systematic uncertainties * Electronic address: auger˙spokespersons@fnal.gov; URL: http:// www.auger.org that are difficult, if not impossible, to assess. To counter these issues, methods using light produced by showers as they cross the atmosphere have been developed. In principle, this allows a calorimetric estimate of the energy. Pioneering work in the USSR in the 1950s [16] led to the use of Cherenkov radiation for this purpose, and this approach has been successfully adopted at the Tunka [17] and Yakutsk [18] arrays. The detection of fluorescence radiation, first achieved in Japan [19] and, slightly later, in the USA [20], has been exploited particularly effectively in the Fly's Eye and HiRes projects to achieve the same objective. The Cherenkov method is less useful at the highest energies as the forward-beaming of the light necessitates the deployment of a large number of detectors while the isotropic emission of the fluorescence radiation enables showers to be observed at distances of 30 km from a single station. For both methods, the on-time is limited to moonless nights, and an accurate understanding of the aerosol content of the atmosphere is needed.
The Pierre Auger Collaboration introduced the concept of a hybrid observatory in which the bulk of the events used for spectrum determination is obtained with an array of detectors deployed on the ground and the integral of the longitudinal profile, measured using a fluorescence detector, is used to calibrate a shower-size estimate made with the ground array. This hybrid approach has led to a substantial improvement in the accuracy of reconstruction of fluorescence events and to a calorimetric estimate of the energy of the primary particles for events recorded during periods when the fluorescence detector cannot be operated. The hybrid approach has also been adopted by the Telescope Array Collaboration [5].
A consistent aim of the Auger Collaboration has been to make the derivation of the energy spectrum as free of assumptions about hadronic physics and the primary composition as possible. The extent to which this has been achieved can be judged from the details set out below. After a brief introduction in Sec. II to relevant features of the Observatory and the data-set, the method of estimation of energy is discussed in Sec. III. In Sec. IV, the approach to deriving the energy spectrum is described, including the procedure for evaluating the exposure and for unfolding the resolution effects, as well as a detailed discussion of the associated uncertainties and of the main spectral features. A search for any dependence of the energy spectrum on declination is discussed in Sec. V, while a comparison with previous works is given in Sec. VI. The results from the measurement of the energy spectrum are summarized in the concluding Sec. VII.

II. THE PIERRE AUGER OBSERVATORY AND THE DATA SETS
A. The Observatory The Pierre Auger Observatory is sited close to the city of Malargüe, Argentina, at a latitude of 35.2 • S with a mean atmospheric overburden of 875 g/cm 2 . A detailed description of the instrument has been published [21], and only brief remarks concerning features relevant to the data discussed in this paper are given.
The surface detector (SD) array comprises about 1600 water-Cherenkov detectors laid out on a 1500 m triangular grid, covering an area of about 3000 km 2 . Each SD has a surface area of 10 m 2 and a height of 1.2 m, holding 12 tonnes of ultra-pure water viewed by 3 ×9" photomultipliers (PMTs). The signals from the PMTs are digitised using 40 MHz 10-bit Flash Analog to Digital Converters (FADCs). Data collection is achieved in real time by searching on-line for temporal and spatial coincidences at a minimum of three locations. When this occurs, FADC data from the PMTs are acquired from which the pulse amplitude and time of detection of signals is obtained. The SD array is operated with a duty cycle close to 100%.
The array is over-looked from four locations, each having six Schmidt telescopes designed to detect fluorescence light emitted from shower excitations of atmospheric nitrogen. In each telescope, a camera with 440 hexagonal PMTs is used to collect light from a 13 m 2 mirror. These instruments, which form the fluorescence detector (FD), are operated on clear nights with low background illumination with an on-time of 15%.
Atmospheric conditions at the site of the Observatory must be known for the reconstruction of the showers. Accordingly, comprehensive monitoring of the atmosphere, particularly of the aerosol content and the cloud cover, is undertaken as described in [21]. Weather stations are located close to the sites of the fluorescence telescopes. Before the Global Data Assimilation system was adopted [22], an extensive series of balloon flights was made to measure the humidity, temperature and pressure in the atmosphere as a function of altitude.

B. The data sets
The data set used for the measurement of the energy spectrum consists of extensive air showers (EAS) recorded by the SD array. EAS detected simultaneously by the SD and the FD play a key role in this work. Dubbed hybrid events, they are pivotal in the determination of the energy of the much more numerous SD events [23]. We use here SD events with zenith angle θ < 60 • , as the reconstruction of showers at larger angles requires a different method due to an asymmetry induced in the distribution of the shower particles by the geomagnetic field and geometrical effects (see [24]). A brief description of the reconstruction of SD and hybrid events is given in [25]: a more detailed description is in [26]. We outline here features relevant to the present analysis.
The reconstruction of the SD events is used to determine the EAS geometry (impact point of the shower axis and arrival direction) as well as a shower-size estimator. To achieve this, the amplitude and the start-time of the signals, recorded at individual SD stations and quantified in terms of their response to a muon travelling vertically and centrally through it (a vertical equivalent muon or VEM), are used. The arrival direction is determined to about 1 • from the relative arrival times of these signals. The impact point and the shower-size estimator are in turn derived by fitting the signal amplitudes to a lateral distribution function (LDF) that decreases monotonically with distance from the shower axis. The shower-size estimator adopted is the signal at 1000 m from the axis, S(1000). For the grid spacing of 1500 m, 1000 m is the optimal distance to minimize the uncertainties of the signal due to the imperfect knowledge of the functional form of the LDF in individual events [27]. The combined statistical and systematic uncertainty decreases from 15% at a shower size of 10 VEM to 5% at the highest shower sizes. The uncertainty on the impact point is of order 50 m. S(1000) is influenced by changes in atmospheric conditions that affect shower development [28], and by the geomagnetic field that impacts on the shower particle-density [29]. Therefore, before using the shower-size estimator in the calibration procedure (Sec. III), corrections of order 2% and 1% for the atmospheric and geomagnetic effects, respectively, are made.
For the analysis in this paper, the SD reconstruction is carried out only for events in which the detector with the highest signal is surrounded by a hexagon of six stations that are fully operational. This requirement not only ensures adequate sampling of the shower but also allows evaluation of the aperture of the SD in a purely geometrical manner in the regime where the array is fully efficient [30]. As shown in Sec. IV, such a regime is attained for events with θ < 60 • at an energy 2.5×10 18 eV. With these selection criteria, the SD data set used below consists of 215,030 events recorded between 1 January 2004 and 31 August 2018.
For hybrid events the reconstruction procedure exploits the amplitude and timing of the signals detected by each PMT in each telescope as well as additional timing information from the SD station with the highest signal. Combining the timing information from FD and SD improves the directional precision to 0.6 • [21]. Hybrid reconstruction provides in addition the longitudinal profile from which the depth of the shower maximum (X max ) and the primary energy are extracted. The light signals in the FD PMTs are converted to the energy deposited in sequential depths in the atmosphere, taking into account the fluorescence and Cherenkov light contributions [31] and their attenuation due to scattering. The longitudinal profile of the energy deposit is reconstructed by means of a fit to a modified Gaisser-Hillas profile [32].
Integration of the longitudinal profile yields a calorimetric measure of the ionisation loss in the atmosphere which is supplemented by the addition of the undetected energy, or "invisible energy", carried into the ground by muons and neutrinos. We denote the sum of these two contributions, our estimate of the energy carried by the incoming primary particle, as E FD . The invisibleenergy correction is estimated with a data-driven analysis and is about 14% at 2.5×10 18 eV falling to about 12% at 10 20 eV [33]. The resolution of E FD is 7.4% at 2.5×10 18 eV and worsens with energy to 8.6% at 6×10 19 eV. It is obtained by taking into account all uncorrelated uncertainties between different showers. In addition to the statistical uncertainty arising from the fit to the longitudinal profile, this resolution includes uncertainties in the detector response, in the models of the state of the atmosphere, and in the expected fluctuations from the invisible energy which, parameterized as a function of the calorimetric energy, is assumed to be identical for any primary of same energy. All the uncorrelated uncertainties are addressed in [34] with further details given in [31]. We note that at higher energies the showers are detected, on average, at larger distances from the FD telescopes because the detection and reconstruction efficiency at larger distances increases with energy. This causes a worsening of the energy resolution because of the interplay between the uncertainty from the aerosols increasing with energy and the uncertainty from photoelectrons decreasing with energy.
The hybrid trigger efficiency, i.e. the probability of detecting a fluorescence event in coincidence with at least one triggered SD station, is 100% at energies greater than 10 18 eV, independent of the mass of the nuclear primaries [35]. The hybrid data set used for the calibration of the SD events comprises 3,338 events with E > 3×10 18 eV collected between 1 January 2004 and 31 December 2017. Other criteria for event selection are detailed in Sec. III.

III. ENERGY ESTIMATION FROM EVENTS RECORDED BY THE SURFACE ARRAY
The energy calibration of the SD shower-size estimator against the energy derived from measurements with the FD is a two-step process. For a cosmic ray of a given energy, the value of S(1000) depends on zenith angle because of the different atmospheric depths crossed by the corresponding shower. As detailed in Sec. III A, we first correct for such an attenuation effect by using the Constant Intensity Cut (CIC) method [36]. The calibration is then made between the corrected shower-size estima-tor, denoted by S 38 , and the energy measured by the FD in hybrid events, E FD : the procedure to obtain the SD energy, E SD , is explained in Sec. III B. The systematic uncertainties associated with the SD energy scale thus obtained are described in Sec. III C. Finally, the estimation of E SD from E FD allows us to derive the resolution, σ SD (E), as well as the bias, b SD (E), down to energies below which the detector is not fully efficient. We explain in Sec. III D the method used to measure b SD (E) and σ SD (E), from which we build the resolution function for the SD to be used for the unfolding of the spectrum. 26 For a fixed energy, S(1000) depends on the zenith angle θ because, once it has passed the depth of shower maximum, a shower is attenuated as it traverses the atmosphere. The intensity of cosmic rays, defined here as the number of events per steradian above some S(1000) threshold, is thus dependent on zenith angle as can be seen in Fig. 1.
Given the highly isotropic flux, the intensity is expected to be θ-independent after correction for the attenuation. Deviations from a constant behavior can thus be interpreted as being due to attenuation alone. Based on this principle, an empirical procedure, the so-called CIC method, is used to determine the attenuation curve as function of θ and therefore a θ-independent showersize estimator (S 38 ). It can be thought of as being the S(1000) that a shower would have produced had it arrived at 38 • , the median angle from the zenith. The small anisotropies in the arrival directions and the zenithal dependence of the resolution on S 38 do not alter the validity of the CIC method in the energy range considered here, as shown in Appendix A.
In practice, a histogram of the data is first built in cos 2 θ to ensure equal exposure; then the events are ordered by S(1000) in each bin. For an intensity high enough to guarantee full efficiency, the set of S(1000) values, each corresponding to the N th largest signal in the associated cos 2 θ bin, provides an empirical estimate of the attenuation curve. Because the mass of each cosmicray particle cannot be determined on an event-by-event basis, the attenuation curve inferred in this way is an effective one, given the different species that contribute at each intensity threshold. The resulting data points are fitted with a third-degree polynomial, S(1000) = S 38 (1 + ax + bx 2 + cx 3 ), where x = cos 2 θ − cos 2 38 • . Fits are shown in the top panel of Fig. 2 for three different intensity thresholds corresponding to I 1 = 2.91×10 4 sr −1 , I 2 = 4.56×10 3 sr −1 and I 3 = 6.46×10 2 sr −1 at 38 • . The attenuation is plotted as a function of sec θ to exθ sec 1 1. hibit the dependence on the thickness of atmosphere traversed. The uncertainties in each data point follow from the number of events above the selected S(1000) values. The N th largest signal in each bin is a realization of a random variable distributed as an order-statistic variable where the total number of ordered events in the cos 2 θ bin is itself a Poisson random variable. Within a precision better than 1%, the standard deviation of the random variable can be approximated through a straightforward Poisson propagation of uncertainties, namely ∆S(N ) The number of bins is adapted to the available number of events for each intensity threshold, from 27 for I 1 so as to guarantee a resolution on the number of events of 1% in each bin, to 8 for I 3 so as to guarantee a resolution of 4%.
The curves shown in Fig. 2 are largely shaped by the electromagnetic contribution to S(1000) which, once the shower development has passed its maximum, decreases with the zenith angle because of attenuation in the increased thickness of atmosphere. The muonic component starts to dominate at large angles, which explains the flattening of the curves. In the bottom panel, the curves are normalized to 1 at 38 • to exhibit the differences for the selected intensity thresholds. Some dependence with the intensity thresholds, and thus with the energy thresholds, is observed at high zenith angles: high-energy showers appear more attenuated than low energy ones. This results from the interplay between the mass composition and the muonic-to-electromagnetic signal ratio at ground level. A comprehensive interpretation of these curves is however not addressed here.
The energy dependence in the CIC curves that is observed is accounted for by introducing an empirical dependence in terms of y = log 10 (S 38 /40 VEM) in the coefficients a, b and c through a second-order polynomial in y. The polynomial coefficients derived are shown in Table I. They relate to S 38 values ranging from 15 VEM to 120 VEM. Outside these bounds, the coefficients are set to their values at 15 and 120 VEM. This is because below 15 VEM, the isotropy is not expected anymore due to the decreasing efficiency, while above 120 VEM, the number of events is low and there is the possibility of localized anisotropies. The shower-size estimator, S 38 , is converted into energy through a calibration with E FD by making use of a subset of SD events, selected as described in Sec. II, which have triggered the FD independently. For the analysis, we apply several selection criteria to guarantee a precise estimation of E FD as well as fiducial cuts to minimise the biases in the mass distribution of the cosmic rays introduced by the field of view of the FD telescopes.
The first set of cuts aims to select time periods during which data-taking and atmospheric conditions are suitable for collecting high-quality data [37]. We require a high-quality calibration of the gains of the PMTs of the FD and that the vertical aerosol optical depth is measured within 1 hour of the time of the event, with its value integrated up to 3 km above the ground being less than 0.1. Moreover, measurements from detectors installed at the Observatory to monitor atmospheric conditions [21] are used to select only those events detected by telescopes without clouds within their fields of view. Next, a set of quality cuts are applied to ensure a precise reconstruction of the energy deposit [37].
We select events with a total track length of at least 200 g/cm 2 , requiring that any gap in the profile of the deposited energy be less than 20% of the total track length and we reject events with an uncertainty in the reconstructed calorimetric energy larger than 20%. We transform the χ 2 into a variable with zero mean and unit variance, z = χ 2 − n dof / √ 2n dof with n dof the number of degrees of freedom, and require that the z values be less than 3. Finally, the fiducial cuts are defined by an appropriate selection of the lower and upper depth boundaries to enclose the bulk of the X max distribution and by requiring that the maximum accepted uncertainty in X max is 40 g/cm 2 and that the minimum viewing angle of light in the telescope is 20 • [37]. This limit is set to reduce contamination by Cherenkov radiation. A final cut is applied to E FD : it must be greater than 3×10 18 eV to ensure that the SD is operating in the regime of full efficiency (see Sec. IV A).
After applying these cuts, a data set of 3,338 hybrid events is available for the calibration process. With the current sensitivity of our X max measurements in this energy range, a constant elongation rate (that is, a single logarithmic dependence of X max with energy) is observed [37]. In this case, a single power law dependence of S 38 with energy is expected from Monte-Carlo simulations. We thus describe the correlation between S 38 and E FD , shown in Fig. 3, by a power law function, where A and B are fitted to data. In this manner the correlation captured through this power-law relationship is fairly averaged over the underlying mass distribution, and thus provides the calibration of the mass-dependent S 38 parameter in terms of energy in an unbiased way over the covered energy range. Due to the limited number of events in the FD data set at the highest energies, deviations from the inferred power law cannot be fully investigated currently. We note however that any indication for a strong change of elongation rate cannot be inferred at the highest energies from our SD-based indirect measurement reported in [15].
[eV]  Figure 3: Correlation between the SD shower-size estimator, S38, and the reconstructed FD energy, EFD, for the selected 3,338 hybrid events used in the fit. The uncertainties indicated by the error bars are described in the text. The solid line is the best fit of the power-law dependence EFD = A S38 B to the data. The reduced deviance of the fit, whose calculation is detailed in Appendix B, is shown in the bottom-right corner.
The correlation fit is carried out using a tailored maximum-likelihood method allowing various effects of experimental origin to be taken into account [38]. The probability density function entering the likelihood procedure, detailed in Appendix B, is built by folding the cosmic-ray flux, observed with the effective aperture of the FD, with the resolution functions of the FD and of the SD. Note that to avoid the need to model accurately the cosmic-ray flux observed through the effective aperture of the telescopes (and thus to rely on mass assumptions), the observed distribution of events passing the cuts described above is used in this probability density function.
The uncertainties in the FD energies are estimated, on an event-by-event basis, by adding in quadrature all uncertainties in the FD energy measurement which are uncorrelated shower-by-shower (see [34] for details). The uncertainties in S 38 are also estimated on an event-byevent basis considering the event-by-event contribution arising from the reconstruction accuracy of S(1000). The error arising from the determination of the zenith angle is negligible. The contribution from shower-to-shower fluctuations to the uncertainty in E SD is parameterized as a relative error in S 38 with 0.13 − 0.08x + 0.03x 2 where x = log 10 (E/eV) − 18.5. It is obtained by subtracting in quadrature the contribution of the uncertainty in S 38 from the SD energy resolution. The latter, as detailed in the following, is measured from data and the resulting shower-to-shower fluctuations are free from any reliance on mass assumption and model simulations.
The best fit parameters are A = (1.86±0.03)×10 17 eV and B = 1.031±0.004 and the correlation coefficient between the parameters is ρ = −0.98. The resulting calibration curve is shown as the red line in Fig. 3. The goodness of the fit is provided by the value of the reduced deviance, namely D/n dof = 3419/3336. The statistical uncertainty on the SD energies obtained propagating the fit errors on A and B is 0.4 % at 3×10 18 eV, increasing up to 1% at the highest energies. The most energetic event used in the calibration is detected at all four fluorescence sites. Its energy is (8.5±0.4)×10 19 eV, obtained from a weighted average of the four calorimetric energies and using the resulting energy to evaluate the invisible energy correction [33]. It has a depth of shower maximum of (763±8) g/cm 2 , which is typical/close to the average for a shower of this energy [37]. The energy estimated from S 38 = 354 VEM is (7.9±0.6)×10 19 eV.

C. ESD: systematic uncertainties
The calibration constants A and B are used to estimate the energy for the bulk of SD events: E SD ≡ AS 38 B . They define the SD energy scale. The uncertainties in the FD energies are estimated, on an event-by-event basis, by adding in quadrature all uncertainties in the FD energy measurement which are correlated shower-by-shower [23].
The contribution from the fluorescence yield is 3.6% and is obtained by propagating the uncertainties in the high-precision measurement performed in the AIRFLY experiment of the absolute yield [39] and of the wavelength spectrum and quenching parameters [40,41]. The uncertainty coming from the characterization of the atmosphere ranges from 3.4% (low energies) to 6.2% (high energies). It is dominated by the uncertainty associated with the aerosols in the atmosphere and includes a minor contribution related to the molecular properties of the atmosphere. The largest correlated uncertainty, associated with the calibration of the FD, amounts to 9.9%. It includes a 9% uncertainty in the absolute calibration of the telescopes and other minor contributions related to the relative response of the telescopes at different wavelengths and relative changes with time of the gain of the PMTs. The uncertainty in the reconstruction of the energy deposit ranges from 6.5% to 5.6% (decreasing with energy) and accounts for the uncertainty associated with the modelling of the light spread away from the image axis and with the extrapolation of the modified Gaisser-Hillas profile beyond the field of view of the telescopes. The uncertainty associated with the invisible energy is 1.5%. The invisible energy is inferred from data through an analysis that exploits the sensitivity of the water-Cherenkov detectors to muons and minimizes the uncertainties related to the assumptions on hadronic interaction models and mass composition [33].
We have performed several tests aimed at assessing the robustness of the analysis that returns the calibration coefficients A and B. The correlation fit was repeated selecting events in three different zenithal ranges. The obtained calibration parameters are reported in Table II. The calibration curves are within one standard deviation of the average one reported above, resulting in energies within 1% of the average ones. Other tests performed using looser selection criteria for the FD events give similar results. By contrast, determining the energy scale in different time periods leads to some deviation of the calibration curves with respect to the average one. Although such variations are partly accounted for in the FD calibration uncertainties, we conservatively propagate these uncertainties into a 5% uncertainty on the SD energy scale.
The total systematic uncertainty in the energy scale is obtained by adding in quadrature all of the uncertainties detailed above, together with the contribution arising from the statistical uncertainty in the calibration parameters. The total is about 14% and it is almost energy independent as a consequence of the energy independence of the uncertainty in the FD calibration, which makes the dominant contribution.

D. ESD: resolution and bias
Our final aim is to estimate the energy spectrum above 2.5×10 18 eV. Still it is important to characterize the energies below this threshold because the finite resolution on the energies induces bin-to-bin migration effects that affect the spectrum. In this energy range, below full efficiency of the SD, systematic effects enter into play on the energy estimate. While the FD quality and fiducial cuts still guarantee the detection of showers without bias towards one particular mass in that energy range, this is no longer the case for the SD due to the higher efficiency of shower detection for heavier primary nuclei [30]. Hence the distribution of S 38 below 3×10 18 eV may no longer be fairly averaged over the underlying mass distribution, and a bias on E SD may result from the extrapolation of the calibration procedure, in addition to the trigger effects that favor positive fluctuations of S 38 at a fixed energy over negative ones. In this section, we determine these quantities, denoted as σ SD (E, θ)/E for the resolution and as b SD (E, θ) for the bias, in a data-driven way. These measurements allow us to characterize the SD resolution function that will be used in several steps of the analysis presented in the next sections. This, denoted as κ(E SD |E; θ), is the conditional p.d.f. for the measured energy E SD given that the true value is E. It is normal-ized such that the event is observed at any reconstructed energy, that is, dE SD κ(E SD |E; θ) = 1. In the energy range of interest, we adopt a Gaussian curve, namely: .  Figure 4: Ratio distribution of the SD energy, ESD, to the FD energy, EFD, from the selected data sample, for three energy ranges. The distributions are all normalized to unity to better underline the difference in their shape. The total number of events for each distribution is 2367, 1261 and 186 from the lower to the higher energy bin, respectively.
The estimation of b SD (E, θ) and σ SD (E, θ) is obtained by analyzing the E SD /E FD histograms as a function of E FD , extending here the E FD range down to 10 18 eV. For Gaussian-distributed E FD and E SD variables, the E SD /E FD variable follows a Gaussian ratio distribution. For a FD resolution function with no bias and a known resolution parameter, the searched b SD (E, θ) and σ SD (E, θ) are then obtained from the data. The overall FD energy resolution is σ FD (E)/E 7.4%. In comparison to the number reported in Sec. II B, σ FD (E)/E is here almost constant over the whole energy range because it takes into account that, at the highest energies, the same shower is detected from different FD sites. In these cases, the energy used in analyses is the mean of the reconstructed energies (weighted by uncertainties) from the two (or more) measurements. This accounts for the improvement in the statistical error.
Examples of measured and fitted distributions of E SD /E FD are shown in Fig. 4 for three energy ranges: the resulting SD energy resolution is σ SD (E)/E = (21.  as a function of E: the resolution is 20% at 2×10 18 eV and tends smoothly to 7% above 2×10 19 eV. Note that no significant zenithal dependence has been observed. The bias parameter b SD (E, θ) is illustrated in Fig. 6 as a function of the zenith angle for four different energy ranges. The net result of the analysis is a bias larger than 10% at 10 18 eV, going smoothly to zero in the regime of full efficiency. Note that the selection effects inherent in the FD field of view induce different samplings of hybrid and SD showers with respect to shower age at a fixed zenith angle and at a fixed energy. These selection cuts also impact the zenithal distribution of the showers. Potentially, the hybrid sample may thus not be a fair sample of the bulk of SD events. This may lead to some misestimation of the SD resolution determined in the data-driven manner presented above. We have checked, using end-to-end Monte-Carlo simulations of the Observatory operating in the hybrid mode, that the particular quality and fiducial cuts used to select the hybrid sample do not introduce significant distortions to the measurements of σ SD (E) shown in deviations of the reconstructed energy histograms remain within 10% (low energies) and 5% (high energies) whatever the assumption on the mass composition. There is thus a considerable benefit in relying on the hybrid measurements,to avoid any reliance on mass assumptions when determining the bias and resolution factors.
From the measurements, a convenient parameterization of the resolution is where the values of the parameters are obtained from a fit to the data: σ 0 = 0.078, σ 1 = 0.16, and E σ = 6.6×10 18 eV. The function and its statistical uncertainty from the fit are shown in Fig. 5. It is worth noting that this parameterization accounts for both the detector resolution and the shower-to-shower fluctuations. Finally, a detailed study of the systematic uncertainties on this parameterization leads to an overall relative uncertainty of about 10% at 10 18 eV and increasing with energy to about 17% at the highest energies. It accounts for the selection effects inherent to the FD field of view previously addressed, the uncertainty in the FD resolution and the statistical uncertainty in the fitted parameterization.
The bias, also parameterized as a function of the energy, includes an additional angular dependence: for log 10 (E/eV) ≤ log 10 (E * /eV) = 18.4, and b SD = 0 otherwise. Here, b 0 = 0.20, b 1 = 0.59 and λ b = 10.0. The parameters are obtained in a two steps process: we first perform a fit to extract the zenith-angle dependence in different energy intervals prior to determining the energy dependence of the parameters. Examples of the results of the fit to the data are shown in Fig. 6. The relative uncertainty in these parameters is estimated to be within 15%, considering the largest uncertainties of the data points displayed in the figure. This is a conservative estimate compared to that obtained from the fit, but this enables us to account for systematic changes that would have occurred had we chosen another functional shape for the parameterization.
The two parameterizations of equations (3) and (4) are sufficient to characterize the Gaussian resolution function of the SD in the energy range discussed here.

IV. DETERMINATION OF THE ENERGY SPECTRUM
In this section, we describe the measurement of the energy spectrum, J(E). Over parts of the energy range, we will describe it using J(E) ∝ E −γ , where γ is the spectral index. In Sec. IV A, we present the initial estimate of the energy spectrum, dubbed the "raw spectrum", after explaining how we determine the SD efficiency, the exposure and the energy threshold for the measurement. In Sec. IV B, we describe the procedure used to correct the raw spectrum for detector effects, which also allows us to infer the spectral characteristics. The study of potential systematic effects is summarised in Sec. IV C, prior to a discussion of the features of the spectrum in Sec. IV D.
A. The raw spectrum An initial estimation of the differential energy spectrum is made by counting the number of observed events, N i , in differential bins (centered at energy E i , with width ∆E i ) and dividing by the exposure of the array, E, The bin sizes, ∆E i , are selected to be of equal size in the logarithm of the energy, such that the width, ∆ log 10 E i = 0.1, corresponds approximately to the energy resolution in the lowest energy bin. The latter is chosen to start at 10 18.4 eV, as this is the energy above which the acceptance of the SD array becomes purely geometric and thus independent of the mass and energy of the primary particle. Consequently, in this regime of full efficiency, the calculation of E reduces to a geometrical problem dependent only on the acceptance angle, surface area and live-time of the array. The studies to determine the energy above which the acceptance saturates are described in detail in [30]. Most notably, we have exploited the events detected in hybrid mode as this has a lower threshold than the SD. Assuming that the detection probabilities of the SD and FD detectors are independent, the SD efficiency as a function of energy and zenith angle, (E, θ), has been estimated from the fraction of hybrid events that also satisfy the SD trigger conditions. Above 10 18 eV, the form of the detection efficiency (which will be used in the unfolding procedure described in Sec. IV B) can be represented by an error function: where p 1 = 0.373 and p 0 (θ) = 18.63 − 3.18 cos 2 θ + 4.38 cos 4 θ − 1.87 cos 6 θ.
For energies above E sat = 2.5×10 18 eV, the detection efficiency becomes larger than 97% and the exposure, E, is then obtained from the integration of the aperture of the array over the observation time [30]. The aperture, A, is in turn obtained as the effective area under zenith angle θ, A 0 cos θ, integrated over the solid angle Ω within which the showers are observed. A 0 is welldefined as a consequence of the hexagonal structure of the layout of the array combined with the confinement criterion described in Sec. II B. Each station that has six adjacent, data-taking neighbors, contributes a cell of area A cell = 1.95 km 2 ; the corresponding aperture for showers with θ ≤ 60 • is A cell = 4.59 km 2 sr. The number of active cells, n cell (t), is monitored second-by-second. The array aperture is then given, second-by-second, by the product of A cell by n cell (t). Finally, the exposure is calculated as the product of the array aperture by the number of live seconds in the period under study, excluding the time intervals during which the operation of the array is not sufficiently stable [30]. This results in a duty cycle larger than 95%.
Between 1 January 2004 and 31 August 2018 an exposure (60,400 ± 1,810) km 2 sr yr was achieved. The resulting raw spectrum, J raw i , is shown in Fig. 7, left panel. The energies in the x-axis correspond to the ones defined by the center of the logarithmic bins (10 18.45 , 10 18.55 , · · · eV). The number of events N i used to derive the flux for each energy bin is also indicated. Upper limits are at the 90% confidence level.
The spectrum looks like a rapidly falling power law in energy with an overall spectral index of about 3. To better display deviations from this function we also show, in the right panel, the same spectrum with the intensity scaled by the cube of the energy: the wellknown ankle and suppression features are clearly visible at ≈ 5×10 18 eV and ≈ 5×10 19 eV, respectively.

B. The unfolded spectrum
The raw spectrum is only an approximate measurement of the energy spectrum, J(E), because of the distortions induced on its shape by the finite energy resolution. This causes events to migrate between energy bins: as the observed spectrum is steep, the migration happens especially from lower to higher energy bins, in a way that depends on the resolution function (see Sec. III D, Eq. (2)). The shape at the lowest energies is in addition affected by the form of the detection efficiency (see Sec. IV A, Eq. (6)) in the range where the array is not fully efficient as events whose true energy is below E sat might be reconstructed with an energy above that.
To derive J(E) we use a bin-by-bin correction approach [42], where we first fold the detector effects into a model of the energy spectrum and then compare the expected spectrum thus obtained with that observed so as to get the unfolding corrections. The detector effects are taken into account through the following relationship, J raw (E SD ; s) = dΩ cos θ dE (E, θ)J(E; s)κ(E SD |E; θ) dΩ cos θ (7) where s is the set of parameters that characterizes the model. The model is used to calculate the number of events in each energy bin, µ i (s) = E ∆Ei dE J(E; s). The bin-to-bin migrations of events, induced by the finite resolution through Eq. (7), is accounted for by calculating the number of events expected between E i and E i + ∆E i , ν i (s), through the introduction of a matrix that depends only on the SD response function obtained from the knowledge of the κ(E SD |E) and (E) functions. To estimate µ i and ν i , we use a likelihood procedure, aimed at deriving the set of parameters s 0 allowing the best match between the observed number of events, N i , and the expected one, ν i . Once the best-fit parameters are derived, the correction factors to be applied to the observed spectrum, c i , are obtained from the estimates of µ i and ν i as c i = µ i /ν i . More details about the likelihood procedure, the elements used to build the matrix and the calculation of the c i coefficients are provided in Appendix C.
Guided by the raw spectrum, we infer the possible functional form for J(E; s) by choosing parametric shapes naturally reproducing the main characteristics visible in Fig. 7. As a first step, we set out to reproduce a rapid change in slope (the ankle) followed by a slow suppression of the intensity at high energies. To do so, we use the 6-parameter function: In addition to the normalization, J 0 , and to the arbitrary reference energy E 0 fixed to 10 18.5 eV, the two parameters γ 1 and γ 2 approach the spectral indices around the energy E 12 , identified with the energy of the ankle. The parameter E s marks the suppression energy around which the spectral index slowly evolves from γ 2 to γ 2 + ∆γ. More precisely, it is the energy at which the flux is one half of the value obtained extrapolating the power law after the ankle. It is worth noting that the rate of change of the spectral index around the ankle is here determined by the parameter ω 12 fixed at 0.05, which is the minimal value adopted to describe the transition given the size of the energy intervals. 1 Unlike a model forcing the change in spectral index to be infinitely sharp, such a choice of transition also makes it possible, subsequently, to test the speed of transition by leaving the parameters free. We have used this function (Eq. (8)) to describe our data for over a decade. However, we find that with the exposure now accumulated, it no longer provides a satisfactory fit, with a deviance D/n dof = 35.6/14. A more careful inspection of Fig. 7 suggests a more complex structure in the region of suppression, with a series of power laws rather than a slow suppression. Consequently, we adopt as a second step a functional form describing a succession of power laws with smooth breaks: with j = i + 1. This functional shape is routinely used to characterize the cosmic-ray spectrum at lower energies (see [43] and references therein). The parameters E 23 and E 34 mark the transition energies between γ 2 and γ 3 , and γ 3 and γ 4 respectively. The values of the ω ij parameters are fixed, as previously, at the minimal value of 0.05. In total, this model has 8 free parameters and leads to a deviance of D/n dof = 17.0/12. That this model better matches the data than the previous one is further evidenced by the likelihood ratio between these models which allows a rejection of Eq. (8) with 3.9σ confidence whose calculation is detailed in Appendix C.
As a third step, we release the parameters ω ij one by one, two by two and all three of them so as to test our sensitivity to the speed of the transitions. Free parameters are only adopted as additions if the improvement to the fit is better than 2σ. Such a procedure is expected to result in a uniform distribution of χ 2 probability for the best-fit models, as exemplified in [44]. For every tested model, the increase in test statistics is insufficient to pass the 2σ threshold.
The adoption of Eq. (9) yields the coefficients c i shown as the black points in Fig. 8 together with their statistical uncertainty. To be complete, we also show with a curve the coefficients calculated in sliding energy windows, to explain the behavior of the c i points. This curve is determined on the one hand by the succession of power laws modeled by J(E, s 0 ), and on the other hand by the  response function. The observed changes in curvature result from the interplay between the changes in spectral indices occurring in fairly narrow energy windows (fixed by the parameters ω ij = 0.05) and the variations in the response function. At high energy, the coefficients tend towards a constant as a consequence of the approximately constancy of the resolution, because in such a regime, the distortions induced by the effects of finite resolution result in a simple multiplicative factor for a spectrum in power law. Overall, the correction factors are observed to be close to 1 over the whole energy range with a mild energy dependence. This is a consequence of the quality of the resolution achieved.
We use the coefficients to correct the observed number of events to obtain the differential intensities as J i = c i J raw i . This is shown in the left panel of Fig. 9. The values of the differential intensities, together the detected and corrected number of events in each energy bin are given in Appendix D. The magnitude of the effect of the forward-folding procedure can be appreciated from the following summary: above 2.5×10 18 eV, where there are 215,030 events in the raw spectrum, there are 201,976 in the unfolded spectrum; the corresponding numbers above 5×10 19 eV and 10 20 eV are 278 and 269, and 15 and 14, respectively. Above 5×10 19 eV (10 20 eV), the integrated intensity of cosmic rays is (4.5 ± 0.3) ×10 −3 km −2 yr −1 sr −1 ( 2.4 +0. 9 −0.6 ×10 −4 km −2 yr −1 sr −1 ). In the right panel of Fig. 9, the fitted function J(E, s 0 ), scaled by E 3 to better appreciate the fine structures, is shown as the solid line overlaid on the data points of the final estimate of the spectrum. The characteristics of the spectrum are given in Table III, with both statistical and systematic uncertainties (for which a comprehensive discussion is given in the next section). These characteristics are further discussed in Sec. IV D. There are several sources of systematic uncertainties which affect the measurement of the energy spectrum, as illustrated in Fig. 10.
The systematic uncertainty in the energy scale gives the largest contribution to the overall uncertainty. As described in Sec. III C, it amounts to about 14% and is obtained by adding in quadrature all the systematic uncertainties in the FD energy estimation and the contribution arising from the statistical uncertainty in the calibration parameters. As the effect is dominated by the uncertainty in the calibration of the FD telescopes, the 14% is almost energy independent. Therefore it has been propagated into the energy spectrum by changing the energy of all events by ±14% and then calculating a new estimation of the raw energy spectrum through Eq. (5) and repeating the forward-folding procedure. When considering the resolution, the bias and the detection efficiency in the parameterization of the response function, the energy scale is shifted by ±14%. The uncertainty in the energy scale translates into an energy-dependent uncertainty in the flux shown by a continuous black line in Fig. 10, top panel. It amounts to 30 to 40% around 2.5×10 18 eV, decreasing to 25% around 10 19 eV, and increasing again to 60% at the highest energies.
A small contribution comes from the unfolding procedure. It stems from different sub-components: (i) the functional form of the energy spectrum assumed, (ii) the uncertainty in the bias and resolution parameterization determined in Sec. III D and (iii) the uncertainty in the detection efficiency determined in Sec. IV A. The impact of contribution (i) has been conservatively evaluated by comparing the output of the unfolding assuming Eq. (8) and Eq. (9) and it is less than 1% at all energies. That of contribution (ii) remains within 2% and is maximal at the highest and lowest energies, while the one of contribution (iii) is estimated propagating the statistical uncertainty in the fit function that parametrizes the detection efficiency (Eq. (6)) and it is within 1% below 4×10 18 eV and negligible above. The statistical uncertainties in the unfolding correction factors also contribute to the total systematic uncertainties in the flux and are taken into account. The overall systematic uncertainties due to unfolding are shown as a gray line in both pan-els of Fig. 10 and are at maximum of 2% at the lowest energies.
A third source is related to the global uncertainty of 3% in the estimation of the integrated SD exposure [30]. This uncertainty, constant with energy, is shown as the blue line in both panels of Fig. 10.
A further component is related to the use of an average functional form for the LDF. The departure of this parameterized LDF from the actual one is source of a systematic uncertainty in S(1000). This can be estimated using a subset of high quality events for which the slope of the LDF [26] can be measured on an event by-event basis. The impact of this systematic uncertainty on the spectrum (shown as a black dotted line in Fig. 10) is around 2% at 2.5×10 18 eV, decreasing to −3% at 10 19 eV, before rising again to 3% above 3×10 19 eV. Other sources of systematic uncertainty have been investigated and are negligible.
We have performed several tests to assess the robustness of the measurement. The spectrum, scaled by E 3 , is shown in top panel of Fig. 11 for three zenith angle intervals. Each interval is of equal size in sin 2 θ such that the exposure is the same, one third of the total one. The ratio of the three spectra to the results of the fit performed in the full field of view presented in Sec. IV B is shown in the bottom panel of the same figure. The three estimates of the spectrum are in statistical agreement. In the region below 2×10 19 eV, where there are large numbers of events, the dependence on zenith angle is below 5%. This is a robust demonstration of the efficacy of our methods.
We have also searched for systematic effects that might be seasonal to test the effectiveness of the corrections applied to S(1000) to account for the influence of the changes in atmospheric temperature and pressure on the shower structure [28], and also searched for temporal effects as the data have been collected over a period of 14 years. Such tests have been performed by keeping the energy calibration curve determined in the full data taking period, as the systematic uncertainty associated with a non-perfect monitoring in time of the calibration of the FD telescopes is included in the overall ±14% uncertainty in the energy scale. The integral intensities above 10 19  The total systematic uncertainty, which is dominated by the uncertainty on the energy scale, is obtained by the quadratic sum of the described contributions and is depicted as a dashed red line in Fig. 10.
The systematic uncertainties on the spectral parameters are also obtained adding in quadrature all the contributions above described, and are shown in Table III. The uncertainties in the energy of the features (E ij ) and in the normalization parameter (J 0 ) are dominated by the uncertainty in the energy scale. On the other hand, those on the spectral indexes are also impacted by the other sources of systematic uncertainties. Bottom panel: relative difference between the spectra in the three zenithal ranges and the fitted spectrum in the full f.o.v.. An artificial shift of ±3.5% is applied to the energies in the x−axis for the spectra obtained with the most and less inclined showers to make easier to identify the different data points.

D. Discussion of the spectral features
The unfolded spectrum shown in Fig. 9 can be described using four power laws as detailed in Table III and equation (9). The well-known features of the ankle and the steepening are very clearly evident. The spectral index, γ 3 , used to describe the new feature identified above 1.3×10 19 eV, differs from the index at lower energies, γ 2 , by ≈ 4σ and from that in the highest energy region, γ 4 , by ≈ 5σ.
The representation of our data, and similar sets of spectral data, using spectral indices is long-established although, of course, it is unlikely that Nature generates exact power laws. Furthermore these quantities are not usually derived from phenomenologically-based predictions. Rather it is customary to compare measurements to such outputs on a point-by-point basis (e.g. [45,46]). Accordingly, the data of Fig. 9 are listed in Table VI.
An alternative manner of presentation of the data is shown in Fig. 12 where spectral indices have been computed over small ranges of energy (each point is computed for 3 bins at low energies growing to 6 at the highest energies). The impact of the unfolding procedure is most clearly seen at the lowest energies (where the energy resolution is less good): the effect of the unfolding procedure is to sharpen the ankle feature. It is also clear from Fig. 12 that slopes are constant only over narrow ranges of energy, one of which embraces the new feature starting just beyond 10 19 eV. Above ≈ 5×10 19 eV, where the spectrum begins to soften sharply, it appears that γ rises steadily up to the highest energies observed. However, as beyond this energy there are only 278 events, an understanding of the detailed behaviour of the slope with energy must await further exposure.

V. THE DECLINATION DEPENDENCE OF THE ENERGY SPECTRUM
In the previous section, the energy spectrum was estimated over the entire field of view, using the local horizon and zenith at the Observatory site to define the local zenithal and azimuth angles (θ, ϕ). Alternatively, we can make use of the fixed equatorial coordinates, right ascension and declination (α, δ), aligned with the equator and poles of the Earth, for the same purpose. The wide range of declinations covered by using events with zenith angles up to 60 • , from δ = −90 • to δ +24.8 • (covering 71% of  Figure 13: Left: Energy spectra in three declination bands of equal exposure. Right: Ratio of the declination-band spectra to that of the full field-of-view. The horizontal lines show the expectation from the observed dipole [47]. An artificial shift of ±5% is applied to the energies in the x−axis of the northernmost/southernmost declination spectra to make it easier to identify the different data points. the sky), allows a search for dependencies of the energy spectrum on declination. We present below the determination of the energy spectrum in three declination bands and discuss the results.
For each declination band under consideration, labelled as k, the energy spectrum is estimated as where N ik and c ik stand for the number of events and the correction factors in the energy bin ∆E i and in the declination band considered k, and E k is the exposure restricted to the declination band k. For this study, the observed part of the sky is divided into declination bands with equal exposure, E k = E/3. The correction factors are inferred from a forward-folding procedure identical to that described in section IV, except that the response matrix is adapted to each declination band (for details see Appendix C). The intervals in declination that guarantee that the exposure of the bands are each E/3 are determined by integrating the directional exposure function, ω(δ), derived in Appendix E, over the declination so as to satisfy δ k δ k−1 dδ cos δ ω(δ) δ3 δ0 dδ cos δ ω(δ) where δ 0 = −π/2 and δ 3 = +24.8 • . Numerically, it is found that δ 1 = −42.5 • and δ 2 = −17.3 • . The resulting spectra (scaled by E 3 ) are shown in the left panel of Fig. 13. For reference, the best fit of the spectrum obtained in section IV B is shown as the black line. No strong dependence of the fluxes on declination is observed.
To examine small differences, a ratio plot is shown in the right panel by taking the energy spectrum observed in the whole field of view as the reference. A weightedaverage over wider energy bins is performed to avoid large statistical fluctuations preventing an accurate visual appreciation. For each energy, the data points are observed to be in statistical agreement with each other. Note that the same conclusions hold when analyzing data in terms of integral intensities, as evidenced for instance in table IV above 8×10 18 eV. Similar statistical agreements are found above other energy thresholds. Hence this analysis provides no evidence for a strong declination dependence of the energy spectrum.
A 4.6% first-harmonic variation in the flux in right ascension has been observed in the energy bins above 8×10 18 eV shown in the right panel of Fig. 13 [47]. It is thus worth relating the data points reported here to these measurements that are interpreted as dipole anisotropies. The technical details to establish these relationships are given in Appendix E.
The energy-dependent lines drawn in Fig. 13-right show the different ratios of intensity expected from the dipolar patterns in each declination band relative to that across the whole field of view. The corresponding data points are observed, within uncertainties, to be in fair agreement with these expectations.
Overall, there is thus no significant variation of the spectrum as a function of the declination in the field of view scrutinized here. A trend for a small declination dependence, with the flux being higher in the Southern hemisphere, is observed consistent with the dipolar patterns reported in [47]. At the highest energies, the event numbers are still too small to identify any increase or decrease of the flux with the declinations in our field of view.  Currently, the Telescope Array (TA) is the leading experiment dedicated to observing UHECRs in the northern hemisphere. As already pointed out, TA is also a hybrid detector making use of a 700 km 2 array of SD scintillators overlooked by fluorescence telescopes located at three sites. Although the techniques for assigning energies to events are similar, there are differences as to how the primary energies are derived, which result in differences in the spectral estimates, as can be appreciated in Fig. 14 where the E 3 -scaled spectrum derived in this work and the one derived by the TA Collaboration [48] are shown.
A useful way to appraise such differences is to make a comparison of the observations at the position of the ankle. Given the lack of anisotropy in this energy range, this spectral feature must be quasi-invariant with respect to direction on the sky. The energy at the ankle measured using the TA data is found to be (4.9 ± 0.1 (stat.))×10 18 eV, with an uncertainty of 21% in the energy scale [49] in good agreement with the one reported here ((5.0 ± 0.1 (stat.) ± 0.8 (sys.))×10 18 eV). Consistency between the two spectra can be obtained in the ankle-energy region up to 10 19 eV by rescaling the energies by +5.2% for Auger and −5.2% for TA. The factors are smaller than the current systematic uncertainties in the energy scale of both experiments. These values encompass the different fluorescence yields adopted by the two Collaborations, the uncertainties in the absolute calibration of the fluorescence telescopes, the influence of the atmospheric transmission used in the reconstruction, the uncertainties in the shower reconstruction, and the uncertainties in the correction factor for the invisible energy. It is worth noting that better agreement can be obtained if the same models are adopted for the fluorescence yield and for the invisible energy correction. Detailed discussions on these matters can be found in [50].
However, even after the rescaling, differences persist above 10 19 eV. At such high energies, anisotropies might increase in size and induce differences in the energy spectra detected in the northern and southern hemispheres. To disentangle possible anisotropy issues from systematic effects, a detailed scrutiny of the spectra in the declination range accessible to both observatories has been carried out [51]. A further empirical, energydependent, systematic shift of +10% (−10%) per decade for Auger (TA) is required to bring the spectra into agreement. A comprehensive search for energy-dependent systematic uncertainties in the energies has resulted in possible non-linearities in this decade amounting to ±3% for Auger and (−0.3 ± 9)% for TA, which are insufficient to explain the observed effect [52]. A joint effort is underway to understand further the sources of the observed differences and to study their impact on the spectral features [53].

VII. SUMMARY
We have presented a measurement of the energy spectrum of cosmic rays for energies above 2.5×10 18 eV based on 215,030 events recorded with zenith angles below 60 • . The corresponding exposure of 60,400 km 2 yr sr, calculated in a purely geometrical manner, is independent of any assumption on unknown hadronic physics or primary mass composition. This measurement relies on estimates of the energies that are similarly independent of such assumptions. This includes the analysis that minimizes the model/mass dependence of the invisible energy estimation as presented in [33]. In the same manner, the flux correction for detector effects is evaluated using a data-driven analysis. Thus the approach adopted differs from that of all other spectrum determinations above 5×10 14 eV where the air-shower phenomenon is used to obtain information.
The measurement reported above is the most precise made hitherto and is dominated by systematic uncertainties except at energies above 5×10 19 eV. The systematic uncertainties have been discussed in detail and it is shown that the dominant one ( 14%) comes from the energy scale assigned using measurements of the energy loss by ionisation in the atmosphere inferred using the fluorescence technique.
In summary, the principal conclusions that can be drawn from the measurement are: 1. The flattening of the spectrum near 5×10 18 eV, the so-called "ankle", is confirmed.
2. The steepening of the spectrum at around 5×10 19 eV is substantiated.
4. No evidence for any dependence of the energy spectrum on declination has been found other than a mild excess from the Southern Hemisphere that is consistent with the anisotropy observed above 8×10 18 eV.
A discussion of the significance of these measurements from astrophysical perspectives can be found in [54].
Appendix A: Expected distribution of cosmic rays in local zenith angle To derive the attenuation curves in a purely datadriven way, we have required the same intensity in equal intervals of sin 2 θ. This requirement relies on the high level of isotropy of the cosmic-ray intensity. In this appendix it is shown that neither the small anisotropies at the energy thresholds of interest, nor the slight zenithal dependencies of the response function of the SD array, alter the constancy of the intensity in terms of sin 2 θ, by more than 0.5% at 3×10 18 eV. This is less than the magnitude of the statistical fluctuations at this energy.
Consider a possible directional dependence for the energy spectrum, J(E, θ(α, δ, t), ϕ(α, δ, t)). Accounting for the energy resolution, the expected event number per steradian and per unit time above any energy threshold E 0 is given by 2 To get the expected number of events in each of the sin 2 θ intervals, it is convenient to consider the left hand side of Eq. (A1) expressed in terms of local sidereal time through the transformation Inclusion of the Dirac function guarantees that the direction α 0 (t) considered throughout the time integration corresponds to the local sidereal time α 0 seen at time t. On inserting Eq. (A1) into Eq. (A2) and carrying out the integration over time, d 2 N/dΩdα 0 becomes where the notation N cell (α 0 ) stands for the total number of active hexagonal cells during the integrated observation time for a flux of cosmic rays from each direction α 0 , N cell (α 0 ) ≡ dt n cell (t)δ(α 0 − α 0 (t)). Due to the long period considered here ( 4,880 sidereal days), an averaging takes place and this function is nearly uniform, N cell (α 0 ) = N 0 cell + δN cell (α 0 ), with δN cell (α 0 )/N 0 cell of the order of a few 10 −3 . The expected d 2 N/d sin 2 θ distribution is then obtained by integrating Eq. (A3) over azimuth ϕ and local sidereal time α 0 .
Characterizing anisotropies, to first order, by a dipole vector with amplitude d and equatorial directions (α d , δ d ), an effective ansatz for the spectrum is then J(E, θ, ϕ, α 0 ) = J 0 (E)(1 + d(α d , δ d ) · n(θ, ϕ)), with n(θ, ϕ) the unit vector on the sphere. The integration over the azimuthal angle in Eq. (A4) selects the coordinate of the dipole vector along the local z−axis defining the zenith angle. The remaining integration over the local sidereal time α 0 cancels the harmonic contribution δN cell (α 0 ) coupled to the isotropic part of the spectrum, and for small anisotropies, the leading order of dN/d sin 2 θ reads as which is the desired expression.
Searches for anisotropies throughout the EeV energy range have been reported previously [47]. We show in Fig. 15 the expected departures from a uniform behaviour for the dN (> E 0 )/d sin 2 θ distribution normalized to its average value, above three different thresholds. The shaded band stands for the statistical fluctuations in each bin of sin 2 θ above the energy threshold under consideration, given the number of events available. Above 3×10 18 eV (left panel), the upper limit previously obtained with a smaller data set between 2×10 18 eV and 4×10 18 eV is used, and the two directions maximizing the departure from a uniform behavior are selected (δ d = ±π/2). Above 8×10 18 eV and 1.6×10 19 eV, dipole parameters reported in [47] are adopted. One can see that the small departures from uniformity all lie within the limits set by the statistical uncertainties, validating the use of the CIC method to derive the attenuation curves at the different energy thresholds.

Appendix B: Energy calibration of S38
To ease the notations in this appendix, we denote here as S the underlying shower-size parameter S 38 , and as S SD its estimator. The probability density function, f 0 (E, S), to detect an hybrid event of underlying energy E and shower size S at ground that would be expected without shower-to-shower fluctuations and with an infinite energy resolution is proportional to where P SD (S) is the SD detection efficiency expressed in terms of S and h(E) is the energy distribution of the events, that is, the underlying energy spectrum multiplied by the effective exposure for the hybrid events. The Dirac function guarantees that the shower size can be modelled with a function S(E) of the true energy. The probability density function for the detection of an hybrid shower that includes the finite energy resolution can be derived by folding f 0 (E, S) with both FD and SD resolution functions, taken as Gaussian distributions with standard deviations σ FD and σ SD , respectively: where E FD is the energy measured by the FD and S SD is the shower size at ground measured by the SD. Here, σ SD accounts for both the shower-to-shower fluctuations of the shower size, σ sh , and for the detector resolution, σ det : The function h(E) is significantly less steep than the energy spectrum because the criteria to select high-quality FD events and to guarantee an unbiased X max distribution are more effective at lower energies. Its estimation through Monte-Carlo simulations to the required accuracy is a difficult task. It is preferable to follow the strategy put forward in [38], that consists of deriving h(E) from the data distribution using the boostrap approximation where the sum runs over the hybrid events. The hybrid events are selected according to the criteria reported in section III C, at energies where the SD is almost fully efficient. This allows us to neglect the effect of P SD (S(E)). Inserting Eq. (B4) into (B2) leads to where σ FDi and σ SDi are the uncertainties on the measurements evaluated on an event-by-event basis. Finally the calibration parameters A and B that enter in f (E FD , S SD ) through the relationship S(E) = (E/A) 1/B are determined maximizing the following likelihood: where the sum with index k runs over the hybrid events selected for the energy calibration. These events have E FDk > 3×10 18 eV and the SD is fully efficient. The sum over i, that defines the probability density function, extends to lower energies to capture the fluctuations of the energy estimators. It is sufficient to select events with energies E FDi > 2.5×10 18 eV, for which the detection efficiency is still very close to 100% (see section IV A) and the power law E = AS B is still valid (see section III B). Given the good FD energy resolution (≈ 7.4%), events below this energy would give a negligible contribution to the likelihood. Once the parameters that best describe the data are determined, the goodness of the fit is estimated through the following deviance definition, where the second term represents an ideal model in which the shower size distribution is perfectly described by the fitted power law.
Appendix C: Details of the forward-folding procedure The response matrix elements, R ij , are the conditional probabilities that the reconstructed energy E SD of an event falls into the bin i given that the true energy E should have been in the bin j: where the zenithal part of the angular integration is performed up to θ max = 60 • . It is used to calculate the number of events which is expected between E i and E i +∆E i , ν i = j R ij µ j , where µ j are the ones expected without detector effects between E j and E j + ∆E j . For sufficiently small bin sizes so that the values of κ and, below full efficiency, are approximately constant, the dependence on J cancels out. The forward-folding fit is thus performed under this approximation (we use ∆ log 10 E = 0.01), allowing the re-calculation of the matrix elements to be avoided at each step of the fit (leading to the definition of a matrix R kl ). We then integrate the expected number of events µ l and ν k (ν k = l R kl µ l ) obtained with the best fit parameters to calculate the number of events in the ∆ log 10 E = 0.1 bins and, from them, the correction coefficients. The R ij matrix calculated in the ∆ log 10 E = 0.1 bins according to Eq. (C1) is reported in Table V. It can be used for testing any model J(E; s) that fits reasonably well the data.
The observed number of events as a function of energy is a single measurement of a random process for which the p.d.f. for observing the set of values N i given a set of expectations ν i follows a multinomial distribution. The total number of events N = i N i being itself a random variable from a Poisson process, the resulting joint p.d.f. for the histogram is the product, over the energy bins considered, of the Poisson probabilities to observe N i events in each bin given an expectation ν i (s). Dropping the constant terms with respect to s, the likelihood function to be minimized then reads: Once the best-fit parameters s are derived, the correction factors c i are then obtained from the estimates of µ i and ν i as c i = µ i /ν i . The goodness-of-fit statistic is provided by its deviance,   which asymptotically behaves as a χ 2 variable with k − s degrees of freedom, where k is the number of measurements (the number of energy bins) and s the number of model parameters [55]. The uncertainties in the energy spectrum that is recovered follow from the covariance matrix of the µ i estimators, which for a Poisson process is given by so that the uncertainties σ Ji scale as √ c i N i /(E ∆E i ). The confidence intervals reported in this paper are then estimated by calculating the 2-sided 16% quantiles of the underlying p.d.f while the 90% confidence-level limits are calculated according to [56].
In section V, the energy spectrum is reported for specific ranges of declinations. The forward-folding procedure used to infer the different spectra is then identical, adapting the response matrix to the declination range under consideration in the following way. The directional raw energy spectrum, J raw (E SD , α, δ; s), is related to the underlying energy spectrum through J raw (E SD , α, δ; s) = 1 ∆t dtdEdΩ cos θ (E, θ)κ(E SD |E, θ)J(E, α, δ; s)δ(Ω − Ω(α, δ, t)), where the Dirac function guarantees that only the local angles (θ, ϕ) pointing to the (α, δ) considered contribute to the integration at time t. By applying the Dirac function, and by using Eq. (E4), the response matrix elements can be shown to be R ij = ∆Ei dE SD ∆Ej dE dh dδ cos δ cos (θ(δ, h)) (E, θ(δ, h))κ(E SD |E, θ(δ, h)) dh dδ cos δ cos (θ(δ, h))H(cos (θ(δ, h)) − cos θ max )∆E j , where the time integration has been substituted for an integration over h = α 0 (t)−α between −π and π, and the Heaviside function, H, guarantees that only the effective zenithal range [0, θ max ] contributes to the integrations. Note that the integration over declination covers only the range under consideration in the numerator, while it covers the whole field of view in the denominator.
The functional shape J(E; s) that best describes the energy spectrum is selected by adopting for the test statistics (TS) the likelihood ratio between an alternative model and a reference one. For each model, the forwardfolding fit is carried out and the corresponding likelihood value is recorded. The TS values are first converted into p-values by integrating the distributions of the TS for the reference model above the value obtained in data. The p-values are then converted into significances assuming 1-sided Gaussian distributions.
The reference model, Eq. (8), is the one that we have used for over a decade. With the new model, Eq. (9), which has two additional parameters to define the feature at 1.3 × 10 19 eV, TS 20. As the two hypotheses are not nested, the likelihood ratio distribution is built by Monte-Carlo to calculate the corresponding p-value. (with the reference model for J(E; s)) by drawing at random a total number of events similar to that of the data. The spectra of these samples were then reconstructed ac-cording to the method presented in Sec. IV.B with both the reference model and the alternative one. For each sample, the TS has been recorded. The distribution obtained is shown in Fig. 16. There are only 10 counts out of 210,000 above TS 20, thus enabling us to reject the reference model at the 3.9σ confidence level.
The test statistic has also been adopted to test our sensitivity to the speed of the transitions. In this case the reference model is Eq. (9), with all paramaters w ij fixed to 0.05, and the alternative one obtained leaving them free in the fit. Since the hypotheses are in this case nested, Wilks theorem applies. For each of the models tested, the increase in TS is less than 2σ. anisotropy measurements.
The directional exposure results from the timeintegration of the directional aperture of an active region of the array, which is considered here as running constantly. To first order, it is well approximated by A cell cos θ, but it is actually slightly larger for showers arriving from the downhill direction due to the small tilt of the Observatory towards the south-east. Although small, it is important to account here for this effect because the directional exposure estimated by neglecting it would distort the cosmic-ray flux in an overall dipolar-shaped way with an amplitude of 0.5% along the declination coordinate. That would consequently bias the energy spectrum estimates in different declination bands. For an angle of incidence (θ, ϕ), the directional aperture per cell is thus given, on average, by A cell (θ, ϕ) 1.95 km 2 (1 + θ tilt tan θ cos (ϕ − ϕ tilt )) cos θ, (E1) with θ tilt = 0.2 • the average inclination to the vertical and ϕ tilt = −30 • the direction in azimuth counterclockwise from east. From this expression, the direc-tional exposure is finally estimated by making use of the time-dependent transformation rules relating the equatorial coordinates (α, δ) to the corresponding local ones (θ(α, δ, t), ϕ(α, δ, t)) at time t, sin θ cos ϕ = − cos δ sin (α 0 (t) − α), (E2) sin θ sin ϕ = cos sin δ − sin cos δ cos (α 0 (t) − α),(E3) cos θ = cos cos δ cos (α 0 (t) − α) + sin sin δ, (E4) where α 0 (t) is the local sidereal time and the latitude of the Observatory. Since the time integration of A cell (θ(α, δ, t), ϕ(α, δ, t)) depends here only on the difference α 0 (t) − α, it can be substituted for an integration over the hour angle h = α 0 (t) − α. Following [57], the constraining θ max (60 • here) translates to an integration over h ranging from −h m to h m with h m = arccos [(cos θ max − sin sin δ)/(cos cos δ)], with the additional constraining that h m = 0 for declinations giving rise to an argument greater than 1 in the arccos function, and h m = π for declinations giving rise to an argument smaller than −1. This leads to the expression ω(δ) ∝ cos cos δ sin h m (δ) + h m (δ) sin sin δ + θ tilt sin ϕ tilt (h m (δ) cos sin δ − sin cos δ sin h m (δ)), which is then used to determine numerically the δ k values in Eq. (11). The normalization of ω is such that E = 2πA 0 ∆t δ3 δ0 dδ cos δ ω(δ). To compare the unfolded spectra obtained for each declination band with those expected from the anisotropy measurements characterized by a vector, d, with amplitude d and equatorial directions (α d , δ d ), we make use of the ansatz dN i dΩ = A 0 ∆t ω(δ) (1 + d(α d , δ d ) · n(α, δ)) ∆Ei dE J 0 (E), considering d constant within ∆E, and with dΩ = cos δdδdα. The choice of the energy bins, ∆E i , follows from that performed in [47] in which the values of d are reported. The integration over right ascension selects the isotropic component and the dipole component along the z−axis defining the declination: Finally, the expected number of events in a given declination band is obtained through a final integration over declination: where the division in three declination bands is explicitly used. The expected intensity is then J ik = N ik /(E/3)/∆E i in each energy bin. The ratios of intensities, r k , depicted by the lines in Fig. 13-right are thus drawn from r = 3 .

(E9)
Note that in these calculations, the dead times of the SD array have not been considered. As discussed in Ap-pendix A, they lead to variations of the event rate of order of a few 10 −3 , imprinting small harmonic dependencies in right ascension to the directional exposure, ω (α, δ) = ω(δ)(1 + δω(α, δ)). Throughout the integrations over 2π in right ascension, these harmonic dependencies cancel exactly when coupled to the isotropic component of the intensity, and give rise to second-order corrections when coupled to the dipolar fraction of amplitude d. The term δω(α, δ) can thus be safely neglected here.