Physics Reach of High-Energy and High-Statistics IceCube Atmospheric Neutrino Data

This paper investigates the physics reach of the IceCube neutrino detector when it will have collected a data set of order one million atmospheric neutrinos with energies in the 0.1 \sim 10^4 TeV range. The paper consists of three parts. We first demonstrate how to simulate the detector performance using relatively simple analytic methods. Because of the high energies of the neutrinos, their oscillations, propagation in the Earth and regeneration due to \tau decay must be treated in a coherent way. We set up the formalism to do this and discuss the implications. In a final section we apply the methods developed to evaluate the potential of IceCube to study new physics beyond neutrino oscillations. Not surprisingly, because of the increased energy and statistics over present experiments, existing bounds on violations of the equivalence principle and of Lorentz invariance can be improved by over two orders of magnitude. The methods developed can be readily applied to other non-conventional physics associated with neutrinos.


I. INTRODUCTION
After the discovery of neutrino oscillations in underground experiments, the observations have been confirmed by experiments using "man-made" neutrinos from accelerators and nuclear reactors [1]. We are entering an era of precision neutrino physics. In this context we discuss the unique potential of IceCube, an experiment that will collect large statistics samples of high energy atmospheric neutrinos. In contrast with its other missions, the beam and its physics exploitation are guaranteed.
With its high statistics data [2] Super-Kamiokande (SK) established beyond doubt that the observed deficit in the µ-like atmospheric events is due to oscillations, a result also supported by the KEK to Kamioka long-baseline neutrino oscillation experiment (K2K) [4] and by the MACRO [5] and Soudan 2 [6] experiments.
It has been recognized that oscillations are not the only possible mechanism for atmospheric ν µ → ν τ flavour transitions [7]. These can also be generated by a variety of nonstandard neutrino interactions characterized by the presence of an unconventional interaction (other than the neutrino mass terms) that mixes neutrino flavours [7]. Examples include violations of the equivalence principle (VEP) [8,9,10], non-standard neutrino interactions with matter [11], neutrino couplings to space-time torsion fields [12], violations of Lorentz invariance (VLI) [13,14] and of CPT symmetry [15,16,17]. From the point of view of neutrino oscillation phenomenology, a critical feature of these scenarios is a departure from the E energy dependence of the conventional oscillation wavelength [18,19]. Although these scenarios no longer explain the data [20,21,22,23,24,25], a combined analysis of the atmospheric neutrino and K2K data can be performed to obtain the best constraints to date on the size of subdominant oscillation effects [26] 1 .
In contrast to the E energy dependence of the conventional oscillation length, new physics predicts neutrino oscillations with wavelengths that are constant or decrease with energy. IceCube, with energy reach in the 0.1 ∼ 10 4 TeV range for atmospheric neutrinos, is the ideal experiment to search for new physics. For most of this energy interval standard ∆m 2 oscillations are suppressed and therefore the observation of an angular distortion of the atmospheric neutrino flux or its energy dependence provide a clear signature for the presence of new physics mixing neutrino flavours.
In this paper we explore the physics that can be probed with the high-statistics highenergy atmospheric data that will be collected by the IceCube detector. In particular we quantify its sensitivity to atmospheric neutrino oscillations driven by new physics effects. The outline is as follows: Our analytic "simulation" of the IceCube detector is described in Sec. II where the expected number of atmospheric neutrino events and their energy 1 Atmospheric neutrinos have also been used to place bounds on other exotic forms of new physics in the neutrino sector such as the possibility of neutrino decay [27,28] or quantum decoherence in the neutrino ensemble [29,30]. distribution are presented. In Sec. III we briefly summarize the formalism for discussing the phenomenolgy of non-standard neutrino oscillations and we derive the evolution equations that describe a high-energy neutrino beam subject to oscillations as well as attenuation and ν regeneration due to τ decay [32,33,34]. In Sec. IV we illustrate the sensitivity of the detector for violations of Lorentz invariance and the equivalence principle.

