Inference of the Mass Composition of Cosmic Rays with energies from $\mathbf{10^{18.5}}$ to $\mathbf{10^{20}}$ eV using the Pierre Auger Observatory and Deep Learning

Abstract

We present measurements of the atmospheric depth of the shower maximum $X_{\mathrm{max}}$, inferred for the first time on an event-by-event level using the Surface Detector of the Pierre Auger Observatory. Using deep learning, we were able to extend measurements of the $X_{\mathrm{max}}$ distributions up to energies of 100 EeV ($10^{20}$ eV), not yet revealed by current measurements, providing new insights into the mass composition of cosmic rays at extreme energies. Gaining a 10-fold increase in statistics compared to the Fluorescence Detector data, we find evidence that the rate of change of the average $X_{\mathrm{max}}$ with the logarithm of energy features three breaks at $6.5\pm 0.6~{}(\mathrm{stat})\pm 1~{}(\mathrm{sys})$ EeV, $11\pm 2~{}(\mathrm{stat})\pm 1~{}(\mathrm{sys})$ EeV, and $31\pm 5~{}(\mathrm{stat})\pm 3~{}(\mathrm{sys})$ EeV, in the vicinity to the three prominent features (ankle, instep, suppression) of the cosmic-ray flux. The energy evolution of the mean and standard deviation of the measured $X_{\mathrm{max}}$ distributions indicates that the mass composition becomes increasingly heavier and purer, thus being incompatible with a large fraction of light nuclei between 50 EeV and 100 EeV.

I Introduction

The arrival directions, energy spectrum, and mass composition are the three important pillars of cosmic ray research. A sound interpretation of the three measurements and their energy dependence, both individually and jointly, is pivotal for a deep understanding of the nature of cosmic rays, including their origin and propagation, and enables the study of astrophysical models. With energies larger than $1\,$EeV ($10^{18}\,$eV), ultra-high-energy cosmic rays (UHECRs) are the most energetic particles ever measured by humankind. One of the lasting puzzles is the origin of the suppression of the cosmic-ray flux observed at around 50 EeV[1, 2, 3, 4]. A precise measurement of the UHECR mass composition can deliver insights into whether the suppression is caused by the interaction of the particles with the cosmic background photons [5, 6], a sign of the maximum energy reached in cosmic accelerators [7], or a combination of both [8, 9]. Due to the low flux at ultra-high energies, the primary composition cannot be measured directly but can only be studied by indirectly analyzing the properties of the induced air showers. Information on the primary mass can be obtained by measuring the atmospheric depth of the shower maximum $X_{\mathrm{max}}$, the depth at which the number of secondary particles reaches its maximum. Investigating the measured $X_{\mathrm{max}}$ distribution, as a function of energy, in terms of its mean and standard deviation (fluctuations), $\langle X_{\mathrm{max}}\rangle$ and $\sigma(X_{\mathrm{max}})$, enables us to study the UHECR mass composition [10]. Heavier particles feature, on average, a smaller $X_{\mathrm{max}}$ since sub-showers are sharing the primary energy, resulting in a maximum higher in the atmosphere. This motivates the investigation of the first moment $\langle X_{\mathrm{max}}\rangle$ of the $X_{\mathrm{max}}$ distribution. Further, the shower-to-shower fluctuations, i.e., the second moment $\sigma(X_{\mathrm{max}})$ of the distribution, is also a mass-sensitive observable. Due to the smaller cross-section and the development of fewer sub-showers, cascades induced by lighter primary particles are subject to larger fluctuations. Furthermore, compared to $\langle X_{\mathrm{max}}\rangle$, the fluctuations $\sigma(X_{\mathrm{max}})$ are sensitive to both the primary mass and the degree of mixing of the primary beam [11] and are almost insensitive to the uncertainties in the hadronic interaction models.

