Mass composition of cosmic rays with energies from [1017.2]eV to [1020]eV using surface and fluorescence detectors of the Pierre Auger Observatory

Ultra-high-energy cosmic rays (UHECRs) are highly energetic particles with []EeV energies, exceeding the capabilities of man-made colliders. They hold information on extreme astrophysical processes that create them and the medium they traverse on their way towards Earth. However, their mass composition at such energies is still unclear, because data interpretation depends on our choice of high energy hadronic interaction models. With its hybrid detection method, the Pierre Auger Observatory has the possibility to detect extensive air showers with an array of surface water-Cherenkov stations (SD) and fluorescence telescopes (FD). We present recent mass composition results from the Pierre Auger Collaboration using observational parameters from SD and FD measurements. Using the full dataset of the Pierre Auger Observatory, implications on composition can be made for energies above [1017.2]eV .


Introduction
With ultra-high-energy cosmic rays reaching collision center-of-mass energies above [40]T eV we are able to observe interactions that are currently unreachable with man-made colliders.The mass composition at such high energies is at present unclear, but holds important information on cosmic ray sources and acceleration processes.For example, mass composition studies are used for the identification of energy spectrum features, such as the transition from galactic to extragalactic sources and the flux suppression at the highest energies.Detection of UHECRs is only possible indirectly through cascades of particles produced in collisions of cosmic rays with atmospheric molecules, resulting in extensive air showers (EAS).Consequently, any information on the initial UHECR depends on our capability to reconstruct its energy, direction, mass and charge.With energies surpassing those of colliders, the collision cross-sections of particles in an EAS can only be extrapolated from existing data, making mass composition studies highly dependent on comparison to simulations.The Pierre Auger Observatory is the largest observatory of UHECRs, located in the Pampa Amarilla near Malargüe, Argentina.It has been continuously taking data since January 2004 and consists of an array of water-Cherenkov surface detectors and a collection of fluorescence telescopes.The surface detector (SD) covers an area of ∼ [3000]km 2 and its 1600 water-Cherenkov stations are arranged in a hexagonal grid with a [1500]m separation.Each station is filled with 12 tonnes of deionized water and three [9]inch photomultiplier tubes (PMTs) observing Cherenkov light that is produced by passing shower particles.Operational time for the SD is nearly 100%, with parts of it selectively removed from analysis during lightning strikes.In 2008, a low-energy upgrade of 61 stations separated by [750]m was built, which further reduces the lower energy limit from [10 18.5 ]eV to [10 17.5 ]eV.Twenty-four fluorescence telescopes located at the periphery of the SD array constitute the fluorescence detector (FD) of the observatory.Measuring fluorescent light ( [300]nm − [450]nm) from deexcitation of nitrogen molecules in the atmosphere, each telescope focuses light onto a 440 PMT pixel camera.The telescopes have a field-of-view between 0 • and 30 • in elevation and have a low-energy limit of [10 17.8 ]eV.A low-energy upgrade called HEAT, built in 2009 near the Coihueco FD building, extends the limit down to [10 17.2 ]eV by covering elevations between 30 • and 60 • above the [750]m SD array.Although its detection is almost calorimetric, the FD is only operational ∼15% of the time, restricted by the amount of ambient light and weather conditions.Using both detection methods, the observatory is able to detect EAS in hybrid mode, giving superior reconstruction of the shower geometry and consequently better energy and zenith angle estimations.

Mass composition studies
A primary UHECR has an energy E far greater than the binding energy of its nucleons.Therefore, the energy can be supposed equally distributed among the nucleons as E = E/A.Compared to a lighter primary, a heavier primary with mass A will therefore interact earlier in the atmosphere and produce a much wider shower front (due to more scattering centers).In order to infer the mass of a primary particle, we must select observational parameters (observables) detected by the SD or FD, which are sensitive to the mass composition of primaries.A range of observables has been identified in [1], including the depth of shower maximum X max from the FD, and lateral distribution function, risetime t 1/2 , number of particles at the ground and muon production depth from the SD.The two latter observables require a good measurement of muons at ground level, which is difficult to separate from the electromagnetic part.An upgrade called AugerPrime will add [4]m 2 scintillator detectors on top of each SD station in order to better separate the muonic and electromagnetic contents of an EAS [2].Some analyses on mass composition have already been made with the asymmetry of risetime and muon production depth and are described in greater detail in [3,4].The observable X max so far has the best discrimination power among different primary particle types, but can only be applied to a smaller dataset consisting of only hybrid events.Risetime, on the other hand, is purely an SD observable and enables the use of the complete dataset of the Pierre Auger Observatory, but includes a mixture of muonic and electromagnetic signals.

Risetime and the Delta method
SD stations measure the energy deposited by shower particles in the form of Cherenkov light.This signal trace holds information on passing particles, with muons producing narrower pulses and electromagnetic particles producing a wider pulse.The risetime t 1/2 is defined as the time it takes for the integrated signal to increase from 10% to 50% of its final magnitude.With each PMT measuring the signal in high-gain and low-gain traces (in case the signal is saturated in the high-gain trace), each triggered station produces up to three values of t 1/2 .With showers induced by heavier primaries developing earlier in the atmosphere than those of lighter primaries, they have a larger muonic content at ground level and thus a shorter risetime.Additionally, risetime depends on the distance from shower axis, zenith angle and Because each event possesses multiple values for risetime, it is important to express them using a single reference value that is mass composition sensitive.This can be achieved with the Delta method [5] by defining a new observable where N is the number of stations triggered by the EAS, ∆ i is the Delta value for each station, t 1/2 is the measured risetime, t bench 1/2 is a benchmark fit of risetimes for a reference energy bin, and σ 1/2 is the average uncertainty on risetime measurements.Dividing data into zenith angle bins and calculating benchmark fits for each of them, removes the dependence of ∆ s on zenith angle and distance from shower axis.

Depth of shower maximum
Fluorescent light that is measured by FD telescopes is mostly caused by interactions of electromagnetic particles with nitrogen molecules, making the depth of shower maximum X max a nearly calorimetric measurement.The longitudinal profile of an EAS describes its development along the shower axis.It reaches its maximum when shower particles have reached their critical energies, which suppresses production of new particles beyond this point.Combining Coihueco FD and HEAT measurements into a HeCo dataset, with a wider field-of-view, gives a better estimation of X max .

Data and analysis
Data from the FD consists of both the standard FD dataset (1 December 2004 -31 December 2015) and the HeCo dataset (1 June 2010 -31 December 2015), with a combined energy range extending above [10 17.2 ]eV.Only hybrid events were used for this analysis, provided measurement conditions were stable, X max was inside the field-of-view of the detector and the resolution was better than [40]g/cm 2 (fiducial and quality cuts).In total, the number of events that passed selection cuts are 25688 and 16778 for standard FD and HeCo, respectively [6].For the SD, the data encompasses both

Estimating X max moments
After the selection process, all events are distributed into nine energy bins for the HeCo and 18 energy bins for the standard FD dataset.Gathered distributions are then corrected for detector resolution and acceptance, and the first two moments of X max are calculated as described in [7].The energy evolution of X max and σ (X max ) describes the average depth of shower maximum and the shower-to-shower fluctuations of the observable, respectively.For comparison, simulations for proton and iron primaries using three different hadronic interaction models are added [6].
Both X max and σ (X max ) are shown in Fig. 1 as a function of energy, with added simulation predictions in order to estimate the mass composition.From simulations, we can see that a constant mass composition has the elongation rate (rate of change of X max ) of ∼ [60]g/cm 2 /decade.There is a clear break at [10 18.33±0.02]eV, with estimated primary mass becoming lighter with increasing energy below it ( [79 ± 1]g/cm 2 /decade) and heavier above it ( [26 ± 2]g/cm 2 /decade).Similarly, fluctuations (Fig. 1, right) up to [10 18.3 ]eV roughly follow the rate for protons and then decrease towards heavier composition with increasing energy.An older set of Pierre Auger Observatory FD measurements have also been compared to data from the Telescope Array (TA).Although their measurement techniques are similar, the treatment of data is different, making direct comparisons difficult.For this reason, the analysis in [8] treats the Pierre Auger Observatory data by using the same reconstruction and analysis approach as for the Middle Drum detector of TA.The final comparison, shown in Fig. 2, equates both methods and results in an adequate agreement between the results of the two observatories.

