The depth of the shower maximum of air showers measured with AERA

. The Auger Engineering Radio Array (AERA), as part of the Pierre Auger Observatory, is an array of radio antennas probing the nature of ultra-high energy cosmic rays at energies around the transition from Galactic to extra-galactic origin. It measures the MHz radio emission of extensive air showers produced by cosmic rays hitting our atmosphere. The elemental composition of cosmic rays is a crucial piece of information in determining what the sources of cosmic rays are and how cosmic rays are accelerated. This composition can be obtained from the mass-sensitive parameter X max , the depth of the shower maximum. We reconstruct X max with a likelihood analysis, by comparing the measured radio footprint on the ground to an ensemble of footprints from Monte-Carlo CORSIKA / CoREAS air shower simulations. We compare our X max reconstruction with ﬂuorescence X max measurements on a per-event basis, a setup unique to the Pierre Auger Observatory, and show the methods to be fully compatible. Furthermore, we extensively validate our reconstruction by identifying and correcting for systematic uncertainties. We determine the resolution of our method as a function of energy and reach a precision better than 15 gcm − 2 at the highest energies. With a bias-free set of around 600 showers, we ﬁnd a light to light-mixed composition at energies between 10 17 . 5 to 10 18 . 8 eV, also in agreement with the Auger ﬂuorescence measurements.


Introduction
The detection of ultra-high energy cosmic rays (UHECRs) is only feasible through indirect detection of extensive air showers. The footprints of the showers, covering extended areas, allow for detection with an array of detectors. Using the captured information the properties of air showers be reconstructed, which can tell us about the cosmic rays and their sources. One key property to differentiate between source models is the particle type, or mass, of the primary cosmic ray. The mass itself however is not a direct observable, but this information can, for example, be extracted from the depth of the shower maximum (X max ). In general, the development of an air shower occurs higher in the atmosphere for a heavy particle than for a light particle, because the heavy particle can be considered a superposition of multiple light particles, each with a fraction of the total energy. These lower-energy particles finish their cascades earlier and hence X max will be shallower for heavier cosmic rays. In this work, we will present results on extracting X max using the radio emission of air showers, compare it with simulations fluorescence measurements, and interpret the results in the context of the mass composition of cosmic rays.

AERA and the Pierre Auger Observatory
The Pierre Auger Observatory [1], located near the town of Malargüe in Argentina, is aimed at detecting UHECRs * e-mail: spokespersons@auger.org * * Full author list at https://www.auger.org/archive/authors_2022_10.html up to the highest energies. It covers an area of 3000 km 2 , making it the largest of its kind in the world. The main components are an array of 1660 water-Cherenkov detectors (WCDs) making up the surface detector (SD) and 27 fluorescence detectors (FD) surrounding and looking over the SD. Located in the west of the surface array is also the Auger Engineering Radio Array (AERA) [2], an array of radio antennas covering roughly 17 km 2 . AERA measures the radio emission produced in air showers at frequencies of 30 to 80 MHz and is sensitive to air showers with energies above roughly 10 17 eV (limited by the radio noise floor) and up to roughly 10 19 eV (exposure-limited by its surface area and current operation time). The combination of radio and various other detection techniques makes the Pierre Auger Observatory an ideal place to study both the properties of cosmic rays and the systematics of their detection techniques.

