Study of cosmic ray events with high muon multiplicity using the ALICE detector at the CERN Large Hadron Collider

ALICE is one of four large experiments at the CERN Large Hadron Collider near Geneva, specially designed to study particle production in ultra-relativistic heavy-ion collisions. Located 52 meters underground with 28 meters of overburden rock, it has also been used to detect muons produced by cosmic ray interactions in the upper atmosphere. In this paper, we present the multiplicity distribution of these atmospheric muons and its comparison with Monte Carlo simulations. This analysis exploits the large size and excellent tracking capability of the ALICE Time Projection Chamber. A special emphasis is given to the study of high multiplicity events containing more than 100 reconstructed muons and corresponding to a muon areal density $\rho_{\mu}>5.9~$m$^{-2}$. Similar events have been studied in previous underground experiments such as ALEPH and DELPHI at LEP. While these experiments were able to reproduce the measured muon multiplicity distribution with Monte Carlo simulations at low and intermediate multiplicities, their simulations failed to describe the frequency of the highest multiplicity events. In this work we show that the high multiplicity events observed in ALICE stem from primary cosmic rays with energies above $10^{16}$ eV and that the frequency of these events can be successfully described by assuming a heavy mass composition of primary cosmic rays in this energy range. The development of the resulting air showers was simulated using the latest version of QGSJET to model hadronic interactions. This observation places significant constraints on alternative, more exotic, production mechanisms for these events.


Introduction
ALICE (A Large Ion Collider Experiment) [1] designed to study Quark-Gluon Plasma (QGP) formation in ultra-relativistic heavy-ion collisions at the CERN Large Hadron Collider (LHC), has also been used to perform studies that are of relevance to astro-particle physics. The use of high-energy physics detectors for cosmic ray physics was pioneered by ALEPH [2], DELPHI [3] and L3 [4] during the Large Electron-Positron (LEP) collider era at CERN. An extension of these earlier studies is now possible at the LHC, where experiments can operate under stable conditions for many years. ALICE undertook a programme of cosmic ray data taking between 2010 and 2013 during pauses in collider operations when there was no beam circulating in the LHC.
Cosmic ray muons are created in Extensive Air Showers (EAS) following the interaction of cosmic ray primaries (protons and heavier nuclei) with nuclei in the upper atmosphere. Primary cosmic rays span a broad energy range, starting at approximately 10 9 eV and extending to more than 10 20 eV. In this study, we find that events containing more than four reconstructed muons in the ALICE Time Projection Chamber (TPC), which we refer to as multi-muon events, stem from primaries with energy E > 10 14 eV. The detection of EAS originating from interactions above this energy, in particular around the energy of the knee in the primary spectrum (E k ∼ 3 × 10 15 eV), has been performed by several large-area arrays at ground level (e.g. [5][6][7]), while deep underground detectors (e.g. [8][9][10]) have studied the high energy muonic component of EAS. The main aims of these experiments were to explore the mass composition and energy spectrum of primary cosmic rays.
The muon multiplicity distribution (MMD) was measured at LEP with the ALEPH detector [11]. This study concluded that the bulk of the data can be successfully described using standard hadronic production mechanisms, but that the highest multiplicity events, containing around 75-150 muons, occur with a frequency which is almost an order of magnitude above expectation, even when assuming that the primary cosmic rays are purely composed of iron nuclei. A similar study was carried out with the DELPHI detector, which also found that Monte Carlo simulations were unable to account for the abundance of high muon multiplicity events [12]. Several proposals have been put forward in the scientific literature to explain this discrepancy. Some authors suggest that hypothetical strangelets form a small percentage of very energetic cosmic rays [13], while others have tried to explain the excess of high muon multiplicity events by the creation of the QGP in interactions involving high mass primary cosmic rays (iron nuclei) with nuclei in the atmosphere [14].
In this paper, we exploit the large size and excellent tracking capability of the ALICE TPC [15] to study the muonic component of EAS. We describe the analysis of the muon multiplicity distribution with particular emphasis on high muon multiplicity events containing more than 100 muons in a single event and corresponding to an areal density ρ µ > 5.9 m −2 . We employ a description of the shower based upon the latest version of QGSJET [16,17], a hadronic interaction model commonly used in EAS simulations.
Details of the environment of ALICE and the detectors used for this analysis are described in the following section, while the selection of the data and the algorithm adopted to reconstruct atmospheric muons are discussed in Section 3. The muon multiplicity distribution and the study of high muon multiplicity events are described in Section 4. The results are presented in Section 5 and, finally, in Section 6 we make some concluding remarks.

