Recent results from the Pierre Auger Observatory on the mass composition and hadronic interactions of ultra-high energy cosmic rays

Recent results on the mass composition and hadronic interactions of ultra-high energy cosmic rays (UHECR) obtained at the Pierre Auger Observatory are reviewed. Different studies based on shower observables recorded with the fluorescence and surface detectors of the Observatory indicate that the mass composition is getting lighter for the energies below ≈ 2 EeV and heavier afterwards. A fully consistent description of the data collected at the Observatory cannot be achieved with any of the up-to-date interaction models.


Introduction
A large variety of data, with unprecedented statistics, has been collected at the Pierre Auger Observatory [1] since the start of the operation in 2004.Analysis of the UHECR mass composition can be performed using different characteristics of the longitudinal and lateral profiles of electromagnetic (EM) and muon shower components with the methods having independent systematic uncertainties.The approaches relying on the studies of the EM component (analysis of the depth of the shower maximum X max ) or using general properties of the extensive air showers (analysis of the correlation between X max and the signal in surface stations or studies of the elongation rates) have smaller uncertainties due to hadronic interaction models compared to the approaches making use of absolute values related to the muon component.In this proceedings the results on the mass composition from the first and second types of the methods are compared and the conclusions on both the evolution of the primary mass with energy and on the drawbacks of the interaction models (EPOS-LHC [2], QGSJetII-04 [3], Sibyll 2.1 [4]) are presented.
The data reviewed in the paper are recorded using the fluorescence (FD) and surface detectors (SD) of the Pierre Auger Observatory, located in the Mendoza province, Argentina.The SD consists of 1660 water-Cherenkov detectors (WCD) placed in a triangular grid with a spacing of 1500 meters over the area of 3000 km 2 .Twenty-seven telescopes of the FD are located at four sites on the perimeter of the SD.The FD operates only during clear moonless nights and its duty cycle is around 15%, while the duty cycle of the SD is ≈100%.Showers triggering both FD and SD are called hybrid events.For them, complementary information from both detectors allows one to reconstruct the shower parameters with higher a e-mail: yushkov.alexey@gmail.comprecision and to perform cross-checks of the operation of the detectors.Three of the analyses presented belowmeasurements of the depth of shower maximum X max , of the correlation between X max and signal in WCDs, and of the muon content in highly inclined events -exploit the advantages of the hybrid design of the Observatory.The measurements of risetime asymmetry and of muon production depth are performed using information from WCDs on particle arriving times sampled at the rate of 40 MHz.
Unless otherwise specified the data used are from the period from January 2004 till December 2012.