Reconstruction and Resolution of X max
The reconstruction of X max with AERA [3,4] makes use of simulation-fitting, based on [5], where the radio signals from a measured air shower are compared to the signals from an ensemble of simulated air showers in order to derive X max . After quality cuts, a set of 2153 air showers with a high-quality reconstruction of shower core position, arrival direction, and energy are selected. For each of these showers, these properties are given as input to the air-shower simulation code in order to create representative simulations. In these simulations X max will vary, both  because of the primary particle that initiated the shower, and because of the stochastic nature of shower development. Hence, we simulate 12 showers induced by iron nuclei and 15 induced by protons to get full coverage of the values of X max that the shower could have realistically had (according to the models). Comparison of the simulations to the measured shower then allows us to extract the bestmatching X max value. For the simulations, we use the airshower simulation code Corsika (v7.7100) [6] where the radio signals are calculated with CoREAS [7]. The highenergy interactions are modelled with QGSJetII-04 [8] and the low-energy ones with UrQMD 1.3cr [9]. Furthermore, for each of the measured air showers that simulations are made for we provide an atmospheric model of the date and time of the event using the Global Data Assimilation System (GDAS) atmospheric model [10].
We then select only the showers with energies above 10 17.5 eV where the trigger of the SD, which we use for the triggering of the radio stations, becomes fully efficient [11]. Crucially, we select from this the showers that are free from selection bias, caused by the radio detection and reconstruction, by evaluating if each of the showers could have been reconstructed had they arrived with any other X max value [3,4]. Showers that do not satisfy this criterion are rejected. After all quality and selection cuts a sample of 594 shower with reconstructed X max is obtained.
The X max reconstruction method also provides an estimation of the uncertainty of the reconstructed X max values, obtained by evaluating the X max reconstruction accuracy of simulated air showers with all noise and detector effects included. Figure 1 shows the median resolution of X max as a function of shower energy (green points), for the sample of 594 showers with energies above 10 17.5 eV. To these points, we fit a function (green line), inspired by the reso- lution in the energy of electromagnetic calorimeters [12] where a = 13.8 ± 6.5 g cm −2 , b = 12.8 ± 2.1 g cm −2 , and c = 11.3 ± 4.5 g cm −2 are the fitted parameters and ⊕ indicates the quadratic sum. We achieve a resolution of X max of better than 15 g cm −2 at the highest energies. At lower energies the resolution decreases due to a decreasing signal-to-noise ratio (leading to relatively lower quality signals) and decreasing shower footprint size (leading to relatively less stations per event). For comparison, the Auger FD resolution of X max (grey lines) [13] is plotted, showing the radio technique to be very competitive.

Comparison to Hybrid Fluorescence Data
For 53 of the air showers in our data set X max has also been measured simultaneously with the FD [14]. This number is limited by the requirement of the fluorescence technique to have dark and clear nights and the relatively large distance between AERA and the FD combined with FD fieldof-view cuts. These 53 showers, however, are extremely valuable. They allow for event-to-event comparison of X max , which gives a handle on understanding systematic uncertainties between the two methods independent of selection effects (unlike comparisons of X max distributions). Fig. 2 shows the differences in X max as measured by FD and AERA. Inset at the top of the figure are the mean and width of the differences. For this calculation, each event is weighted by its combined uncertainty on X max : w = (σ 2  . Mean (left) and width (right) of X max for our measurements (black markers), Auger FD measurements [13] (gray markers), and pure proton and iron nuclei composition from air shower simulations using three different hadronic interaction models (red and blue lines) [13]. Statistical uncertainties are shown with error bars and systematic uncertainties (see also Fig. 4) with caps. 5.7 g cm −2 ). The mean difference of −3.9 ± 11.2 g cm −2 shows there is no significant bias between the two methods and sets a limit on systematic uncertainties between the reconstructions in this energy range (predominantly at the lower energy range of AERA, between 10 17.5 and 10 18.0 eV). Furthermore, the compatibility provides an independent check on the fluorescence technique and gives a new handle on investigating shower physics. This new constraint on the X max scale might prove useful in the context of the muon deficit in air shower simulation codes with respect to measurements [15].

Moments of the X max Distribution
We show the two central moments of the X max distribution in Fig. 3 for the set of 594 showers with high-quality X max (black markers). For the interpretation of these results, we show the predictions from air shower simulations with three different hadronic interaction models for a composition of only protons or only iron nuclei (red and blue lines, respectively). AERA finds a relatively light composition in the energy range of 10 17.5 to 10 18.8 eV, in agreement with the Auger FD measurements [13].

