The energy spectrum of cosmic rays beyond the turn-down around $10^{17}$ eV as measured with the surface detector of the Pierre Auger Observatory

We present a measurement of the cosmic-ray spectrum above 100\,PeV using the part of the surface detector of the Pierre Auger Observatory that has a spacing of 750~m. An inflection of the spectrum is observed, confirming the presence of the so-called \emph{second-knee} feature. The spectrum is then combined with that of the 1500\,m array to produce a single measurement of the flux, linking this spectral feature with the three additional breaks at the highest energies. The combined spectrum, with an energy scale set calorimetrically via fluorescence telescopes and using a single detector type, results in the most statistically and systematically precise measurement of spectral breaks yet obtained. These measurements are critical for furthering our understanding of the highest energy cosmic rays.


Introduction
The steepening of the energy spectrum of cosmic rays (CRs) at around 10 15.5 eV, first reported in [1], is referred to as the "knee" feature. A widespread view for the origin of this bending is that it corresponds to the energy beyond which the efficiency of the accelerators of the bulk of Galactic CRs is steadily exhausted. The contribution of light elements to the all-particle spectrum, largely dominant at GeV energies, remains important up to the knee energy after which the heavier elements gradually take over up to a few 10 17 eV [2][3][4][5][6]. This fits with the long-standing model that the outer shock boundaries of expanding supernova remnants are the Galactic CR accelerators, see e.g. [7] for a review. Hydrogen is indeed the most abundant element in the interstellar medium e-mail: spokespersons@auger.org that the shock waves sweep out, and particles are accelerated by diffusing in the moving magnetic heterogeneities in shocks accordingly to their rigidity. That the CR composition gets heavier for two decades in energy above the knee energy could thus reflect that heavier elements, although subdominant below the knee, are accelerated to higher energies, until the iron component falls off steeply at a point of turndown around 10 16.9 eV. Such a bending has been observed in several experiments at a similar energy, referred to as the "second knee" or "iron knee" [8][9][10][11]. The recent observations of gamma rays of a few 10 14 eV from decaying neutral pions, both from a direction coincident with a giant molecular cloud [12] and from the Galactic plane [13], provide evidence for CRs indeed accelerated to energies of several 10 15 eV, and above, in the Galaxy. A dozen of sources emitting gamma rays up to 10 15 eV have even been reported [14], and the production could be of hadronic origin in at least one of them [15]. However, the nature of the sources and the mechanisms by which they accelerate CRs remain in general undecided. In particular, that particles can be effectively accelerated to the rigidity of the second knee in supernova remnants is still under debate, see e.g. [16].
Above 10 17 eV, the spectrum steepens in the interval leading up to the "ankle" energy, ∼5×10 18 eV, at which point it hardens once again. The inflection in this energy range is not as sharp as suggested by the energy limits reached in the Galactic sources to accelerate iron nuclei beyond the ironknee energy [17]. Questions arise, then, on how to make up the all-particle spectrum until the ankle energy. The hardening around 10 17.3 eV in the light-particle spectrum reported in [18] is suggestive of an extragalactic contribution to the all-particle spectrum steadily increasing. It has even been argued that an additional component is necessary to account for the extended gradual fall-off of the spectrum and for the mass composition in the iron-knee-to-ankle region, be it of Galactic [17] or extragalactic origin [19]. While the concept that the Galactic-to-extragalactic transition occurs somewhere between 10 17 eV and a few 10 18 eV is well-accredited, a full understanding of how it occurs is hence lacking. The approximately power-law shape of the spectrum in this energy range may mask a complex superposition of different components and phenomena, the disentanglement of which rests on the measurements of the all-particle energy spectrum, and of the abundances of the different elements as a function of energy, both of them challenging from an experimental point of view. On the one hand, the energy range of interest is accessible only through indirect measurements of CRs via the extensive air showers that they produce in the atmosphere. Therefore, the determination of the properties of the CRs, especially their mass and energy, is prone to systematic effects. On the other hand, different experiments, different instruments and different techniques of analysis are used to cover this energy range, so that a unique view of the CRs is only possible by combining measurements the matching of which inevitably implies additional systematic effects.
The aim of this paper is to present a measurement of the CR spectrum from 10 17 eV up to the highest observed energies, based on the data collected with the surface-detector array of the Pierre Auger Observatory. The Observatory is located in the Mendoza Province of Argentina at an altitude of 1400 m above sea level at a latitude of 35.2 • S, so that the mean atmospheric overburden is 875 g/cm 2 . Extensive air showers induced by CR-interactions in the atmosphere are observed via a hybrid detection using a fluorescence detector (FD) and a surface detector (SD).
The FD consists of five telescopes at four sites which look out over the surface array, see Fig. 1. Four of the telescopes (shown in blue) cover an elevation range from 0 • to 30 • while the fifth, the High Elevation Auger Telescopes (HEAT), covers an elevation range from 30 • to 58 • (shown in red). Each telescope is used to collect the light emitted from air molecules excited by charged particles. After first selecting the UV band with appropriate filters (310 to 390 nm), the light is reflected off a spherical mirror onto a camera of 22×20 hexagonal, 45.6 mm, photo-multiplier tubes (PMTs). In this way, the longitudinal development of the particle cascades can be studied and the energy contained within the electromagnetic sub-showers can be measured in a calorimetric way. Thus the FD can be used to set an energy scale for the Observatory that is calorimetric and so is independent of simulations of shower development.
The SD, the data of which are the focus of this paper, consists of two nested hexagonal arrays of water Cherenkov detectors (WCDs). The layout, shown in Fig. 1, includes the SD-1500, with detectors spread apart by 1500 m and totaling approximately 3000 km 2 of effective area. The detectors of the SD-750 are instead spread out by 750 m, yielding an effective area of 24 km 2 . SD-750 and SD-1500 include identical WCDs, cylindrical tanks of pure water with a 10 m 2 base and a height of 1.2 m. Three 9" PMTs are mounted to the top of each tank and view the water volume. When relativistic secondaries enter the water, Cherenkov radiation is emitted, reflected via a Tyvek lining into the PMTs, and digitized using 40 MHz 10-bit Flash Analog to Digital Converters (FADCs). Each WCD along with its digitizing electronics, communication hardware, GPS, etc., is referred to as a station.
Using data collected over 15 years with the SD-1500, we recently reported the measurement of the CR energy spectrum in the range covering the region of the ankle up to the highest energies [20,21]. In this paper we extend these measurements down to 10 17 eV using data from the SD-750: not only is the detection technique consistent but the same methods are used to treat the data and build he spectrum. The paper is organized as follows: we first explain how, with the SD-750 array, the surface array is sensitive to primaries down to 10 17 eV in Section 2; in Section 3, we describe how we reconstruct the showers up to determining the energy; we illustrate in Section 4 the approach used to derive the energy spectrum from SD-750; finally, after combining the spectra measured by SD-750 and SD-1500, we present the spectrum measured using the Auger Observatory from 10 17 eV upwards in Section 5 and discuss it in the context of other measurements in Section 6.

