Determination of the invisible energy of extensive air showers from the data collected at Pierre Auger Observatory

In order to get the primary energy of cosmic rays from their extensive air showers using the fluorescence detection technique, the invisible energy should be added to the measured calorimetric energy. The invisible energy is the energy carried away by particles that do not deposit all their energy in the atmosphere. It has traditionally been calculated using Monte Carlo simulations that are dependent on the assumed primary particle mass and on model predictions for neutrino and muon production. In this work the invisible energy is obtained directly from events detected by the Pierre Auger Observatory. The method applied is based on the correlation of the measurements of the muon number at the ground with the invisible energy of the showers. By using it, the systematic uncertainties related to the unknown mass composition and to the high energy hadronic interaction models are significantly reduced, improving in this way the estimation of the energy scale of the Observatory.


Introduction
Above 10 15 eV cosmic rays are detected indirectly through the extensive air showers (EAS) they produce in the atmosphere. Most of the cosmic ray energy is carried by electromagnetic particles of the EAS, which can be detected by their secondary electromagnetic signatures, e.g. radio, Cherenkov or fluorescence light. In the case of the fluorescence detection, the fluorescence radiation emitted by the nitrogen molecules of air excited by the charged particles of the EAS is produced in proportion to the energy dissipation, allowing a reconstruction of the longitudinal profile of the energy deposit (dE/dX) of the shower as a function of the atmospheric depth X. The atmosphere is used as a calorimeter and the integral (dE/dX) dX, called the calorimetric energy of the shower, E cal , is measured.
E cal underestimates the total shower energy (E 0 ) because neutrinos do not suffer electromagnetic interactions and high energy muons reach ground level after releasing only a portion of their energy into the atmosphere. Thus, an estimation of the primary energy E 0 with the fluorescence detection technique is obtained by adding to E cal a correction to account for the invisible energy (E inv ) carried by the particles that do not dissipate all their energy in the atmosphere. E inv amounts to about 10% -20% of E 0 .
E inv can be calculated directly from the energy deposited in the atmosphere by the different components of simulated air showers [1]. There are large differences in the values of the ratio of E inv to E 0 as a function of E cal for different hadronic models and primary masses, as could be * e-mail: mariazzi@fisica.unlp.edu.ar * * e-mail: auger\protect\_spokespersons@fnal.gov seen in Fig. 1. The estimation of E inv is affected by the irreducible uncertainties associated with the models describing the hadronic interactions and also by the mass composition of cosmic rays.
The models to get E inv can be improved further using the primary mass composition estimated with the fluorescence detectors [2] so that the spread between the predictions is significantly reduced for a given mass. However, the uncertainties associated with the hadronic interaction models are difficult to estimate and are ultimately unknown [3]. Even after the updates with LHC data, the models still fail to describe several properties of the shower development related to muons [15], and this can introduce unpredictable biases in the E inv estimation.
Thus the strategy followed in this work is to estimate E inv using the correlations that exist between E inv and shower observables measured at the Pierre Auger Observatory [4], correlations that to a large extent are not sensitive to the hadronic interaction models and primary mass composition.
The Pierre Auger Observatory [4] is a hybrid observatory, because the measurements are done combining the data of a Surface Detector (SD) that is sensitive to the muon content of EAS and a Fluorescence Detector (FD). The SD consists of 1660 water-Cherenkov detectors (WCDs) arranged on a hexagonal grid of 1.5 km spacing extending over a total area of ∼ 3000 km 2 . The FD consists of 24 telescopes placed in four sites located along the perimeter of the Observatory that overlook the atmosphere above the surface array. The FD operates during clear and moonless nights with a duty cycle of about 14% [5]. . Average invisible energy fraction as a function of E cal calculated with Monte Carlo simulations using the hadronic interaction models tuned with LHC data (right) and the models developed before the LHC data were available (left). The predictions for proton and iron primaries are shown with solid and dashed lines, respectively. The simulations were performed using the CORSIKA code for the models EPOS