Moments of X max
For the calculation of the two central moments of the X max distribution, we account for various sources of systematic uncertainty. The contributions are shown in Fig. 4 as a function of energy. There are several small contributions from the models that we use to simulate our showers; the choice of hadronic interaction model has been shown to vary X max by 5 g cm −2 [16] and the GDAS atmospheric model has a systematic uncertainty of 2.2 g cm −2 [16,17]. Furthermore, our SD energy measurements have a systematic uncertainty, which we propagate to X max , resulting in 2.9 g cm −2 (and 0.3 g cm −2 for σ(X max ) ).
Next, we account for the effects of our X max reconstruction method (red bars in Fig. 4). The reconstruction bias of X max for our method depends on X max itself, hence we quantify, by using the reconstruction of X max with simulated showers, how X max and σ(X max ) vary for compositions ranging from pure proton to pure iron nuclei. The resulting range in X max we quote as, admittedly conservative, systematic uncertainty. We note that when comparing our X max results to a particular composition measurement (e.g., from a different detector) the AERA systematic uncertainties would be much reduced. For example, for the composition as measured by the Auger FD the total AERA systematic uncertainties would reduce to about ±10 g cm −2 at all energies.
For the acceptance cuts in Sec. 3 we determined whether AERA would have been able to detect and reconstruct our measured showers if they would have had a different X max . To be more precise, we required that at least 90% of a Gumbel X max distribution [18,19] for both proton and iron nuclei showers would have been reconstructed. This calculation was performed using the simulation set we have for each measurement [3]. A requirement of 100% would unnecessarily cut many events because they might not have been detected if they had a very unlikely X max value. To still account for this effect we calculate the effect of the non-unity acceptance on the two central moments. This results in less than 4 g cm −2 on X max and 5 g cm −2 on σ(X max ) (blue bars in Fig. 4). We also account for systematic uncertainties from the calculation of the central moments. Because of the relatively low number of events per energy bin these calculations can be sensitive to single outliers. We quantify this with two contributions: the first accounts for events with very small uncertainties (i.e., events where the uncertainty on the X max uncertainty estimation is an under-fluctuation). This can have an especially significant effect on the second moment of X max where one has to remove the detector resolution from the width of the X max distribution. Hence, we perform the calculations of the moments with lower limits for the X max uncertainty, ranging between 5 and 15 g cm −2 , and quote the range of the moments as a systematic uncertainty (see orange bars in Fig. 4). The second effect is for the few events with very large uncertainty (larger than 80 g cm −2 ). Here we repeat the calculations with an event cut on the X max uncertainty at 80 g cm −2 and quote the difference as systematic uncertainty (cyan bars). This latter effect is minor since there are only a few events with higher uncertainties.
Finally, to check for any residual biases that have not been accounted for by the previous contributions we check for any dependencies in X max as a function of geometry parameters such as arrival direction (zenith and azimuth angle), shower core position, and angle to the geomagnetic field. These represent the effects on acceptance from, respectively, the asymmetric directional sensitivity of the radio antennas, the irregular spacing of radio antennas, and the dependence of the radio emission strength on the angle between the magnetic field and the shower axis. It is conceivable that these effects could lead to part of the showers never being seen at all (even though we checked with the acceptance that each of our showers that we do see would have been seen at any X max ). We have checked the trends of X max as a function of these geometry parameters and find a maximal residual bias of 3 to 10 g cm −2 (depending on shower energy) is allowed within statistical uncertainties. To be conservative we include this as a systematic uncertainty (see green bars in Fig. 4).

Conclusions
We have presented the results of the measurement of X max with AERA at the Pierre Auger Observatory. We observe a rather light composition at energies between 10 17.5 and 10 18.8 eV, the energy range where the origin of cosmic rays is expected to transition from Galactic to extragalactic. With our results, we also show compatibility with Auger fluorescence measurements of the central moments of X max with a set of 594 air showers and show the compatibility of X max measurements of individual showers with a set of 53 events measured simultaneously with both detectors. This is strong independent support for the validity of both techniques. Furthermore, we have shown our method is able to reach a precision of X max reconstruction better than 15 g cm −2 at high energies, competitive with other techniques.