The ALICE experiment
ALICE is located at Point 2 of the LHC accelerator complex, approximately 450 m above sea level in a cavern 52 m underground with 28 m of overburden rock. The rock absorbs all of the electromagnetic and hadronic components of the observed EAS, so that only muons with an energy E, at the surface, larger than 16 GeV reach the detectors [18]. The geometry of ALICE is typical of a collider experiment. A large solenoidal magnet forms a central barrel that houses several detectors, including a large, cylindrical TPC. Outside the solenoid, and on one end, there is a single-arm, forward spectrometer, which was not used in this analysis. A complete description of the apparatus is given in [1].
The ALICE TPC is the largest detector of its type ever built. It was used to reconstruct the trajectory of cosmic ray muons passing through the active volume of the detector, which comprises a cylindrical gas volume divided into two halves by a central membrane. The TPC has an inner radius of 80 cm, an outer radius of 280 cm and a total length of 500 cm along the LHC beam direction. At each end of the cylindrical volume there are multi-wire proportional chambers with pad readout. For the purpose of detecting cosmic ray muons, the total area of the detector due to its horizontal cylindrical geometry is approximately 26 m 2 . However, after placing a cut on the minimum length required to reconstruct a cosmic ray track the maximum effective area reduces to approximately 17 m 2 . The apparent area of the detector also varies with the zenith angle of the incident muons. Track selection is discussed in more detail in Section 3. An example of a single atmospheric muon event is shown in Fig. 1. Three detector subsystems were used to provide dedicated triggers for this study: ACORDE (Alice COsmic Ray DEtector) [19], SPD (Silicon Pixel Detector) [20] and TOF (Time Of Flight detector) [21].
ACORDE is an array of 60 scintillator modules located on the three upper faces of the octagonal yoke of the solenoid, covering 10% of its surface area. A trigger was formed by the coincidence of signals in two different modules (a two-fold coincidence), although the trigger can also be configured to select events when a single module fires or when more than two modules fire.
The SPD is part of the Inner Tracking System located inside the inner field cage of the TPC. It is composed of two layers of silicon pixel modules located at a distance of 39 mm and 76 mm from the LHC beam axis, respectively. The layers have an active length of 28.3 cm, centred upon the nominal interaction point of the LHC beams. The SPD was incorporated into the trigger by requiring a coincidence between signals in the top and bottom halves of the outermost layer.
The TOF is a cylindrical array of multi-gap resistive-plate chambers that completely surrounds the outer radius of the TPC. The TOF trigger requires a signal in a pad corresponding to a cluster of readout channels covering an area of 500 cm 2 in the upper part of the detector and another signal in a pad in the opposite lower part forming a back-to-back coincidence with respect to the central axis of the detector. The configuration of the pads involved in the trigger can be changed via software. In some periods of data taking, this flexibility has been exploited to require a signal in an upper pad and in the opposing pad plus the two adjacent pads forming a back-to-back ±1 coincidence.
Cosmic ray data were acquired with a combination (logical OR) of at least two out of the three trigger conditions (ACORDE, SPD and TOF) depending on the run period. The trigger efficiency was studied with a detailed Monte Carlo simulation, which is discussed in Section 4. Most events were classified as either single muon events or multi-muon events, with a small percentage of "interaction" events where very energetic muons have interacted with the iron yoke of the magnet producing a shower of particles that pass through the TPC.