Identification of Showers with the SD-750: From the Trigger to the Data Set
The implementation of an additional set of station-level trigger algorithms in mid-2013 is particularly relevant for the operation of the SD-750. Their inclusion in this work extends the energy range over which the SD-750 triggers with > 98% probability from 10 17.2 eV down to 10 17 eV.
To identify showers, a hierarchical set of triggers is used which range in scope from the individual station-level up to the selection of events and the rejection of random coincidences. The trigger chain, extensively described in [22], has been used since the start of the data taking of the SD-1500, and was successively adopted for the SD-750. In short, station-level triggers are first formed at each WCD. They are then combined with those from other detectors and examined for spatial and temporal correlations, leading to an array trigger, which initiates data acquisition. After that, a similar hierarchical selection of physics events out of the combinatorial background is ultimately made.
We describe in this section the design of the triggers (Section 2.1). We then illustrate their effect on the data, at the level of the amplitude of detected signals (Section 2.2) and on the timing of detected signals in connection with the event selection (Section 2.3). Finally we describe the energy at which acceptance is 100% (Section 2.4). A more detailed description of the trigger algorithms can be found in Appendix A.

The Electromagnetic Triggers
Using the station-level triggers, the digitized waveforms are constantly monitored in each detector for patterns consistent with what would be expected as a result of air-shower secondary particles (primarily electrons and photons of 10 MeV on average, and GeV muons) entering the water volume . The typical morphologies include large signals, not necessarily spread in time, such as those close to the shower core, or sequences of small signals spread in time, such as those nearby the core in low-energy showers, or far from the core in high-energy ones. Atmospheric muons, hitting the WCDs at a rate of 3 kHz, are the primary background. The output from the PMTs has only a small dependence on the muon energy. The electromagnetic and hadronic background, while also present, yields a total signal that is usually less than that of a muon. Consequently, the atmospheric muons are the primary impediment to developing a station-level trigger for small signal sizes without contaminating the sampling of an air shower with spurious muons.
The response of an individual WCD to secondary particles has been studied using unbiased FADC waveforms and dedicated studies of signals from muons [23].
Originally, two triggers were implemented into the station firmware, called threshold (TH), more adept to detect muons, and time-over-threshold (ToT), more suited to identify the electromagnetic component. Both of these have settings which require the signal to be higher in amplitude or longer than what is observed for a muon traveling vertically through the water volume. As such, they have the inherent limitation of being insensitive to signals which are smaller than (or equal to) that of a single muon, thus prohibiting the measurement of pure electromagnetic signals, which are generally smaller.
To bolster the sensitivity of the array to such small signals, two additional triggers were designed. The first, timeover-threshold-deconvolved (ToTd), first removes the typical exponential decay created by Cherenkov light inside the water volume, after which the ToT algorithm is applied. The second, multiplicity-of-positive-steps (MoPS), is designed to select small, non-smooth signals, a result of many electromagnetic particles entering the water over a longer period of time than a typical muon pulse. This is done by counting the number of instances in the waveform where consecutive bins are increasing in amplitude. Both of the trigger algorithms are described in detail in Appendix A.
The implementation of the ToTd and MoPS (the rate of which is around 0.3 Hz, compared to 0.6 Hz of ToT and 20 Hz of TH) did not require any modification in the logic of the array trigger, which calls for a coincidence of three or more SD stations that pass any combination of the triggers described above with compact spacing, spatially and temporally [22]. We note that in spite of the low rate of the ToTd and MoPS relative to TH and ToT, the array rate more than doubled after their implementation. This, as will be shown in the following, is due to the extension of measurements to the more abundant, smaller signals.

Effect of ToTd and MoPS on Signals Amplitudes
The ToTd and MoPS triggers extend the range over which signals can be observed at individual stations into the region which is dominated by the background muons that are created in relatively low energy air showers. By remaining insensitive to muon-like signals, these two triggers increase the sensitivity of the SD to the low-energy parts of the showers that have previously been below the trigger threshold.
The effects of the additional triggers can be seen in the distribution of the observed signal sizes. An example of such a distribution, based on one month of air-shower data, is shown in Fig. 2. The signal sizes are shown in the calibration unit of one vertical equivalent muon (VEM), the total deposited charge of a muon traversing vertically through the water volume [22]. For the stations passing only the ToT and TH triggers (shown in solid black), the distribution of  deposited signals is the convolution of three effects, the uniformity of the array, the decreasing density of particles as a function of perpendicular distance to the shower axis (henceforth referred to as the axial distance), and the shape of the CR spectrum resulting in the negative slope above 7 VEM. Furthermore there is a decreasing efficiency of the ToT and TH at small signal sizes. The range of additional signals that are now detectable via the ToTd and MoPS triggers are shown in dashed red. As expected, ToTd and MoPS triggers increase the probability of the SD to detect small amplitude signals, namely between 0.3 and 5 VEM. That the high-signal tail of this distribution ends near 10 VEM is consistent with a previous study [24] that estimated that the ToT+TH triggers were fully efficient above this value. The additional sensitivity to small air-shower signals also increases the multiplicity of triggered stations per event. This increase is characterized in Fig. 3, which shows the number of additional triggered stations per event as a function of the number of stations that pass the TH and ToT triggers, after removing spuriously triggered stations. The median increase of multiplicity in each horizontal bin is shown by the black circles and indicates a typical increase of one station per event.

Effects of ToTd and MoPS on Signal Timing
The increased responsiveness of the ToTd and MoPS algorithms to smaller signals, specifically due to the electromagnetic component, has an effect also on the observed timing of the signals. In general, the electromagnetic signals are expected to be delayed with respect to the earliest part of the shower which is muon-rich, the delay increasing with axial distance. Further, in large events, stations that pass these triggers tend to be on the edge of the showers, where the front is thicker, thus increasing the variance of the arrival times. Such effects can be seen through the distribution of the start times for stations that pass the ToTd and MoPS triggers. The residuals of the pulse start times with respect to a plane front fit of the three stations with the largest signals in the event are shown in Figure 4 for different trigger types. The entries shown in blue correspond to stations that passed the ToT algorithm, the ones in green to stations that pass the TH trigger (but not the ToT trigger), and those in red to stations that pass the ToTd and/or MoPS triggers, only. For each of the trigger types, there is a clear peak near zero, which reflects the approximately planar shower front close to the core. Stations that pass the TH condition, but not the ToT one, tend to capture isolated muons, including background muons arriving randomly in time. This explains the vertical offset, flat and constant, in the green curve. In turn, the lack of such a baseline shift in the blue and red distributions gives evidence that the ToT, TOTd and MoPS algorithms reject background muons effectively. This is particularly successful for the ToTd and MoPS that accept very small signals, of approximately 1 VEM in size. One can see that these distributions have different shapes and that, in particular, the start time distributions of signals that pass the ToTd and MoPS have much longer tails than those of the TOT triggers, including a second distribution beginning around 1.5 µs possibly due to heavily delayed electromagnetic particles. The extended time portion of showers accessed by the ToTd and MoPS triggers has implications on the procedure used to select physical events from the triggered ones [22]. In this process, non-accidental events, as well as non-accidental stations, are disentangled on the basis of their timing. First, we identify the combination of three stations where they form a triangle, in which at least two legs are 750 m long, and where they have the largest summed signal among all such possible configurations. These stations make up the event seed and the arrival times of the signals are fit to a plane front. Additional stations are then kept if their temporal residual, Δ , is within a fixed window, low < Δ < high . Motivated by the differing time distributions, updated low and high values were calculated based on which trigger algorithm was satisfied. Using the distributions of timing residuals, shown in Fig. 4, the baseline was first subtracted. Then the limits of the window, low and high , were chosen such that the middle 99% of the distribution was kept. The trigger-wise limits are summarized in Table 1.