Measurement of the depth of shower maximum X max [5, 6]
Direct determination of X max and the calorimetric measurement of the energy of the primary particle are based on the detection of fluorescence light emitted by the molecules of atmospheric nitrogen due to their interactions with the shower particles.The parameters of the atmosphere, acting as a calorimeter, are thoroughly monitored at the Observatory and only events passing a pre-selection on good atmospheric conditions (cloud coverage, aerosol density) are accepted for the analysis.These events are subject to further selection ensuring unbiased measurement of the X max distributions including specific requirements on the field of view of the FD telescopes and the quality cuts on the profile.For full details on the selection and reconstruction procedures the reader is referred to [5,6].
In June 2010 three High Elevation Auger Telescopes (HEAT) started to take data in addition to the standard 24 telescopes of the FD.The HEAT in combination with the standard telescopes at the Coihueco site (called further HeCo) have a field of view 2   way nearby showers with energies down to 10 17 eV can be observed.To make the data set recorded by HeCo complementary and independent, events (1377 in total) with energies < 10 18.3 eV observed by HeCo during the period 01/06/2010 -15/08/2012 are excluded from the data of the standard FD.
The final data set contains 5490 HeCo events with energies 10 17.0 −10 18.3 eV and 18382 events from the standard FD with energies > 10 17.8 eV.The X max resolution is around 30 g/cm 2 at 10 17.0 eV and improves with the energy to 15 g/cm 2 for E > 10 18.5 eV.The systematic uncertainty is below (−15; + 10) g/cm 2 for most of the energies.
The comparison of the first two central moments of the measured X max distributions with predictions of different hadronic interaction models for protons and iron nuclei is given in Fig. 1.The observed elongation rate (ER) is described well with a broken line, the energy of the break is lg(E 0 /eV) ≈ 18.3 (E 0 ≈ 2 EeV) and the slopes below and above the break are equal to ≈ 85 and ≈ 26 [g/cm 2 / decade] correspondingly.The ER in simulations for all hadronic interaction models and all constant compositions lies in the range of 54−64 [g/cm 2 /decade], thus the observed ERs indicate that the primary mass is decreasing below the break around 2 EeV and is increasing for the higher energies.Using the superposition model (see e.g.[7,8]) the mean logarithmic mass can be found from the measured X max using the relation X max ≈ X p max − ln A D 10 / ln (10), where X p max is the value for protons.The evolution of ln A with energy is shown in Fig. 2 (top row).The observed transition from lighter to heavier ISVHECRI 2016 composition implies that the primary beam is not pure and represents a mix of different nuclei.The degree of mixing (the variance of the logarithmic masses σ 2 ln A ) can be estimated using the measured width of the X max distributions [9]: σ 2 (X max ) ≈ σ 2 sh + σ 2 ln A (D 10 / ln(10)) 2 , where σ sh is the shower-to-shower X max fluctuations for a fixed mass.The results on σ 2 ln A , presented in Fig. 2 (bottom row), show different degrees of mixing in the primary beam depending on the hadronic interaction model used.Though for QGSJetII-04the variance of the masses takes negative values for some of the energies, the discrepancy is not larger than two standard statistical deviations and is non-significant if one takes into account the systematic errors on X max measurements and on the FD energy scale.

Measurement of the correlation between X max and S(1000) [10]
The spread of the masses in the primary beam σ (ln A) = ln 2 A − ln A 2 can be estimated from the measurement of the correlation between X max and the signal in WCDs at 1000 m from the shower axis S(1000) [10,11].In simulations this correlation is close to zero or is slightly positive for pure compositions.Mixed compositions show a negative correlation: heavier primary nuclei have shallower X max ( X max ∼ − ln A) and larger S(1000)(due to a larger number of muons N µ ∼ A 1−β , β 0.9 [12]), thus in a sample with a mixed composition the shallower showers on average will have larger signals.Such behavior is a general property of extensive air showers and thus the results of the correlation analysis are weakly sensitive to the choice of a particular interaction model.
In this analysis, to avoid a decorrelation due to the spreads of energies and zenith angles, X max and S(1000) are scaled to a reference energy of 10 EeV and S(1000) is additionally scaled to zenith angle of 38 • .The scaled variables are denoted further as X * max and S * 38 .The correlation between X * max and S * 38 is evaluated using a ranking correlation coefficient r G proposed in [13].The conclusions do not change when other correlation coefficients are used.
To have a data set with an unbiased mass composition the event selection is the same as in the X max analysis [6].In addition, to have a reliable reconstruction of S(1000), it is required that the SD station with the highest signal should have at least 5 active neighbour stations.For the energy range lg(E/eV) = 18.5−19.0and the zenith angle range 0 • −65 • , the final data set contains n = 1376 events.The statistical uncertainty can be approximated [10] by r G 0.9/ √ n, thus r G (data) = 0.024.The systematic error on r G is small 0.01: invariance of the ranking correlation coefficients to any transformations leaving ranks of the events unchanged makes the correlation analysis almost insensitive to the systematic uncertainties in X max and S(1000).
In Fig. 3 the distributions of X * max and S * 38 in data and in simulations for protons and iron with EPOS-LHCare shown together with the corresponding values of the correlation coefficients.The correlation coefficients for protons for all models are non-negative: r G (EPOS-LHC) = 0.00, r G (QGSJetII-04) = 0.08 and r G (Sibyll 2.1) = 0.06.For heavier nuclei the correlation coefficients are larger than those for protons and for all p − He mixtures r G 0 as well.Thus, the negative correlation found in data r G (X * max , S * 38 ) = −0.125 ± 0.024 (stat) cannot be reproduced using any pure composition and the data can be explained only by mixed compositions containing primary nuclei heavier than helium A > 4. The conclusions do not change when the key parameters of the interaction models (cross-section, multiplicity, elasticity, pion charge ratio) are modified by a factor ≈ 1.5.
The spread of the primary masses σ (ln A) can be estimated from Fig. 4 where the correlation value found in data is compared to the values in simulated mixtures with different fractions of protons, helium, oxygen and iron nuclei.The correlation r G (X * max , S * 38 )gets more negative for the larger spreads of the masses.For all three models (results for Sibyll 2.1, not shown in Fig. 4, are almost the same as for QGSJetII-04) the spread of masses corresponding to the correlation found in data, lies in the range 1.0 σ (ln A) 1.7.

