Neutrino events within muon bundles at neutrino telescopes

The atmospheric neutrino flux includes a component from the prompt decay of charmed hadrons that becomes significant only at $E\ge 10$ TeV. At these energies, however, the diffuse flux of cosmic neutrinos discovered by IceCube seems to be larger than the atmospheric one. Here we study the possibility to detect a neutrino interaction in down-going atmospheric events at km$^3$ telescopes. The neutrino signal will always appear together with a muon bundle that reveals its atmospheric origin and, generically, it implies an increase in the detector activity with the slant depth. We propose a simple algorithm that could separate these events from regular muon bundles.


Introduction
The flux of atmospheric leptons, both muons and neutrinos, is sensitive to the multiplicity and the inelasticity in proton-air, pion-air and gamma-air collisions, probing a forward kinematical region and a high energy regime that are difficult to access at colliders. It is apparent that an accurate description of these hadronic collisions is essential to connect the energy and composition of primary cosmic rays (CRs) with the data at neutrino telescopes and air-shower observatories.
One of the possibilities that has received a lot of attention throughout the years [1][2][3][4][5][6][7] is the production of atmospheric charm. Pions of energy above 30 GeV become less effective producing leptons in the air, as their decay length grows longer than their interaction length. This softens the high-energy spectrum of atmospheric neutrinos, changing their power law from approximately E −2.7 to about E −3.7 [8]. Charmed hadrons, on the other hand, are less frequent inside air showers, but they have a much shorter lifetime than pions and kaons. At energies up to the PeV scale D mesons and Λ c baryons always decay before they can lose energy, so their relative contribution to the atmospheric lepton flux increases with E. It is expected that, depending on the zenith inclination, * at energies around 100 TeV [4,9] this charm component may dominate the atmospheric lepton flux.
Moreover, any estimate of the neutrino flux from charm decays cannot avoid a significant degree of uncertainty. The reason is easy to understand. The primary CR flux is very steep, and secondary hadrons will be produced according to the same E −2.7 power law. Consider then an atmospheric D meson of energy E. We may wonder what is the most likely energy of its parent hadron. The D may come from a hadron of energy just 10 times larger (i.e., the D took a fraction x = 0.1 of the collision energy), but also from a parent 1000 times more energetic (x = 10 −3 ). Of course, a collision with x = 0.1 is more unlikely than one with x = 10 −3 , but this may be compensated by the fact that hadrons of energy 10 3 E are much more rare than those of just 10 E. It turns out that a few collisions where the charmed hadrons take a large fraction of the collision energy could increase very substantially their production power law in the atmosphere and thus the flux of neutrinos resulting from their decay.
Perturbative QCD calculations [10] focus on transverse charm and are able to reproduce very accurately the LHC data, but they do not include non-perturbative effects that may be important at forward rapidities. In particular, the factorization theorem used in these calculations implies that the fragmentation of the charm quarks produced in the collision should be independent from the initial state. Fixed target experiments like E791, however, contradict this scheme [11]. In π − collisions with Carbon and Platinum targets at 500 GeV they observed forward events of large x where the c goes into a D 0 or thec into a D − much more likely than into a D + or aD 0 , respectively. These leading charm hadrons appearing in the fragmentation region share a valence quark with the incident pion, suggesting a process of coalescence during hadronization. Another possibility that may be difficult to probe at colliders is that of diffractive charm. One may think, for example, of a 10 TeV proton scattering off an air nucleus with a diffractive mass m * p ≈ 5 GeV and then going into a final Λ cD pair that carries all (or most of) the initial energy. A 1% component of intrinsic charm [12] in protons and pions could favor these processes and imply that the forward charm [13][14][15] contribution to the atmospheric neutrino flux completely dominates over the perturbative one. † Unfortunately, the search for atmospheric neutrinos from charm decays has been so far unsuccessful. However, IceCube observed in 2013 [17,18] a diffuse flux of cosmic neutrinos that at E > 30 TeV is several times larger than the total atmospheric flux. In the 30-500 TeV region its spectral index α seems similar to what we may expect from atmospheric charm (α ≈ 2.7), whereas at PeV energies the cosmic flux becomes harder (α ≈ 2.0-2.3). Although this flux is a great discovery, it makes the possibility to detect neutrinos from charm even more difficult. In upgoing or near-horizontal events both fluxes are indistinguishable [19], as they are expected with the same angular distribution and imply a similar ratio of shower to track (with a muon after the interaction) events. Actually, the best fit obtained by IceCube from the data on high energy events is no neutrinos from charm at all. Obviously, their analysis is performed trying to minimize the atmospheric background, i.e., cutting any events where muons enter the detector from a down-going direction.
Here we will explore the opposite possibility. We will focus on down-going events, where the neutrino signal appears together with a muon bundle that, in turn, guarantees its atmospheric origin. Arguably, this is what will be needed to determine the prompt neutrino flux. Other approaches (spectrum and lateral separation of large p T muons [20][21][22]) focus on down-going events as well. Our analysis will involve two main aspects that we study by using the air shower simulator CORSIKA [23]: (i) the relation † For a CR spectral index α ∼ 2.7, the flux of atmospheric charm is proportional to the 1.7-moment (Z) of the yield in hadronic collisions [16]. Therefore, D mesons produced with a 50% of the energy of the projectile in 0.01% of hadronic collisions (Z = 3.0 × 10 −5 ) would contribute to the prompt flux 4 times more than D's produced with 0.1% of the energy in 100% of hadronic collisions (Z = 7.9 × 10 −6 ). between a neutrino of given energy and the energy of its parent air shower and (ii) the characterization of muon bundles from CR primaries of any energy and composition (in Sections 2 and 3). Then we will analyze the longitudinal energy depositions through the ice or water in a down-going event with or without a ν interaction; these depositions determine the detector activity at km 3 observatories like IceCube or KM3NeT [24]. Finally, we propose an algorithm based on four observables that could be used to separate events with an atmospheric neutrino interaction from events with just stochastic energy depositions of a muon bundle (in Section 4).