II. SIMULATION OF MUON EVENT RATES IN ICECUBE
In a high energy neutrino telescope muon neutrinos are detected via their charged current (CC) interactions in the matter surrounding the detector. Such interactions produce muons which reach the detector. High energy muons have very large average range resulting in an effective volume of the detector significantly larger than the instrumented volume.
In our semianalytical calculation we will obtain the expected number of ν µ induced events from dEν d cos θ is the differential muon neutrino neutrino flux in the vicinity of the detector after evolution in the Earth matter (see next section for details). We use as input the neutrino fluxes from Honda [35] which we extrapolate to match at higher energies the fluxes from Volkova [36]. At high energy prompt neutrinos from charm decay are important. In order to estimate the uncertainty associated with the poorly known charm meson production cross sections at the relevant energies, we compute the expected number of events for two different models of charm production: the recombination quark parton model (RQPM) developed by Bugaev et al [37] and the model of Thunman et al (TIG) [38] that predicts a smaller rate.
is the differential CC interaction cross section producing a muon of energy E 0 µ and n T is the number density of nucleons in the matter surrounding the detector and T is the exposure time of the detector. Equivalently, muon events arise fromν µ interactions that are evaluated by an equation similar to Eq.(2). After production with energy E 0 µ , the muon ranges out in the rock and in the ice surrounding the detector and looses energy. We denote by F (E 0 µ , E fin µ , l) the function that describes the energy spectrum of the muons arriving at the detector. Thus F (E 0 µ , E fin µ , l) represents the probability that a muon produced with energy E 0 µ arrives at the detector with energy E fin µ after traveling a distance l. We compute the function F (E 0 µ , E fin µ , l) by propagating the muons to the detector taking into account energy losses due to ionization, bremsstrahlung, e + e − pair production and nuclear interactions according to Ref. [39]. We include in F (E 0 µ , E fin µ , l) the possibility of fluctuations around the average muon energy loss (using the average energy loss would equalize l to the average muon range distance). Thus in our calculation we keep E 0 µ , E fin µ , and, l as independent variables. For simplicity we use n T and F (E 0 µ , E fin µ , l) in ice and we account for the effect of the rock bed below the ice in the form of an additional angular dependence of the effective area for upward going events (see Eq. (10 below)).
The details of the detector are encoded in the effective area A 0 ef f . We use the following phenomenological parametrization of the A 0 ef f to simulate the response of the IceCube detector after events that are not neutrinos have been rejected (this is achieved by quality cuts referred to as "level 2" cuts in Ref. [40]) In A 0 (E f in µ ) we include the energy dependence of the effective area due to trigger requirements (see Ref. [40]). We find good agreement with the results for the experiment MC simulation if we introduce a simple straight line dependence on log 10 where A 0 is an overall normalization constant which is fixed to reproduce the expected number of events in the absence of oscillations: 91000 events/yr after level 2 cuts for conventional atmospheric neutrinos. We next have to "simulate" cuts introduced in Ref. [40] in the muon tracklength l min and the number of optical modules reporting signals in an event N CH,min . R(l min ) represents the smearing in the minimum track length cut, l min = 300 m, due to the uncertainty in the track length reconstruction which we parametrize by a Gaussian with σ l = 50 m. The angular dependence of the effective area for downgoing events (θ < 80 • ) is determined by the level 2 cut on the minimum number of channels N CH > N CH,min (cos θ) = 150 + 250 cos θ. We translate this requirement in an E f in µ -dependent angular constraint as where N CH E f in µ is the average channel multiplicity produced by a muon which reaches the detector with energy E f in µ and σ 2 N CH is the spread on the distribution. Using Fig. 7 in Ref. [40] we obtain the parametrization Finally, we can account for the presence of the rock bed below the detector by introducing a phenomenological angular dependence of the effective area for upward going muons independent of the muon energy. We show in Fig. 1 the effective area A ef f , defined as the ratio of the number of upgoing muon events, with/without the inclusion of A 0 (E f in µ ) × R(l min ) and the level 2 cuts on l min , and compare our results to the detector simulations after cuts from Fig.5 of Ref. [40]. Our calculation correctly reproduces the energy threshold of the effective area.   In Fig. 3 we show the expected spectrum of events in the absence of oscillations after level 2 cuts as a function of the muon energy at the detector, E f in µ (full line). For comparison we also show the spectrum as a function of the muon energy before ranging E 0 µ . From the figure we read that in 10 years of operation IceCube will collect more than 7 × 10 5 atmospheric neutrino events with energies E f in µ > 100 GeV. These events are generated by neutrinos with large enough energy for the standard ∆m 2 oscillations to be very much suppressed so they should behave as flavour eigenstates. This high-statistics high-energy event sample offers a unique opportunity to test new physics mechanisms for leptonic flavour mixing as we discuss next.