Event reconstruction and data selection
The TPC tracking algorithm [22] was designed to reconstruct tracks produced in the interaction region of the two LHC beams. It finds tracks by working inwards from the outer radius of the detector where, during collider operation, the track density is lowest. The present analysis used the same tracking algorithm but removed any requirement that tracks should pass through a central interaction point. However, the tracking algorithm has not been optimised for very inclined (quasi horizontal) tracks. Therefore, to avoid reconstruction inaccuracies associated with the most inclined showers, we restricted the zenith angle of all events to the range 0 • < θ < 50 • .
As a consequence of reconstructing tracks from the outer radius of the TPC inwards, cosmic ray muons are typically reconstructed as two separate tracks in the upper and lower halves of the TPC as shown in Fig. 1. We refer to these tracks as up and down tracks. Following this first pass of the reconstruction a new algorithm was applied to match each up track with its corresponding down track to reconstruct the full trajectory of the muons and to eliminate double counting. Starting with single muon events (producing two TPC tracks), where the matching of tracks is straightforward, the reconstruction has been tuned to handle events containing hundreds of muons. High multiplicity Monte Carlo events have been used to optimise the matching performance.
Each TPC track can be reconstructed with up to 159 individual space points. In order to maximise the detector acceptance for this analysis, tracks were required to have a minimum of 50 space points and, in events where the magnetic field was on, a momentum greater than 0.5 GeV/c to eliminate all possible background from electrons and positrons. In multi-muon events, accepted tracks were required to be approximately parallel since atmospheric muons coming from the same EAS arrive almost parallel at ground level. The parallelism cut involves forming the scalar product of the direction of the analysed track t a with a reference track t r , requiring that t a · t r = cos(∆Ψ) > 0.990 to accept the analysed track. The reference track was chosen to give the largest number of tracks satisfying the parallelism cut. This requirement introduces an additional momentum cut due to the bending of muon tracks in the magnetic field. The momentum cut is a function of the azimuth angle of the muon track and varies between 1 and 2 GeV/c. Finally, each up track was matched to the nearest down track if the distance of closest approach between them at the horizontal mid plane of the TPC was d xz < 3 cm. This value was chosen to be large enough to maximise the matching efficiency in high multiplicity Monte Carlo events, while keeping combinatorial background to a minimum.
A muon reconstructed with two TPC tracks (up and down) is called a "matched muon". When a TPC track fulfils all the criteria to be a muon track: number of space points, momentum and parallelism, but does not have a corresponding track within d xz < 3 cm in the opposite side of the TPC, this track is still accepted as a muon candidate but flagged as a "single-track muon". Most single-track muons are found to cross the TPC near its ends where part of the muon trajectory falls outside the detector.
To quantify the performance of the tracking and matching algorithms, we studied the multiplicity dependence of the reconstruction efficiency using Monte Carlo simulated events. We generated 1000 events for 20 discrete values of the muon multiplicity, varying between 1 and 300, which were then reconstructed using the same algorithms applied to real events. In each event, muons were generated parallel to each other like in EAS and cross the whole TPC volume. Fig. 2 shows the mean values (MEAN) and root-mean-square (RMS) of the relative difference between the number of generated and reconstructed muons, as a function of the number of generated muons. The root-mean-square represents the resolution on the number of reconstructed muons and is typically less than 4%, while for the highest multiplicities it is around 2%. The mean value is less than 1% up to N µ ≈ 50, increasing to 5% at high muon multiplicities (N µ ≈ 300).  To illustrate the similarity of the data and the Monte Carlo simulation, Fig. 3 shows the ratio of the number of muons reconstructed as single tracks (either up or down tracks) to the total number of reconstructed muons (both single and matched tracks) for different multiplicities. The ratio obtained from the data is compared with the ratios obtained from simulated samples of pure proton primary cosmic rays and pure iron primaries. Over the range of intermediate muon multiplicities shown, the ratio varies between 0.2 and 0.4 with good agreement between data and simulations. There is no significant difference between the simulated proton and iron samples. Only multi-muon events are discussed further in this paper. We define multi-muon events as those events with more than four reconstructed muons in the TPC (N µ > 4). In total, we collected a sample of 7487 multi-muon events.

Analysis of the data and simulation
To obtain the MMD we have corrected the measured distribution for the efficiency of the trigger. The correction was calculated from a Monte Carlo simulation that is described later in this section. Given the complementary coverage of the TOF barrel to the TPC, the TOF trigger was mainly responsible for selecting events in the low-to-intermediate range of muon multiplicities (7 ≤ N µ ≤ 70). The efficiency of the TOF trigger as a function of the muon multiplicity is shown in Fig. 4. The efficiency is lower at low muon multiplicity due to the back-to-back coincidence requirement of the TOF trigger. The efficiency of the ACORDE trigger has a similar, increasing trend with the muon multiplicity. The multiplicities at which the two triggers reach full (100%) efficiency are N µ > 10 (TOF) and N µ > 15 (ACORDE). Given the much smaller area of the SPD in comparison with the TPC, the efficiency of the SPD trigger is significantly lower than both ACORDE and TOF. It makes only a minor contribution to the MMD in the low-to-intermediate range of muon multiplicities.
The MMD obtained from the whole data sample and corrected for trigger efficiency is shown in Fig. 5.
Values for the systematic uncertainty in the number of events as a function of multiplicity have been estimated by varying the parameters of the track reconstruction and matching algorithms. We find a smooth distribution up to a muon multiplicity of around 70 and then 5 events with a muon multiplicity greater than 100. We define the events with N µ > 100 high muon multiplicity (HMM) events. Given the nature and topology of high multiplicity events, all trigger conditions contributed to this sample with The difficulty in describing EAS, and consequently the number of muons reaching ground level, mainly arises from uncertainties in the properties of multi-particle production in hadron-air interactions. These interactions are often described phenomenologically within Monte Carlo event generators. Model parameters, such as total and inelastic hadron-proton cross sections, inelastic scattering slopes and diffractive structure functions, are constrained by measurements obtained from accelerator experiments.
In this analysis we have adopted the CORSIKA [23] event generator incorporating QGSJET [16] for the hadronic interaction model to simulate the generation and development of EAS. CORSIKA version 6990 incorporating QGSJET II-03 has been used to study the MMD distribution and HMM events; CORSIKA version 7350 incorporating QGSJET II-04 has been used to check and confirm the results for HMM events. The significant differences between the two versions of QGSJET are the inclusion of Pomeron loops in the formalism of QGSJET II-04 and a retuning of the model parameters using early LHC data for the first time [24]. Most relevant to the present study is that pion exchange is assumed to dominate forward neutral hadron production in the QGSJET II-04, which has been shown to enhance the production of ρ 0 mesons resulting in an enhancement of the muon content of EAS by about 20% [25].
In previous studies of cosmic ray muon events at LEP, QGSJET 01 was used to model hadronic interactions. Apart from the way in which nonlinear effects are modelled, another significant difference between this earlier version of the model and QGSJET II-03/04 is the deeper shower maximum, X max , used in the later versions. This results in a steeper lateral muon distribution and an associated increase of the muon density close to the core of the shower, which can also have an impact on the observed rate of HMM events.
When generating cosmic ray events, the core of each shower was scattered randomly at ground level over an area covering 205 × 205 m 2 centred upon the nominal LHC beam crossing point. This area was chosen to minimise the number of events to be generated without creating any bias on the final results. We found that, when the core was located outside this area, only a very small number of events gave rise to muons crossing the TPC and these events were always of low multiplicity (N µ < 4). Therefore, neglecting these events does not affect the results reported in this paper.
To have a fast and flexible method of estimating several important parameters and observables involved in the analysis, we started with a simplified Monte Carlo simulation. This simulation did not explicitly model interactions in the rock above the experiment. Instead, the trajectories of the muons arriving at the surface were simply extrapolated as straight lines to the depth of ALICE while imposing an energy cut E µ > 16 GeV/ cos(θ ), where θ is the zenith angle of the muon. All muons passing this cut and crossing an area of 17 m 2 , corresponding to the horizontal cross-sectional area of the TPC, were considered to be detected.
To understand the complete sample of the recorded data, including the origin of low muon multiplicity events, we generated events initiated by the interaction of proton and iron ( 56 Fe) primaries with energies E > 10 12 eV. This revealed that most single muon events stem from primaries in the energy range 10 12 < E < 10 13 eV, while primaries in the energy range 10 13 < E < 10 14 eV produce muon multiplicities typically in the range from 1 to 4, independent of the mass of the primary cosmic rays. Primaries with energies below 10 14 eV therefore produce a negligible contribution to multi-muon events (N µ > 4) that are of interest in this study. Consequently, only energies E > 10 14 eV were considered in the full simulation.
The first step in the analysis was to attempt to reproduce the measured MMD in the low-intermediate range of multiplicity (7 ≤ N µ ≤ 70). Samples of proton and iron primary cosmic rays were generated in the energy range 10 14 < E < 10 18 eV and with zenith angles in the interval 0 • < θ < 50 • . The composition of cosmic rays in this energy range is a mixture of many species of nuclei in a ratio that is not well-known and which varies with energy. To simplify the analysis and interpretation of the data we have modelled the primary cosmic ray flux using a pure proton sample, representing a composition dominated by light nuclei, and a pure iron sample, representing a composition dominated by heavy nuclei. In relation to the MMD, the proton sample provides a lower limit on the number of events for a given multiplicity, while the iron sample provides an upper limit. A typical power law energy spectrum, E −γ , has been adopted with a spectral index γ = 2.7 ± 0.03 for energies below the knee (E k = 3 × 10 15 eV) and γ k = 3.0 ± 0.03 for energies above the knee. The total (all particle) flux of cosmic rays has been calculated by summing the individual fluxes of the main chemical elements at 1 TeV [26] where measurements are most precise. The flux was estimated to be F(1 TeV) = 0.225± 0.005 (m 2 s sr TeV) −1 .
All events generated with energies E > 10 14 eV were subsequently considered for a complete analysis using a detailed simulation taking into account all possible interactions in matter surrounding the experiment. In each event, all muons were extrapolated to the horizontal mid-plane of the experiment and flagged if they hit an enlarged area of 36 m 2 centred upon the TPC with no restriction on the energy of the muons. All flagged muons were recorded along with their position and momentum at ground level and used as input to the ALICE simulation framework. In this framework, the ALICE experimental hall and the environment above and around the apparatus as well as all the detectors are accurately described. Flagged muons were propagated through this environment with GEANT3 [27]. Any muon that crossed the detector apparatus was treated by a detector response simulation that produced pseudo-raw data, which was then processed with the same reconstruction code that was applied to real data, including the TPC tracking algorithm and the track matching algorithm developed for this analysis.

The muon multiplicity distribution
We generated simulated events equivalent to 30.8 days live time to permit direct comparison with the data without the need to apply an arbitrary normalisation factor. A comparison of the trigger corrected, measured MMD with the simulations is shown in Fig. 6. For ease of comparison, the points obtained with the simulations were fitted with a power-law function to obtain the curves for proton and iron.
At lower multiplicities, corresponding to lower primary energies, we find that the data approach the proton curve, which represents a light ion composition of the primary cosmic ray flux, while higher multiplicity data lie closer to the iron curve, representing a heavier composition. The limited statistics in the range N µ > 30 does not allow for a precise, quantitative study of the composition but suggests that the average mass of the primary cosmic ray flux increases with increasing energy, a finding consistent with several previous experiments [28][29][30][31].
The errors in Fig. 6 are shown separately (statistical and systematic) for data, while for Monte Carlo they are the quadrature sum of the statistical and systematic uncertainties. The systematic errors in the simulations take into account uncertainties in the flux of cosmic rays at 1 TeV, the slope of the energy spectrum below and above the knee, the description of the rock above the experiment and the uncertainty in the the number of days of data taking (detector live time). The largest contribution to the systematic error is due to the uncertainty in the spectral index below the knee (γ = 2.7 ± 0.03), which results in an uncertainty of approximately 15% in the MMD. The error in the description of the rock above the experiment corresponds to an uncertainty in the energy threshold of the muons reaching the detector, which results in a systematic error of approximately 4%. Each of the other uncertainties gives a contribution of around 2% to the systematic error. For muon multiplicities N µ > 30, statistical uncertainties are dominant.
Following success in describing the magnitude and shape of the MMD over this intermediate range of multiplicities (7 ≤ N µ ≤ 70) we have used the same simulation framework to study the frequency of HMM events. Since these are particularly rare events, a very high statistics sample of simulated HMM events was required to permit a meaningful quantitative comparison.

High muon multiplicity events
Taking the dataset as a whole, corresponding to 30.8 days and a mixture of running conditions, we find 5 HMM events with muon multiplicities N µ > 100 (as can be seen in Fig. 5) giving a rate of 1.9 × 10 −6 Hz. Each of these events were examined closely to exclude the possibility of "interaction" events. The highest multiplicity event reconstructed in the TPC was found to contain 276 muons, which corresponds to a muon areal density of 18.1 m −2 . For illustration, a display of this event is shown in Fig. 7. The zenithal and azimuthal angular distributions of the muons from the same HMM event are shown in Fig. 8, while the spatial distribution of matched and single-track muons at the TPC mid plane is shown in Fig. 9. We note that the majority of single-track muons are reconstructed near the ends of the TPC where muons may enter or leave the active volume without producing a track either the upper or lower halves of the detector.
One of the aims with this study is to compare the rate of HMM events obtained from simulations to the measured rate. To limit the effect of fluctuations in the number of simulated HMM events, we have simulated a live time equivalent to one year with CORSIKA 6990 using QGSJET II-03 for the hadronic  interaction model. The simplified Monte Carlo used as a first step of the analysis demonstrated that only primaries with energy E > 10 16 eV contribute to these events. Therefore, only events in the range of primary energy 10 16 < E < 10 18 eV have been generated to achieve an equivalent of 365 days exposure for both proton and iron primaries.
The rate of HMM events obtained with the Monte Carlo can be compared with the observed rate. Since we have simulated samples of HMM events corresponding to one year live time, the statistical uncertainty in the simulated rate will be lower than that in the measured rate. Results obtained for the number of HMM events expected in one year from both the simplified Monte Carlo and the full simulation (the first of five statistically independent simulations) are shown in the first row of Table 1. Comparison of the results demonstrates that the detailed modelling of the underground environment has about a 30% effect on the number of HMM events. Due to the small numbers of HMM events we reused the same simulated EAS sample to perform four additional simulations by randomly assigning the core of each shower over the usual surface level area of 205 × 205 m 2 . Given that the acceptance of the TPC is almost 3000 times smaller, this ensures that the samples are statistically independent. A summary of the results obtained for all five simulations is presented in Table 1 for both CORSIKA 6990 with QGSJET II-03 and CORSIKA 7350 with QGSJET II-04.
Final values for the HMM event rate for proton and iron primaries were calculated by taking the average value obtained from the five simulations, while the statistical uncertainty was estimated from the standard deviation of the 5 values from the mean. Table 2 summarises the mean number of HMM events expected in one year for each primary ion calculated from the full simulation.
There are two major contributions to the systematic uncertainty on the number of HMM events. The first contribution stems from the muon reconstruction algorithm. To estimate its contribution we took the first simulated sample, corresponding to 365 days of data taking, for each element and each CORSIKA code version and redetermined the number of HMM events using different tunes of the track selection and matching algorithms. The second contribution stems from the uncertainties of the parameters used in the simulations, as discussed in section 4.  The systematic uncertainties have been added in quadrature to the statistical uncertainty in the final comparison of the observed rate of HMM events with that obtained from the Monte Carlo simulations.