Neutrinos and their parent cosmic ray
Let us start with the following question. Suppose we observe an atmospheric neutrino of energy E ν = 10 TeV entering a km 3 telescope from a zenith inclination θ z = 45 • . What is the energy of its parent CR? Obviously, this neutrino may have been produced by a CR of any energy E > E ν , so the actual question is: What is the probability distribution of the parent energy? The answer will depend on two basic quantities: the yield of neutrinos of energy E ν produced per proton air shower of energy E, and the primary CR spectrum and composition at E > E ν .
We may express the neutrino yield per proton shower as f pν (x, E), where x = E ν /E is the fraction of the shower energy taken by the neutrino. Notice that x takes values between 0 and 1, that the integral of f pν (x, E) between these two values gives the total number of neutrinos produced inside the shower, and that if instead we integrate x f pν (x, E) we will get the fraction of the shower energy carried by all these neutrinos.
We have used CORSIKA with SIBYLL 2.3C [9] as the hadronic interaction model to deduce the ν e and ν µ yields from proton primaries of E = (10 3 , 10 4 , . . . , 10 8 ) GeV, and we have obtained a simple fit that performs well in this energy interval (see the details in appendix A). In Fig. 1 we plot the total yields at three different energies from θ z = 45 • (we provide the zenith angle dependence in the appendix) together with our fit. The plots show that lower energy showers are more likely to include a neutrino carrying a large fraction x of the shower energy. These yields must be understood as the sum of two contributions: conventional neutrinos from pion and kaon decays plus neutrinos from the decay of charmed hadrons. The lower plot expresses the relative contribution of these two components for an average 10 6 GeV proton shower. We find that CORSIKA gives the same ν e and ν µ yields from charm and an almost perfect scaling (i.e., the charm contribution in this plot does not depend on E; we will neglect the ≈ 2% ν τ component from D s decays). The plot also shows, for example, that in the 10 6 GeV proton shower charm decays dominate the production of ν e 's of E ν > 7 TeV, or that a ν µ of E ν > 100 TeV inside the same shower is still 4 times more likely conventional than from charm.
From these yields in proton showers we can easily estimate the ones for other primaries, like He or Fe. In particular, assuming that a nucleus of mass number A and energy E is the superposition of A nucleons of energy E/A, we obtain As mentioned above, the second key element to relate an atmospheric neutrino with for ν e (red) and ν µ (blue) both for p (solid) and He (dashed) primaries.
its parent shower is the primary CR flux. At energies below E knee = 10 6.5 GeV we will assume that it is dominated by proton and He nuclei with slightly different spectral indices [25], and a similar number of protons and He nuclei at E ≈ 10TeV. Beyond the CR knee, up to E ankle = 10 9.5 GeV, the composition is uncertain, while the total flux becomes Φ = 330 (E/GeV) −3.0 [particle (GeV cm 2 s sr) −1 ]. Throughout our analysis we will consider the limiting cases with a pure proton or a pure Fe composition at E > E knee and will take a central case where the composition is assumed to be protons and He nuclei in the proportion estimated at E = E knee .
Our results are summarized in Fig. 2. At 10 TeV and θ z = 45 • the atmospheric neutrino flux Φ ν includes a 6% of ν e and a 94% of ν µ . For a pure proton composition above E knee , 66.3% (67.5%) of Φ νµ (Φ νe ) comes from proton showers, whereas this percentage decreases to 61.6% (61.3%) for a pure Fe composition. When the parent is a proton (solid lines in Fig. 2 Right), the fraction of energy taken by the neutrino is distributed according to where Φ p ν is the contribution from proton primaries to Φ ν . An analogous expression describes the fraction of energy taken by neutrinos coming from a He primary (dashed lines in the same figure). We obtain that a ν µ (blue lines) carries in average 13% of the shower energy when the primary is a proton or 3% when it is a He nucleus. For electron neutrinos (red lines) these average fractions are a bit smaller: 9% and 2%, respectively.
Our results may seem somewhat surprising. It is apparent that most of the neutrinos produced by a CR primary of energy E will carry a very small fraction of the shower energy (see Fig. 1), however, the rare events where the neutrino takes a large fraction of this energy dominate Φ ν . The steep fall of the CR flux with the energy suppresses the contribution of neutrinos with a small x, i.e., inside very energetic showers. We find, for example, that when the primary is a proton 75% of muon neutrinos of E ν = 10TeV come from showers of E < 232 TeV, and that the ratio x = E ν /E grows even larger at lower neutrino energies (e.g., at E ν = 1 TeV, E < 21 TeV). These results, fully compatible with the ones in [9], imply that most neutrino events take place inside relatively weak muon bundles.