Effect of the ToTd and MoPS on the energy above which acceptance is fully-efficient
Most relevant to the measurement of the spectrum is the determination of the energy threshold above which the SD-750 becomes fully efficient. To derive this, events observed by the FD were used to characterize this quantity as a function of energy and zenith angle. The FD reconstruction requires only a single station be triggered to yield a robust determination of the shower trajectory. Using the FD events with energies above 10 16.8 eV, the lateral trigger probability (LTP), the chance that a shower will produce a given SD trigger 16  The detection efficiency of the SD-750 for air showers with < 40 • is shown for the original (dashed red) and expanded (solid blue) station-level trigger sets with bands indicating the systematic uncertainties. The trigger efficiency was determined using data above 10 16.8 eV and is extrapolated below this energy (shown in gray). as a function of axial radius, was calculated for all trigger types. The LTP was then parameterized as a function of the observed air-shower zenith angle and energy. It is important to note that because the LTP is derived using observed air showers as a function of energy, this calculation reflects the efficiency as a function of energy based on the true underlying mass distribution of primary particles. Further details of this method can be found in [25].
The SD-750 trigger efficiency was then determined via a study in which isotropic arrival directions and random core positions were simulated for fixed energies between 10 16.5 and 10 18 eV. Each station on the array was randomly triggered using the probability given by the LTP. The set of stations that triggered were then checked against the compactness criteria of the array-level triggers, as described in [22]. The resulting detection probability for showers with zenith angles < 40 • is shown as a solid blue line in Fig. 5 as a function of energy. The detection efficiency becomes almost unity (> 98%) at around 10 17 eV. For comparison, we show in the same figure, in dashed red, the detection efficiency curve for the original set of station-triggers, TH and ToT, in which the full efficiency is attained at a larger energy, i.e., around 10 17.2 eV.
A description for the detection efficiency, ( ), below 10 17 eV, will be important for unfolding the detector effects close to the threshold energy (see Section 4). This quantity was fit using the results of the LTP simulations with < 40 • The energy-cut corresponding to the full-efficiency threshold increases with zenith angle, due to the increasing attenuation of the electromagnetic component with slant depth. The zenith angle 40 • was chosen as a balance to have good statistical precision and a low energy threshold. and is well-parameterized by where erf ( ) is the error function, = 16.4 ± 0.1 and = 0.261 ± 0.007. For events used in this analysis, there is an additional requirement regarding the containment of the core within the array: only events in which the detector with the highest signal is surrounded by a hexagon of six stations that are fully operational are used. This criterion not only ensures adequate sampling of the shower but also allows the aperture of the SD-750 to be evaluated in a purely geometrical manner [22]. With these requirements, the SD-750 data set used below consists of about 560,000 events with < 40 • and > 10 17 eV recorded between 1 January 2014 and 31 August 2018. The minimum energy cut is motivated by the lowest energy to which we can cross-calibrate with adequate statistics the energy scale of the SD with that of the FD (see Section 3.3). The corresponding exposure, E, after removal of time periods when the array was unstable (<2% of the total) is E = (105 ± 4) km 2 sr yr.

Energy Measurements with the SD-750
In this section, the method for the estimation of the airshower energy is detailed together with the resulting energy resolution of the SD-750 array. The measurement of the actual shower size is first described in Section 3.1 after which the corrections for attenuation effects are presented in Section 3.2. The energy calibration of the shower size after correction for attenuation is presented in Section 3.3. The energy resolution function is finally derived in Section 3.4.

Estimation of the Shower Size
The general strategy for the reconstruction of air showers using the SD-750 array is similar to that used for the SD-1500 array which is detailed extensively in [26]. In this process, the arrival direction is obtained using the start times of signals, assuming either a plane or a curved shower front, as the degrees of freedom allow. The lateral distribution of the signal is then fitted to an empirically-chosen function to infer the size of the air shower, which is used as a surrogate for the primary energy. The reconstruction algorithm thus produces an estimate of the arrival direction and the size of the air shower via a log-likelihood minimization. This is primarily due to the instabilities in the wireless communications systems as well as periods where large fractions of the array were not functioning.
The lateral fall-off of the signal, ( ), with increasing distance, , to the shower axis in the shower plane is modeled with a lateral distribution function (LDF). The stochastic variations in the location and character of the leading interaction in the atmosphere result in shower-to-shower fluctuations of the longitudinal development that propagate onto fluctuations of the lateral profile, sampled at a fixed depth. Showers induced by identical primaries at the same energy and at the same incoming angle can thus be sampled at the ground level at a different stage of development. The LDF is consequently a quantity that varies on an event-by-event basis. However, the limited degrees of freedom, as well as the sparse sampling of the air-shower particles reaching the ground, prevent the reconstruction of all the parameters of the LDF for individual events. Instead, an average LDF, ( ) , is used in the reconstruction to infer the expected signal, ( opt ), that would be detected by a station located at a reference distance from the shower axis, opt [27,28]. This reference distance is chosen so as to minimize the fluctuations of the shower size, down to 7% in our case. The observed distribution of signals is then adjusted to ( ) by scaling the normalization, ( opt ), in the fitting procedure.
The reference distance, or optimal distance, opt , has been determined on an event-by-event basis by fitting the measured signals to different hypotheses for the fall-off of the LDF with distance to the core as in [28]. Via a fit of many power-law-like functions, the dispersion of signal expectations has been observed to be minimal at opt 450 m, which is primarily constrained by the geometry of the array. The expected signal at 450 m from the core, (450), has thus been chosen to define the shower-size estimate.
The functional shape chosen for the average LDF is a parabola in a log-log representation of ( ) as a function of the distance to the shower core, where = ln( /(450 m)), and and are two structure parameters. The overall steepness of the fall-off of the signal from the core is governed by , while the concave deviation from a power-law function is given by . The values of and have been obtained in a data-driven manner, by using a set of air-shower events with more than three stations, none of which have a saturated signal. The zenith angle and the shower size are used to trace the age dependence of the structure parameters based on the following parameterization in terms of the reduced variables = sec − 1.27 and = ln (450) − 5: For any specific set of values p = { , }, the reconstruction is then applied to calculate the following 2 -like quantity, globally to all events: The sum over stations is restricted to those with observed signals larger than 5 VEM to minimize the impact of upward fluctuations of the station signals far from the core and hence to avoid biases from trigger effects, and to stations more than 150 m away from the core. The uncertainty , is proportional to √︁ , [26]. tot is the total number of stations in all such events. The best-fit { , } values are collected in Table 2.