Azimuthal asymmetry in the risetime of signals in WCDs [14]
Muons produce short and spiky time traces in WCDs, while for EM particles the distribution of the arrival times is larger and smoother.Since the relative contribution of muons to the SD signal increases with the growth of the mass of the primary nuclei, the risetime, t 1/2 , defined as the time of the increase of the total integrated signal in WCDs from 10% to 50%, is a mass-sensitive observable.Due to changes in the muon/EM signal ratio and in the thickness of the shower front with distance to the shower axis r and zenith angle θ , risetime is a function of r and θ .Risetime grows almost linearly t 1/2 (r ) ∝ r for a large range of r , thus t 1/2 /r is used to cancel out the distance dependence and to include WCDs from all distances in the analysis.
In inclined showers signal risetime in individual WCDs depends on their azimuth angle ζ in the shower plane: particles hitting 'early' stations (|ζ | < π/2) cross smaller atmosphere thickness and have smaller angle to the shower axis compared to the particles hitting 'late' stations (|ζ | > π/2).With the growth of zenith angle fewer EM particles reach late detectors due to absorption by the additional atmosphere (the signals thus become more muon-rich) and risetime in the late detectors becomes shorter compared to that in early detectors.As shown in Fig. 5 (top) the asymmetry in risetimes reaches its maximum at a certain zenith angle and decreases for larger angles where the muon component starts to give a dominant contribution to the signals.
The dependence of the risetime on azimuth angle can be parameterized as t  The behavior of the asymmetry parameter with ln(sec(θ )) at lg(E/eV) = 19.1 in stations within the range r = 500 − 1000 m is shown in Fig. 5 (bottom) together with a fit with a Gaussian function.The value of sec(θ ) max at which the asymmetry parameter reaches its maximum provides information on the mass of the primary particle.The behavior of sec(θ ) max in the SD data is studied using events with the energies corresponding to 100% SD efficiency E > 3 × 10 18 eV for the zenith angle range of θ = 30 • − 62 • .For an accurate reconstruction of t 1/2 only the stations with the 100% trigger probability (signal > 10 VEM) are used.The stations should not have saturated traces and the reliable determination of t 1/2 should be possible, these requirements limit the range of the WCDs distances to r = 500−2000 m.For the data recorded from January 2004 to October 2014 a total of 191534 signal traces from 54584 events are selected.The evolution of sec(θ ) max with energy is analyzed in two ranges of the ISVHECRI 2016 distances to the shower axis r = 500−1000 m and r = 1000−2000 m.Systematic error on sec(θ ) max is ≈ 0.014.
The comparison of sec(θ ) max in the data and simulations for protons and iron nuclei is shown in Fig. 6.The values for data are bracketed by the QGSJetII-04and EPOS-LHCpredictions for protons and iron.There are indications of an increase of the primary mass with energy, in agreement with other results reported in this proceedings, but uncertainties from the interaction models on the absolute value of sec(θ ) max do not allow one to make robust conclusions of average mass of the primary nuclei.For QGSJetII-04such conclusions have an unphysical dependence on the range of the distances to the shower axis favouring lighter composition at larger distances r = 1000−2000 m.
From Fig. 7, where the values of sec(θ ) max converted to ln A are presented, one can see that the average mass from the risetime asymmetry analysis lies between the values obtained from the measurements of X max and the muon production depth (MPD, discussed in the next section).This might be due to the different role played by the muon and EM components in each of the analyses: for X max determination the EM component is of major importance, MPD analysis relies on the reconstruction of properties of muons, and the sec(θ ) max behavior results from a complicated interplay between both muon and EM components.