Leading muon and muon bundle
Atmospheric muon neutrinos will always be produced together with a µ ± of similar energy. Since these neutrinos carry a significant fraction of the shower energy, it follows that there will be a leading muon of energy well above the average muon energy in the bundle at the core of the air shower. This leading muon will be absent in ν e events.
It is straightforward to parametrize its energy distribution using the CORSIKA simulations described in the previous Section. Suppose that a proton shower of energy E produces a ν µ of energy E ν = xE with x > 10 −3 ; let us define the energy of the leading muon as E µ ≡ e αµ xE (i.e., E µ = E ν for α µ = 0). We find that the distribution of α can be fitted with ‡ with E given in GeV. The energy distribution of the muon accompanying the neutrino produced inside a shower of energy E ν /x is theñ  This distribution will be independent from the zenith inclination of the primary but not its composition. When the primary is a nucleus of mass number A, the distribution is obtained just by changing E → E/A and x → Ax in the expression above.
Suppose that a 200 TeV proton shower produces a 10 TeV ν µ ; we find that the leading muon has in this case an average energy of E µ = 7.5 TeV, and that with a 50% probability E µ < 5.4 TeV. We find remarkable that, although in average muons take more energy than neutrinos in meson decays, the leading muon inside a shower with a very energetic (x > 10 −3 ) neutrino carries a smaller fraction of energy. Our result reflects that the neutrinos emitted forward in the meson decay contribute to Φ ν more than the ones emitted backwards, and in the first case the muon takes a smaller fraction of the meson energy.
To study the possibility to detect atmospheric neutrino interactions in down-going events, a precise characterization of the muon bundle in the core of the air shower is also essential: we use CORSIKA to obtain a fit of their number and energy distribution. In a proton shower of energy E from θ z = 45 • the muons of E µ ≥ 500 GeV are distributed according to (see Fig. 3) (all energies in GeV) up to E µ ≈ 0.2E, with a 30% dispersion with respect to this central value. The total number of muons N p µ (E) is then obtained by integrating this expression. Again, we will approximate the bundle in the shower started by a nucleus as the sum of A proton bundles of energy E/A. As for the zenith angle dependence, it can be approximated by the same factor that multiplies the conventional yield in the expression (A.2) given in the appendix. Our results on the spectrum and number of muons in a bundle are consistent with the ones discussed in [21].
Once the muons penetrate the ice or water, they will lose energy through four basic processes: ionization, pair production, bremsstrahlung and photohadronic interactions. We will use the differential cross sections dσ/dν for these processes in [26], where ν is the fraction of the muon energy deposited in these collisions with Hydrogen and Oxygen nuclei. To simulate the propagation of each individual muon we define steps of 25 m and separate soft collisions that imply a continuous energy loss from harder stochastic processes. In the first type we include both ionization and radiative collisions of ν < 10 −2.5 .
In Fig. 4 we provide examples of the propagation of muons and of muon bundles through several km of ice, together with the average energy deposited per 100 meters at different depths for bundles from proton primaries of 10 5 , 10 6 and 10 7 GeV. The average in the plot is obtained for 10 4 showers of each energy; it approximately scales like E 0.81 with the shower energy and like exp − X 1100 m.w.e. with the slant depth.