Results
In Table 3 we present the results of this analysis where we compare the rate of simulated HMM events with the measured rate. We note that the pure iron sample simulated with CORSIKA 7350 and QGSJET II-04 produces a HMM event rate in close agreement with the measured value. The equivalent rate obtained with CORSIKA 6990 and QGSJET II-03 is lower, although still consistent with the measured rate.
The difference between the two simulations comes primarily from the hadronic model used to generate the EAS. It is more difficult to reconcile the measured rate of HMM events with the simulated rate obtained using proton primaries, independent of the version of the model. However, the large uncertainty in the measured rate prevents us from drawing a firm conclusion about the origin of these events, although heavy nuclei appear to be the most likely candidates. Therefore, an explanation of HMM events in terms of a heavy primary cosmic ray composition at high energy and EAS described by conventional hadronic mechanisms appears to be compatible with our observations. This is consistent with the fact that they stem from primaries with energies E > 10 16 eV, where recent measurements [32,33] suggest that the composition of the primary cosmic ray spectrum is dominated by heavier elements.
Finally, we have investigated the distribution of simulated EAS core positions at the location of ALICE for each of the HMM events simulated with iron primaries using CORSIKA 7350 and QGSJET II-04 in Table 1, equivalent to 5 years of data taking. The distribution is shown in Fig. 10, where the colour of each point indicates the energy associated with the primary cosmic ray so as to give a visual representation of the correlation between the distance of the core from the centre of ALICE at surface level and the energy of the primary cosmic ray. We note that the shower cores of all HMM events fall within an area of approximately 140 × 140 m 2 centred upon ALICE, which is located at the origin in Fig. 10. The average distance of the shower core from the centre of ALICE for all events is 19 m and the RMS value of the distribution is 16 m. Primaries with an energy E > 3 × 10 17 eV, corresponding to the highest energy interval studied in this analysis, produce larger showers that may give rise to HMM events when the shower core falls farther from the location of ALICE. In this case, the mean of the shower core distribution from the centre of ALICE is 37 m and the RMS value of the distribution is 18 m.

Summary
In the period 2010 to 2013, ALICE acquired 30.8 days of dedicated cosmic ray data recording approximately 22.6 million events containing at least one reconstructed muon. Comparison of the measured muon multiplicity distribution with an equivalent sample of Monte Carlo events suggests a mixed-ion primary cosmic ray composition with an average mass that increases with energy. This observation is in agreement with most experiments working in the energy range of the knee. Following the successful description of the magnitude of the MMD in the low-to-intermediate range of muon multiplicities we used the same simulation framework to study the frequency of HMM events. High muon multiplicity events were observed in the past by experiments at LEP but without satisfactory explanation. Similar high multiplicity events have been observed in this study with ALICE. Over the 30.8 days of data taking reported in this paper, 5 events with more than 100 muons and zenith angles less than 50 • have been recorded. We have found that the observed rate of HMM events is consistent with the rate predicted by CORSIKA 7350 using QGSJET II-04 to model the development of the resulting air shower, assuming a pure iron composition for the primary cosmic rays. Only primary cosmic rays with an energy E > 10 16 eV were found to give rise to HMM events. This observation is compatible with a knee in the cosmic ray energy distribution around 3 × 10 15 eV due to the light component followed by a spectral steepening, the onset of which depends on the atomic number (Z) of the primary.
The expected rate of HMM events is sensitive to assumptions made about the dominant hadronic production mechanisms in air shower development. The latest version of QGSJET differs from earlier versions in its treatment of forward neutral meson production resulting in a higher muon yield and has been retuned taking into account early LHC results on hadron production in 7 TeV proton-proton collisions. This is the first time that the rate of HMM events, observed at the relatively shallow depth of ALICE, has been satisfactorily reproduced using a conventional hadronic model for the description of extensive air showers; an observation that places significant constraints on alternative, more exotic, production mechanisms.
Compared to the previous studies at LEP, there are two distinguishing aspects of this work that have led to these new insights into the origin of HMM events. The first has been the ability to generate large samples of very energetic cosmic rays, allowing for a more reliable estimate of the expected rate of these events. The second, and more important, aspect has been the recent advances in the hadronic description of EAS. This is a continually evolving field. We note that in a preparatory study [18] carried out by ALICE in 2004, using an older version of CORSIKA (version 6031), no HMM events were observed in the MMD distribution simulated for 30 days of data taking with a pure iron primary cosmic ray composition. In the present work, Table 3 gives a quantitative comparison of the rate of HMM events predicted by two more recent versions of CORSIKA and QGSJET, illustrating the evolution of the hadronic description of EAS in recent years. Only in the latest version of the model there has been a significant increase in the rate of HMM events that better approaches the rate observed in this study.