Using fluorescence telescopes, $X_{\mathrm{max}}$ can be directly reconstructed by observing the longitudinal shower development. Nevertheless, due to the observations being confined to dark and moonless nights, the duty cycle is limited. In contrast, sparse surface-detector arrays have a duty cycle close to $100\%$ and sample the secondary shower particles at the ground. Thus, they cannot directly observe $X_{\mathrm{max}}$, making its reconstruction challenging. However, information about the shower development and $X_{\mathrm{max}}$ is contained in the lateral number density and distribution of arrival times of particles reaching the ground. By studying the risetimes of the time-dependent signals, conclusions on the average composition have already been drawn in the past [12]. However, to infer the UHECR mass composition beyond mere $\langle X_{\mathrm{max}}\rangle$ measurements, more sophisticated methods are needed to fully exploit the complex data. The advent of deep learning [13, 14] provides new analysis techniques for large and complex data sets. First approaches have already been successfully applied to LHC data [15] and physics in general [16]. The recent progress offers supplementary and improved reconstruction algorithms for neutrino [17, 18, 19] and cosmic-ray observatories [20]. This includes the deep-learning-based reconstruction of $X_{\mathrm{max}}$ [21, 22, 23, 24] and muon signals [25] using the temporal structure of signals measured by the Surface Detector of the Pierre Auger Observatory.

In this work, we use this novel reconstruction technique to study the mass composition of UHECRs in terms of $\langle X_{\mathrm{max}}\rangle$ and $\sigma(X_{\mathrm{max}})$ in the energy range from 3 to 100 EeV. With about $50{,}000$ events, this is, at the time, the most comprehensive study of the UHECR mass composition and the first measurement of $\sigma(X_{\mathrm{max}})$ beyond 50 EeV. A comprehensive discussion of the analysis, including the technical details of the novel $X_{\mathrm{max}}$ reconstruction, is given in an accompanying publication [26].

II Methodology

In the past two decades, our understanding of UHECRs has grown enormously due to the construction of the Pierre Auger Observatory [27] and the Telescope Array Project [28]. The Pierre Auger Observatory is the world’s largest UHECR instrument and features a hybrid detection technique, being a combination of surface detectors and fluorescence telescopes, to measure cosmic-ray-induced air showers. In total, 1660 water-Cherenkov detectors, spanning 3000 km², are arranged in a triangular 1500-meter-grid and form the Surface Detector (SD) — the centerpiece of the Observatory with a duty cycle close to 100%. The SD is overlooked by 24 telescopes located at four sites that form the Fluorescence Detector (FD) of the Observatory. Additionally, three high-elevation telescopes overlook an infilled array of 61 stations with 750 m spacing that enable measurements below 3 EeV. The FD detects the emitted fluorescence light from air molecules excited by the particle shower penetrating the atmosphere. The requirement for dark and moonless nights limits the duty cycle of the FD to about $15\%$.

The typical size of an air-shower footprint with $E>10$ EeV amounts to tens of $\mathrm{km}^{2}$, and it usually triggers more than ten stations of the SD. Each of these stations is equipped with three photomultiplier tubes (PMTs) that record the time-dependent responses to shower particles digitized and sampled in steps of $25$ ns. The resulting three traces are then calibrated in units of VEM (vertical equivalent muons), i.e., the average signal produced by muons traversing the detector vertically, provided by an in-situ calibration on a minute timescale using atmospheric muons. Several station-level measurements characterize each event in our analysis: the arrival time of the first particles at the respective station and, for each PMT, a trace of $3~{}\upmu$s time length (120 time steps) containing the signal.

In this work, we use two different data sets: a hybrid data set, featuring both an FD and SD reconstruction used to calibrate the reconstruction algorithm to the $X_{\mathrm{max}}$ scale of the FD, and the full SD data set for performing the high-statistics measurement of $\langle X_{\mathrm{max}}\rangle$ and $\sigma(X_{\mathrm{max}})$. To guarantee high data quality, we apply a dedicated event selection. We only select events with an energy $E_{\mathrm{SD}}>3$ EeV to ensure full trigger efficiency [29]. Additionally, we require a zenith angle $\theta<60^{\circ}$, a hexagon of working stations around the station with the largest signal [30]. Furthermore, a fiducial SD cut is applied [26] to ensure an unbiased $X_{\mathrm{max}}$ measurement, accepting only events inside a zenith-angle range where the absolute $X_{\mathrm{max}}$ reconstruction bias of the deep-learning-based algorithm is smaller than $10$ $\mathrm{g\,cm^{-2}}$. This fiducial zenith-angle range depends on energy and was derived from simulations [26]. The SD data set comprises events collected between 1 January 2004 and 31 August 2018. After selection, the data set consists of $48{,}824$ events. For the calibration of the novel algorithm, hybrid events featuring both an FD and SD reconstruction are used. Hence, in addition to the previous cuts, a strict FD selection is applied. We accept only events with good atmospheric conditions and small uncertainties on the observed shower profile. In particular, we reject events with $X_{\mathrm{max}}$ reconstructed outside the telescope field of view. Since $X_{\mathrm{max}}$ depends on the primary mass, such a cut can cause a selection bias. To avoid such a bias, we apply a fiducial cut that ensures uniform acceptance for most of the $X_{\mathrm{max}}$ distribution [10]. $1642$ hybrid events remain after selection.

