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

TL;DR

This paper presents a new method to directly measure the invisible energy in cosmic ray air showers using data from the Pierre Auger Observatory, reducing uncertainties in energy estimation.

Contribution

It introduces a direct measurement approach for invisible energy based on muon counts, minimizing reliance on model-dependent simulations.

Findings

Reduced systematic uncertainties in energy estimates

Improved accuracy in cosmic ray energy scale

Method applicable to extensive air shower data

Abstract

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.

Equations

$${ϵ}_c^{\pi }=\frac{E_0}{N^n},$$

$$E_{\mathrm{inv}}={ϵ}_c^{\pi }\text{leavevmode}{N}_{\mu },$$

$$E_{\mathrm{inv}}={ϵ}_c^{\pi }{\left(\frac{E_0}{{ϵ}_c^{\pi }}\right)}^{\beta },$$

$$E_{\mathrm{inv}}^A={ϵ}_c^{\pi }{\left(\frac{E_0}{{ϵ}_c^{\pi }}\right)}^{\beta }{A}^{1-\beta }.$$

$${\rho }_{\mu }(\vec{r})=N_{19}\text{leavevmode}{\rho }_{\mu ,19}(\vec{r};\theta ,\varphi ),$$

$$E_{\mathrm{inv}}=C\text{leavevmode}{\left(R_{\mu }^{\mathrm{MC}}\right)}^{\delta },$$

$$E_{\mathrm{inv}}={ϵ}_c^{\pi }{\beta }_0{\left(\frac{E_0}{{ϵ}_c^{\pi }}\right)}^{\beta }.$$

$$E_0={\gamma }_0(\Delta X){\left[S(1000)\right]}^{\gamma },$$

$$E_{\mathrm{inv}}$$

$$=$$

$$A(\Delta X)$$

$$B$$

$$E_{\mathrm{inv}}=A(\Delta X)\text{leavevmode}{\left[S(1000)\right]}^B\text{leavevmode}{\left(\frac{{\tilde{\gamma }}_0(\Delta X)}{{\gamma }_0(\Delta X)}\right)}^{\beta }\text{leavevmode}\frac{{\tilde{\beta }}_0}{{\beta }_0},$$

$$E_{\mathrm{inv}}=a{\left(\frac{E_{\mathrm{cal}}}{10^{18}\text{leavevmode}\mathrm{eV}}\right)}^b.$$

Full text

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.

1 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 $\int\left(dE/dX\right)dX$, called the calorimetric energy of the shower, $E_{\rm cal}$, is measured.

$E_{\rm 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_{\rm cal}$ a correction to account for the invisible energy ($E_{\rm inv}$) carried by the particles that do not dissipate all their energy in the atmosphere. $E_{\rm inv}$ amounts to about 10% - 20% of $E_{0}$.

$E_{\rm inv}$ can be calculated directly from the energy deposited in the atmosphere by the different components of simulated air showers Barbosa . There are large differences in the values of the ratio of $E_{\rm inv}$ to $E_{0}$ as a function of $E_{\mathrm{cal}}$ for different hadronic models and primary masses, as could be seen in Fig. 1. The estimation of $E_{\rm 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_{\mathrm{inv}}$ can be improved further using the primary mass composition estimated with the fluorescence detectors Auger-mass-comp 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 Pierog-ICRC2017 . Even after the updates with LHC data, the models still fail to describe several properties of the shower development related to muons Auger-HadInt-UHECR18 , and this can introduce unpredictable biases in the $E_{\mathrm{inv}}$ estimation.

Thus the strategy followed in this work is to estimate $E_{\rm inv}$ using the correlations that exist between $E_{\mathrm{inv}}$ and shower observables measured at the Pierre Auger Observatory Auger , correlations that to a large extent are not sensitive to the hadronic interaction models and primary mass composition.

The Pierre Auger Observatory Auger 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 $\sim 3000\leavevmode\nobreak\ {\rm 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% AugerFD .

2 Phenomenology of the invisible energy

The Heitler model and its extension to hadronic cascades Heitler ; Matthews 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_{\mathrm{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:

[TABLE]

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_{\mathrm{inv}}$ is proportional to the number of muons ($N_{\mu}$) reaching ground level.

[TABLE]

This expression will be the guiding thread to estimate $E_{\mathrm{inv}}$ with a measurement of $N_{\mu}$ in inclined showers.

Another important feature of the model is the power law of $E_{\mathrm{inv}}$:

[TABLE]

where $\beta=\ln(\frac{2}{3}N)/\ln N$. Air shower simulations predict values of $\beta$ in the range from 0.88 to 0.92 K-U-review . $\beta$ also fixes the $E_{\mathrm{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:

[TABLE]

This relationship will be the guiding thread to estimate $E_{\mathrm{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_{\mathrm{inv}}$, which is calculated following the method described in Barbosa .

The correlation between $E_{\mathrm{inv}}$ and $N_{\mu}$ has been studied simulating showers with different primary masses and hadronic interactions models using the CORSIKA code (CORSIKA, , & references therein), as it is shown in Fig. 2. In spite of the very large spread in the predictions of $N_{\mu}$ and $E_{\mathrm{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_{\mathrm{inv}}$ from the measurements of $N_{\mu}$.

3 Estimation of $E_{\mathrm{inv}}$ using Auger data

Two different reconstruction techniques are used for the SD events: one for the so-called vertical showers with zenith angles $\theta<60^{\circ}$ SDExp , and one for the inclined showers with $\theta>60^{\circ}$ HASRec . WCDs are sensitive to the electromagnetic and hadronic components of a shower.

The most straightforward way to estimate $E_{\rm 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_{\rm inv}$, as seen in Sec. 2 (Eq. (2)). The muon number cannot be directly measured for vertical events. However, $E_{\rm inv}$ can be obtained from the energy estimator using the power law relationship between $E_{\rm inv}$ and $E_{0}$ (see Eq. (4)).

3.1 $E_{\rm inv}$rom inclined showers

The reconstruction of inclined events HASRec 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 $(\theta,\phi)$ 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 $\vec{r}$ is then parameterised with

[TABLE]

where $\rho_{\mu,19}(\vec{r};\theta,\phi)$ 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_{\mu}^{\rm MC}$ for each simulated event. $R_{\mu}^{\rm 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_{\mu}^{\rm MC}$ is within 5% for several hadronic interaction models and primaries AugerMuonSize . 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_{\mu}$, represents an unbiased estimator of the total number of muons at ground level.

The correlation between $E_{\rm 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 (CORSIKA, , & references therein) and the other with AIRES using the model QGSJet01 (AIRES, , & references therein).

For each simulated event, we calculate the values of $E_{\rm inv}$ and of the muon number at ground level $R_{\mu}^{\rm MC}$.

For all the samples of simulated events, the correlation between $E_{\rm inv}$ and $R_{\mu}^{\rm MC}$ is well described by a power-law

[TABLE]

where the values of the parameters $C$ and $\delta$ are obtained from a fit to the events. Examples of the correlation between $E_{\rm inv}$ and $R_{\mu}^{\rm 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_{\rm inv}$ in the data from the measurement of $R_{\mu}$ that, as seen before, is the unbiased estimator of $R_{\mu}^{\rm MC}$. Since the mass composition of the data is not precisely known, the estimation of $E_{\rm inv}$ is obtained using the parameterisation of $E_{\rm inv}$ as a function of $R_{\mu}$ for a mixture of 50% protons and 50% iron. This is done taking the average of the two $E_{\rm 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 events111 events simulated with the detector response and $R_{\mu}$ reconstructed with the same algorithm used for the data. For each event, we compute $E_{\rm inv}$ from $R_{\mu}$ using the estimation for the mixed proton and iron composition, and we compare it with the true value of $E_{\rm inv}$. The average values of the residuals as a function of the true value of $E_{\rm 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 $\pm 10\%$ which is an indication of the overall systematic uncertainty in $E_{\rm inv}$ estimation, which is dominated by the model and mass dependence of the values of $C$ and $\delta$.

3.2 $E_{\rm inv}$) from vertical showers

As seen in Sec. 2, $E_{\mathrm{inv}}$ is a power law function of $E_{0}$

[TABLE]

$\beta_{0}$, equal to $A^{1-\beta}$ in the extended Heitler model Matthews (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 Auger , by correcting for the shower attenuation using the constant intensity cut method CIC . To estimate $E_{\rm inv}$ from $S(1000)$, we use the functional form

[TABLE]

where $\Delta X=X-X_{\rm max}$ is the atmospheric slant depth between ground level at the Auger site and the depth of the shower maximum development, and $\gamma_{0}(\Delta X)$ is related to the attenuation of $S(1000)$ with $\Delta X$.

Combining Eq. (7) and (8) one obtains

[TABLE]

where

[TABLE]

The parameter $B$ and those defining the function $A(\Delta X)$ are determined using Monte Carlo simulations. Using the QGSJetII-03 hadronic interaction model, we find $\beta=0.925$ and $\gamma=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 $\Delta X$ one can obtain an event-by-event estimate of $E_{\rm inv}$. The function $A(\Delta X)$ is calculated using events simulated with the QGSJetII-03 hadronic interaction model for a mixed composition of 50% protons and 50% iron. $A(\Delta 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 $\gamma_{0}(\Delta X)$ among the simulations used to parametrise $A(\Delta X)$, and the ones used to simulate the events. Note that the function $\gamma_{0}(\Delta X)$ includes the conversion factor needed to obtain $E_{0}$ from $S(1000)$ which is strongly model dependent.

A better estimation of $E_{\mathrm{inv}}$ can be obtained taking into account these differences using the following equation

[TABLE]

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(\Delta X)$, respectively. $\beta$ is fixed to 0.925. The functions $\gamma_{0}$ are obtained from Eq. (8) using $E_{0}$ and $S(1000)$. The ratio $\tilde{\beta}_{0}/\beta_{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_{\mathrm{inv}}$ using the improved parameterisation of Eq. (13) are shown in Fig. 4 (right). The true value of $E_{\rm inv}$ can be recovered within a few % for all models and primaries. Note also how we improve the estimation of $E_{\mathrm{inv}}$ for QGSJetII-03, despite the primary mass composition used to parametrise $A(\Delta X)$ being different to that of the simulated events used to test the analysis method.

3.3 Parameterisation of $E_{\rm inv}$ as a function of $E_{\rm cal}$

The analysis methods described in Sec. 3.1 and 3.2 allow us to obtain an event-by-event estimation of $E_{\rm 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\times 10^{18}$ eV for the inclined HASRec and $3\times 10^{18}$ eV for the vertical events SDExp , the trigger is biased towards events with a higher number of muons, and thus higher $E_{\rm inv}$ and consequently larger systematic uncertainties. In order to get an estimation of $E_{\mathrm{inv}}$ useful for all FD events, including the ones with energies below the full SD trigger efficiency, the event-by-event estimation of $E_{\rm inv}$ is parameterised as a function of $E_{\rm 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 Auger-EnSc a parameterisation was obtained.

The correlation between $E_{\rm inv}$ and $E_{\rm cal}$ is well approximated by a power law relationship

[TABLE]

The data and the fitted function are shown in Fig. 5.

For a quantitative comparison of the two data-driven estimations of $E_{\mathrm{inv}}$ one has to take into account $E_{\mathrm{inv}}$ zenith angle dependence. Since the majority of the events are below $60^{\circ}$, the $E_{\rm inv}$ parameterisation from the inclined data set, that is on average 5% than the vertical one, has been corrected. The two data-driven $E_{\mathrm{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_{\rm max}$ measurements Auger-mass . This is due to the muon deficit AugerMuonSize as models fail to describe the properties of shower development related to muons and therefore to $E_{\rm inv}$.

It is worth noting that the two estimates are partially correlated since they both use the measurement of $N_{\mu}$. However, they are affected by different systematics.

The estimations of $E_{\mathrm{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_{\rm cal}^{A}\simeq 2\times 10^{18}\rm{eV}$ measured by Auger Auger-mass ; Auger-mass-comp-ICRC17 . The function is obtained by extrapolating the parameterisation obtained from data down to $E_{\rm cal}^{A}$ and, below this energy, using a model inspired function that matches the parameterisation at $E_{\rm cal}^{A}$. For the latter, we use the function of Eq. (4) in which the mean composition as a function of energy is taken from the Auger FD measurements Auger-mass-comp-ICRC17 together with a value of $\beta=0.9$ that reproduces the simulations at lower energies. The extrapolation of $E_{\mathrm{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_{\mathrm{inv}}$ using the data collected by the AMIGA muon detectors AMIGA installed at the Observatory and using the 750m-spacing sub-array of WCDs Auger .

4 Conclusions

A data driven estimation of $E_{\rm inv}$ of cosmic ray showers detected by the Pierre Auger Observatory, was presented. Two analysis methods for inclined ($60^{\circ}<\theta<80^{\circ}$) and vertical ($\theta<60^{\circ}$) events were developed. $E_{\rm inv}$ has been parameterised as a function of $E_{\rm 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_{\rm inv}$ results are considerably higher than the predictions given by simulations. This is a consequence of the muon deficit in models AugerMuonSize , 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 Auger-mass ; Auger-mass-comp-ICRC17 . This is due to the sensitivity of $N_{\mu}$ 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_{\mathrm{inv}}$ as a function of $E_{\mathrm{cal}}$.

The measurement of $N_{\mu}$ makes the analysis of $E_{\mathrm{inv}}$ from inclined showers rather straightforward and intrinsically better than the analysis used for vertical events.