Neutrino events within a bundle
Our objective is to establish criteria to separate muon bundles that include a neutrino interaction from those bundles that do not. These criteria or cuts should be very efficient eliminating plain bundles while selecting a significant fraction of the events with a neutrino interaction.
For each event we define four basic parameters related to observables that can be measured with different degree of precision in telescopes like IceCube or KM3NeT: 1. X A : age of the track, i.e., the slant depth from the ground to the point of entry in the detector. X A depends on the inclination and the coordinates of the event.
2. E max : maximum energy deposition within a 100 meter interval along the track crossing the detector.
3. E − : total energy deposited in the detector before the maximum deposition E max divided by the number of 100 meter intervals. Our unit length is set at 100 meters, the typical separation between strings at km 3 telescopes.
4. E + : total energy in the detector after E max divided by the number of 100 meter intervals.
The number of 100 meter intervals before and after the maximum deposition will depend, like X A , on the inclination and coordinates of each event. We will define cuts in terms of the ratios E max /E − and E + /E − .
Let us first consider charged current (CC) ν e events, with all the neutrino energy deposited in a single 100 m interval. A typical 1 TeV event will come together with the weak muon bundle of a E ≤ 20 TeV shower, able to reach the telescope only from vertical directions. These events would imply a value of E max ≥ 30E − . A similar deposition E max could as well be produced by a muon that reaches the telescope with an energy of, for example, 2 TeV. However (i) such muons deposit around 50 GeV in each 100 m interval previous to E max , and (ii) they usually appear inside more energetic showers, together with other muons that also contribute to E − and reduce the value of E max /E − . In addition, this type of depositions subtracts a significant fraction of energy to the muon, implying a drop in the signal after E max . Notice that this effect would be absent when the energy deposition is caused by the neutrino. Requiring that E + ≥ 0.9E − we would make sure that those events do not pass the cut.
If we increase the energy by a factor of 10 and target 10 TeV CC ν e events, two competing effects are noticeable. On one hand, stochastic energy depositions grow linearly with the energy of a muon, while the growth of its continuous energy loss is a bit slower. § This first effect suggests that we should increase the minimum value of E max /E − required to select neutrino events. On the other hand, however, the scaling also implies stronger muon bundles giving a more sustained and regular deposition: if a 50 TeV muon reaches the detector after crossing a depth X A , it is likely that other less energetic muons in the same bundle will reach as well. Although both effects tend to cancel, we conclude that a more effective cut to select ν e events must depend on E max : As we have already mentioned, the cut is a condition on the ratio between the energy deposited after and before the maximum deposition. If E + /E − < 0.9 then there is a chance that E max has been caused by a very energetic muon, whereas a revival of the signal by a factor of 1.5 after E max , E + /E − > 1.5, is only expected in ν µ events (see below). The E + < 0.02E − possibility in the search for ν e events is added to include neutrinos interacting after all muons in the bundle have stopped. Finally, we will also require that the track intersecting the detector must have a minimum length of 500 m, with at least 200 m before the maximum energy deposition (i.e., E − is obtained as the average over at least two 100 meter intervals).
The characterization of CC ν µ events is equally simple. The two main differences with the case just discussed are that (i) the ν µ will deposit in the interaction point only a fraction of its energy and (ii) it will create a muon of similar energy. Again, it is essential that the main contribution to the atmospheric neutrino flux comes from primaries of energy (per nucleon) just 5-20 times larger. A typical event will consist of a ν µ together with a leading muon and a bundle: first the propagation (a large enough age X A ) weakens the muon track entering the detector, then there is a significant energy deposition (E max E − ) followed by a track that is revived by the final muon (E + > E − ). We can use This condition is fully effective when the track inside the detector includes at least two 100 m length intervals before and after E max . § Notice that this includes ionization but also radiative processes of ν < 10 −2.5 . Distribution of E max /E − for 6 × 10 5 bundle "events" (see text) from 10 4 proton showers of 10 PeV. We have separated the events with E max < 1 TeV. Right. Distribution of E + /E − for the 42 events that pass the first cut. We include the distribution when we reduce (dashes) or increase (dots) in a 10% the E max /E − cut.
We have simulated and analyzed a sample of 10 4 muon bundles from proton and He showers of energy between 10 4 and 10 8 GeV at X A > 1500 m.w.e. Our procedure has been the following. First we generate the bundle. Then we make a Monte Carlo simulation of its propagation through the ice or water, where it defines a track of energy depositions. We take the track after a slant depth of 1500 m.w.e. and divide it in 500 or 900 meter intervals ("short" and "long" events): each one of these intervals might be intersecting the telescope and define an event. For each event, we take 100 meter segments, we determine the total energy deposition in each segment, we find E max , E − and E + and finally we apply the cuts.
We find no segment with 5 or 9 length intervals (i.e., a 500-900 meter track inside the detector) that passes the cuts established above and gives a false positive. In Fig. 5 (left) we plot the distribution of E max /E − for E shower = 10 PeV, which is the energy with the largest fraction of events passing the first cut (all energies give qualitatively similar results). In the plot we separate the events with E max ≤ 1 TeV, which have the cut at E max /E − ≈ 30 (at higher values of E max the cut is closer to 60). When we apply this first cut plus the requirement of at least two length intervals before E max , 42 out of the 6 × 10 6 events survive. Then we apply the cut on E + /E − . In Fig. 5 (right) we show how this variable is distributed among the 42 events: none of them passes the second requirement to be classified as a ν e (0.8 < E + /E − < 1.5) or a ν µ (1.5 < E + /E − ) event.  Figure 6: Energy depositions in 100 meter intervals of ice or water for a 10 6 GeV proton (left) or a He (right) shower. Cyan and blue dots indicate, respectively, the depositions before and after the ν µ CC interaction (see text). We include in red dots the depositions of the muon bundle without the leading muon; the ν e energy deposition defining E max in this case has been omitted.
If we relaxed in a 10% the cut on E max /E − , 71 events would pass it and 2 of them would give a false positive (they both would be declared ν e events, see Fig. 5). In contrast, a 10% increase in the minimum value of E max /E − implies that only 21 events pass the first cut and all of them are clearly excluded by their value of E + /E − . The events more likely to give a false positive appear when a single muon is produced with a large fraction (above 1%) of the shower energy [28].
As for the real neutrino events, when we include an arbitrary neutrino interaction that passes the cuts we find that the prescription separating ν µ from ν e events is very efficient. In particular, we find that ν µ CC interactions are never taken as a ν e event, whereas the opposite case (ν e interactions confused with a ν µ CC event) has a frequency below 10%.
Let us illustrate these results with a couple of examples. In Fig. 6 we provide the energy depositions produced by a 10 6 GeV proton (left) or a He (right) shower. The red dots correspond to a typical muon bundle for a ν e event (with no leading muon); the proton shower in the plot includes 10 muons of energy between 500 GeV and 2.1 TeV, whereas the He shower generates 17 muons of energy between 500 GeV and 3.6 TeV. In both cases, we see that any ν e CC event of energy above 10 TeV (omitted in the plot) at a slant depth X A ≥ 2000 m.w.e. would pass the cut for ν e events defined in Eq. (4.1).
In a contained ν µ event (cyan and blue dots in the same figure) we expect a leading muon after the interaction: we have added to the muon bundle a 30 TeV leading ν µ plus a 14 TeV muon (proton shower) or a 14 TeV ν µ plus a 10 TeV muon (He shower). The ν µ experiences then a CC interaction of inelasticity 0.46 and 0.22, respectively, in both cases at 4100 m.w.e. If the interaction occurred inside the detector, the revival of the track produced by the final muon would imply that the event passes the ν µ cut. If the same deposition E max were produced by a an isolated 10 TeV muon, instead of stronger the track afterwards would have become significantly weaker.