Figure 1: Application of the DNN to hybrid data. Correlation between fluorescence observations of the FD and DNN predictions using SD data after calibration. The 1642 events show a clear correlation of $\rho=0.7$ and a bias $<1~{}\mathrm{g\,cm^{-2}}$.

Information on the shower development and the primary particle mass is encoded in the temporal structure of the recorded SD signals, i.e., the signal traces and arrival times [12, 31]. The novel $X_{\mathrm{max}}$ reconstruction applied in this work is based on a deep neural network (DNN) tailored to the particular situation of the SD. Two connected sub-networks form the DNN architecture to exploit the patterns of different shower components in the time-resolved particle density. For example, muons usually produce signal spikes, while the signals from each electron, positron, and photon are individually smaller and are spread out in time because of multiple scattering [32]. In the first part, the signal traces are analyzed using long-short term memory (LSTM) layers [33] on a station-by-station level. The next part of the DNN is used to exploit the spatial distribution of the signal footprint induced on the SD grid by combining the outputs of the first part. Therefore, we utilize hexagonal convolutions [34] to account for the structure of the detector grid. We additionally use residual connections that improve the training behavior and performance of deep networks [35, 36]. Finally, the DNN predicts a single event-wise value for $X_{\mathrm{max}}$ [26].

(a)

(b)

Figure 2: Energy evolution of (a) the mean depth of shower maximum $\langle X_{\mathrm{max}}\rangle$ and (b) the fluctuations of shower maximum $\sigma(X_{\mathrm{max}})$ as determined using the FD reconstruction (grey open squares) [37] and the SD-based DNN predictions (black circles). Red (blue) lines indicate expectations for a pure proton (iron) composition for various hadronic models. The number of events in each bin is indicated in panel (b).

The DNN was trained using the simulated detector responses [38] of $400{,}000$ showers induced by proton, helium, oxygen, and iron with energies from 1 to 160 EeV following an $E^{-1}$ spectrum. The showers were simulated with CORSIKA [39] using the EPOS-LHC interaction model. Realistic operation conditions of the detector are imitated [21] during training and act as a form of data augmentation [14]. Whilst training, the DNN learns to infer an event-by-event $X_{\mathrm{max}}$ value following the simulated $X_{\mathrm{max}}$ distribution using the characteristic signal patterns of representative cosmic ray masses. For more details on the algorithm, its training, and its design, we refer to Ref. [21].

During the Observatory lifetime, the shapes of the recorded SD signals have changed due to the aging of the PMTs, the electronics, the reflective Tyvek liner, and potential changes in the water transparency. The degrading of the signal shapes is monitored using the $A/P$ (area-over-peak) observable, which is the ratio of pulse integral and pulse amplitude for a signal produced by a single muon. It is modeled as constant in the utilized simulations but decreases with time in the array [40, 41], causing a decrease in the $X_{\mathrm{max}}$ predictions when applied to measured data. The predictions and their aging dependence on $A/P$ can be described and corrected via a parameterization [21, 26]. In addition, we calibrate the DNN to remove a dependence of $X_{\mathrm{max}}$ on the distance of the station with the largest signal to the reconstructed shower core. Furthermore, we examine seasonal and diurnal variations of the $X_{\mathrm{max}}$ reconstruction and correct it by fitting an oscillating function. The amplitudes found are below $2~{}\mathrm{g\,cm^{-2}}$.

Finally, we use hybrid data to calibrate the SD-based DNN algorithm to the FD $X_{\mathrm{max}}$ scale. Since UHECRs feature energies several orders of magnitude above what can be reached with human-built accelerators, air-shower simulations need to make use of extrapolations of accelerator data. The phenomenological modeling and the extrapolations used differ for each hadronic interaction model. As fluorescence telescopes directly observe $X_{\mathrm{max}}$, they offer the possibility of removing the dependence of the SD-based algorithm on the particular interaction model and significantly reduce the systematic uncertainties of the $\langle X_{\mathrm{max}}\rangle$ measurement.