Correction of Attenuation Effects
There are two significant observational effects that impact the precision of the estimation of the shower size. Both of these effects are primarily a result of the variable slant depth that a shower must traverse before being detected with the SD. Since the mean atmospheric overburden is 875 g/cm 2 at the location of the Observatory, nearly all observed showers in the energy range considered in this analysis have already reached their maximum size and have started to attenuate [29]. Thus, an increase in the slant depth of a shower results in a more attenuated cascade at the ground, directly impacting the observed shower size.
The first observational effect is related to the changing weather at the Observatory. Fluctuations in the air pressure equate to changes in the local overburden and thus showers observed during periods of relatively high pressure result in an underestimated shower size. Similarly, the variations in the air density directly change the Molière radius which directly affects the spread of the shower particles. The increased lateral spread of the secondaries, or equivalently, the decrease in the density of particles on the ground, also leads to a systematically underestimated shower size. Both the airdensity and pressure have typical daily and yearly cycles that imprint similar cycles upon the estimation of the shower size.
The relationship between these two atmospheric parameters and the estimated shower sizes has been studied using events detected with the SD [30]. From this relationship, a model was constructed to scale the observed value of (450) to what would have been measured had the shower been instead observed at a time with the daily and yearly average atmosphere. When applying this correction to individual air showers, the measurements from the weather stations located at the FD sites are used. The values of (450) are scaled up or down according to these measurements, resulting in a shift of at most a few percent. The shower size is eventually the proxy of the air-shower energy, which is calibrated with events detected with the FD (see Section 3.3). Since the FD operates only at night when, in particular, the air density is relatively low, the scaling of (450) to a daily and yearly average atmosphere corrects for a 0.5% shift in the assigned energies.
The second observational effect is geometric, wherein showers arriving at larger zenith angles have to go through more atmosphere before reaching the SD. To correct for this effect, the Constant Intensity Cut (CIC) method [31] is used. The CIC method relies on the assumption that cosmic rays arrive isotropically, which is consistent with observations in the energy range considered [32]. The intensity is thus expected to be independent of arrival direction after correcting for the attenuation. Deviations from a constant behavior can thus be interpreted as being due to attenuation alone. Based on this property, the CIC method allows us to determine the attenuation curve as function of the zenith angle and therefore to infer a zenith-independent shower-size estimator.
We empirically chose a functional form which describes the relative amount of attenuation of the air shower, The scaling of this function is normalized to the attenuation of a shower arriving at 35 • by choosing = sin 2 35 • − sin 2 . For a given air shower, the observed shower size can be scaled using Eq. (6) to get the equivalent signal of a shower arriving with the reference zenith angle, 35 , via the relationship (450) = 35 CIC ( ). Isotropy implies that d /d sin 2 is constant. Thus, the shape of CIC ( ) is determined by finding the parameters and for which the CDF of events above (450) > cut CIC ( ) is linear in sin 2 using an Anderson-Darling test [33]. The parameter cut defines the size of a shower with = 35 • at which the CIC tuning is performed, the choice of which is described below.
Since the attenuation that a shower undergoes before being detected is related to the depth of shower maximum and the particle content, the shape of CIC ( ) is dependent on both the energy and the average mass of the primary particles at that energy. Further, this implies that a single choice of cut could introduce a mass and/or energy bias. Thus, Eq. (6) was extended to allow the polynomial coefficients, ∈ { , }, to be functions of (450) via ( (450)) = 0 + 1 + 2  where = lg( (450)/VEM). The function CIC ( , (450)) was tuned using an unbinned likelihood.
The fit was performed so as to guarantee equal intensity of the integral spectra using eight threshold values of cut between 10 and 70 VEM, evenly spaced in log-scale. These values were chosen to avoid triggering biases on the low end and the dwindling statistics on the high end. The best fit parameters are given in Table 3. The resulting 2D distribution of the number of events, in equal bins of sin 2 and lg 35 , is shown in Fig. 6, bottom panel. It is apparent that the number of events above any sin 2 value is equalized for any constant line for lg 35 0.7. The magnitude of the CIC correction is (−27 ± 4)% for vertical showers (depending on (450)) and +15% for a zenith angle of 40 • . The conversion of the shower size, corrected for attenuation, is based on a special set of showers, called golden hybrid events, which can be reconstructed independently by the FD and by the SD. The FD allows for a calorimetric estimate of the primary energy except for the contribution carried away by particles that reach the ground. The amount of this so-called invisible energy, 20% at 10 17 eV and 15% at 10 18 eV, has been evaluated using simulations [34] tuned to measurements at 10 18.3 eV so as to correct for the discrepancy in the muon content of simulated and observed showers [35]. The empirical relationship between the FD energy measurements, FD , and the corrected SD shower size, 35 , allows for the propagation of the FD energy scale to the SD events.