Summary and discussion
The determination of the atmospheric neutrino flux at energies from about 1 TeV up to several 100 TeV is essential both in the search for atmospheric charm and for a precise characterization of the high energy diffuse flux recently discovered by IceCube. However, this atmospheric flux is difficult to access with ν telescopes, as at E ≈ 10 TeV it seems to be 5-10 times weaker than the astrophysical one. Any possibility to disentangle these two components in the flux and search for neutrinos from charm decays requires the detection of ν interactions in down-going events, where the presence of additional muons will reveal the atmospheric origin.
Here we have explored that type of events. Our analysis focuses on the muon bundle produced in the core of the air shower together with the neutrino. In particular, we have studied the energy depositions as the bundle propagates in ice or water. The longitudinal pattern of depositions would translate into a particular signal in a km 3 telescope. Our objective has been to show that this pattern could be different enough when it includes a neutrino interaction.
Our first observation has been that most atmospheric neutrinos are produced inside air showers that are just ten times more energetic. As a consequence, its relative effect on the signal associated to the muon bundle tends to be very large. The typical topology is a weak signal entering the detector, followed by a large energy deposition, and finally a stronger signal in case of a CC ν µ interaction or a weak one in a ν e or NC interaction. Generically, neutrino events imply a signal that increases with the slant depth inside the telescope, while muon bundles tend to imply the opposite effect.
We have defined cuts based on the ratios E max /E − and E max /E + (see Section 4) that seem to exclude muon bundles of any energy. A muon can certainly have a stochastic deposition of half its energy, but not without leaving a trace both before (E − ) and after (E + ) this E max . In 10 4 simulations of muon bundles, we find that when the ratio E max /E − is very large then the signal E + is significantly weaker (relative to E − ) than in a CC ν µ or a ν e event (E + < 0.8 E − ). The bundle events that are closest to the cuts include one single muon carrying a significant fraction of the shower energy that deposits a large fraction of its energy when the rest of the bundle is already weak. Actually, the search for this type of muon events could be of interest by itself [21,27,28] and seems also possible.
Our results should be considered just a first step in the search for neutrino interactions in down-going events at ν telescopes. We show that there are basic physics criteria that could separate these events from plain muon bundles. Of course, to determine in detail whether or not an analysis along these lines could give positive results in actual observations would depend on the experimental conditions (volume, energy resolution, triggers, etc.) at each observatory. Neutrino telescopes have been built to look for high energy sources and avoid the atmospheric background. However, they have also pursued other more unlikely but equally interesting objectives: IceCube has been able to define a strategy to target transient event of energy as low as 1-10 GeV [29], to look for high p T muons [20,22], to determine the atmospheric muon flux at E µ ≥ 10 TeV [21] or to reconstruct starting muon tracks [30]. A more precise characterization of the atmospheric neutrino flux at 1-100 TeV seems a very interesting objective as well.

A Neutrino yields
In our parametrization of the yields in proton showers we have separated the conventional f conv pν (x, E) and the prompt f charm pν (x, E) contributions. The yields refer to the sum of ν i +ν i , with i = e, µ, and we express them in terms of four energy and flavor dependent parameters (A, B, C, D) as: where m µ is the muon mass. From the CORSIKA-SIBYLL 2.3C simulation (10 4 showers of each energy with E min = 10 −3 E shower and 50 showers with E min = 1 GeV) we deduce the value of the 4 parameters for each flavor at six different proton energies, and then we interpolate (linearly in log E) inside each energy interval.
For the conventional yield at θ z = 45 • we obtain the values given in Table 1. The angular dependence (see Fig. 7) may be described in terms of the zenith angle of the line sight at h = 30 km, θ * (θ z ), defined in [8]: with R ⊕ the radius of the Earth. We fit For the neutrinos from charm decays we obtain similar ν µ and ν e yields and no energy dependency in the 4 parameters: