The energy spectrum of cosmic rays beyond the turn-down around 1017 eV as measured with the surface detector of the Pierre Auger Observatory

Eur. Phys. J. C (2021) 81:966 https://doi.org/10.1140/epjc/s10052-021-09700-w Regular Article - Theoretical Physics The energy spectrum of cosmic rays beyond the turn-down around 1017 eV as measured with the surface detector of the Pierre Auger Observatory Pierre Auger Collaboration The Pierre Auger Observatory, Av. San Martín Norte 306, 5613 Malargüe, Mendoza, Argentina; URL: http://www.auger.org Received: 10 May 2021 / Accepted: 29 September 2021 © The Author(s) 2021 Abstract We present a measurement of the cosmic-ray spectrumabove 100PeVusingthepartofthesurfacedetector 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 second-knee feature. The spectrum is then combined with that of the 1500m 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. 1 Introduction The steepening of the energy spectrum of cosmic rays (CRs) at around 1015.5eV, 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, remainsimportantupto the kneeenergyafterwhichtheheavier elements gradually take over up to a few 1017 eV [2–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 that the shock waves sweep out, and particles are accelerSupplementary Information The online version contains supplementary material available at https://doi.org/10.1140/epjc/ s10052-021-09700-w. e-mail: [email protected] ated 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 ≃1016.9eV. Such a bending has been observed in several experiments at a similar energy, referred to as the “second knee” or “iron knee” [8–11]. The recent observations of gamma rays of a few 1014 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 1015 eV, and above, in the Galaxy. A dozen of sources emitting gamma rays up to 1015 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 1017 eV, the spectrum steepens in the interval leading up to the “ankle” energy, ∼5×1018 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 1017.3eV 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]. WhiletheconceptthattheGalactic-to-extragalactictransition occurs somewhere between 1017 eV and a few 1018 eV is 0123456789().: V,-vol 123 966 Page 2 of 25 Eur. Phys. J. C (2021) 81:966 Fig. 1 The layout of the SD and FD of the Pierre Auger Observatory are shown above. The respective fields of view of the five FD sites are shown in blue and orange. The 1600 SD locations which make up the SD-1500 are shown in black while the stations which belong only to the SD-750 and the boarder of this sub-array are highlighted in cyan well-accredited,afullunderstandingofhowitoccursishence 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. TheaimofthispaperistopresentameasurementoftheCR spectrum from 1017 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 1400m above sea level at a latitude of 35.2◦S, so that the mean atmospheric overburden is 875g/cm2. 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–390nm), the light is reflected off a spherical mirror onto a camera of 22×20 hexagonal, 45.6mm, 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 1500m and totaling approximately 3000km2of effective area. The detectors of the SD-750 are instead spread out by 750m, yielding an effective area of 24km2. SD-750 and SD-1500 include identical WCDs, cylindrical tanks of pure water with a 10m2base and a height of 1.2m. 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, reflectedviaaTyvekliningintothePMTs,anddigitizedusing 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 1017 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 1017 eV in Sect. 2; in Sect. 3, we describe how we reconstruct the showers up to determining the energy; we illustrate in Sect. 4the 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 usingtheAugerObservatory from 1017 eVupwardsin Sect.5 and discuss it in the context of other measurements in Sect. 6. 123 Eur. Phys. J. C (2021) 81:966 Page 3 of 25 966 2 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 1017.2eV down to 1017 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 SD1500, 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 (Sect. 2.1). We then illustrate their effect on the data, at the level of the amplitude of detected signals (Sect. 2.2) and on the timing of detected signals in connection with the event selection (Sect. 2.3). Finally we describe the energy at which acceptance is 100% (Sect. 2.4). A more detailed description of the trigger algorithms can be found in Appendix A. 2.1 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.1 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 3kHz, 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. 1The 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.3Hz, compared to 0.6Hz of ToT and 20Hz 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. 2.2 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 whichis dominatedbythebackgroundmuonsthatarecreated inrelativelylowenergyairshowers.Byremaininginsensitive 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 123 966 Page 4 of 25 Eur. Phys. J. C (2021) 81:966 Fig. 2 Distribution of the signal sizes at individual stations which pass the TH and ToT triggers (solid black) and signals which pass only the ToTd and/or MoPS triggers (dashed red) Fig. 3 The increase in station multiplicity when including the ToTd and MoPS triggers versus the original multiplicity with only ToT and TH. The black circles show the median increase in that multiplicity bin 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 ≃7VEM. 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 5VEM. That the high-signal tail of this distribution ends near 10VEM 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 Fig. 4 Distributions of start times with respect to a plane front for stations that pass the ToT and TH algorithms, in blue and in green, respectively. The signals due to ToTd and MoPS are shown in red. Positive residuals correspond to a delay with respect to the plane wave expectation 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. 2.3 Effects of ToTd and MoPS on signal timing The increased responsiveness of the ToTd and MoPS algorithms tosmaller 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 Fig. 4for 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 123 Eur. Phys. J. C (2021) 81:966 Page 5 of 25 966 Table 1 Temporal window limits tlow and thigh used to remove stations from an event, for each station-level trigger algorithm Trigger type tlow (ns) thigh (ns) ToT −397 1454 ToTd −468 2285 MoPS −477 2883 TH −485 1379 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 1VEM 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, weidentifythecombinationofthreestationswheretheyform atriangle,inwhichatleasttwolegsare750mlong,andwhere 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, t, is within a fixed window, tlow <t<thigh. Motivated by the differing time distributions, updated tlow and thigh 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, tlow and thigh, were chosen such that the middle 99% of the distribution was kept. The trigger-wise limits are summarized in Table 1. 2.4 Effect of the ToTd and MoPS on the energy above which acceptance is fully-efficient Mostrelevanttothemeasurementofthespectrumisthedetermination 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 1016.8eV, the lateral trigger probability (LTP), Fig. 5 The detection efficiency of the SD-750 for air showers with θ<40◦is shownfortheoriginal(dashedred) and expanded(solidblue) station-level trigger sets with bands indicating the systematic uncertainties. The trigger efficiency was determined using data above 1016.8eV and is extrapolated below this energy (shown in gray) the chance that a shower will produce a given SD trigger 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 1016.5 and1018 eV.Eachstationonthearraywasrandomlytriggered 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 shownasasolidbluelineinFig.5asafunctionofenergy.The detectionefficiencybecomesalmostunity(>98%)ataround 1017 eV.2For 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 1017.2eV. A description for the detection efficiency, (E), below 1017 eV, will be important for unfolding the detector effects close to the threshold energy (see Sect. 4). This quantity was fit using the results of the LTP simulations with θ<40◦and 2The 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. 123 966 Page 6 of 25 Eur. Phys. J. C (2021) 81:966 is well-parameterized by (E)=1 21+erf lg(E/eV)−μ σ,(1) where erf(x)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 E>1017 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 Sect. 3.3). The corresponding exposure, E,afterremovalof time periods when the array was unstable3(<2% of the total) is E=(105 ±4)km2sryr. 3 Energy measurements with the SD-750 Inthissection,themethodfortheestimationoftheair-shower 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 Sect. 3.1 after which the corrections for attenuation effects are presented in Sect. 3.2. The energy calibration of the shower size after correction for attenuation is presented in Sect. 3.3. The energy resolution function is finally derived in Sect. 3.4. 3.1 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 SD1500 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. 3This 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, S(r), with increasing distance, r, to the shower axis in the shower plane is modeledwithalateral distributionfunction (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,preventthe reconstruction of all the parametersofthe LDF for individual events. Instead, an average LDF, S(r), is used in the reconstruction to infer the expected signal, S(ropt), that would be detected by a station located at a reference distance from the shower axis, ropt [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 S(r)by scaling the normalization, S(ropt), in the fitting procedure. The reference distance, or optimal distance, ropt, 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 ropt ≃450m, which is primarily constrained by the geometry of the array. The expected signal at 450m from the core, S(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 S(r)as a function of the distance to the shower core, lnS(r)=ln S(450)+βρ+γρ 2,(2) where ρ=ln(r/(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 t=sec θ−1.27 and u=ln S(450)−5: β=(β0+β1t+β2t2)(1+β3u), (3) γ=γ0+γ1u.(4) For any specific set of values p={βi,γ i}, the reconstruction is then applied to calculate the following χ2-like quantity, 123 Eur. Phys. J. C (2021) 81:966 Page 7 of 25 966 Table 2 Best-fit {βi,γi} values defining the structure parameters of the LDF Parameter Value β02.95 ±0.02 β1−1.0±0.2 β20.7±0.2 β30.02 ±0.01 γ00.26 ±0.09 γ1−0.02 ±0.01 globally to all events: Q2(p)=1 Ntot Nevents  k=1 Nk  j=1 (Sk,j−S(rj,p))2 σ2 k,j .(5) The sum over Nkstations is restricted to those with observed signals larger than 5VEM 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 150m away from the core. The uncertainty σk,jis proportional to Sk,j[26]. Ntot is the total number of stations in all such events. The best-fit {βi,γi} values are collected in Table 2. 3.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 875g/cm2at 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 inamoreattenuatedcascadeattheground,directlyimpacting 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 anunderestimatedshowersize.Similarly,thevariationsinthe 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 air-density 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 S(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 S(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 Sect. 3.3). Since the FD operates only at night when, in particular, the air density is relatively low, the scaling of S(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, fCIC(θ) =1+ax +bx2.(6) The scaling of this function is normalized to the attenuation ofa showerarrivingat35◦by choosing x=sin235◦−sin2θ. Foragivenair shower, theobservedshowersizecanbescaled using Eq. (6) to get the equivalent signal of a shower arriving with the reference zenith angle, S35, via the relationship S(450)=S35 fCIC(θ). Isotropy implies that dN/dsin 2θis constant. Thus, the shape of fCIC(θ) is determined by finding the parameters aand bfor which the CDF of events above S(450)> Scut fCIC(θ) is linear in sin2θusing an Anderson-Darling test [33]. The parameter Scut defines the size of a shower with θ=35◦at which the CIC tuning is performed, the choice of which is described below. Sincetheattenuationthatashowerundergoesbeforebeing detected is related to the depth of shower maximum and the particle content, the shape of fCIC(θ) 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 Scut could introduce a mass and/or energy bias. Thus, Eq. (6)was extended to allow the polynomial coefficients, k∈{a,b}, 123 966 Page 8 of 25 Eur. Phys. J. C (2021) 81:966 Fig. 6 Top: histogram of reconstructed shower sizes and zenith angles. The solid black line represents the shape of fCIC at 10VEM. Bottom: same distribution but as a function of corrected shower size, S35,and zenith angle. The dashed black line indicates the mapping of the solid black line in the top figure after inverting the effects of the CIC correction Table 3 The energy dependence of the CIC parameters (Eq. (6)) are given below k0k1k2 a2.42 −0.886 0.268 b−4.56 5.61 −2.47 to be functions of S(450)via k(S(450)) =k0+k1y+k2y2 where y=lg(S(450)/VEM). The function fCIC(θ, S(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 Scut between 10 and 70VEM, 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 sin2θand lg S35,is shown in Fig. 6, bottom panel. It is apparent that the number of events above any sin2θvalue is equalized for any constant line for lg S35 0.7. The magnitude of the CIC correction is (−27 ±4)% for vertical showers (depending on S(450)) and +15% for a zenith angle of 40◦. Fig. 7 Correlation between the SD shower-size estimator, S35,andthe reconstructed FD energy, EFD, for the selected hybrid events 3.3 Energy calibration of the shower size 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 1017 eV and ≃15% at 1018 eV, has been evaluated using simulations [34] tuned to measurementsat1018.3eVsoastocorrectforthediscrepancy in the muon content of simulated and observed showers [35]. The empirical relationship between the FD energy measurements, EFD, and the corrected SD shower size, S35, allows for the propagation of the FD energy scale to the SD events. FD events were selected based on quality and fiducial criteria aimed at guaranteeing a precise estimation of EFD 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 EFD ≥1017 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 S35 and EFD. Fourty-five 123 Eur. Phys. J. C (2021) 81:966 Page 9 of 25 966 events in the energy range between 1016.5eV and 1017 eV are included in the likelihood as described in [37]. As S35 depends on the mass composition of the primary particles, the relation between S35 and EFD, 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 Xmaxsuggest that this composition trend follows a logarithmic evolution up to an energy of 1018.3eV, 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, ESD =ASB 35,(7) which is expected from Monte-Carlo simulations in the case of a single logarithmic dependence of Xmax with energy. The energy of an event with S35 =1VEM arriving at the reference angle, A, and the logarithmic slope, B, 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 Aand Bwhile accounting for the clustering of events in the range 1017.4to 1017.7eV observed in Fig. 7due 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 cosmic-ray 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 cosmicray intensity observed through the effective aperture of the telescopes(and thustoreduce relianceonmass assumptions), the observed distribution of events passing the cuts described above is used. The FD energy resolution, σFD(E)/EFD,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(S35)/S35,is,on the other hand, comprised of two terms, the detector sampling fluctuations, σdet(S35), and the shower-to-shower fluctuations,σsh(S35).Theformerisobtainedfrom 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 S35 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 derivation of Aand Bfollows a two-step procedure. A first iteration of the fit is carried out by using an 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 (≃1017 and ≃1019 eV, respectively) Systematic Uncertainty (%) Absolute fluorescence yield 3.6 Atmosphere and scattering 2–6 FD Calibration 10 Longitudinal profile reconstruction 7–5.5 Invisible energy 3–1.5 educated guess for σsh(S35), as expected from Monte-Carlo simulationsforamass-compositionscenariocompatiblewith data [29]. The total resolution σSD(S35)/S35 is then extracted from data as explained next in Sect. 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 A=(13.2±0.3)PeV and B=1.002 ±0.006. The goodness of the fit is supported by the χ2/NDOF =2120/1978 (p=0.013). We use these values of Aand Bto calibrate the shower sizes in terms of energiesbydefining the SD estimator of energies, ESD, according to Eq. (7). The SD energy scale is set by the calibration procedure and thus it inherits the Aand Bcalibration-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. 3.4 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, ESD and EFD, respectively. Unlike for the SD, the FD directly provides a view of the shower development so a total energy resolution, σFD(E), can be estimated for each of the golden hybrid events. Using the known σFD(E), the resolution of SD can be determined by studying the distribution of the ratio of the two energy measurements. 123 966 Page 16 of 25 Eur. Phys. J. C (2021) 81:966 6 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 1017 eV. The measurement includes a total exposure of 105km2sryr 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 1016.9eV 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 1017 eV and thus that it is their decline that largely defines the shape of the all-particle spectrum. The hypothesis is also supported by a preliminary study of the distributions of the depths of the shower maximum, Xmax, 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 protonsto Fe-nuclei,hasbeenderived. Fromallthree models, a fall-off of the Fe component above 1017 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 Z, namely Emax(Z)≃ZEproton max , to explain the knee structures. The measurements of the all-particle flux from various experiments [9–11,44–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 respective 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×1019 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, Xmax values will be determined down to 1017 eV using the three HEAT telescopes. On the other hand, the determination of the muon component of EAS above 1017 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 Z, namely Emax(Z)≃ZEproton max . In addition, we will be able to extend the measurement of the energy spectrum below 1017 eV with a denser array of 433m-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. Acknowledgements The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort from the technical and administrative staff in Malargüe. We are very grateful to the following agencies and organizations for financial support: Argentina – Comisión Nacional de Energía Atómica; Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT); Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET); Gobierno de la Provincia de Mendoza; Municipalidad de Malargüe; NDM Holdings and Valle Las Leñas; in gratitude for their continuing cooperation over land access; Australia – the Australian Research Council; Belgium – Fonds de la Recherche Scientifique (FNRS); Research Foundation Flanders (FWO); Brazil – Conselho Nacional de Desenvolvimento Científico e Tecnológico(CNPq);FinanciadoradeEstudoseProjetos(FINEP);Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ); São Paulo Research Foundation (FAPESP) Grants No. 2019/101512, No. 2010/07359-6 and No. 1999/05404-3; Ministério da Ciência, Tecnologia, Inovações e Comunicações (MCTIC); Czech Republic – Grant No. MSMT CR LTT18004, LM2015038, LM2018102, CZ.02.1.01/0.0/0.0/16_013/0001402, CZ.02.1.01/0.0/0.0/18_046/ 0016010 and CZ.02.1.01/0.0/0.0/17_049/0008422; France – Centre de Calcul IN2P3/CNRS; Centre National de la Recherche Scientifique(CNRS);ConseilRégionalIle-de-France;DépartementPhysique Nucléaire et Corpusculaire (PNC-IN2P3/CNRS); Département Sciences de l’Univers (SDU-INSU/CNRS); Institut Lagrange de Paris 123 Eur. Phys. J. C (2021) 81:966 Page 17 of 25 966 (ILP) Grant No. LABEX ANR-10-LABX-63 within the Investissements d’Avenir Programme Grant No. ANR-11-IDEX-0004-02; Germany – Bundesministerium für Bildung und Forschung (BMBF); Deutsche Forschungsgemeinschaft (DFG); Finanzministerium BadenWürttemberg; Helmholtz Alliance for Astroparticle Physics (HAP); Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF); Ministerium für Innovation, Wissenschaft und Forschung des Landes Nordrhein-Westfalen; Ministerium für Wissenschaft, Forschung und Kunst des Landes Baden-Württemberg; Italy – Istituto Nazionale di Fisica Nucleare (INFN); Istituto Nazionale di Astrofisica (INAF); Ministero dell’Istruzione, dell’Universitá e della Ricerca (MIUR); CETEMPS Center of Excellence; Ministero degli Affari Esteri (MAE); México – Consejo Nacional de Ciencia y Tecnología (CONACYT) No. 167733; Universidad Nacional Autónoma de México (UNAM); PAPIIT DGAPA-UNAM; The Netherlands – Ministry of Education, Culture and Science; Netherlands Organisation for Scientific Research (NWO); Dutch national e-infrastructure with the support of SURF Cooperative; Poland -Ministry of Science and Higher Education, grant No. DIR/WK/2018/11; National Science Centre, Grants No. 2013/08/M/ST9/00322, No. 2016/23/B/ST9/01635 and No. HARMONIA 5–2013/10/M/ST9/00062, UMO-2016/22/M/ST9/00198; Portugal – Portuguese national funds and FEDER funds within Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia (COMPETE); Romania – Romanian Ministry of Education and Research, the Program Nucleu within MCI (PN19150201/16N/2019 and PN19060102) and project PN-III-P11.2-PCCDI-2017-0839/19PCCDI/2018 within PNCDI III; Slovenia – Slovenian Research Agency, grants P1-0031, P1-0385, I0-0033, N10111; Spain – Ministerio de Economía, Industria y Competitividad (FPA2017-85114-P and PID2019-104676GB-C32), Xunta de Galicia (ED431C 2017/07), Junta de Andalucía (SOMM17/6104/UGR, P18-FR-4314) Feder Funds, RENATA Red Nacional Temática de Astropartículas (FPA2015-68783-REDT) and María de Maeztu Unit of Excellence (MDM-2016-0692); USA – Department of Energy, Contracts No. DE-AC02-07CH11359, No. DE-FR02-04ER41300, No. DEFG02-99ER41107 and No. DE-SC0011689; National Science Foundation, Grant No. 0450696; The Grainger Foundation; Marie CurieIRSES/EPLANET; European Particle Physics Latin American Network; University of Delaware Research Foundation (UDRF) – 2019; and UNESCO. Data Availability Statement This manuscript has data included as electronic supplementary material. The online version of this article contains supplementary material, which is available to authorized users. Open Access This article is licensed under a Creative Commons Attribution4.0 InternationalLicense,which permitsuse,sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP3. Appendix A: The electromagnetic trigger algorithms 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 ≃150ns (≃6 ADC bins) pulse with an amplitude of ≃1IVEM, where IVEM 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. TheToTd trigger uses the typical decay time of Cherenkov light inside the water volume, τ=67ns, 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 ≈60ns decay constant, is compressed into one or two time bins. The exponential tail of the signal is deconvolved using Di=Si−Si−1e−t/τ 1−e−t/τ ,(A.1) where Siis the signal in the i-th time-bin and t=25ns is the ADC bin-width. For an exponential decay with the mean decaytime, the deconvolvedvalues, Di,wouldbe zero. Howeverforanexponentialdecaywithstatisticalnoise that is proportional to √Si,theset{Di}would exponentially decrease withanincreaseddecaylength τ=2τ. After performing the deconvolutioninEq.(A.1),thetriggerissatisfiedif≥13ADC bins (≥325ns) are above 0.2IVEM, in coincidence between two of the three PMTs, within a sliding 3µs (120bin) 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.Only11binsareabove0.2IVEM 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 apositive step.5For each positive step, the total vertical increase, j, must be above that of typical noise, and below the characteristic amplitude of a vertical muon, namely 3<j≤31. If more than four of the positive-step instances 5For example, four bins with Si≤Si+1≤Si+2≤Si+3is considered one positive step, not three positive steps. 123 966 Page 18 of 25 Eur. Phys. J. C (2021) 81:966 Fig. 19 Top: Example waveform which passes the ToTd algorithm. Middle: The deconvolution (Eq. (A.1)) of the first waveform along with the threshold to pass the algorithm (dashed red line). Bottom: Example waveform which passes the MoPS algorithm 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 bin migration is corrected to produce the unfolded spectrum.Themagnitudeofthecorrectionfactor,asdescribed by Eq. (12), is shown in Fig. 20 along with the statistical uncertainty band. The energy spectrum of the SD-750 array isreportedinTable8andthecorrelationmatrixofthespectral parameters at the nominal energy scale in Table 9(statistical uncertainties). Finally, the combined energy spectrum is Fig. 20 The scaling factor that has been applied to the raw spectrum to produce the unfolded spectrum (see Eq. (12)) and the statistical uncertainty Table 8 SD-750 spectrum data. The correlations between systematic uncertainties are provided in the Supplementary material lg(E/eV)NJ±σstat±σsyst km2yr sr eV 17.05 217094 6.568 +0.015 +2.0 −0.015 −1.8×10−14 17.15 132828 3.302 +0.010 +1.0 −0.010 −0.9×10−14 17.25 79931 1.625 +0.006 +0.5 −0.006 −0.5×10−14 17.35 47509 7.860 +0.038 +2.5 −0.038 −2.3×10−15 17.45 27889 3.738 +0.023 +1.2 −0.023 −1.1×10−15 17.55 16407 1.775 +0.014 +0.6 −0.014 −0.5×10−15 17.65 9695 8.44 +0.09 +2.9 −0.09 −2.5×10−16 17.75 5653 3.95 +0.05 +1.5 −0.05 −1.1×10−16 17.85 3317 1.86 +0.03 +0.7 −0.03 −0.5×10−16 17.95 1990 8.91 +0.20 +3.6 −0.20 −2.6×10−17 18.05 1158 4.14 +0.12 +1.8 −0.12 −1.2×10−17 18.15 651 1.85 +0.07 +0.9 −0.07 −0.5×10−17 18.25 367 8.35 +0.43 +4.1 −0.45 −2.4×10−18 18.35 235 4.26 +0.27 +2.2 −0.28 −1.2×10−18 18.45 139 2.01 +0.17 +1.1 −0.17 −0.6×10−18 18.55 79 9.0+1.0+5.2 −1.0−2.5×10−19 18.65 45 4.1+0.7+2.4 −0.6−1.2×10−19 18.75 31 2.3+0.4+1.3 −0.4−0.6×10−19 18.85 29 1.7+0.3+0.9 −0.3−0.5×10−19 19.10 36 2.8+0.5+1.6 −0.5−0.8×10−20 19.40 7 5.7+2.4+3.2 −2.4−1.6×10−21 reported in Table 10 and the correlation matrix of the spectral parameters at the nominal energy scale in Table 11 (statistical uncertainties). 123 Eur. Phys. J. C (2021) 81:966 Page 19 of 25 966 Table 9 Elements of the correlation matrix (statistical uncertainties) of the spectral parameters describing the SD-750 energy spectrum at the nominal energy scale J0γ1E12 γ2ω01 J010.978 −0.067 −0.120 0.998 γ11−0.094 −0.109 0.967 E12 1−0.814 −0.059 γ21−0.123 ω01 1 Table 10 Combined SD spectrum data. The correlations between systematic uncertainties are provided in the Supplementary material lg(E/eV)J±σstat±σsyst km2yr sr eV 17.05 6.341 +0.015 +2.1 −0.015 −1.9×10−14 17.15 3.191 +0.010 +1.1 −0.010 −0.9×10−14 17.25 1.577 +0.006 +0.5 −0.006 −0.5×10−14 17.35 7.643 +0.039 +2.6 −0.039 −2.3×10−15 17.45 3.650 +0.024 +1.3 −0.024 −1.1×10−15 17.55 1.739 +0.015 0.6 −0.015 −0.5×10−15 17.65 8.32 +0.09 +3.0 −0.09 −2.4×10−16 17.75 3.90 +0.05 +1.4 −0.05 −1.1×10−16 17.85 1.85 +0.03 +0.7 −0.03 −0.5×10−16 17.95 8.87 +0.20 +3.3 −0.20 −2.5×10−17 18.05 4.14 +0.12 +1.6 −0.12 −1.2×10−17 18.15 1.90 +0.07 +0.7 −0.07 −0.5×10−17 Table 10 continued lg(E/eV)J±σstat±σsyst km2yr sr eV 18.25 8.47 +0.43 +3.3 −0.44 −2.4×10−18 18.35 4.17 +0.28 +1.7 −0.27 −1.2×10−18 18.45 1.929 +0.007 +0.7 −0.007 −0.5×10−18 18.55 9.041 +0.042 +2.9 −0.042 −2.0×10−19 18.65 4.294 +0.026 +1.3 −0.026 −0.9×10−19 18.75 2.167 +0.016 +0.6 −0.016 −0.4×10−19 18.85 1.226 +0.011 +0.3 −0.011 −0.2×10−19 18.95 6.82 +0.08 +1.6 −0.08 −1.3×10−20 19.05 3.79 +0.05 +0.9 −0.05 −0.7×10−20 19.15 2.07 +0.03 +0.5 −0.03 −0.4×10−20 19.25 1.04 +0.02 +0.2 −0.02 −0.2×10−20 19.35 0.53 +0.01 +1.6 −0.01 −1.3×10−20 19.45 2.49 +0.08 +0.9 −0.08 −0.7×10−21 19.55 1.25 +0.05 +0.5 −0.05 −0.3×10−21 19.65 5.99 +0.32 +2.4 −0.32 −1.8×10−22 19.75 1.95 +0.17 +0.9 −0.17 −0.7×10−22 19.85 8.1+1.0+4.0 −0.9−2.8×10−23 19.95 1.8+0.5+1.0 −0.4−0.7×10−23 20.05 5.5+2.5+3.3 −1.8−2.2×10−24 20.15 2.9+1.7+1.9 −1.2−1.2×10−24 123 966 Page 20 of 25 Eur. Phys. J. C (2021) 81:966 Table 11 Elements of the correlation matrix (statistical uncertainties) of the spectral parameters describing the combined SD energy spectrum J0γ1E12 γ2E23 γ3E34 γ4ω01 J01−0.470 0.552 −0.357 −0.383 −0.095 −0.033 0.035 −0.258 γ11−0.585 0.524 0.877 0.358 0.075 −0.085 0.966 E12 1−0.896 −0.425 0.192 0.119 0.110 −0.493 γ210.455 −0.385 −0.217 −0.154 0.475 E23 10.252 −0.063 −0.174 0.858 γ310.474 0.136 0.366 E34 10.805 0.075 γ41−0.078 ω01 1 References 1. G.V. Kulikov, G.B. Khristiansen, On the size spectrum of extensive air showers. J. Exp. Theor. Phys. 35, 635 (1958) 2. HEGRA Collaboration, Energy spectrum and chemical compositionofcosmicraysbetween0.3 PeVand 10eVdeterminedfromthe Cherenkov light and charged particle distributions in air showers. Astron. Astrophys. 359, 682 (2000). arXiv:astro-ph/9908202 3. J.W. Fowler, L.F. Fortson, C.C.H. Jui, D.B. Kieda, R.A. Ong, C.L. Pryke et al., A measurement of the cosmic ray spectrum and composition at the knee. Astropart. Phys. 15, 49 (2001). https://doi.org/ 10.1016/S0927-6505(00)00139-0.arXiv:astro-ph/0003190 4. EAS-TOP Collaboration, The cosmic ray primary composition in the“knee”regionthroughtheEASelectromagneticandmuonmeasurements at EAS-TOP. Astropart. Phys. 21, 583 (2004). https:// doi.org/10.1016/j.astropartphys.2004.04.005 5. MACRO, EAS-TOP Collaboration, The primary cosmic ray composition between 1015 and 1016 eV from extensive air showers electromagnetic and TeV muon data. Astropart. Phys. 20, 641 (2004). https://doi.org/10.1016/j.astropartphys.2003.10.004. arXiv:astro-ph/0305325 6. A.P.Garyaka, R.M. Martirosov,S.V. Ter-Antonyan,N.Nikolskaya, Y.A. Gallant, L. Jones et al., Rigidity-dependent cosmic ray energy spectra in the knee region obtained with the GAMMA experiment. Astropart. Phys. 28, 169 (2007). https://doi.org/10.1016/j. astropartphys.2007.04.004.arXiv:0704.3200 7. P. Blasi, The origin of galactic cosmic rays. Astron. Astrophys. Rev. 21, 70 (2013). https://doi.org/10.1007/s00159-013-0070-7. arXiv:1311.7346 8. KASCADE-Grande Collaboration, The spectrum of high-energy cosmic rays measured with KASCADE-Grande. Astropart. Phys. 36, 183 (2012). https://doi.org/10.1016/j.astropartphys.2012.05. 023.arXiv:1206.3834 9. IceCube Collaboration, Cosmic ray spectrum and composition from PeV to EeV using 3 years of data from IceTop and IceCube. Phys. Rev. D 100, 082002 (2019). https://doi.org/10.1103/ PhysRevD.100.082002.arXiv:1906.04317 10. Telescope Array Collaboration, The cosmic-ray energy spectrum between 2 PeV and 2 EeV Observed with the TALE detector in monocular mode. Astrophys. J. 865, 74 (2018). https://doi.org/10. 3847/1538-4357/aada05.arXiv:1803.01288 11. N.M. Budnev et al., The primary cosmic-ray energy spectrum measured with the Tunka-133 array. Astropart. Phys. 117, 102406 (2020). https://doi.org/10.1016/j.astropartphys.2019.102406 12. A. Albert et al., Evidence of 200 TeV photons from HAWC. Astrophys. J. Lett. 907, L30 (2021). https://doi.org/10.3847/2041-8213/ abd77b.arXiv:2012.15275 13. Tibet ASgamma Collaboration, First detection of sub-PeV diffuse gamma rays from the galactic disk: evidence for ubiquitous galactic cosmic rays beyond PeV energies. Phys. Rev. Lett. 126, 141101 (2021). https://doi.org/10.1103/PhysRevLett.126.141101. arXiv:2104.05181 14. LHAASO Collaboration, Ultrahigh-energy photons up to 1.4 petaelectronvolts from 12 gamma-ray galactic sources. Nature 594, 33–36 (2021). https://doi.org/10.1038/s41586-021-03498-z 15. LHAASO Collaboration, Discovery of the Ultra-high energy gamma-ray source LHAASO J2108+5157. arXiv:2106.09865 16. P. Cristofari, P. Blasi, E. Amato, The low rate of Galactic pevatrons. Astropart. Phys. 123, 102492 (2020). https://doi.org/10.1016/j. astropartphys.2020.102492.arXiv:2007.04294 17. A.M. Hillas, Can diffusive shock acceleration in supernova remnants account for high-energy galactic cosmic rays? J. Phys. G 31, R95 (2005). https://doi.org/10.1088/0954-3899/31/5/R02 18. KASCADE-Grande Collaboration, KASCADE-Grande measurements of energy spectra for elemental groups of cosmic rays. Astropart. Phys. 47 (2013) 54. https://doi.org/10.1016/j. astropartphys.2013.06.004.arXiv:1306.6283 19. R. Aloisio, V. Berezinsky, P. Blasi, Ultra high energy cosmic rays: implications of Auger data for source spectra and chemical composition. JCAP 10, 020 (2014). https://doi.org/10.1088/1475-7516/ 2014/10/020.arXiv:1312.7459 20. Pierre Auger Collaboration, Features of the energy spectrum of cosmic rays above 2.5×1018 eV using the Pierre Auger Observatory. Phys. Rev. Lett. 125, 121106 (2020). https://doi.org/10.1103/ PhysRevLett.125.121106.arXiv:2008.06488 21. PierreAugerCollaboration,Measurementofthecosmic-rayenergy spectrum above 2.5×1018 eV using the Pierre Auger Observatory. Phys. Rev. D 102, 062005 (2020). https://doi.org/10.1103/ PhysRevD.102.062005.arXiv:2008.06486 22. Pierre Auger Collaboration, Trigger and Aperture of the Surface Detector Array of the Pierre Auger Observatory. Nucl. Instrum. Meth. A 613, 29 (2010). https://doi.org/10.1016/j.nima.2009.11. 018.arXiv:1111.6764 23. Pierre Auger Collaboration, Calibration of the surface array of the PierreAuger Observatory. Nucl.Instrum. Meth.A 568, 839 (2006). https://doi.org/10.1016/j.nima.2006.07.066.arXiv:2102.01656 24. Pierre Auger Collaboration, The Pierre Auger Cosmic Ray Observatory. Nucl. Instrum. Meth. A 798, 172 (2015). https://doi.org/10. 1016/j.nima.2015.06.058.arXiv:1502.01323 25. Pierre Auger Collaboration, The lateral trigger probability function for the ultra-high energy cosmic ray showers detected by the Pierre Auger Observatory. Astropart. Phys. 35, 266 (2011). https://doi. org/10.1016/j.astropartphys.2011.08.001.arXiv:1111.6645 26. Pierre Auger Collaboration, Reconstruction of events recorded with the surface detector of the pierre auger observatory. JINST 15, P10021 (2020). https://doi.org/10.1088/1748-0221/15/ 10/P10021.arXiv:2007.09035 123 Eur. Phys. J. C (2021) 81:966 Page 21 of 25 966 27. A.M. Hillas, Derivation of the EAS spectrum. Acta Phys. Acad. Sci. Hung. 29, 355 (1970) 28. D.W. Newton, J. Knapp, A.A. Watson, The optimum distance at which to determine the size of a giant air shower. Astropart. Phys. 26, 414 (2007). https://doi.org/10.1016/j.astropartphys.2006.08. 003.arXiv:astro-ph/0608118 29. J. Bellido (Pierre Auger Collaboration), Depth of maximum of airshower profiles at the Pierre Auger Observatory: measurements above 1017.2eV and composition implications. PoS ICRC2017, 506 (2017). https://doi.org/10.22323/1.301.0506 30. Pierre Auger Collaboration, Impact of atmospheric effects on the energy reconstruction of air showers observed by the surface detectors of the Pierre Auger Observatory. JINST 12, P02006 (2017). https://doi.org/10.1088/1748-0221/12/02/ p02006.arXiv:1702.02835 31. J. Hersil, I. Escobar, D. Scott, G. Clark, S. Olbert, Observations of extensive air showers near the maximum of their longitudinal development. Phys. Rev. Lett. 6, 22 (1961). https://doi.org/ 10.1103/PhysRevLett.6.22 32. Pierre Auger Collaboration, Cosmic-ray anisotropies in right ascension measured by the Pierre Auger Observatory. Astrophys. J. 891, 142 (2020). https://doi.org/10.3847/1538-4357/ab7236. arXiv:2002.06172 33. T.W. Anderson, D.A. Darling, A test of goodness of fit. J. Am. Stat. Assoc. 49, 765 (1954) 34. Pierre Auger Collaboration, Data-driven estimation of the invisible energy of cosmic ray showers with the Pierre Auger Observatory. Phys. Rev. D 100, 082003 (2019). https://doi.org/10.1103/ PhysRevD.100.082003.arXiv:1901.08040 35. Pierre Auger Collaboration, Direct measurement of the muonic content of extensive air showers between 2 ×1017 and 2 ×1018 eV at the Pierre Auger Observatory. Eur. Phys. J. C 80, 751 (2020). https://doi.org/10.1140/epjc/s10052-020-8055-y 36. Pierre Auger Collaboration, Depth of maximum of air-shower profiles at the Pierre Auger Observatory. II. Composition implications. Phys. Rev. D 90, 122006 (2014). https://doi.org/10.1103/ PhysRevD.90.122006.arXiv:1409.5083 37. H.P. Dembinski, B. Kégl, I.C. Mari¸s, M. Roth, D. Veberiˇc, A likelihood method to cross-calibrate air-shower detectors. Astropart. Phys. 73, 44 (2016). https://doi.org/10.1016/j.astropartphys.2015. 08.001.arXiv:1503.09027 38. Pierre Auger Collaboration, The energy scale of the Pierre Auger Observatory. PoS ICRC2019, 231 (2020). https://doi.org/ 10.22323/1.358.0231 39. D.Hecketal.,CORSIKA:aMonteCarlocodetosimulateextensive air showers. Report fzka 6019 (1998) 40. T. Pierog, I. Karpenko, J.M. Katzy, E. Yatsenko, K. Werner, EPOS LHC: test of collective hadronization with data measured at the CERN Large Hadron Collider. Phys. Rev. C 92, 034906 (2015). https://doi.org/10.1103/PhysRevC.92.034906.arXiv:1306.0121 41. H.P. Dembinski, R. Engel, A. Fedynitch, T. Gaisser, F. Riehn, T. Stanev, Data-driven model of the cosmic-ray flux and mass composition from 10 GeV to 1011 GeV. PoS ICRC2017, 533 (2017). https://doi.org/10.22323/1.301.0533 42. M. Bonamente, Distribution of the C statistic with applications to the sample mean of Poisson data. J. Appl. Stat. 47, 2044 (2020). https://doi.org/10.1080/02664763.2019.1704703. arXiv:1912.05444 43. B. Dawson (Pierre Auger Collaboration), The Energy Scale of the Pierre Auger Observatory. PoS ICRC2019, 231 (2019). https:// doi.org/10.22323/1.358.0231 44. M. Nagano, M. Teshima, Y. Matsubara, H. Dai, T. Hara, N. Hayashida et al., Energy spectrum of primary cosmic rays above 1017 eV determined from the extensive air shower experiment at Akeno. J. Phys. G 18, 423 (1992). https://doi.org/10.1088/ 0954-3899/18/2/022 45. S. Ter-Antonyan, Sharp knee phenomenon of primary cosmic ray energy spectrum. Phys. Rev. D 89, 123003 (2014). https://doi.org/ 10.1103/PhysRevD.89.123003.arXiv:1405.5472 46. KASCADE-Grande Collaboration, KASCADE-Grande energy spectrum of cosmic rays interpreted with post-LHC hadronic interaction models. PoS ICRC2015, 359 (2016). https://doi.org/10. 22323/1.236.0359 47. E.N. Gudkova, N.M. Nesterova, Results of the further analysis of data from the Tien Shan array in the energy spectrum of primary cosmic rays in the energy Range of 2 ×1013 -3×1017 eV. Phys. Atomic Nuclei 83, 629 (2020). arXiv:2010.04236 48. TIBETIII Collaboration,TheAll-particle spectrumofprimarycosmicraysinthewideenergyrangefrom 1014 eVto1017 eVobserved with the Tibet-III air-shower array. Astrophys. J. 678, 1165 (2008). https://doi.org/10.1086/529514.arXiv:0801.1803 49. S.P. Knurenko, Z.E. Petrov, R. Sidorov, I.Y. Sleptsov, S.K. Starostin, G.G. Struchkov, Cosmic ray spectrum in the energy range 1015–1018 eV and the second knee according to the small Cherenkov setup at the Yakutsk EAS array, Proc. of 33rd ICRC (2013). arXiv:1310.1978 50. O.A. Gress, T.I. Gress, E.E. Korosteleva, L.A. Kuzmichev, B.K. Lubsandorzhiev, L.V. Pan’kov et al., The study of primary cosmic rays energy spectrum and mass composition in the energy range 0.5–50 PeV with TUNKA Eas Cherenkov array. Nucl. Phys. BProc. Suppl. 75, 299 (1999) 51. Pierre Auger Collaboration, Depth of maximum of air-shower profiles at the Pierre Auger Observatory: measurements above 1017.2eV and composition implications. PoS ICRC2017, 506 (2018). https://doi.org/10.22323/1.301.0506 52. S. Ostapchenko, QGSJET-II: physics, recent improvements, and resultsfor airshowers,inEPJ Web of Conferences,vol.52,p.02001 (EDP Sciences, 2013) 53. F. Riehn, H.P. Dembinski, R. Engel, A. Fedynitch, T. Gaisser, T. Stanev, The hadronic interaction model Sibyll 2.3c and Feynman scaling. PoS ICRC2017, 301 (2017). https://doi.org/10.22323/1. 301.0301.arXiv:1709.07227 54. O. Deligny, (Pierre Auger and Telescope Array Collaborations), The energy spectrum of ultra-high energy cosmic rays measured at the Pierre Auger Observatory and at the Telescope Array. PoS ICRC2019, 234 (2019). https://doi.org/10.22323/1.358.0234 55. V. Novotny (Pierre Auger Collaboration), Measurement of the spectrum of cosmic rays above 1016.5eV with Cherenkovdominatedeventsatthe PierreAugerObservatory.PoS ICRC2019, 374 (2019). https://doi.org/10.22323/1.358.0374 123 966 Page 22 of 25 Eur. Phys. J. C (2021) 81:966 Pierre Auger Collaboration P. Abreu71, M. Aglietta51,53,J.M.Albury 12, I. Allekotte1,A.Almela 8,11, J. Alvarez-Muñiz78,R.AlvesBatista 79, G. A. Anastasi51,62, L. Anchordoqui86, B. Andrada8, S. Andringa71,C.Aramo 49, P. R. Araújo Ferreira41, J. C. Arteaga Velázquez66,H.Asorey 8, P. Assis71, G. Avila10, A. M. Badescu74, A. Bakalova31, A. Balaceanu72, F. Barbato44,45, R. J. Barreira Luz71, K. H. Becker37, J. A. Bellido12,68, C. Berat35, M. E. Bertaina51,62,X.Bertou 1, P. L. Biermann95, P. Billoir34, V. Binet6,K.Bismark 8,38, T. Bister41, J. Biteau36, J. Blazek31,C.Bleve 35, M. Boháˇcová31, D. Boncioli45,56, C. Bonifazi25, L. Bonneau Arbeletche20, N. Borodai69, A.M.Botti 8, J. Brack97, T. Bretz41, P. G. Brichetto Orchera8, F. L. Briechle41, P. Buchholz43, A. Bueno77, S. Buitink14, M. Buscemi46, M. Büsken8,38, K. S. Caballero-Mora65, L. Caccianiga48,58, F. Canfora79,80, I. Caracas37, J.M.Carceller 77, R. Caruso46,57, A. Castellina51,53, F. Catalani18, G. Cataldi47, L. Cazon71, M. Cerda9, J. A. Chinellato21, J. Chudoba31, L. Chytka32,R.W.Clay 12, A. C. Cobos Cerutti7, R. Colalillo49,59, A. Coleman92 , M. R. Coluccia47, R. Conceição71, A. Condorelli44,45, G. Consolati48,54, F. Contreras10, F. Convenga47,55, D. Correia dos Santos27, C. E. Covault84, S. Dasso3,5, K. Daumiller40,B.R.Dawson 12, J.A.Day 12, R.M.deAlmeida 27, J. de Jesús8,40, S.J.deJong 79,80, G. De Mauro79,80, J. R. T. de Mello Neto25,26, I. De Mitri44,45, J. de Oliveira17, D. de Oliveira Franco21,F.dePalma 47,55, V. de Souza19, E. De Vito47,55, M. del Río10, O. Deligny33, A. Di Matteo51, C. Dobrigkeit21, J.C.D’Olivo 67, L. M. Domingues Mendes71, R. C. dos Anjos24, D. dos Santos27,M.T.Dova 4,J.Ebr 31, R. Engel38,40, I. Epicoco47,55, M. Erdmann41, C. O. Escobar94, A. Etchegoyen8,11, H. Falcke79,80,81,J.Farmer 91,G.Farrar 89, A. C. Fauth21, N. Fazzini94, F. Feldbusch39, F. Fenu51,53,B.Fick 88, J. M. Figueira8, A. Filipˇciˇc75,76, T. Fitoussi40, T. Fodran79, M. M. Freire6, T. Fujii91,98, A. Fuster8,11, C. Galea79, C. Galelli48,58, B. García7, A. L. Garcia Vegas41, H. Gemmeke39, F. Gesualdi8,40, A. Gherghel-Lascu72,P.L.Ghia 33, U. Giaccari79, M. Giammarchi48, J. Glombitza41, F. Gobbi9, F. Gollan8, G. Golup1, M. Gómez Berisso1, P. F. Gómez Vitale10, J. P. Gongora10, J. M. González1, N. González13, I. Goos1,40,D.Góra 69, A. Gorgi51,53, M. Gottowik37, T. D. Grubb12, F. Guarino49,59, G. P. Guedes22, E. Guido51,62, S. Hahn8,40, P. Hamal31, M. R. Hampel8, P. Hansen4, D. Harari1,V.M.Harvey 12, A. Haungs40, T. Hebbeker41, D. Heck40, G. C. Hill12,C.Hojvat 94, J. R. Hörandel79,80,P.Horvath 32, M. Hrabovský32, T. Huege14,40, A. Insolia46,57,P.G.Isar 73, P. Janecek31, J. A. Johnsen85, J. Jurysek31, A. Kääpä37, K. H. Kampert37, N. Karastathis40, B. Keilhauer40,J.Kemp 41, A. Khakurdikar79, V. V. Kizakke Covilakam8,40, H. O. Klages40, M. Kleifges39, J. Kleinfeller9, M. Köpke38, N. Kunka39, B. L. Lago16,R. G. Lang19, N. Langner41,M. A. Leigui de Oliveira23,V. Lenok40,A. Letessier-Selvon34,I. Lhenry-Yvon33, D. Lo Presti46,57, L. Lopes71, R. López63,L.Lu 93, Q. Luce38, J. P. Lundquist75, A. Machado Payeras21, G. Mancarella47,55, D. Mandat31, B. C. Manning12, J. Manshanden42, P. Mantsch94, S. Marafico33, A. G. Mariazzi4,I.C.Mari¸s13, G. Marsella46,60, D. Martello47,55, S. Martinelli8,40, H. Martinez19, O. Martínez Bravo63, M. Mastrodicasa45,56, H. J. Mathes40, J. Matthews87, G. Matthiae50,61, E. Mayotte37, P. O. Mazur94, G. Medina-Tanco67,D.Melo 8, A. Menshikov39, K.-D. Merenda85, S. Michal32, M. I. Micheletti6, L. Miramonti48,58, D. Mockler13,38, S. Mollerach1, F. Montanet35, C. Morello51,53,M.Mostafá 90, A.L.Müller 8, M.A.Muller 21,K.Mulrey 14, R. Mussa51, M. Muzio89, W. M. Namasaka37, A. Nasr-Esfahani37, L. Nellen67, M. Niculescu-Oglinzanu72, M. Niechciol43, D. Nitz88, D. Nosek30, V. Novotny30, L. Nožka32, A. Nucita47,55, L. A. Núñez29, M. Palatka31, J. Pallotta2, P. Papenbreer37, G. Parente78, A. Parra63, J. Pawlowsky37, M. Pech31, F. Pedreira78,J.P¸ekala69, R. Pelayo64, J. Peña-Rodriguez29, E. E. Pereira Martins8,38,J. Perez Armand20,C. Pérez Bertolli8,40,M. Perlin8,40,L. Perrone47,55,S. Petrera44,45,T.Pierog40, M. Pimenta71, V. Pirronello46,57, M. Platino8, B. Pont79, M. Pothast79,80, P. Privitera91, M. Prouza31, A. Puyleart88, S. Querchfeld37, J. Rautenberg37, D. Ravignani8, M. Reininghaus8,40,J.Ridky 31, F. Riehn71, M. Risse43,V.Rizi 45,56, W. Rodrigues de Carvalho20, J. Rodriguez Rojo10, M. J. Roncoroni8,M.Roth 40, E. Roulet1, A. C. Rovero5, P. Ruehl43, S. J. Saffi12, A. Saftoiu72, F. Salamida45,56, H. Salazar63, G. Salina50, J. D. Sanabria Gomez29, F. Sánchez8, E. M. Santos20, E. Santos31, F. Sarazin85, R. Sarmento71, C. Sarmiento-Cano8, R. Sato10,P.Savina 33,47,55, C. M. Schäfer40, V. Scherini47, H. Schieler40, M. Schimassek8,38, M. Schimp37, F. Schlüter8,40, D. Schmidt38, O. Scholten14,83, P. Schovánek31, F. G. Schröder40,92, S. Schröder37, J. Schulte41, A. Schulz38, S. J. Sciutto4, M. Scornavacche8,40,A.Segreto 46,52, S. Sehgal37, R.C.Shellard 15,G.Sigl 42, G. Silli8,40,O.Sima 72,99,R.Šmída 91, P. Sommers90, J. F. Soriano86, J. Souchard35, R. Squartini9, M. Stadelmaier8,40, D. Stanca72, S. Staniˇc75, J. Stasielak69, P. Stassi35, A. Streich8,38, M. Suárez-Durán13, T. Sudholz12, T. Suomijärvi36, A. D. Supanitsky8, Z. Szadkowski70, A. Tapia28, C. Taricco51,62, C. Timmermans79,80, O. Tkachenko40, P. Tobiska31, C. J. Todero Peixoto18,B.Tomé 71, Z. Torrès35,A.Travaini 9, P. Travnicek31, C. Trimarelli45,56, M. Tueros4, R. Ulrich40, M. Unger40, L. Vaclavek32, M. Vacula32, J. F. Valdés Galicia67, L. Valore49,59, E. Varela63, A. Vásquez-Ramírez29, D. Veberiˇc40, C. Ventura26,I.D.VergaraQuispe 4, V. Verzi50, J. Vicha31, J. Vink82, S. Vorobiov75, H. Wahlberg4, C. Watanabe25, A.A.Watson 96, M. Weber39, A. Weindl40, L. Wiencke85, 123 Eur. Phys. J. C (2021) 81:966 Page 23 of 25 966 H. Wilczy´nski69, M. Wirtz41, D. Wittkowski37, B. Wundheiler8, A. Yushkov31, O. Zapparrata13,E.Zas 78, D. Zavrtanik75,76, M. Zavrtanik75,76, L. Zehrer75 1Centro Atómico Bariloche and Instituto Balseiro (CNEA-UNCuyo-CONICET), San Carlos de Bariloche, Argentina 2Centro de Investigaciones en Láseres y Aplicaciones, CITEDEF and CONICET, Villa Martelli, Argentina 3Departamento de Física and Departamento de Ciencias de la Atmósfera y los Océanos, FCEyN, Universidad de Buenos Aires and CONICET, Buenos Aires, Argentina 4IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina 5Instituto de Astronomía y Física del Espacio (IAFE, CONICET-UBA), Buenos Aires, Argentina 6Instituto de Física de Rosario (IFIR)-CONICET/U.N.R. and Facultad de Ciencias Bioquímicas y Farmacéuticas U.N.R., Rosario, Argentina 7Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), Universidad Tecnológica Nacional-Facultad Regional Mendoza (CONICET/CNEA), Mendoza, Argentina 8Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), Buenos Aires, Argentina 9Observatorio Pierre Auger, Malargüe, Argentina 10 Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe, Argentina 11 Universidad Tecnológica Nacional-Facultad Regional Buenos Aires, Buenos Aires, Argentina 12 University of Adelaide, Adelaide, SA, Australia 13 Université Libre de Bruxelles (ULB), Brussels, Belgium 14 Vrije Universiteit Brussels, Brussels, Belgium 15 Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil 16 Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Nova Friburgo, Brazil 17 Instituto Federal de Educação, Ciência e Tecnologia do Rio de Janeiro (IFRJ), Rio de Janeiro, Brazil 18 Escola de Engenharia de Lorena, Universidade de São Paulo, Lorena, SP, Brazil 19 Instituto de Física de São Carlos, Universidade de São Paulo, São Carlos, SP, Brazil 20 Instituto de Física, Universidade de São Paulo, São Paulo, SP, Brazil 21 Universidade Estadual de Campinas, IFGW, Campinas, SP, Brazil 22 Universidade Estadual de Feira de Santana, Feira de Santana, Brazil 23 Universidade Federal do ABC, Santo André, SP, Brazil 24 Universidade Federal do Paraná, Setor Palotina, Palotina, Brazil 25 Instituto de Física, Universidade Federal do Rio de Janeiro, Rio de Janeiro, RJ, Brazil 26 Observatório do Valongo, Universidade Federal do Rio de Janeiro (UFRJ), Rio de Janeiro, RJ, Brazil 27 Universidade Federal Fluminense, EEIMVR, Volta Redonda, RJ, Brazil 28 Universidad de Medellín, Medellín, Colombia 29 Universidad Industrial de Santander, Bucaramanga, Colombia 30 Institute of Particle and Nuclear Physics, Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic 31 Institute of Physics of the Czech Academy of Sciences, Prague, Czech Republic 32 Palacky University, RCPTM, Olomouc, Czech Republic 33 CNRS/IN2P3, IJCLab, Université Paris-Saclay, Orsay, France 34 Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Sorbonne Université, Université de Paris, CNRS-IN2P3, Paris, France 35 Univ. Grenoble Alpes, CNRS, Grenoble Institute of Engineering Univ. Grenoble Alpes, LPSC-IN2P3, 38000 Grenoble, France 36 Université Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France 37 Department of Physics, Bergische Universität Wuppertal, Wuppertal, Germany 38 Institute for Experimental Particle Physics, Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany 39 Institut für Prozessdatenverarbeitung und Elektronik, Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany 40 Institute for Astroparticle Physics, Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany 41 III. Physikalisches Institut A, RWTH Aachen University, Aachen, Germany 42 II. Institut für Theoretische Physik, Universität Hamburg, Hamburg, Germany 43 Department Physik-Experimentelle Teilchenphysik, Universität Siegen, Siegen, Germany 44 Gran Sasso Science Institute, L’Aquila, Italy 123 966 Page 24 of 25 Eur. Phys. J. C (2021) 81:966 45 INFN Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy 46 INFN, Sezione di Catania, Catania, Italy 47 INFN, Sezione di Lecce, Lecce, Italy 48 INFN, Sezione di Milano, Milan, Italy 49 INFN, Sezione di Napoli, Naples, Italy 50 INFN, Sezione di Roma “Tor Vergata”, Rome, Italy 51 INFN, Sezione di Torino, Turin, Italy 52 Istituto di Astrofisica Spaziale e Fisica Cosmica di Palermo (INAF), Palermo, Italy 53 Osservatorio Astrofisico di Torino (INAF), Turin, Italy 54 Dipartimento di Scienze e Tecnologie Aerospaziali, Politecnico di Milano, Milan, Italy 55 Dipartimento di Matematica e Fisica “E. De Giorgi”, Università del Salento, Lecce, Italy 56 Dipartimento di Scienze Fisiche e Chimiche, Università dell’Aquila, L’Aquila, Italy 57 Dipartimento di Fisica e Astronomia, Università di Catania, Catania, Italy 58 Dipartimento di Fisica, Università di Milano, Milan, Italy 59 Dipartimento di Fisica “Ettore Pancini”, Università di Napoli “Federico II”, Naples, Italy 60 Dipartimento di Fisica e Chimica ”E. Segrè”, Università di Palermo, Palermo, Italy 61 Dipartimento di Fisica, Università di Roma “Tor Vergata”, Rome, Italy 62 Dipartimento di Fisica, Università Torino, Turin, Italy 63 Benemérita Universidad Autónoma de Puebla, Puebla, Mexico 64 Unidad Profesional Interdisciplinaria en Ingeniería y Tecnologías Avanzadas del Instituto Politécnico Nacional (UPIITA-IPN), Mexico, D.F., Mexico 65 Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas, Mexico 66 Universidad Michoacana de San Nicolás de Hidalgo, Morelia, Michoacán, México 67 Universidad Nacional Autónoma de México, Mexico, D.F., Mexico 68 Facultad de Ciencias Naturales y Formales, Universidad Nacional de San Agustin de Arequipa, Arequipa, Peru 69 Institute of Nuclear Physics PAN, Krakow, Poland 70 Faculty of High-Energy Astrophysics, University of Łód´z, Łód´z, Poland 71 Laboratório de Instrumentação e Física Experimental de Partículas – LIP and Instituto Superior Técnico-IST, Universidade de Lisboa-UL, Lisbon, Portugal 72 “Horia Hulubei” National Institute for Physics and Nuclear Engineering, Bucharest-Magurele, Romania 73 Institute of Space Science, Bucharest-Magurele, Romania 74 University Politehnica of Bucharest, Bucharest, Romania 75 Center for Astrophysics and Cosmology (CAC), University of Nova Gorica, Nova Gorica, Slovenia 76 Experimental Particle Physics Department, J. Stefan Institute, Ljubljana, Slovenia 77 Universidad de Granada and C.A.F.P.E., Granada, Spain 78 Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de Santiago de Compostela, Santiago de Compostela, Spain 79 IMAPP, Radboud University Nijmegen, Nijmegen, The Netherlands 80 Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science Park, Amsterdam, The Netherlands 81 Stichting Astronomisch Onderzoek in Nederland (ASTRON), Dwingeloo, The Netherlands 82 Faculty of Science, Universiteit van Amsterdam, Amsterdam, The Netherlands 83 Kapteyn Astronomical Institute, University of Groningen, Groningen, The Netherlands 84 Case Western Reserve University, Cleveland, OH, USA 85 Colorado School of Mines, Golden, CO, USA 86 Department of Physics and Astronomy, Lehman College, City University of New York, Bronx, NY, USA 87 Louisiana State University, Baton Rouge, LA, USA 88 Michigan Technological University, Houghton, MI, USA 89 New York University, New York, NY, USA 90 Pennsylvania State University, University Park, PA, USA 91 University of Chicago, Enrico Fermi Institute, Chicago, IL, USA 92 Department of Physics and Astronomy, Bartol Research Institute, University of Delaware, Newark, DE, USA 93 Department of Physics and WIPAC, University of Wisconsin-Madison, Madison, WI, USA 94 Now at Fermi National Accelerator Laboratory, Fermilab, Batavia, IL, USA 123 Eur. Phys. J. C (2021) 81:966 Page 25 of 25 966 95 Now at Max-Planck-Institut für Radioastronomie, Bonn, Germany 96 Now at School of Physics and Astronomy, University of Leeds, Leeds, UK 97 Now at Colorado State University, Fort Collins, CO, USA 98 Now at Hakubi Center for Advanced Research and Graduate School of Science, Kyoto University, Kyoto, Japan 99 Now at University of Bucharest, Physics Department, Bucharest, Romania 123