Composition from high $p_\mathrm{T}$ muons in IceCube

Cosmic rays with energies up to $10^{11}\,\mathrm{GeV}$ enter the atmosphere and produce showers of secondary particles. Inside these showers muons with high transverse momentum ($p_\mathrm{T} \gtrsim 2\,\mathrm{GeV}$) are produced from the decay of heavy hadrons, or from high $p_\mathrm{T}$ pions and kaons very early in the shower development. These isolated muons can have large transverse separations from the shower core up to several hundred meters, together with the muon bundle forming a double or triple track signature in IceCube. The separation from the core is a measure of the transverse momentum of the muon's parent particle. Assuming the validity of perturbative quantum chromodynamics (pQCD) the muon lateral distribution depends on the composition of the incident nuclei, thus the composition of high energy cosmic rays can be determined from muon separation measurements. Vice versa these muons can help to understand uncertainties due to phenomenological models as well as test pQCD predictions of high energy interactions involving heavy nuclei. After introducing the physics scenario of high $p_\mathrm{T}$ muons in kilometer-scale neutrino telescopes we will review results from IceCube in its 59-string configuration as starting point and discuss recent studies on composition using laterally separated muons in the final detector configuration.


Introduction
The primary cosmic ray spectrum has been studied for many decades and incident cosmic rays up to energies of about 10 11 GeV were first observed in the early 1960's [1]. At energies around the knee of the spectrum (∼ 3 PeV) direct measurements by balloon and satellite detectors have shown that cosmic rays include nuclei from hydrogen to iron as expected from supernovae [2]. However, at energies above the PeV-range these experiments have very limited statistics and measurements of cosmic rays are only possible by observing extensive air showers with groundbased experiments. In this case studies of the cosmic ray mass composition rely on indirect observables such as the ratio of the measured electromagnetic energy to the number of muons [3,4] or the atmospheric depth X max , where an air shower reaches the maximum number of particles [5]. These indirect measurements strongly depend on phenomenological models and simulations that relate the observations to the mass composition. Thus, the results will be sensitive to the assumed high energy hadronic interaction models [6].
High energy ( 1 TeV) muons are produced early in the shower development and can therefore be a direct probe of the initial interaction [7]. These muons are produced by the decay of secondary pions and kaons (conventional muons) or from the decay of hadrons containing heavy quarks, mostly charm (prompt muons) [8]. The bulk of muons produced inside the shower are conventional but a e-mail: soldin@uni-wuppertal.de b http://icecube.wisc.edu/collaboration/authors/current at energies above 100 TeV prompt muons are expected to dominate [9].
Inside an air shower high energy muons can be produced that have large transverse momentum (p T ) imparted to them by their parents. These muons will separate from the shower core while traveling to the ground, forming laterally separated muons (LS muons) with separations up to several hundred meters from the core. The resulting lateral separation is a direct measure of the p T of the muon's parent. Experimentally a transition from soft to hard interactions, that can be described in the context of perturbative quantum chromodynamics (pQCD), i.e. for values of p T 2 GeV, should be visible in the p T spectrum and thereby in the lateral separation distribution. The p T spectrum is expected to fall off exponentially with a transition to a power law at roughly 2 GeV [10]. Moreover, assuming the validity of pQCD calculations the muon p T distributions depend on the incident nuclei [11]. Thus, the lateral separation distribution of high energy muons provides a complementary approach to study the cosmic ray mass composition and may help to understand the uncertainties due to phenomenological models as well as test pQCD predictions at high energies and low Bjorken-x.
In IceCube high energy muons are detected with a 1 km 3 array of optical sensors buried at depths between 1450 and 2450 m in the Antarctic ice at the South Pole [12]. Using 5160 digital optical modules (DOMs) on 86 vertical strings Cherenkov radiation produced by charged particles traversing through the ice is detected, such as high energy muons. The surface detector IceTop consists of 162 tanks of ice, each instrumented with two DOMs, that detect cosmic ray air showers at the surface [13]. Ice-Cube and IceTop therefore provide unique opportunities to study high energy muons from cosmic ray air showers.