Estimating evolution of ∆ s with energy
Calculated ∆ s values for each event are distributed into ten energy bins for the [750]m and 14 energy bins for the [1500]m array.Both datasets are assigned a benchmark energy bin, for which benchmark fits are defined, where the value of ∆ s is by definition zero.∆ s as a function of energy is shown in Fig. 3, with added simulation predictions.Compared to For comparison, simulations for proton and iron primaries using two different hadronic interaction models are added.Benchmarking energy bins are [ [10 17.7 ]eV, [10 17.8 ]eV] on the left and [ [10 19.1 ]eV, [10 19.2 ]eV] on the right figure [5].
simulations, the composition shows a tendency of becoming lighter with increasing energy for the [750]m array and becoming heavier for the [1500]m array.The break occurs at around [10 18.3 ]eV.

Results and composition implications
In order to estimate the mass composition of primary UHECRs, we must convert ∆ s and X max into average logarithmic mass ln A , and σ (X max ) into σ (ln A).For FD measurements, a parameterization described in [9] where X max p and σ 2 sh are mean values for proton showers.The mean primary mass and spread of masses using three different hadronic interaction models are shown in Fig. 4. The break in the composition is still present at [10 18.33±0.02]eV, while variances become predominantly smaller with increasing energy above [10 18.3 ]eV.In some regions variances also have an unphysical negative value, which indicates that models predict a broader spread of masses than what is expected from data.Assuming the validity of the superposition model, the average logarithmic mass for SD measurements is calculated as where ∆ s p is the mean value for proton showers and ∆ s Fe the mean value for iron showers.The mean primary mass for two different hadronic interaction models and a comparison to X max measurements are shown in Fig. 5.The trend in ln A and the break in the estimated and QGSJetII-04 (right) [5].ln A determined from FD measurements (gray) are added for comparison [10].
composition are comparable to those obtained from FD measurements, but the values indicate a heavier composition over the complete energy range.This difference is caused by models not describing the muon component well enough.Calculating Pearson's correlation of ∆ s and X max enables the transformation of ∆ s into a calibrated value of X max extracted from SD data only.X max calculated in this way is transformed into the average logarithmic mass and shown in Fig. 6.Comparison shows a good agreement between SD and FD measurements, with SD measurements being statistically stronger (almost double the number events).The and QGSJetII-04 (right).Results from the Delta method have been calibrated to X max [5] and compared to FD measurements (gray) [10].
highest energy bins have been fragmented further for SD and show a possible reduction of the trend towards heavier composition.
As an additional analysis approach, the distributions of X max have been parameterized as described in [11] and fit for fractions of p, He, N and Fe in each energy bin.The final fits are shown in Fig. 7 for three different high-energy hadronic interaction models.Fractions The p-value serves as a goodness of fit estimator [6].

Figure 1 .
Figure 1.Mean (left) and standard deviation (right) of the distribution of X max as a function of energy.For comparison, simulations for proton and iron primaries using three different hadronic interaction models are added[6].

Figure 2 .
Figure 2. Comparison of X max data as a function of energy for the TA and Pierre Auger Observatory folded with the same reconstruction and detector effects [8].

Figure 5 .
Figure 5. ln A as a function of energy for two different hadronic interaction models: EPOS-LHC (left) and QGSJetII-04 (right)[5].ln A determined from FD measurements (gray) are added for comparison[10].

Figure 6 .
Figure 6.ln A as a function of energy for two different hadronic interaction models: EPOS-LHC (left) and QGSJetII-04 (right).Results from the Delta method have been calibrated to X max[5] and compared to FD measurements (gray)[10].

Figure 7 .
Figure 7. Mass fractions for p, He, N and Fe applied to parameterized X max distributions.The thick parts of the error bars show statistical uncertainties, while systematics are shown with thinner error bars.The p-value serves as a goodness of fit estimator[6].