Energy Calibration of the Shower Size
FD events were selected based on quality and fiducial criteria aimed at guaranteeing a precise estimation of FD as well as at minimizing any acceptance biases towards light or heavy mass primaries introduced by the field of view of the FD telescopes. The cuts used for the energy calibration are similar to those described in [29,36]. They include the selection of data when the detectors are properly operational and the atmosphere properties like clouds coverage and the vertical aerosol depth are suitable for a good determination of the air-shower profile. A further quality selection includes requirements on the uncertainties of the energy assignment (less than 12%) and of the reconstruction of the depth at the maximum of the air-shower development (less than 40 g cm −2 ). A possible bias due to a selection dependency on the primary mass is avoided by using an energy dependent fiducial volume determined from data as in [29].
Restricting the data set to events with FD ≥ 10 17 eV, (to ensure that the SD is operating in the regime of full efficiency) there are 1980 golden-hybrid events available to establish the relationship between 35 and FD . Fourty-five events in the energy range between 10 16.5 eV and 10 17 eV are included in the likelihood as described in [37]. As 35 depends on the mass composition of the primary particles, the relation between 35 and FD , shown in Fig. 7, accounts for the trend of the composition change with energy inherently as the underlying mass distribution is directly sampled by the FD. Measurements of max suggest that this composition trend follows a logarithmic evolution up to an energy of 10 18.3 eV, beyond which the number of events available for this analysis is too small to affect the results in any way [36]. So we choose a power-law type relationship, which is expected from Monte-Carlo simulations in the case of a single logarithmic dependence of max with energy. The energy of an event with 35 = 1 VEM arriving at the reference angle, , and the logarithmic slope, , are fitted to the data by means of a maximum likelihood method which models the distribution of golden-hybrid events in the plane of energies and shower sizes. The use of these events allows us to infer and while accounting for the clustering of events in the range 10 17.4 to 10 17.7 eV observed in Fig. 7 due to the fall-off of the energy spectrum combined with the restrictive golden-hybrid acceptance for low-energy, dim showers. A comprehensive derivation of the likelihood function can be found in [37]. The probability density function entering the likelihood procedure, detailed in [37], is built by folding the cosmicray intensity, as observed through 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 intensity observed through the effective aperture of the telescopes (and thus to reduce reliance on mass assumptions), the observed distribution of events passing the cuts described above is used. The FD energy resolution, FD ( )/ FD , is typically between 6% and 8% [38]. It results from the statistical uncertainty arising from the fit to the longitudinal profile, the uncertainties in the detector response, the uncertainties in the models of the state of the atmosphere, and the uncertainties in the expected fluctuations from the invisible energy. The SD shower-size resolution, SD ( 35 )/ 35 , is, on the other hand, comprised of two terms, the detector sampling fluctuations, det ( 35 ), and the shower-to-shower fluctuations, sh ( 35 ). The former is obtained from the sum of the squares of the uncertainties from the reconstructed shower size and zenith angle, and from the attenuation-correction terms that make up the 35 assignment. The latter stem from the stochastic nature of both the depth of first interaction of the primary and the subsequent development of the particle cascade. This contribution thus depends on the CR mass composition and on the hadronic interactions in air showers. For this reason, the Table 4 The systematic uncertainties on the FD energy scale are given below. Lines with multiple entries represent the values at the low and high end of the considered energy range ( 10 17 and 10 19 eV, respectively).

Systematic Uncertainty
Absolute fluorescence yield 3.6% Atmosphere and scattering 2 to 6% FD Calibration 10% Longitudinal profile reconstruction 7 to 5.5% Invisible energy 3 to 1.5% derivation of and follows a two-step procedure. A first iteration of the fit is carried out by using an educated guess for sh ( 35 ), as expected from Monte-Carlo simulations for a mass-composition scenario compatible with data [29]. The total resolution SD ( 35 )/ 35 is then extracted from data as explained next in Section 3.4 and used in a second iteration.
The resulting relationship is shown as the red line in Fig. 7 with best-fit parameters such that = (13.2 ± 0.3) PeV and = 1.002±0.006. The goodness of the fit is supported by the 2 /NDOF = 2120/1978 ( = 0.013). We use these values of and to calibrate the shower sizes in terms of energies by defining the SD estimator of energies, SD , according to Eq. (7). The SD energy scale is set by the calibration procedure and thus it inherits the and calibration-parameters uncertainties and the FD energy-scale uncertainties, listed in Table 4. The systematic uncertainty, after addition in quadrature, of the energy scale is about 14% and is almost energy independent. The energy independence is a consequence of the 10% uncertainty of the FD calibration, which is the dominant contribution.

Resolution Function of the SD-750 Array
The SD resolution as a function of energy is needed in several steps of the analysis. In the regime of full efficiency, it can be considered as a Gaussian function centered on the true energy, the width of which reflects the statistical uncertainty associated with the detection and reconstruction processes on one hand, and the stochastic development of the particle cascade on the other hand. The combination of the two can be estimated for the golden hybrid events, thus allowing us to account for the contribution of the shower-to-shower fluctuations in a data-driven way.
Each event observed by the SD and FD results in two independent measurements of the air-shower energy, SD and FD , respectively. Unlike for the SD, the FD directly provides a view of the shower development so a total energy resolution, FD ( ), can be estimated for each of the golden hybrid events. Using the known FD ( ), the resolution of SD can be determined by studying the distribution of the ratio of the two energy measurements.   The total SD energy resolution, as calculated using the golden hybrid events (red circles) is shown in bins with equal statistics. The parameterization of the resolution is shown by the solid blue line and the corresponding 68% confidence interval in dashed lines. The energy resolution, calculated using mass-weighted MC air showers (gray squares), is shown as a verification of the method.
For two independent, Gaussian-distributed random variables, and , their ratio, = / , produces a ratio distribution that depends on the means ( , ) and standard deviations ( , ) of the two variables, PDF( ; , , , ). Likewise, the ratio of the two energy measurements, = SD / FD , follows such a distribution to first order. Because the FD sets the energy scale of the Observatory, there is inherently no bias in the energy measurements with respect to its own scale and thus, on average, FD ( ) = 1. Using the golden hybrid data set, the ratio distribution was fit in an unbinned likelihood analysis, PDF( ; SD ( ), 1, SD ( ), FD ( )).
An example of the measured energy-ratio distributions is shown in Fig. 8 with the fitted curve overlaid on the data points. Carrying out the fit in different energy bins, the SD resolution, shown by the red points in Fig. 9, is represented by, The corresponding curve is overlaid in blue, bracketed by the 68% confidence region.
To measure the spectrum above the 10 17 eV threshold, the knowledge of the resolution function, which induces binto-bin migration of events, and of the detection efficiency are also required for energies below this threshold. As a verification, particularly in the energy region where Eq. (8) is extrapolated, a Monte-Carlo analysis was performed. A set of 325,000 CORSIKA [39] air showers were used, consisting of proton, helium, oxygen, and iron primaries with energies above 10 16 eV. EPOS-LHC [40] was used as the hadronic interaction model. The air showers were run through the full SD simulation and reconstruction algorithms. The events were weighted based on the primary mass according to the Global Spline Fit (GSF) model [41] to account for the changing mass-evolution near the second knee and ankle. The reconstructed values of (450) were corrected by applying the energy-dependent CIC method to obtain values for 35 and these values were then calibrated against the Monte-Carlo energies. During the calibration, a further weighting was performed based on the energy distribution of golden hybrid events to account for the hybrid detection efficiency. Following the calibration procedure, each MC event was assigned an energy in the FD energy scale (i.e. MC → 35 → FD ).
The SD energy resolution was calculated using the massweighted simulations and is shown in gray squares in Fig. 9. Indeed, the simulated and measured SD resolutions show a similar trend and agree to within the uncertainties, supporting the golden hybrid method.
In the energy region at-and-below 10 17 eV, systematic effects also enter into play on the energy estimate. An energydependent offset, a bias, is thus expected in the resolution function for several reasons: 1. The application of the trigger below threshold, combined with the finite energy resolution, cause an overestimate of the shower size, on average, which is then propagated to the energy assignment. 2. The linear relationship assumed in Eq. (7) cannot account for a possible sudden change in the evolution of the masscomposition with energy. Such a change would require a broken power law for the energy calibration relationship. 3. In the energy range where the SD is not fully efficient, the SD efficiency is larger for light primary nuclei, thus preventing a fair sampling of 35 values over the underlying mass distribution.
Because there is an insufficient number of FD events which pass the fiducial cuts below 10 17 eV, the bias was Bias parameterization Fig. 10 The bias of the energy assignment for the SD-750 was studied using Monte Carlo simulations, weighted according to the GSF model [41]. The ratio of the assigned and expected values as a function of energy are shown (red circles) along with the parameterization (blue line) given in Eq. (9). Table 5 Best-fit parameters for the relative energy bias of the SD-750, SD ( ), given in Eq. (9). characterized, using the same air-shower simulations as used for the resolution cross-check. The remaining relative energy bias is shown in Fig. 10. The ratio between the reconstructed and expected values are shown as the red points as a function of FD . A larger bias of 20% is seen at low energies, where upward fluctuations are necessarily selected by the triggering conditions. In the range considered for the energy spectrum, > 10 17 eV, the bias is 3% or less. To complete the description of the SD resolution function, the relative bias was fit to an empirical function,