Phenomenology of the invisible energy
The Heitler model and its extension to hadronic cascades [9,10] provide a qualitative description of EAS which is suitable enough to serve as a guiding thread in the next sections, where the starting points of the data-driven approaches to estimate E inv will be inspired by some of the expressions outlined below.
In the model, only pions are produced in the hadronic interactions, all with the same energy and the same particle multiplicity (N). The neutral pions decay almost immediately into two photons, generating an electromagnetic cascade. Charged pions interact hadronically until the average energy of the charged pions is decreased to such a level that their time-dilated decay length becomes smaller than their hadronic interaction length . This energy is referred to as the pion critical energy: where E 0 is the primary particle energy and n is the number of interactions suffered by the charged pions. One important feature of the model is that E inv is proportional to the number of muons (N µ ) reaching ground level.
This expression will be the guiding thread to estimate E inv with a measurement of N µ in inclined showers. Another important feature of the model is the power law of E inv : where β = ln( 2 3 N)/ ln N. Air shower simulations predict values of β in the range from 0.88 to 0.92 [11]. β also fixes the E inv dependence on the mass number A of the primary. Neglecting collective effects in the first interactions so that the cascade is the superposition of A cascades initiated by primary protons of energy E 0 /A one has: This relationship will be the guiding thread to estimate E inv from vertical showers measurements. Monte Carlo simulations take into account all the complex phenomena occurring throughout the EAS development giving more quantitative predictions of E inv , which is calculated following the method described in [1].
The correlation between E inv and N µ has been studied simulating showers with different primary masses and hadronic interactions models using the CORSIKA code [12, & references therein], as it is shown in Fig. 2. In spite of the very large spread in the predictions of N µ and E inv , the correlation is good and is similar for all models and primaries, suggesting that it is possible to obtain a robust estimation of E inv from the measurements of N µ .

Estimation of E inv using Auger data
Two different reconstruction techniques are used for the SD events: one for the so-called vertical showers with zenith angles θ < 60 • [18], and one for the inclined showers with θ > 60 • [6]. WCDs are sensitive to the electromagnetic and hadronic components of a shower.
The most straightforward way to estimate E inv is to use inclined showers, in which the electromagnetic component of the shower is largely absorbed and it is possible to measure the total number of muons arriving at ground level which is an observable expected to be proportional to E inv , as seen in Sec. 2 (Eq. (2)). The muon number cannot be directly measured for vertical events. However, E inv can be obtained from the energy estimator using the power law relationship between E inv and E 0 (see Eq. (4)).

E inv rom inclined showers
The reconstruction of inclined events [6] is based on the fact that the muon number distribution at ground level can be described by a density scaling factor that depends on E 0 and primary mass, and by a lateral shape that, for a given arrival direction (θ, φ) of the shower, is consistently reproduced by different hadronic interaction models and depends only weakly on E 0 and primary mass.
The muon number density as a function of the position at ground r is then parameterised with where ρ µ,19 ( r; θ, φ) is a reference distribution conventionally calculated for primary protons at 10 19 eV using the hadronic interaction model QGSJetII-03, and the scale factor N 19 represents the shower size relative to the normalization of the reference distribution. The performance of the reconstruction is validated on simulated events. The reconstructed value of N 19 is compared with its true value R MC µ for each simulated event. R MC µ is defined as the ratio of the total number of muons at ground level to the total number of muons in the reference model. The relative deviation of N 19 from R MC µ is within 5% for several hadronic interaction models and primaries [16]. A bias correction is then applied to N 19 in order to reduce the residuals to within 3% of the most recent models tuned with LHC data. In this way, the corrected value of N 19 , which in the following is called R µ , represents an unbiased estimator of the total number of muons at ground level.
The correlation between E inv and the total number of muons at ground level is studied with two data sets: one simulated with CORSIKA using the hadronic interaction models EPOS LHC and QGSJET II-04 [ For each simulated event, we calculate the values of E inv and of the muon number at ground level R MC µ . For all the samples of simulated events, the correlation between E inv and R MC µ is well described by a power-law where the values of the parameters C and δ are obtained from a fit to the events. Examples of the correlation between E inv and R MC µ are shown for in Fig. 3 (left), where the lines show the fitted power law relationships.
The relationship of Eq. (6) is used to estimate E inv in the data from the measurement of R µ that, as seen before, is the unbiased estimator of R MC µ . Since the mass composition of the data is not precisely known, the estimation of E inv is obtained using the parameterisation of E inv as a function of R µ for a mixture of 50% protons and 50% iron. This is done taking the average of the two E inv estimations in Fig. 3 (left) obtained for proton and iron primaries using the EPOS LHC hadronic interaction model.
The performance of the analysis is studied on fully simulated events 1 . For each event, we compute E inv from R µ using the estimation for the mixed proton and iron composition, and we compare it with the true value of E inv . The average values of the residuals as a function of the true value of E inv are shown in Fig. 3 (right) for proton and iron primaries for EPOS LHC, QGSJET II-04 and QGSJet01 hadronic interaction models. The residuals are within ±10% which is an indication of the overall systematic uncertainty in E inv estimation, which is dominated by the model and mass dependence of the values of C and δ.

E inv ) from vertical showers
As seen in Sec. 2, E inv is a power law function of E 0 β 0 , equal to A 1−β in the extended Heitler model [10] (see Eq. (4)), is a parameter introduced in order to account for the large variations in the predictions of the number of muons that are obtained using different hadronic interaction models once E 0 and primary mass are fixed.
In the reconstruction of vertical events, E 0 is estimated from S (1000), the signal at 1000 m from the core [4], by correcting for the shower attenuation using the constant intensity cut method [17]. To estimate E inv from S (1000), we use the functional form  where ∆X = X − X max is the atmospheric slant depth between ground level at the Auger site and the depth of the shower maximum development, and γ 0 (∆X) is related to the attenuation of S (1000) with ∆X.
Combining Eq. (7) and (8) one obtains where The parameter B and those defining the function A(∆X) are determined using Monte Carlo simulations. Using the QGSJetII-03 hadronic interaction model, we find β = 0.925 and γ = 1.0594, so that their product is B = 0.98. Different interaction models yield the same value of B to within 2%. This value will be used from now on, so that with Eq. (10) and the measurements of S (1000) and ∆X one can obtain an event-by-event estimate of E inv . The function A(∆X) is calculated using events simulated with the QGSJetII-03 hadronic interaction model for a mixed composition of 50% protons and 50% iron. A(∆X) is parameterised with the fourth-degree polynomial. The performance of the analysis is tested with proton and iron events simulated with the hadronic interaction models QGSJetII-03, and EPOS 1.99. The average value of the residuals as a function of E 0 , shown in Fig. 4 (left), are between −5% and 20%. The spread in the residuals is mainly due to the difference in the predictions of the number of muons and of the attenuation function γ 0 (∆X) among the simulations used to parametrise A(∆X), and the ones used to simulate the events. Note that the function γ 0 (∆X) includes the conversion factor needed to obtain E 0 from S (1000) which is strongly model dependent.
A better estimation of E inv can be obtained taking into account these differences using the following equation where the quantities with and without the accent tilde are calculated for the data sample that we are analysing and for the one used to parametrise A(∆X), respectively. β is fixed to 0.925. The functions γ 0 are obtained from Eq. (8) using E 0 and S (1000). The ratioβ 0 /β 0 is estimated from the ratio of the number of muons at ground level for the two data sets, information that is available in the CORSIKA events. The residuals in E inv using the improved parameterisation of Eq. (13) are shown in Fig. 4 (right). The true value of E inv can be recovered within a few % for all models and primaries. Note also how we improve the estimation of E inv for QGSJetII-03, despite the primary mass composition used to parametrise A(∆X) being different to that of the simulated events used to test the analysis method.