Muon lateral separation
An isolated high p T muon together with a muon bundle forming the most compact region of the shower, produce a double track signature in IceCube. Figure 1 shows a typical simulated double track event in IceCube. The muon's lateral separation from the shower core is related to the muon's p T via where H int is the primary interaction height, θ is the zenith angle and E µ is the muon's energy at generation which is well approximated by the energy at the Earth's surface. For the relevant muon energies and separations the effect of Earth's magnetic field and multiple scattering is negligible. Hence, the initial p T is the dominant effect producing the separation, although the separation distribution will be biased toward events produced at high altitudes, high p T , and low energy [10]. The horizontal spacing between the IceCube strings of 125 m will be a rough estimate of the minimum resolvable track separation. For a 1 TeV vertical muon produced at an altitude of 50 km this corresponds to a p T of about 2.5 GeV. More inclined arrival directions will decrease this threshold. Thus in this work we will consider transverse momenta of p T ≥ 2 GeV. In high energy primary interactions the parents of LS muons will most likely be produced in high p T particle jets. Due to energy and momentum conservation high p T jets are typically produced back-to-back in the center of mass frame forming a dijet event [14]. If in both jets high p T particles decay into muons, two isolated muons and the central muon bundle can produce a triple track signature in IceCube where the two isolated muons are located roughly on opposite sites of the bundle. Although the expected flux will be low, this would be a distinctive signature and a direct probe of hadronic high energy cosmic ray interactions.
A first analysis of laterally separated muons in IceCube was performed using one year of data in the 59-string configuration [10]. In this analysis high energy data is used where events are required to generate at least 630 photoelectrons corresponding to a primary energy threshold of about 1 TeV [15]. After high energy filtering 6.4 × 10 7 events remain in the data sample and additional reconstructions are applied.

Reconstruction
Double track signatures have a unique topology: a muon bundle accompanied by a LS muon. Thus a dedicated reconstruction algorithm based on a two track hypothesis is needed to reconstruct the muon bundle and a LS muon with the same timing and direction. A fast first guess reconstruction is used that is based on a linear relationship between arrival times of the Cherenkov light at the DOMs and the wavefront which gives the approximate direction of the muon bundle. The hits are rotated into a plane perpendicular to the first guess track and sorted into two sets of hits using the k-means clustering algorithm [16] that sorts the hits according to the closeness of the mean of the cluster. The larger cluster of hits is expected to belong to the muon bundle and is reconstructed using a maximumlikelihood function that accounts for the arrival time of the Cherenkov photons and the scattering of light in ice [17]. Hits that belong to the LS muon are identified by their arrival time relative to the muon bundle: only hits from the smaller cluster that arrive 100 ns earlier than the expectation for light from the muon bundle are considered as LS muon hits. Afterwards this method is iterated to improve the accuracy of the reconstruction, especially to eliminate the few cases where the initial clustering algorithm didn't perform well. The LS muon hits are then reconstructed seperately using the previous maximum likelihood function.

Event selection
Only events are kept where both the reconstruction of muon bundle track as well as the LS muon track reconstruction succeeded, i.e. the likelihood fits found a minimum. The main background of this analysis are cosmic ray air showers without a laterally separated muon that can be divided into two distinct types of background: single air showers and multiple independent coincident showers.
Background from multiple coincident showers is reduced by requiring that the LS muon and muon bundle track directions agree within 5 • and that they arrive within ±450 ns from each other. An irreducible background remains from showers that arrive from the same direction and at the same time. This background is estimated with off-time data, events from the same direction but arriving between 450 and 1350 ns from each other.
Single shower events can be eliminated using the topological properties of single and double tracks. Therefore a reduced likelihood function made for a single track hypothesis is used where single showers are wellreconstructed while double tracks are not. Requiring an output of the (negative) reduced likelihood larger than 7.5 reduces the single shower background significantly. The closer two tracks are, the more difficult it is to separate them. Hence, events with track separations d T < 135 m are removed. To remove surviving events that still show poorly reconstructed directions, we require LS muon hits on more than two strings and the LS muon track to trigger more than 8 DOMs.