Parameter Value Uncertainty
The corresponding best fit parameters (blue line in Fig. 10) are given in Table 5.

Measurement of the Energy Spectrum
To build the energy spectrum from the reconstructed energy distribution, we need to correct the raw spectrum, obtained as raw = /(EΔ ), for the bin-to-bin migrations of events due to the finite accuracy with which the energies are assigned. The energy bins are chosen to be regularly sized in decimal logarithm, Δ lg = 0.1, commensurate with the energy resolution. The level of migration is driven by the resolution function, the detection efficiency in the energy range just below the threshold energy, and the steepness of the spectrum. To correct for these effects, we use the binby-bin correction approach presented in [21]. It consists of folding the detector effects into a proposed spectrum function, ( , k), with free parameters, k, such that the result describes the set of the observed number of events . The set of expectations, , is obtained as (k) = (k), where the coefficients (reported in a matrix format in the Supplementary material) describe the bin-to-bin migrations, and where are the expectations in the case of an ideal detector obtained by integrating the proposed spectrum over and +Δ scaled by E. The optimal set of free parameters, k, is inferred by minimizing a log-likelihood function built from the Poisson probabilities to observe events when (k) are expected.
To choose the proposed function, we plot in Fig. 11 the residuals (red dots) of the SD-750 raw spectrum with respect to a reference function, ref ( ), that fits the SD-1500 spectrum below the ankle energy down to the SD-1500 threshold energy, 10 18.4 eV. A re-binning was applied at and above 10 19 eV to avoid too large statistical fluctuations.
The reference function in this energy range, as reported in [21], is  a transition is also observed, with much lower sensitivity, using data from the SD-750 array. Below 10 18.7 eV and down to 10 17.4 eV, one can see a shift of the raw SD-750 spectrum compared to ref ( ). This is expected from a combination of primarily the resolution effects to be unfolded and of a possible mismatch, within the energy-dependent budget of uncorrelated uncertainties, of the SD-1500 and SD-750 SD energy scales. Below 10 17.4 eV, a slight roll-off begins. Overall, these residuals are suggestive of a power-law function to describe the data leading up to the ankle energy where the spectrum hardens, with a gradually changing spectral index over the lowest energies studied. Consequently, the proposed function is chosen as three power laws with transitions occurring over adjustable energy ranges, ( , k) = 0 10 17 eV with = + 1. The normalization factor 0 , the three spectral indices , and the transition parameter 01 constitute the free parameters in k. The transition parameter 12 , constrained with much more sensitivity using data from the SD-1500, is fixed at 12 = 0.05 [21].
Combining all the ingredients at our disposal, we obtain the final estimate of the spectrum, , unfolded for the effects of the response of the detector and shown in Fig. 12. It is obtained as where the and coefficients are estimated using the bestfit parametersk. Their ratios define the bin-by-bin corrections used to produce the unfolded spectrum. The correction applied extends from 0.84 at 10 17 eV to 0.99 around the ankle (see Appendix B). The best-fit spectral parameters are Table 6 Best-fit values of the spectral parameters (Eq. (11)). The parameter 12 is fixed to the value constrained in [21]. Note that the parameters 0 and 01 correspond to features below the measured energy region and are treated only as aspects of the unfolding fixed to their best-fit values to infer the uncertainties of the measured spectral parameters.

Parameter
Value ± stat ± syst 0 /(km 2 yr sr eV) (1.09 ± 0.04 ± 0.28) ×10 −13 01 0.49 ± 0.07 ± 0.34 1 3.34 ± 0.02 ± 0.09 12  reported in Table 6, while the statistical correlations between the parameters are detailed in Appendix B ( Table 9). The goodness-of-fit of the forward-folding procedure is attested by the deviance of 15.9, which, if considered to follow the C statistics [42], can be compared to the expectation of 16.2 ± 5.6 to yield a -value of 0.50.
The fitting function is shown in Fig. 13, superimposed to the spectrum scaled by 2.6 , allowing one to better appreciate its characteristics, from the turn-over at around 10 17 eV up to a few 10 19 eV, thus including the ankle. The turnover is observed with a very large exposure, unprecedented at such energies. However, as indicated by the magnitude of the transition parameter, 01 0.49, the change of the spectral index occurs over an extended Δ lg 0.5 energy range, so that the spectral index 0 cannot be observed but only indirectly inferred. Also, the value of the energy break, 01 1.24×10 17 eV, turns out to be close to the threshold en-Note that the -value for a proposed function which does not include a transition from 0 to 1 can be rejected with more than 20 confidence. ergy. These two facts thus imply that, while a spectral break is found beyond any doubt, it cannot wholly be characterised, as only the higher energy portion is actually observed. Consequently, the fit values describing 01 and 0 are not to be considered as true measurements but as necessary parameters in the fit function, the statistical resolutions of which are on the order of 35%. Once we infer their best-fit values, we use these values as "external parameters" to estimate the uncertainties of the other spectral parameters. This procedure gives rise to an increase of the systematic uncertainties, but is necessary as 01 and 0 are not directly observed. Beyond the smooth turn-over around 01 , the intensity can be described by a power-law shape as ( ) ∝ − 1 , up to 12 = (3.9 ± 0.8) ×10 18 eV, the ankle energy, the value of which is within 1.4 of that found with the much larger exposure of the SD-1500 measurement of the spectrum, namely (5.0 ± 0.1)×10 18 eV. Also the value of 1 = 3.34 ± 0.02 is within 1.8 of that obtained with the SD-1500 between 10 18.4 and 10 18.7 eV (3.29 ± 0.02).
The characteristics of the measured spectrum can also be studied by looking at the evolution of the spectral index as a function of energy, ( ). Rather than relying on the empirically chosen unfolding function, this slope parameter can be directly fit using the values calculated in ( ). Power-law fits were performed for a sliding window of width Δ lg = 0.3. The resulting estimations of the so obtained spectral indexes are shown in Fig. 14. The values of the spectral index fits present a consistent picture of the evolution. Beginning at the lowest energies shown, ( ) increases first quite rapidly, finally approaching a value of 3.3 leading up to the ankle asymptotically. Unsurprisingly, this is the value found for The systematic uncertainties that affect the measurement of the spectrum are dominated by the overall uncertainty of the energy scale, detailed in [43], and is, itself, dominated by the absolute calibration of the fluorescence telescopes (10%). The total uncertainty in the energy scale is / = 14%. Once propagated, the steepness of the spectrum as a function of energy amplifies this uncertainty, roughly as / = ( 1 − 1) / , resulting in a total flux uncertainty of / 35%. However, for a more exact calculation of the uncertainty, the energies of the individual events were shifted by ±14% and the unfolding procedure was repeated. The result is shown as dashed red lines in Fig. 15.
Beyond that of the energy scale, the additional uncertainties are subdominant but are important to understand as they have energy dependence and some are uncorrelated with other flux measurements made at the Observatory. Such knowledge is particularly important for the combination of the two SD spectra presented later in Section 5. The most relevant of these energy-dependent uncertainties is associated with the procedure of the forward-folding itself. The uncertainties in the resolution function and in the detection efficiency all contribute a component to the overall unfolding uncertainty. The forward-folding process was hence repeated by shifting, within the statistical uncertainties, the parameterizations of the energy resolution (Eq. (8)) and efficiency parameterization, and by bracketing the bias with the pure proton/iron mass primaries below full efficiency. The impact of the resolution uncertainties on the unfolding procedure is the larger, in particular at the highest energies. On the other hand, the energy bias and reduced efficiency below 10 17 eV only impacts the first few bins. These various components are summed in quadrature and are shown by the dotted blue line in Fig. 15. These influences are clearly seen to impact the spectrum by <4%. The last significant uncertainty in the flux is related to the calculation of the geometric exposure of the array. This quantity has been previously studied and is 4% for the SD-750 which directly translates to a 4% energy-independent shift in the flux [24].
The resulting systematic uncertainties of the spectral parameters are given in Table 6. For completeness, beyond the summary information provided by the spectrum parameterization, the correlation matrix of the energy spectrum is given in the Supplementary material. It is obtained by repeating the analysis on a large number of data sets, sampling randomly the systematic uncertainties listed above.