Parameterisation of E inv as a function of E cal
The analysis methods described in Sec. 3.1 and 3.2 allow us to obtain an event-by-event estimation of E inv from the data collected by the Pierre Auger Observatory. The analysis is limited to events sufficiently energetic to ensure a full SD trigger efficiency. At energies lower than 4 × 10 18 eV for the inclined [6] and 3 × 10 18 eV for the vertical events [18], the trigger is biased towards events with a higher number of muons, and thus higher E inv and consequently larger systematic uncertainties. In order to get an estimation of E inv useful for all FD events, including the ones with energies below the full SD trigger efficiency, the event-by-event estimation of E inv is parameterised as a function of E cal above the full trigger efficiency, with the function being extrapolated to lower energies.
Analysing hybrid events collected from 1 January 2004 to 31 December 2015, selected with the same cuts used for the energy calibration of the SD energy estimators [7] a parameterisation was obtained.
The correlation between E inv and E cal is well approximated by a power law relationship The data and the fitted function are shown in Fig. 5. For a quantitative comparison of the two data-driven estimations of E inv one has to take into account E inv zenith angle dependence. Since the majority of the events are below 60 • , the E inv parameterisation from the inclined data set, that is on average 5% than the vertical one, has been corrected. The two data-driven E inv estimations are compared in Fig. 6 (left) together with the theoretical predictions for post-LHC hadronic interaction models. They are still larger than the predictions for iron primaries, in contradiction with the mean mass obtained using X max measurements [8]. This is due to the muon deficit [16] as models fail to describe the properties of shower development related to muons and therefore to E inv .
It is worth noting that the two estimates are partially correlated since they both use the measurement of N µ . However, they are affected by different systematics.
The estimations of E inv obtained above the energy of full SD trigger efficiency can be extrapolated to lower energies taking into account the change in the mean mass composition evolution with energy at E A cal ≃ 2 × 10 18 eV measured by Auger [8,20]. The function is obtained by extrapolating the parameterisation obtained from data down to E A cal and, below this energy, using a model inspired function that matches the parameterisation at E A cal . For the latter, we use the function of Eq. (4) in which the mean composition as a function of energy is taken from Figure 6. Left:E inv for inclined and vertical events compared with predictions from simulations. Systematic uncertainty are shown with the shaded bands. The estimate for inclined events is extrapolated to low energies. Right:Auger data-driven estimation of E inv compared with the parameterisations for protons, iron and mixed composition reported in [1] and the one in use by Telescope Array [22].
the Auger FD measurements [20] together with a value of β = 0.9 that reproduces the simulations at lower energies. The extrapolation of E inv obtained from the inclined events, shown with the black dashed line in Fig. 6 (left), will be replaced in the near future with a more accurate estimation of E inv using the data collected by the AMIGA muon detectors [21] installed at the Observatory and using the 750m-spacing sub-array of WCDs [4].

Conclusions
A data driven estimation of E inv of cosmic ray showers detected by the Pierre Auger Observatory, was presented. Two analysis methods for inclined (60 • < θ < 80 • ) and vertical (θ < 60 • ) events were developed. E inv has been parameterised as a function of E cal and extrapolated to energies below the SD full trigger efficiency. The two estimations agree at a level well within the systematic uncertainties, that are estimated to be of the order of 10%−15%.
E inv results are considerably higher than the predictions given by simulations. This is a consequence of the muon deficit in models [16], a deficit due to the failure of the hadronic interaction models to describe the properties of shower development related to muons. Moreover, the results are consistent with the evolution of the mass with energy as measured by Auger [8,20]. This is due to the sensitivity of N µ to the primary mass and, at lower energy, due to the use of the mean mass composition to find the functional form that describes E inv as a function of E cal .
The measurement of N µ makes the analysis of E inv from inclined showers rather straightforward and intrinsically better than the analysis used for vertical events.