III. PROPAGATION IN MATTER OF HIGH ENERGY OSCILLATING NEUTRINOS
As described in the introduction, new physics (NP) scenarios can result in lepton flavour mixing in addition to "standard" ∆m 2 oscillations. We concentrate on ν µ -ν τ flavour mixing mechanisms for which the propagation of neutrinos (+) and antineutrinos (−) is governed by the following Hamiltonian [16]: where ∆m 2 is the mass-squared difference between the two neutrino mass eigenstates, σ ± n accounts for a possible relative sign of the NP effects between neutrinos and antineutrinos show the spectrum as a function of the initial muon energy E 0 µ . and ∆δ n parametrizes the size of the NP terms. The matrices U θ and U ξn,±ηn are given by: U ξn,±ηn = cos ξ n sin ξ n e ±iηn − sin ξ n e ∓iηn cos ξ n ; (12) by η n we denote the possible non-vanishing relative phases. If NP strength is constant along the neutrino trajectory the oscillation probabilities take the form [16]: with sin 2 2Θ = 1 R 2 sin 2 2θ + R 2 n sin 2 2ξ n + 2R n sin 2θ sin 2ξ n cos η n , R = 1 + R 2 n + 2R n (cos 2θ cos 2ξ n + sin 2θ sin 2ξ n cos η n ) , where, for simplicity, we have assumed scenarios with one NP source characterized by a unique ∆δ n .
Eq. (11) describes, for example, flavour mixing due to new tensor-like interactions for which n = 1 leading to a contribution to the oscillation wavelength inversely proportional to the neutrino energy. This is the case for ν µ 's and ν τ 's of different masses in the presence of violation of the equivalence principle due to non-universal coupling of the neutrinos, γ 1 = γ 2 (ν 1 and ν 2 being related to ν µ and ν τ by a rotation ξ vep ), to the local gravitational potential φ [8,9,31] For constant potential φ, this mechanism is phenomenologically equivalent to the breakdown of Lorentz invariance resulting from different asymptotic values of the velocity of the neutrinos, c 1 = c 2 , with ν 1 and ν 2 being related to ν µ and ν τ by a rotation ξ vli [13,14].
For vector-like interactions, n = 0, the oscillation wavelength is energy-independent. This may arise, for instance, from a non-universal coupling of the neutrinos, k 1 = k 2 (ν 1 and ν 2 is related to the ν µ and ν τ by a rotation ξ Q ), to a space-time torsion field Q [12]. Violation of CPT resulting from Lorentz-violating effects such as the operator,ν α L b αβ µ γ µ ν β L , also leads to an energy independent contribution to the oscillation wavelength [15,16,17] which is a function of the eigenvalues of the Lorentz violating CPT-odd operator, b i , and the rotation angle, ξ CPT , between the corresponding neutrino eigenstates ν i and the flavour eigenstates ν α .
The flavour oscillations of atmospheric ν µ 's in this scenarios is described by Eq. (13) with the identification: for VEP (17) where for the first three scenarios σ + = σ − while for the CPT violating case σ + = −σ − . At present the strongest limits on NP neutrino oscillations arise from the non-observation of departure from the ∆m 2 oscillation behaviour in atmospheric neutrinos at SK and the confirmation of ν µ oscillations with the same oscillation parameters from K2K. In Eqs. (17)(18)(19)(20) we quote the 3σ bounds from the up-to-date combined analysis of SK and K2K data performed in Ref. [26].
For most of the neutrino energies considered here, the standard ∆m 2 oscillations are suppressed and the NP effect is directly observed. As a consequence, the results will be independent of the phase η n and we can chose the NP parameters in the range The Hamiltonian of Eq. (11) describes the coherent evolution of the ν µ -ν τ ensemble for any neutrino energy. High-energy neutrinos propagating in the Earth can also interact inelastically with the Earth matter either by charged current and neutral current (NC) and as a consequence the neutrino flux is attenuated. This attenuation is qualitatively and quantitatively different for ν τ 's and ν µ 's. Muon neutrinos are absorbed by CC interactions while tau neutrinos are regenerated because they produce a τ that decays into another tau neutrino before losing energy [32]. As a consequence, for each ν τ lost in CC interactions, another ν τ appears (degraded in energy) from the τ decay and the Earth never becomes opaque to ν ′ τ s. Furthermore, as pointed out in Ref. [33], a new secondary flux ofν µ 's is also generated in the leptonic decay τ → µν µ ν τ .
Attenuation and regeneration effects of incoherent neutrino fluxes can be consistently described by a set of coupled partial integro-differential cascade equations (see for example [34] and references therein). In this way, for example, the observed ν µ and oscillation-induced ν τ fluxes (and the associated event rates in a high energy neutrino telescope) from astrophysical sources has been evaluated. Alternatively, these effects can be accounted for in a Monte Carlo simulation of the neutrino propagation in matter [32,33,41]. Whatever the technique used, because of the long distance traveled by the neutrinos from the source, the oscillations average out and the neutrinos arriving at the Earth can be treated as an incoherent superposition of mass eigenstates.
For atmospheric neutrinos this is not the case because oscillation, attenuation, and regeneration effects occur simultaneously when the neutrino beam travels across the Earth's matter. For the phenomenological analysis of conventional neutrino oscillations this fact can be ignored because the neutrino energies covered by current experiments are low enough for attenuation and regeneration effects to be negligible. Especially for non-standard scenario oscillations, future experiments probe high-energy neutrinos for which the attenuation and regeneration effects have to be accounted for simultaneously. In order to do so, we modify the neutrino flavour oscillation equations and couple them to the τ evolution equations as described next.
We find it convenient to use the density matrix formalism to describe neutrino flavour oscillations. The evolution of the neutrino ensemble is determined by the Liouville equation for the density matrix where H is given by Eq. (11). The survival probability in Eq. (13) is given by P µµ (t) = Tr[Π νµ ρ(t)], where Π νµ = ν µ ⊗ ν µ is the ν µ state projector, and with initial condition ρ(0) = Π νµ . An equivalent equation can be written for the antineutrino density matrix. For the case of oscillations between two neutrino states the hermitian operators ρ, H and the flavour projectors Π νµ and Π ντ can be expanded in the basis formed by the unit matrix and the three Pauli matrices σ i . In particular we can write and the evolution of the neutrino ensemble is determined by a precession-like equation of the three-vector p(t) In this formalism attenuation effects due to CC and NC interactions can be introduced by relaxing the condition Tr(ρ) = 1. In this case and where we have explicitly exhibited the energy dependence and n T (x) is the number density of nucleons at the point x = ct. σ α CC (E) is the cross section for CC interaction, ν α + N → l α + X, and σ NC (E) is the cross section for ν α + N → ν α + X which is flavour independent. Thus we obtain four equations that describe the evolution of the neutrino system because one has to take into account both the flavour precession of the vector p(E, t) as well as the neutrino intensity attenuation encrypted in the evolution of p 0 (E, t).
ν τ regeneration and neutrino energy degradation can be accounted for by coupling these equations to the shower equations for the τ flux, F τ (E τ , t) (we denote by F the differential fluxes dφ/(dE d cos θ)). For convenience we define the neutrino flux density matrix is the initial neutrino flux. The equa-tions can be written as: λ τ dec (E τ , x) = γ τ c τ τ . τ τ is the τ lifetime and γ τ = E τ /m τ is its gamma factor. We have calculated the CC and NC distributions using the MRST-g parton distributions [42], and we have taken the τ decay distribution dN dec (Eτ ,Eν) dEν from Ref. [34] and dN dec (Ēτ ,Eν ) dEν from Ref. [43]. The third term in Eq. (29) represents the neutrino regeneration by NC interactions and the fourth term represents the contribution from ν τ regeneration, ν τ → τ − → ν τ , describing the energy degradation in the process. The secondary ν µ flux fromν τ regeneration,ν τ → τ + →ν τ µ + ν µ , is described by the last term where we denote by over-bar the energies and fluxes of the τ + . Br µ = 0.18 is the branching ratio for this decay. In Eq. (30) the first term gives the loss of taus due to decay and the last term gives the τ generation due to CC ν τ interactions. In writing these equations we have neglected the tau energy loss, which is only relevant at much higher energies.
An equivalent set of equations can be written for the antineutrino flux density matrix and the for the τ + flux. Both sets of equations are coupled due to the secondary neutrino flux term.
We solve this set of ten coupled evolution equations that describe propagation through the Earth numerically using the matter density profile of the Preliminary Reference Earth Model [44] and obtain the neutrino fluxes in the vicinity of the detector from In Fig. 4 we illustrate the interplay between the different terms in Eqs. (29) and (30). The figure covers the example of VLI-induced oscillations with δc/c = 10 −27 and max-imal ξ vli mixing. The upper panels show the final ν µ and ν τ fluxes for vertically upgoing neutrinos after traveling the full length of the Earth for the initial conditions The figure illustrates that the attenuation in the Earth suppresses the neutrino fluxes at higher energies. The effect of the attenuation in the absence of oscillations is given by the dotted thin line in the left panel. Even in the presence of ocillations this effect can be well described by an overall exponential suppression [39,43] both for ν µ 's and the oscillated ν τ 's. In other words, we closely reproduce the curve for "oscillation + attenuation" simply by multiplying the initial flux by the oscillation probability and an exponential damping factor: where X(θ) is the column density of the Earth.
The main effect of energy degradation by NC interactions (the third term in Eq. (29)) that is not accounted for in the approximation of Eq.(33) is the increase of the flux in the oscillation minima (the flux does not vanish in the minimum) because higher energy neutrinos end up with lower energy as a consequence of the NC interactions. The difference between the dash-dotted line and the dashed line is due to the interplay between the ν τ regeneration effect (fourth term in Eq. (29)) and the flavour oscillations. As a consequence of the first effect, we see in the right upper panel that the ν τ flux is enhanced because of the regeneration of higher energy ν τ 's, ν τ (E) → τ − → ν τ (E ′ < E), that originated from the oscillation of higher energies ν µ 's. In turn this excess of ν τ 's produces an excess of ν µ 's after oscillation which is seen as the difference between the dashed curve and the dash-dotted curve in the left upper panel. Finally the secondary effect ofν τ regeneration (last term in Eq. (29)),ν τ (E) → τ + → µ +ν τ ν µ (E ′ < E), results into the larger ν µ flux (seen in the left upper panel as the difference between the dashed and the thick full lines). This, in turn, gives an enhancement in the ν τ flux after oscillations as seen in the right upper panel.
The lower panels show the final ν µ and ν τ fluxes for an atmospheric-like energy spectrum In this case all regeneration effects are suppressed. Regeneration effects result in the degradation of the neutrino energy and the more steeply falling the neutrino energy spectrum, the smaller the contribution to the total flux. Therefore, in this case, the final fluxes can be relatively well described by the approximation in Eq.(33).

IV. EXAMPLE OF PHYSICS REACH: VLI-INDUCED OSCILLATIONS
Neutrino oscillations introduced by NP effects result in an energy dependent distortion of the zenith angle distribution of atmospheric muon events. We quantify this effect in IceCube by evaluating the expected angular and E f in µ distributions in the detector using Eq. (2) in  conjunction with ν µ (andν µ ) fluxes obtained after evolution in the Earth for different sets of NP oscillation parameters. Together with ν µ -induced muon events, oscillations also generate µ events from the CC interactions of the ν τ flux which reaches the detector producing a τ that subsequently decays as τ → µν µ ν τ and produces a µ in the detector. Using the techniques discussed in the previous sections we compute the number of ν τ -induced muon events as can be found in Ref. [43]. Equivalently we compute the number ofν τinduced muon events.
For illustration we concentrate on oscillations resulting from VLI that lead to an oscillation wavelength inversely proportional to the neutrino energy. The results can be directly applied to oscillations due to VEP.
We show in Fig. 5 the zenith angle distributions for muon induced events for different values of the VLI parameter δc/c and maximal mixing ξ vli = π/4 for different threshold energy E f in µ > E threshold normalized to the expectations for pure ∆m 2 oscillations. The full lines include both the ν µ -induced events (Eq.(2)) and ν τ -induced events (Eq.(35)) while the last ones are not included in the dashed curves. We see that for a given value of δc/c there is a range of energy for which the angular distortion is maximal. Above that energy, the oscillations average out and result in a constant suppression of the number of events. Inclusion of the ν τ -induced events events leads to an overall increase of the event rate but slightly reduces the angular distortion (see also Fig. 6) as a consequence of the "anti-oscillations" of the ν τ 's as compared to the ν µ 's.
In order to quantify the energy-dependent angular distortion we define the vertical-tohorizontal double ratio where by E f in,i µ we denote integration in an energy bin of width 0.2 log 10 (E f in,i µ ) using that IceCube measures energy to 20% in log 10 E for muons.
In what follows we will use the double ratio in Eq. (35) as the observable to determine the sensitivity of IceCube to NP-induced oscillations. We have chosen a double ratio to eliminate uncertainties associated with the overall normalization of the atmospheric fluxes at high energies.
Also, in the definition of the double ratio we have conservatively included only events well below the horizon cos θ < −0.2 to avoid the possible contamination from missreconstructed atmospheric muons which can still survive after level 2 cuts in the angular bins closer to the horizon [40].
In Fig. 6 we plot the expected value of this ratio for different values of δc/c. As mentioned above, IceCube measures energy to 20% in log 10 E for muons. Accordingly, we have divided   the data in 16 E f in µ bins: 15 bins between 10 2 and 10 5 GeV and one containing all events above 10 5 GeV. In the figure the full lines include both the ν µ -induced events (Eq.(2)) and ν τ -induced events (Eq.(35)) while the last ones are not included in the dashed curves. As described above, the net result of including the ν τ -induced events is a slight decrease of the maximum expected value of the double ratio. The data points in the figure show the expected statistical error corresponding to the observation of no NP effects in 10 years of IceCube.
In order to estimate the expected sensitivity we assume that no NP effect is observed and define a simple χ 2 function as where σ stat,i is computed from the expected number of events in the absence of NP effects (see Table I).
We show in Fig. 7 the sensitivity limits in the [δc/c, ξ vli ]-plane at 90, 95, 99 and 3 σ CL obtained from the condition χ 2 (δc/c, ξ vli ) < χ 2 max (CL, 2dof). In order to estimate the uncertainty associated with the poorly known prompt neutrino fluxes we show in the figure the results obtained using the RQPM model (filled regions) and the TIG model (full lines).    The difference is about 50% in the strongest bound on δc/c. The figure illustrates the improvement on the present bounds by more than two orders of magnitude even within the context of this very conservative analysis. The loss of sensitivity at large δc/c is a consequence of the use of a double ratio as an observable. Such an observable is insensitive to NP effects if δc/c is large enough for the oscillations to be always averaged leading only to an overall suppression.
When data becomes available a more realistic analysis is likely to lead to a further improvement of the sensitivity.

V. SUMMARY
In this paper we have investigated the physics that can be probed with the high-statistics high-energy data on atmospheric neutrinos which will be collected by the IceCube detector. In order to do so first we have developed a semianalytical simulation of the detector performance. In particular, we present in Eqs. for the response of the IceCube detector after events that are not neutrinos have been rejected using the quality cuts referred to as level 2 cuts. We conclude that in 10 years of operation IceCube will collect more than 700 thousand atmospheric neutrino events with energies E f in µ > 100 GeV which offer a unique opportunity to test new physics mechanisms for leptonic flavour mixing which are not suppressed at high energy. In general these effects are expected to induce an energy dependent angular distortion of the events.
Next, because of the relatively high energy of the neutrino sample, NP induced flavour oscillations, propagation in the Earth, regeneration of neutrinos due to τ decay must be treated in a consistent way. In Sec. III we have presented the corresponding evolution equations. We conclude that for steeply falling neutrino energy spectra, such as the atmospheric neutrino one, the dominant effect together with flavour oscillations is the attenuation of the oscillation amplitude due to inelastic CC and NC interactions of the neutrinos in the Earth in conjunction with the production ν τ -induced muon events due to the chain ν τ → τ → µν µ ν τ in the vicinity of the detector. ν τ -induced muon events can increase the event sample by at most O(10%).
Finally we have applied these results to realistically evaluate the reach of IceCube in studying physics beyond conventional neutrino oscillations induced by violation of Lorentz invariance and/or the equivalence principle. In Fig. 7 we show how even with a very conservative analysis the range of testable sizes of these effects can be easily extended by more than two orders of magnitude. The methods developed are readily applicable to probe speculations on other non-conventional physics associated with neutrinos.