The Combined SD-750 and SD-1500 Energy Spectrum
The spectrum obtained in Section 4 extends down to 10 17 eV and at the high-energy end overlaps with the one recently reported in [21] using the SD-1500 array. The two spectra are superimposed in Fig. 16. Beyond the overall consistency observed between the two measurements, a combination of them is desirable to gather the information in a single energy spectrum above 10 17 eV obtained with data from both the SD-750 and the SD-1500 of the Pierre Auger Observatory. We present below such a combination considering adjustable re-scaling factors in exposures, E, and SD energy scales, SD , within uncorrelated uncertainties. The combination is carried out using the same bin-bybin correction approach as in Section 4. The joint likelihood function, L (s, E, SD ), is built from the product of the individual Poissonian likelihoods pertaining to the two SD measurements, L 750 and L 1500 . These two individual likeli-hoods share the same proposed function, with = + 1 and 0 = 10 18.5 eV. As in [21], the transition parameters 12 , 23 and 34 are fixed to 0.05. In this way, the same parameters s are used during the minimisation process to calculate the set of expectations (s, E, SD ) of the two arrays. For each array, a change of the associated exposure E → E + E impacts the coefficients accordingly, while a change in energy scale SD → SD + SD impacts as well the observed number of events in each bin. Additional likelihood factors, L E and L SD , are thus required to control the changes of the exposure and of the energy-scale within their uncorrelated uncertainties. The likelihood factors described below account for E and SD changes associated with the SD-750 only. We have checked that allowing additional free parameters, such as the E corresponding to the SD-1500, does not improve the deviance of the best fit by more than one unit, and thus their introduction is not supported by the data. Both likelihood factors are described by Gaussian distributions with a spread given by the uncertainty pertaining to the exposure and to the energy-scale. The joint likelihood function reads then as L (s, E, SD ) = L 750 × L 1500 × L E × L SD .
The allowed change of exposure, E, is guided by the systematic uncertainties in the SD-750 exposure, E /E = 4%. Hence, the constraining term for any change in the SD-750 exposure reads, dropping constant terms, as Likewise, uncertainties in and , and , translate into uncertainties in the SD-750 energy scale. Statistical contributions stem from the energy calibration of 35 , which are by essence uncorrelated to those of the SD-1500. Other uncorrelated contributions of the systematic uncertainties from the FD energy scales propagated to the SD-1500 and SD-750 could enter into play. The magnitude of such systematics, syst , is difficult to quantify. By testing several values for syst , we have checked, however, that such contributions have a negligible impact on the combined spectrum. Hence, the constraining term for any change in energy scale can be considered to stem from statistical uncertainties only and reads as / eV E  13)) is overlaid in red along with the one sigma error band in gray. Table 7 Best-fit values of the combined spectral parameters (Eq. (13)). The parameter 12 , 23 and 34 are fixed to the value constrained in [21]. Note that the parameters 0 and 01 correspond to features below the measured energy region and should be treated only as aspects of the combination.