Measurements of the muon production depth X µ max [15]
Muons in the atmosphere travel along almost straight lines and the production height z of each muon along the shower axis can be reconstructed from a quite simple geometry shown in Fig. 8 [15]: z ≈ 1/2 r 2 /c(t − t ε ) − c(t − t ε ) + − z π , here t εmean kinematic delay; z π -pion decay length.After conversion of heights z to the mass units X µ = ∞ z ρ(z )dz one gets a distribution of the MPD for an event and the depth of the maximum muon production X µ max is determined from the fit with a Gaisser-Hillas function as shown in Fig. 8.To reduce the EM contamination of the traces in the WCDs only inclined showers with zenith angles θ = [55 • ; 65 • ] are used in the current analysis.Additionally, a cut on the primary energies E > 20 EeV is applied to have sufficient number of muons for the reconstruction of the profile, and a cut on the station distances to the core 1700 m < r < 4000 m ensuring at the same time similar proton and iron selection efficiencies and good X µ max resolution.The final data set contains 481 events, the systematic uncertainty on X µ max is around 17 g/cm 2 , the resolution is 100 (80) g/cm 2 at 10 19.3 eV for p (Fe) and 50 g/cm 2 at 10 20.0 eV.
The evolution of the MPD in data and the predictions for protons and iron, obtained with EPOS-LHCand QGSJetII-04, are shown in Fig. 9. Comparison of the measured ER D µ 10 = d X µ max /d lg(E/eV) = −25 ± 22 (stat.)± 21 (sys.)[g/cm 2 /decade] (Fig. 9) with the values for protons (35.9 ± 1.2 [g/cm 2 /decade]) and iron (48.0 ± 1.2 [g/cm 2 /decade]) indicates the transition from lighter to heavier composition with increasing energy.Due to large uncertainties on the absolute values of X µ max the reconstructed ln A strongly depends on the interaction model.For EPOS-LHCthe mass composition is heavier than iron (cf.Figs. 7, 9).The deeper X µ max for EPOS-LHCis due to a larger elasticity in protonair and pion-air interactions, which results from a better, compared to QGSJetII-04, description of rapidity gap distributions measured at the LHC [16].Thus the better agreement of ln A from X µ max and X max measurements for QGSJetII-04 (Fig. 7) might be rather accidental.