Results
After applying all selection criteria 34, 754 events remain in the data, where the expected number of random coincident showers obtained from off-time data is 456. Figure  2 shows the lateral separation distribution of the remaining events as well as the irreducible background estimated from off-time data. Also shown are model predictions from Sibyll 2.1 [18], QGSJET01c [19], and DPMJET 2.55 [20] generated using the CORSIKA [21] simulation program where the simulated cosmic ray spectrum is based on the Hörandel polygonato model [22] with primary energies between 600 GeV and 10 11 GeV.
To test the expected transition from an exponential to a power law in the lateral distribution the distribution was translated to sea level after subtracting the background and then fitted with an exponential function that transitions  Figure 3. The lateral data distribution at sea level and the best fit parameters for equation (2). The exponential part is plotted as dotted red line and the power law as dashed blue line [10]. into a power law with A, B, C, and n allowed to vary. The resulting fit and the corresponding best fit parameters are shown in Figure  3. The fit has a χ 2 /DOF of 30.8/19 with a probability of 4% being a good fit for this distribution while a fit to a purely exponential function has a χ 2 /DOF of 60.5/21 with a probability of 0.001%. Hence, an exponential that transitions to a power law is favoured, as expected from pQCD.
While different hadronic interaction models seem to describe the track separation distribution of the experimental data fairly well, the arrival angle distribution shows a significant disagreement between simulation and data. Figure 4 shows the cosine of the zenith angle of the muon bundle for events that survive all selection criteria. While Sibyll and QGSJET overpredict the event rate at high zenith angles and underpredict the rate for more vertical  Figure 4. Zenith angle distribution of the reconstructed muon bundle track after applying all selection criteria for IC59 data and simulations, as well as irreducible background estimated from off-time data [10].
EPJ Web of Conferences events, DPMJET seems to describe the data better, at least at high zenith angles. A Kolmogorov-Smirnov test of the simulated distributions against the data finds that none of the simulations are in good agreement. The highest correlation is found for DPMJET with a probability of only 5 × 10 −12 that it has been drawn from the same distribution as the data [10].
Besides different approaches of p T modeling between the interaction models it is interesting to note that DPM-JET is the only model that includes a hard component of heavy hadrons including charm while agreeing best with data. However, the zenith angle distribution of prompt muons for DPMJET is rather flat, thus the impact on the zenith angle discrepancy is expected to be small. Moreover, further studies have shown that the difference in zenith angle could arise from a larger fraction of muons from kaon decays. DPMJET for example shows a higher fraction of muons produced by kaons at high zenith angles compared to QGSJET [10]. The discrepancy may also be related to the primary mass composition, since the interaction models show a different behavior in zenith angle for different primary compositions. This indicates that an interplay between kaon/pion treatment, charm abundance, and composition could explain the zenith angle discrepancy.

Composition from high p T muons
pQCD predictions of p T spectra depend on the mass composition of the initial cosmic rays [11]. Nuclei with energy E prim and atomic number A have the same parton distributions as A nucleons, each with energy E prim divided by A. Since the lateral separation of muons is a direct measure of the initial p T it will also be sensitive to the cosmic ray composition. As pointed out, the zenith angle discrepancy, for example, may be related to differences in mass composition. A Kolmogorov-Smirnov test indicates that protononly distributions show the best match with experimental data from IC59 but no final conclusions have been drawn.
To improve the situation, more sophisticated studies on the influence of primary mass composition on high p T muons are pursued. Therefore the first interaction of the cosmic ray with the atmosphere 1 for different primary energies and mass compositions is simulated using the CRMC package [23] as interface to get access to different hadronic interaction models. Generating 5 × 10 4 collisions for each initial energy of log 10 (E prim /GeV) = 4.0, 4.1, . . . , 7.9, 8.0 allows obtaining transverse momentum and energy spectra of particles that may decay into muons: pions (π ± ), kaons (K ± , K 0 ), and heavy hadrons (prompt: and their antiparticles). The resulting p T spectra are fitted in the range 2 GeV ≤ p T ≤ 12 GeV with a power law of the form where a and b are allowed to vary. For a given primary energy only a weak correlation between the secondary particles's p T and its energy is assumed and the energy spectra of each component can also be obtained from the simulations. Figures 5 and 6 show the resulting p T spectrum fits obtained from EPOS-LHC [24], QGSJET II-4 [25], Sibyll 2.1 [18], and HIJING 1.3 [26] for primary energies of 10 5 GeV and 10 6 GeV respectively, which is the typical energy range of cosmic ray events in IceCube after high energy filtering. The distributions are shown for different primary mass compositions: proton-only (red) and irononly (blue). Note that HIJING is the only hadronic interaction model that includes a prompt component. Although different interaction models give very different p T spectra the iron distributions show a significantly steeper spectrum compared to the proton spectra for all models and energies considered here, as expected.
It is interesting to note that especially in iron collisions EPOS produces a significantly larger number of high p T particles per collision with respect to other models while Sibyll seems to underestimate the number of high p T pions and kaons. Moreover, different interaction models produce very different fractions of high p T pions and kaons: Sibyll seems to generate the lowest number of kaons with high transverse momentum compared to other models while EPOS and QGSJET show the largest fraction of kaons. As pointed out before, this may affect the zenith angle distributions and could explain the observed zenith angle discrepancy in IC59. Also notice that the prompt spectra obtained from HIJING show a significantly flatter slope compared to the pion and kaon contributions, as expected from pQCD. Hence, studies of high p T muons seem to be also sensitive to the prompt muon flux, especially at energies where prompt muons are expected to dominate ( 100 TeV).
To obtain p T and energy spectra of secondary pions, kaons, and the prompt component for a given primary energy, each generated p T spectrum is fitted with a power law of the form (3) and a spline fit is applied to the resulting spectral indices b. The spectral indices and corresponding spline fits obtained from EPOS-LHC and HIJING 1.3 simulations are shown in Figure 7. As expected, large differences between the spectral indices of proton and iron distributions can be observed depending on the initial primary energy and the interaction model used. Since the differences between proton and iron spectra are larger than the model uncertainties and because muons will approximately have the same transverse momentum as their parent particles the difference between different primary mass compositions will be visible in the lateral separation distribution of muons.
Using IceCube/IceTop coincident events it is possible to reconstruct double track events in IceCube as shown in Sec. 2.1 to obtain the lateral separation and to use IceTop energy reconstructions [13] to estimate the primary energy of the events. This allows measurement of the spectral index of the LS muon spectrum in bins of the primary energy and thus studies on cosmic ray mass composition based on high p T muons. Moreover, dedicated simulations allow testing of different hadronic interaction models and pQCD predictions at high energies and low Bjorken-x.
For an incident cosmic ray spectrum of the form dN/dE ∝ E −γ where the spectral index γ changes from 2.7 to 3.0 at 3 PeV (the knee) the transverse momentum for a certain primary energy is generated using the spectral index obtained from the spline fits in Figure 7. The resulting p T distributions from 10 7 collisions generated with EPOS and HIJING for primary energies in the range of 10 4 GeV to 10 8 GeV are shown in Figure 8. Following an E −γ primary spectrum 95% of the events are generated with initial energies of log 10 (E prim /GeV) < 4.77 using similar spectral indices (see Figure 7). Thus, the resulting spectra can be approximated by power law distributions and are once again fitted using a power law with the best fit parameters shown in Table 1. The iron distributions show a significantly steeper spectrum and at high p T the expected rates for proton primaries are roughly ten times higher than for iron (lower panel).
Using this method we can simulate LS muon events by adding high p T muons from pion, kaon, and heavy hadron decays to pre-simulated CORSIKA air showers. Transverse momentum and energy at generation are obtained from the parent particle spectrum according to the primary energy. Since the fraction of primary collisions producing at least one high p T particle is known from the simulations, the flux of laterally separated muons can then be normalized to the known muon flux at the ground [9,27]. After detector simulation, triggering, and filtering these simulations allow an optimization of the double track reconstruction for IC86 and the development of selection criteria for a mass composition analysis based on high p T muons.
An analysis of laterally separated muons using one year of IceCube data in its final 86-string configuration is in preparation where data from the in-ice detector is used to measure the lateral separation of muons and Ice-Cube/IceTop coincident events allow studies on the cos-  T applied to EPOS and HIJING simulations for proton (red) and iron (blue) primaries respectively and corresponding spline fits (lines). Assuming an E −γ primary spectrum 95% of the events are generated with initial energies log 10 (E prim /GeV) < 4.77 (dashed gray line). mic ray mass composition. Furthermore, this analysis will include a dedicated search for triple track signatures introduced in Sec. 2. Therefore we modified the reconstruction algorithm described in Sec. 2.1 based on a triple track hypothesis with a clustering into three sets of hits, where the largest of the clusters is defined as the bundle track and the other two sets of hits as LS muon tracks respectively. Although the flux will be low, this will extend high p T muon analyses to very unique signal events with known origin that will therefore be a direct probe of the underlying primary cosmic ray interaction.