Parameter
Value ± stat ± syst 0 / (km 2 yr sr eV) ( is shown in Fig. 17. Here, the observed number of events 750 in each bin is calculated at the re-scaled energies, while the effective exposure, E eff , is the shifted one of the SD-750 in the energy range where ,1500 = 0, the one of the SD-1500 in the energy range where ,750 = 0, and the sum E 750 + E + E 1500 in the overlapping energy range. The set of spectral parameters are collected in Table 7, while the corresponding correlation matrix is reported in Appendix B (Table 11) for E, and fixed to their best-fit values. The change in exposure is E/E = +1.4%, while the one in energy scale follows from / = −2.5% and / = +0.8%. The goodness-of-fit is evidenced by a deviance of 37.2 for an expected value of 32 ± 8. We also note that the parameters describing the spectral shape are in agreement with those of the two individual spectra from the SD arrays.
The impact of the systematic uncertainties, dominated by those in the energy scale, on the spectral parameters are reported in Table 7. For completeness, beyond the summary information provided by the spectrum parameterization, the correlation matrix of the energy spectrum itself is also given in the Supplementary material.  [44], GAMMA [45], IceTop [9], KASCADE-Grande [46], TALE [10], Tien Shan [47], Tibet-III [48], Tunka-133 [11], Yakutsk [49]. The experiments that set their energy scale using calorimetric observations are indicated by solid colored markers while those with an energy scale based entirely on simulations are shown by gray markers.

Discussion
We have presented here a measurement of the CR spectrum in the energy range between the second knee and the ankle, which is covered with high statistics by the SD-750, including 560,000 events with zenith angles up to 40 • and energies above 10 17 eV. The measurement includes a total exposure of 105 km 2 sr yr and an energy scale set by calorimetric observations from the FD telescopes. We note a significant change in the spectral index and with a width that is much broader than that of the ankle feature.
Such a change has been observed by a number of other experiments, and via various detection methods. Most notably, the nature of this feature was linked to a softening of the heavy-mass primaries beginning at 10 16.9 eV by the KASCADE-Grande experiment, leading to the moniker iron knee [8]. Additional analyses by the Tunka-133 [50] and IceCube [9] collaborations have given further evidence that high-mass particles are dominant near 10 17 eV and thus that it is their decline that largely defines the shape of the allparticle spectrum. The hypothesis is also supported by a preliminary study of the distributions of the depths of the shower maximum, max , measured at the Auger Observatory [36,51]. These have been parametrized according to the hadronic models EPOS-LHC [40], QGSJetII-04 [52] and Sibyll2.3 [53]. From these parametrizations, the evolution over energy of the fractions of different mass groups, from protons to Fe-nuclei, has been derived. From all three models, a fall-off of the Fe component above 10 17 eV is inferred. The consistency of all these observations strongly supports a scenario of Galactic CRs characterised by a rigidity-dependent maximum acceleration energy for particles with charge , namely max ( ) proton max , to explain the knee structures. The measurements of the all-particle flux from various experiments [9][10][11][44][45][46][47][48][49] in the energy region surrounding the second knee are shown in Fig. 18. Experiments which set their energy scale using calorimetric measurements are plotted using colored markers (Auger SD-750, TA TALE, TUNKA-133, Yakutsk) while the measurements shown in gray markers represent MC-based energy assignments. The spread between various experiments is statistically significant. However, all these measurements are consistent with the SD-750 spectrum within the 14% energy scale systematic uncertainty. Understanding the nature of the off-sets in the energy scales is beyond the scope of this paper. However, we note that the TALE spectrum agrees rather well with the SD-750 spectrum, offset by 5 to 6% in energy. The agreement is notable given that at-and-above the ankle, an energy scale off-set of around 11% is required to bring the spectral measurements with SD-1500 of the Auger Observatory and the SD of the Telescope Array into agreement [54].
Additionally, we have presented a robust method to combine energy spectra. Using the result from the SD-750 and a previously reported measurement using the SD-1500, a unified SD spectrum was calculated by combining the respec-tive observed fluxes, energy resolutions, and exposures. The result has partial coverage of the second knee and full coverage of the ankle, an additional inflection at 1.4×10 19 eV, and the suppression. This procedure is applied to spectra inferred from a single detector type (i.e. water-Cherenkov detectors), but can be used for the combination of any spectral measurements for which the uncorrelated uncertainties can be estimated.
The impressive regularity of the all-particle spectrum observed in the energy region between the second knee and the ankle can hide an underlying intertwining of different astrophysical phenomena, which might be exposed by looking at the spectrum of different primary elements. In the future, further measurements will allow separation of the intensities due to the different components. On the one hand, max values will be determined down to 10 17 eV using the three HEAT telescopes. On the other hand, the determination of the muon component of EAS above 10 17 eV will be possible using the new array of underground muon detectors [35], co-located with the SD-750. This will help us in studying whether the origin of the second knee stems from, for instance, the steep fall-off of an iron component, as expected for Galactic CRs characterized by a rigidity-dependent maximum acceleration energy for particles with charge , namely max ( ) proton max . In addition, we will be able to extend the measurement of the energy spectrum below 10 17 eV with a denser array of 433 m-spaced detectors and with the analysis of the Cherenkov light in FD events [55]. The extension will allow us to lower the threshold and to further explore the second-knee region in more detail.
The ToTd and MoPS triggers were designed to be insensitive to atmospheric muons such that they enable the detection of small electromagnetic signals from air showers. The typical morphology of a waveform from a ∼GeV muon is a 150 ns ( 6 ADC bins) pulse with an amplitude of 1 VEM , where VEM is the maximum amplitude of a signal created by a muon that traverses the water volume vertically [23]. Thus, the ToTd and MoPS algorithms are used to look for signals that do not fit this criteria.
The two additional triggers build upon the ToT trigger in two ways, applying more sophisticated analyses to the signal waveform. They are aimed at further suppressing the muon background so as to enhance the sensitivity to pure electromagnetic signals, which are generally smaller.
The ToTd trigger uses the typical decay time of Cherenkov light inside the water volume, = 67 ns, to deconvolve the exponential tail of the pulses before applying the ToT condition. This has the effect of reducing the influence of muons in the trigger, since the typical signal from a muon, with fast rise time and ≈60 ns decay constant, is compressed into one or two time bins. The exponential tail of the signal is deconvolved using where is the signal in the -th time-bin and Δ = 25 ns is the ADC bin-width. For an exponential decay with the mean decay time, the deconvolved values, , would be zero. However for an exponential decay with statistical noise that is proportional to √ , the set { } would exponentially decrease with an increased decay length = 2 . After performing the deconvolution in Eq. (A.1), the trigger is satisfied if ≥13 ADC bins (≥325 ns) are above 0.2 VEM , in coincidence between two of the three PMTs, within a sliding 3 µs (120 bin) time window. An example of a waveform which passes the ToTd trigger and its deconvolution are shown in the top two plots of Fig. 19. Only 11 bins are above 0.2 VEM in the original waveform such that it cannot pass the traditional TOT algorithm. However the deconvolution has the 13 bins required to be above the threshold.
The second, MoPS, counts the number of instances, in a sliding 3 µs window, in which there is a monotonic increase of the signal amplitude. Each such instance of successive increases in the digitized waveform is what we define as a positive step. For each positive step, the total vertical increase, , must be above that of typical noise, and below the characteristic amplitude of a vertical muon, namely 3 < ≤ 31. If more than four of the positive-step instances fall within this range, the trigger condition is satisfied. An example of a waveform which passes the MoPS trigger is shown in the bottom plot of Fig. 19.

Appendix B: Spectrum Data
We report in this appendix several data of interest. Note that more can be found in the Supplemental Material in electronic format. The scaling factor that has been applied to the raw spectrum to produce the unfolded spectrum (see Eq. (12)) and the statistical uncertainty.
The bin migration is corrected to produce the unfolded spectrum. The magnitude of the correction factor, as described by Eq. (12), is shown in Fig. 20 along with the statistical uncertainty band. The energy spectrum of the SD-750 array is reported in Table 8 and the correlation matrix of the spectral parameters at the nominal energy scale in Table 9 (statistical uncertainties). Finally, the combined energy spectrum is reported in Table 10 and the correlation matrix of the For example, four bins with ≤ +1 ≤ +2 ≤ +3 is considered one positive step, not three positive steps. spectral parameters at the nominal energy scale in Table 11 (statistical uncertainties).  Table 9 Elements of the correlation matrix (statistical uncertainties) of the spectral parameters describing the SD-750 energy spectrum at the nominal energy scale.   Table 11 Elements of the correlation matrix (statistical uncertainties) of the spectral parameters describing the combined SD energy spectrum.