Measurement of the fluctuations in the number of muons in extensive air showers with the Pierre Auger Observatory

We present the first measurement of the fluctuations in the number of muons in extensive air showers produced by ultra-high energy cosmic rays. We find that the measured fluctuations are in good agreement with predictions from air shower simulations. This observation provides new insights into the origin of the previously reported deficit of muons in air shower simulations and constrains models of hadronic interactions at ultra-high energies. Our measurement is compatible with the muon deficit originating from small deviations in the predictions from hadronic interaction models of particle production that accumulate as the showers develop.

We present the first measurement of the fluctuations in the number of muons in extensive air showers produced by ultra-high energy cosmic rays.We find that the measured fluctuations are in good agreement with predictions from air shower simulations.This observation provides new insights into the origin of the previously reported deficit of muons in air shower simulations and constrains models of hadronic interactions at ultra-high energies.Our measurement is compatible with the muon deficit originating from small deviations in the predictions from hadronic interaction models of particle production that accumulate as the showers develop.

INTRODUCTION
Ultra High Energy Cosmic Rays (UHECRs) are particles coming from outer space, with energies exceeding 10 18 eV.They provide the only experimental opportunity to explore particle physics beyond energies reachable by Earth-based accelerators, which go up to cosmic ray energies of 9 × 10 16 eV.
The Pierre Auger Observatory [1] detects extensive air showers that are initiated by the UHECRs colliding with the nuclei in the atmosphere.Information about UHE-CRs is extracted using simulations based on hadronic interaction models which rely on extrapolations of accelerator measurements to unexplored regions of phase space, most notably the forward and highest-energy region.In addition, accelerator experiments at the highest energies either probe the interactions between protons or of protons with heavy nuclei, while most interactions within air showers are between pions and light nuclei.
A further challenge is that the UHECR mass has to be measured despite not being yet completely decoupled from the hadronic uncertainties.The observable with the least dependence on hadronic interactions is the atmospheric depth at which the longitudinal development of the electromagnetic (EM) component of the shower reaches the maximum number of particles, namely X max [2].
In hadronic cascades the energy of each interacting particle is distributed among the secondaries, mostly pions.Neutral pions rapidly decay into two photons feeding a practically decoupled electromagnetic cascade (other resonances decaying into π 0 's, electrons and or photons also contribute).Charged pions (and other long lived mesons like kaons) tend to further interact until their individual energies are below a critical value, below which they are more likely to decay.Muons, which are products of hadronic decays, are thus predominantly produced in the final shower stages.In sufficiently inclined showers, the pure EM component is absorbed in the atmosphere and the particles that reach the ground (muons and muon decay products) directly sample the muon content [3,4], reflecting the hadronic component of the shower.
Air showers are mainly detected at the Pierre Auger Observatory by the Surface Detector (SD), an array of water-Cherenkov detector stations, and the Fluorescence Detector (FD) consisting of 24 fluorescence telescopes.By selecting the sub-sample of events reconstructed with both the SD and the FD and with zenith angles exceeding 62 • , both the muon content and the energy of the shower are simultaneously measured.
The results obtained indicate that all the simulations underestimate the number of muons in the showers [5,6].These analyses come with the caveat that they cannot distinguish a muon rescaling from a shift in the absolute energy scale of the FD measurement.However, muon content and energy scale were disentangled in a complementary technique based on showers with zenith angles below 60 • .Using the longitudinal profile of the shower in the atmosphere obtained with the FD and the signals at the ground measured with the SD, it was shown that the muonic component still has to be scaled up to match observed data, while no rescaling of the EM component and the FD energy is required [7].The measurements with the FD also show that both the position of the shower maximum in the atmosphere (X max ) and the entire shape of the EM shower are well described by the simulations [8,9].At lower energies, down to ∼ 10 17. 3 eV, in a measurement using the sub-array of buried scintillators of the Pierre Auger Observatory, a direct count of the muons independent of EM contamination was obtained, which also shows that simulations produce too few muons [10].There is much evidence that all the simulations underpredict the average number of muons in the showers: a comprehensive study of muon number measurements made with different experiments has shown that the muon deficit in simulations starts around ∼ 10 16 eV and steadily increases with energy.Depending on model and experiment, the deficit at ∼ 10 20 eV ranges between tens of percent up to a factor of two [11].
The increased statistics obtained at the Pierre Auger Observatory allows us to now take a further step and explore fluctuations in the number of muons between showers, hereinafter referred to as physical fluctuations.The ratio of the physical fluctuations to the average number of muons (relative fluctuations) has been shown to be mostly dominated by the first interaction, rather than the lower energy interactions deeper in the shower development [12,13].Here, we exploit the sensitivity of fluctuations to the first-interaction to explore hadronic interactions well above the energies achievable in accelerator experiments.

METHODOLOGY
Our analysis here is based on the set of inclined air showers (62 • < θ < 80 • ) that are reconstructed both with the SD and with the FD between 1 January 2004 and 31 December 2017.For each event we obtain independent measurements of the muon content (with the SD) and the calorimetric energy (with the FD).To ensure the showers can be reconstructed with small uncertainties, we select only events with at least four triggered stations in the SD array and we further require that all the stations surrounding the impact point of the shower on the ground are operational at the time of the event.Only events with good atmospheric conditions, (few clouds and a low aerosol content) are accepted in order to guarantee a good energy reconstruction with the FD.In addition it is required that the entire shower profile and in particular X max is within the field of view of our telescopes.Since heavy primaries penetrate the atmosphere less than light ones, the acceptance with this selection would be mass dependent.To avoid this bias we constrain the field of view to the region where all values of X max are accepted.Further details are given in [5,14].These selection criteria result in a total number of events of 786.
In addition only events with energy larger than 4 × 10 18 eV, which ensures full trigger efficiency of the SD [3], are used to extract the fluctuations (281 events).
The number of muons is reconstructed by fitting a 2-D model of the lateral profile of the muon density at the ground to the observed signals in the SD array.The free parameters of the fit are the zenith and azimuth angles of the shower, the impact point of the shower on the ground (shower core position) and a normalization factor with respect to a reference muon density profile in simulated proton showers at 10 19 eV [3].There exists a residual pure EM component in showers with low zenith angles and stations very close to the shower core position (at 400 m and 64 • it is ∼ 6%) which has been subtracted using a parameterisation [4].The dimensionless normalization factor we obtain from the fit is then transformed to the dimensionless quantity, R µ , which is given by the integrated number of muons at the ground divided by a reference given by the average number of muons in simulated proton showers at 10 19 eV and the given zenith angle.At 10 19 eV and an inclination of 60 • , R µ = 1 corresponds to 2.148 × 10 7 muons.For more details see [5].
In the following we refer to R µ as the number of muons, for short.
The calorimetric energy of the air showers, E cal , is reconstructed by integrating the longitudinal shower profiles observed with the FD [9,15].The total energy of the shower is then obtained by adding the average energy carried away by muons and neutrinos, the so-called invisible energy, E = E cal + E inv .At 10 19 eV E inv accounts for 14% of the total energy in air showers [16][17][18][19][20].
In Fig. 1 the muon number, R µ , is shown as a function of the measured energy.Markers on the top of the frame define the bins in energy for which we will extract the fluctuations, with the numbers of events in each bin shown above.The bins are chosen such that the number of events in each is similar.Based on models of air shower development and given the gradual change of the composition in this energy range (single logarithmic dependence on energy) [8,[21][22][23], the number of muons is related to the primary energy by a single power-law which can be fitted following a procedure described in the text below.The best-fit parameters are given at the beginning of the next section.The scattering in the data has three sources: experimental uncertainties in the energy s E and in the muon number s µ from event reconstruction (both represented by the error bars), and the physical fluctuations in the muon number denoted as σ.Given Eq. ( 1) the variance of the muon number is In this work, we adopt a method based on maximizing the likelihood of a probability distribution function (PDF).The PDF incorporates the various contributions to the fluctuations, treating each energy bin independently while also accounting for the effect of the migration of events between bins [5,24].
The model assumes that measurements of E and R µ follow Gaussian distributions centered at the true value, with widths given by the detector resolution s E and s µ , which are the uncertainties obtained in each individual event reconstruction [3,25].
Physical fluctuations are also assumed to follow a Gaussian distribution of width σ.Simulations have shown this is an acceptable approximation given the event number in each bin.
The total PDF is obtained through the convolution of the detector response and the physical fluctuations with the probability distribution of the hybrid events measured at the Pierre Auger Observatory.The loglikelihood function is then given by ln L(a, b, σ1 , .., σ6 ) = . ( The probability of hybrid events, h(E), (product of the energy spectrum of CRs and the efficiency of detection), can be obtained from the data itself as explained in [24] and [10,26].The RHS of Eq. ( 2) depends on the parameters a, b via Eq.(1).To obtain the energy dependence of the fluctuations we parameterize σ by six independent values such that σ(E) = σk • R µ (E) where the constants σk are the relative fluctuations in the k-th energy bin with limits where k runs from one to six.In Eq. ( 2), k = 0 corresponds to the contributions from the interval [0, E thr ] where the SD is not fully efficient.
The fluctuations here are assumed to take the value of the first fitted bin, σ0 ≡ σ1 .The sum over the index i in Eq. ( 2) (the usual sum over the log-likelihoods of events), includes only events above the energy threshold of 4 ×10 18 eV.The function C(E) is the normalization factor from the double Gaussian.The result of the fit for the parameters a and b are shown in Fig. 1.The fluctuations are shown in Fig. 2. The distribution of the number of muons and the PDF in the individual energy bins can be found in the Supplemental Material [27].
The dominant systematic uncertainties of σ come from the uncertainties in the resolutions s E and s µ .For s µ we estimate the uncertainty using simulations and data.In simulations the uncertainty was estimated by the spread in a sample of simulated showers where each shower is reconstructed multiple times, each time changing only the impact point at the ground.For data we reconstruct the same event multiple times, leaving out the signals from one of the detector stations.The average relative resolution, s µ /R µ , and its systematic uncertainty is thus (10 ± 3) % at 10 19 eV.
We verified the values of s E by studying the difference in the energy reconstruction of events measured independently by two or more FD stations.The width of the distribution of these energy-differences is found to be compatible with s E .We therefore take the statistical one-σ uncertainties of this cross check as a conservative upper limit of the systematic uncertainty of s E [28].The average relative energy resolution, s E /E , is about (8.4 ± 2.9) % at 10 19 eV.
We have further confirmed that there are no significant contributions to the fluctuations from differences between Any residual electromagnetic component in the signal would affect the lower zenith angles more.We therefore split the event sample at the median zenith angle (66 • ) and compare the resulting fluctuations.We find no significant difference between the more and the less inclined sample.
In another test we do find a small modulation of R µ with the azimuth angle (< 1%), which we correct for.This modulation is related to the approximations used in the reconstruction, which deal with the azimuthal asymmetry of the muon densities at the ground due to the Earth's magnetic field [3].Finally we have run an endto-end validation of the whole analysis method described in this article on samples of simulated proton, helium, oxygen and iron showers.
Due to the almost linear relation between R µ and E, the systematic uncertainty on σ due to the uncertainty of the absolute energy scale of 14 % [25] practically cancels out in the relative fluctuations.The systematic uncertainty in the absolute scale of R µ of 11 % [5] drops out for the same reason.The systematic effects for the bin around 10 19 eV are summarized in Tab.I.Over all energies the systematic uncertainties are below 8%.
The measured relative fluctuations as a function of the energy are shown in Fig. 2. We note that the measurement falls within the range that is expected from current hadronic interaction models for pure proton and pure iron primaries [30][31][32][33][34][35][36][37][38].To estimate the effect of a mixed composition, we take the fractions of the four mass components (proton, helium, nitrogen and iron) derived from the X max measurements [8,39,40] and, using Measured relative fluctuations in the number of muons as a function of the energy and the predictions from three interaction models for proton (red) and iron (blue) showers.The gray band represents the expectations from the measured mass composition interpreted with the interaction models.The statistical uncertainty in the measurement is represented by the error bars.The total systematic uncertainty is indicated by the square brackets.
the simulations of the pure primaries, calculate the corresponding fluctuations in the number of muons.The gray band in Fig. 2 encompasses the predicted σ/ R µ of the three interaction models Qgsjet II-04, Epos-lhc and Sibyll 2.3d given the inferred composition mix for each [41].
In Fig. 3 the effects of different composition scenarios on both the fluctuations and the average number of muons can be shown by drawing, at a fixed primary energy of 10 19 eV, the relative fluctuations σ/ R µ against the average number of muons R µ .Given any one of the interaction models, any particular mixture of the four components p, He, N and Fe falls somewhere within one of the areas enclosed by the corresponding colored lines.The points of pure composition in this contour are labeled accordingly.For each model the expected values for σ/ R µ and R µ given the composition mixture obtained from the X max measurements [8] is indicated within each contour by the correspondingly colored star marker.The shaded areas surrounding the star markers indicate the statistical and systematic uncertainties inherited from the X max measurements [42].Finally our measurement with statistical and systematic uncertainty is shown by the black marker.
Within the uncertainty none of the predictions from the interaction models and the X max -composition (star markers) are consistent with our measurement.The predictions from the interaction models Qgsjet II-04, Epos-lhc and Sibyll 2.3d can be reconciled with our measurement by an increase in the average number of Data (black, with error bars) compared to models for the fluctuations and the average number of muons for showers with a primary energy of 10 19 eV.Fluctuations are evaluated in the energy range from 10 18.97 eV and 10 19.15 eV.The statistical uncertainty is represented by the error bars.The total systematic uncertainty is indicated by the square brackets.The expectation from the interaction models for any mixture of the four components p, He, N, Fe is illustrated by the colored contours.The values preferred by the mixture derived from the Xmax measurements are indicated by the star symbols.The shaded areas show the regions allowed by the statistical and systematic uncertainties of the Xmax measurement [42].muons of 43%, 35%, and 26%, respectively.For the fluctuations, no rescaling is necessary for any model.
Taken together, the average value and fluctuations of the muon flux constrain the way hadronic interaction models should be changed to agree with air shower data.
To see this we briefly discuss the origin of the fluctuations.
The average number of muons in a proton shower of energy E has been shown in simulations to scale as N * µ = C E β where β 0.9 [12,13,22,23].
If we assume all the secondaries from the first interaction produce muons following the same relation as given for protons above, we obtain the number of muons in the shower as where index j runs over m secondary particles which reinteract hadronically and x j = E j /E is the fraction of energy fed to the hadronic shower by each [43].In this expression the fluctuations in N µ are induced by α 1 in the first generation which fluctuates because the multiplicity m and the energies x j of the secondaries fluctuate [13].We can continue this reasoning for the subsequent generations to obtain here the subindex i runs over n generations, until the cascade stops.We note that for the calculation of α 2 , in the second generation, there are m particles contributing.Assuming the distributions of the α's for each one are similar, when adding up the muons produced by each, the fluctuations produced by one are statistically likely to be compensated by another.In other words the α 2 distribution is narrower by a factor ∼ 1/ √ m.The deeper the generation, the sharper the corresponding α i is expected to be.As a result, the dominant part of the fluctuations comes from the first interaction.This has also been observed with simulations.The model can be generalised for primary nuclei with mass A using the superposition model, and fixing the number of participants to A protons, which reduce the different contributions to the fluctuations by a factor ∼ 1/ √ A.
There are two options to increase the average number of muons in air showers.One is to increase α in a specific generation, notably the first where the energy is the highest and exotic phenomena could conceivably play a role, i.e. α 1 → α 1 + δα 1 .Note that if only the first generation is modified (implying some sort of threshold effect for new physics) the increase in N µ is linear with the modification.There are several examples in the literature where this approach has been used assuming different mechanisms [44][45][46][47][48].For the fluctuations the change depends on the model.Alternatively, the number of muons can be increased by introducing small deviations in the hadronic energy fraction, δα, in all generations.Accumulated along a number n of generations, these small deviations build up as N µ ∝ (α + δα) n .For instance, a 5% deviation per generation converts into ∼ 30% deviation after six generations [49].On the other hand, a change of 5% in the fluctuations of α is not amplified in the muon fluctuations because of the suppression in later generations.This approach characterises the increase in the number of muons in the current hadronic interaction models with regard to previous models.It is also compatible with the increase of the discrepancy in the average number of muons across a wide range of energies reported in [11].
The present analysis finding that fluctuations are consistent with model predictions means that the increase in muon number may be a small effect accumulating over many generations, or a very particular modification of the first interaction that changes N µ without changing the fluctuations [50].

SUMMARY
We have presented for the first time a measurement of the fluctuations in the number of muons in inclined air showers, as a function of the UHECR primary energy.Within the current uncertainties, the relative fluctuations show no discrepancy with respect to the expectation from current high-energy hadronic interaction models and the composition taken from X max measurements.This agreement between models and data for the fluctuations, combined with the significant deficit in the predicted total number of muons, points to the origin of the models' muon deficit being a small deficit at every stage of the shower which accumulates along the shower development, rather than a discrepancy in the first interaction.Adjustments to models to address the current muon deficit, must therefore not alter the predicted relative fluctuations.
The Pierre Auger Observatory is currently undergoing an upgrade which includes the deployment of scintillators on top of the SD stations [51] to help disentangle the muonic and electromagnetic content of the showers, as well as an array of radio antennas [52].It has been shown that radio arrays can provide an estimate of the calorimetric energy [53] and therefore it will soon be possible to perform an analysis similar to the one presented here with much larger statistics using hybrid events measured by the high-duty-cycle radio and surface detector arrays [52].

σ
FIG. 1. Number of muons as a function of the measured energy.The black line is the fitted Rµ = a (E/(10 19 eV)) b .Markers on the top of the frame define the bins in which the fluctuations are evaluated.The numbers give the events in each bin.The effect of the uncertainty of the absolute energy scale is indicated by σsys(E).
FIG. 2.Measured relative fluctuations in the number of muons as a function of the energy and the predictions from three interaction models for proton (red) and iron (blue) showers.The gray band represents the expectations from the measured mass composition interpreted with the interaction models.The statistical uncertainty in the measurement is represented by the error bars.The total systematic uncertainty is indicated by the square brackets.

TABLE I .
Contributions to the systematic uncertainty in the relative fluctuations around 10 19 eV (10 18.97 eV to 10 19.15 eV).