Muon content of highly inclined showers [17]
Using highly inclined events (zenith angles > 62 • ) an almost direct measurement of the muon shower content can be performed.The approach consists in production of the reference muon density lateral profiles ρ µ,19 ( − → r ; θ, φ) with a shape weakly depending on energy, on the type of the primary nuclei, on the hadronic interaction model and on the software used for shower simulations.For a given muon density ρ µ , a muon scale factor N 19 with respect to ρ µ,19 ( − → r ; θ, φ) is then introduced: ρ µ = N 19 ρ µ,19 ( − → r ; θ, φ).By construction N 19 does not depend on the zenith angle.Finally one compares scaling factors in data and in Monte-Carlo (MC) simulations R µ = N data 19 /N MC 19 .The production of the reference profiles ρ µ,19 ( − → r ; θ, φ) for different zenith θ and azimuth φ angles is done using QGSJetII-03 protons at 10 19 eV.Hybrid events with energies above 10 18.6 eV (≈ 4 EeV), corresponding to 100% trigger efficiency of the SD, are accepted for the analysis, the zenith angles are limited to the [62 • ; 80 • ] range.In total 174 events remain after the selection.
The muon scales R µ in data and simulations for protons and iron are shown in Fig. 10.Similarly to the measurements of X max , risetime asymmetry and MPD, the slope d R µ /d lg(E/eV) in data differs from the MC predictions for constant compositions and is compatible with a transition to a heavier mass composition.The absolute values of R µ in data correspond to masses close to or heavier than iron, thus being in contradiction with the results on ln A from the analysis of X max .Similar conclusions were obtained on the muon content of the vertical showers [18].To quantify the discrepancy, the mean logarithmic mass from the analysis of X max distributions can be converted into the mean logarithmic muon scale ln R µ [17] thus leading to a one-to-one relation between ln R µ and X max .The results presented in Fig. 11 show, that the muon content ln R µ in data is (30−80)% larger compared to simulations for primary masses ln A corresponding to the measured X .max Nevertheless, the minimum difference between data and EPOS-LHCis only 1.4σ sys .

Discussion
The results of studies of the UHECR mass composition at the Pierre Auger Observatory, presented in these proceedings, indicate that ln A is not constant for all energies above 10 17 eV.Specifically, the ER of X max (d X max /d lg(E/eV)) shows that the composition is getting lighter below 2 EeV and then heavier for higher energies.Conversion of the moments of X max distributions to the moments of the logarithmic mass produces results that are consistent with a mixed composition within  systematic uncertainties when EPOS-LHCis used.For QGSJetII-04the conclusions are less certain and even an unphysical variance of logarithmic masses (σ 2 ln A < 0) is observed for E 3 EeV, though the values differ from zero by less than 2σ stat .
Analysis of the correlation between X max and S(1000) provides a strong evidence on the presence of several types of nuclei near the ankle (≈ 5 EeV) in the primary spectrum with the spread of the masses σ (ln A)>1 for all interaction models used (EPOS-LHC, QGSJetII-04, Sibyll 2.1).The data are described well by mixed  compositions containing nuclei heavier than helium A > 4. The results are robust with respect to modifications of the key parameters of the interaction models.In [19] the X max distributions measured at the Pierre Auger Observatory were interpreted in terms of fractions of protons, helium, nitrogen and iron.For QGSJetII-04and Sibyll 2.1the composition was found to be close to 0.5 p − 0.5 He (σ (ln A) ≈ 0.69) for energies lg(E/eV) = 18.5−19.0.The correlation coefficient for 0.5 p − 0.5 He mix for these models is r G ≈ 0.08 and is in contradiction with the value found in data, indicating thus that QGSJetII-04and Sibyll 2.1do not describe X max distributions and (X max , S(1000)) correlation in a self-consistent way.For EPOS-LHCthe fits of X max distributions produce a mixed composition including protons, helium and nitrogen with σ (ln A) ≈ 1.2 and r G ≈ −0.094 in a better agreement with the correlation analysis.
The ERs for the azimuthal asymmetry in the risetime and for the MPD differ from the ERs for constant

Figure 1 .
Figure 1.The mean and the standard deviation of the measured X max distributions compared to the simulations with EPOS-LHC, QGSJetII-04and Sibyll 2.1for protons and iron.

Figure 2 .
Figure 2. The evolution of the mean (top) and of the variance (bottom) of ln A determined from the measurements of X max .

Figure 8 .Figure 9 .
Figure 8. Top: the geometry used for the estimation of the muon production height z.Bottom: the muon production depth (MPD) profile of a real event with the energy of 33 ± 1 EeV fitted with the Gaisser-Hillas function (red line).

Figure 10 .
Figure 10.The muon scales R µ in data and in the MC simulations with EPOS-LHCand QGSJetII-04for protons and iron.

Figure 11 .
Figure 11.Comparison of ln R µ at 10 19 eV in data to the values in simulations for protons, iron and for the mixes with ln A from the measurements of X max .