Conclusions
Laterally separated muons can be used as a direct probe of high energy cosmic ray interactions because high energy particles with a high transverse momentum are typically produced in the first interaction. It was shown that the lateral separation of muons can be measured in Ice-Cube using dedicated reconstruction methods and selection criteria. Data from the 59-string configuration indicates a transition in the lateral separation distribution from an exponential to a power law, as expected from pQCD. Moreover, simulations with most recent hadronic interaction models have shown that LS muons produced by the decay of high p T particles are sensitive to the mass composition of the initial cosmic ray. Depending on the interaction model used the proton and iron spectra show significant differences. At high p T the expected rates for proton primaries are roughly ten times higher than for iron which is larger than the uncertainties between different interaction models and will therefore be visible in the lateral separation distribution. Using IceCube/IceTop coincident events allows measurement of the primary mass composition as a function of the primary energy. Hence, the analysis of LS muons from cosmic ray air showers using data model a b χ 2 /DOF proton EPOS-LHC 1.1 · 10 8 −5.6 142/97 HIJING 1.3 9.7 · 10 7 −5.4 193/97 iron EPOS-LHC 7.8 · 10 8 −7.8 575/97 HIJING 1.3 2.1 · 10 9 −9.0 1238/95 Table 1. Best fit parameters of power law fits shown in Figure 7. The estimated statistical errors are in the order of less than 1%.
from IceCube/IceTop in its 86-string configuration will be sensitive to the primary mass composition and represents a complementary approach of composition analysis with respect to other cosmic ray experiments. Furthermore, studies on high p T muons, together with dedicated simulations, provide useful information on high energy interactions of heavy nuclei at low Bjorken-x and can help to understand uncertainties in phenomenological models as well as test pQCD predictions.