Neutrino Detection with Inclined Air Showers

The possibilities of detecting high energy neutrinos through inclined showers produced in the atmosphere are addressed with an emphasis on the detection of air showers by arrays of particle detectors. Rates of inclined showers produced by both down-going neutrino interactions and by up-coming $\tau$ decays from earth-skimming neutrinos as a function of shower energy are calculated with analytical methods using two sample neutrino fluxes with different spectral indices. The relative contributions from different flavors and charged, neutral current and resonant interactions are compared for down-going neutrinos interacting in the atmosphere. No detailed description of detectors is attempted but rough energy thresholds are implemented to establish the ranges of energies which are more suitable for neutrino detection through inclined showers. Down-going and up-coming rates are compared.


Introduction
The observation of Ultra High Energy Cosmic Rays (UHE CR) of energy ∼ 10 20 eV [1,2,3,4] has stimulated much theoretical and experimental activity in the field of Astroparticle Physics. These cosmic rays are are expected to be suppressed by interactions with the Cosmic Microwave Background (CMB) radiation, what is known as the Greisen Zatsepin Kuzmin (GZK) cutoff [5,6]. Although many mysteries remain to be solved [7], we know it is virtually impossible to produce these energetic particles without associated fluxes of gamma rays and neutrinos from pion decays. The gamma rays pair produce in the background photon fields and the electrons loose energy through synchrotron emission in a "cosmic cascading" process which degrades their energy down to the 100-MeV region. Neutrinos travel undisturbed carrying a "footprint" of the production model. It is remarkable that the diffuse gamma rays in the 100-MeV region, UHE neutrinos and UHE CR are deeply related. Their precise measurements, crucial for unraveling the mysterious origin of the highest energy particles in the Universe, must be considered part of the same priority.
High energy Cosmic Ray production mechanisms are basically of two types, namely acceleration models and "top-down" scenarios [8]. Neutrino fluxes in the EeV range are difficult to avoid but their fluxes are uncertain [9]. In models that accelerate protons and nuclei, pions are believed to be produced by cosmic ray interactions with matter or radiation at the source. In "top-down" scenarios protons and neutrons are produced from quark and gluon fragmentation, a mechanism which is known to produce about 30 times more pions than nucleons. Furthermore protons and nuclei also produce pions in their unavoidable interactions responsible for the GZK cutoff which decay to produce the cosmological neutrinos [10]. Cosmological neutrinos could dominate in some acceleration models but in top-down scenarios they should be well below the neutrinos from the decays of the fragmented pions. The UHE proton to neutrino flux ratio is thought to carry important information concerning the origin of UHE CR.
High energy neutrinos produce extensive air showers and can be detected by the very same detectors that measure the cosmic ray spectrum. The main challenge lies in separating showers initiated by neutrinos by those induced by regular cosmic rays. It was already suggested in the 1960's that this could be done at high zenith angles [11] because the atmosphere slant depth is quite large. Neutrinos, having very small cross sections, can interact at any point along their trajectories while most of the air showers induced by protons, nuclei or photons are absorbed before reaching ground level. The signature for neutrino events is thus inclined showers that interact deep in the atmosphere.
Inclined showers were first observed in the 60's by several groups [12,13]. The observed rate for shower sizes between 10 3 -10 5 particles detected in the Tokyo (INS) array [12] is however consistent with hard processes induced by energetic muons that interact deep into the atmosphere [14,15,16,17]. For shower sizes above 10 5 particles, AKENO has published an upper bound on muon poor air showers at zenith angles greater than 60 0 [18]. These observations have been very useful for establishing bounds on models for high energy diffuse neutrino fluxes at the Earth [19,20]. Inclined showers were also detected in Haverah Park (HP), a 12 km 2 air shower array in the UK. The largest events detected by Haverah Park were actually coming at very large zenith angles [13]. The inclined data set was recently studied and shown to be consistent with hadronic origin [21]. When the Auger project was conceived as the largest and most accurate air shower detector [22], it became clear that it would have a competitive acceptance for inclined showers induced by neutrinos compared to devoted neutrino experiments in planning or construction [23,24]. Other experiments have been proposed to search for deep inclined showers in the atmosphere from satellites at orbital altitudes that allow large volumes of atmosphere to be monitored to increase the chances of detecting neutrinos [25,26] or from behind mountains [27,28,29].
Although neutrinos of τ flavor are heavily suppressed at production, neutrino flavor oscillations [30], a maximal θ 23 mixing and a 2:1 ratio of ν µ : ν e at production, lead to approximately equal fluxes of all three flavors after propagation over cosmological distances [31]. When the ν τ undergoes a charged current interaction it produces a τ which, having a short lifetime, typically decays in flight. The process can induce two showers separated by a distance gap that is on average proportional to the τ energy. These events were soon discussed in the context of deep underground neutrino telescopes and referred to as "double bang events" [32]. More recently double bang events have been explored for higher energy neutrinos in the atmosphere [33]. As a ν τ propagates through the Earth, charged current interactions followed by τ decays (which always include a lower energy ν τ ) effectively shift the neutrino energy "regenerating" the ν τ flux [34]. For the higher energy neutrinos these processes can occur several times. For neutrinos traversing most of the Earth the emerging flux has a characteristic "pile-up" at an energy region around the PeV scale, except for those entering with high zenith angles.
It was remarked in 1999 [35] that the perspectives of detecting ν τ fluxes can increase dramatically for neutrinos that enter the Earth surface just below the horizon traversing a relatively small earth matter depth. These neutrinos can undergo a charged current interaction to produce a τ that actually exits the Earth and decays in the atmosphere producing an upcoming air shower. These "Earth skimming" neutrinos can be very effective in producing τ 's because of the long τ range at high energies (∼ 10 km) which sets the scale of the effective volume. This has turned into a promising field that has been studied by a number of authors [36,37,38,39,40,41,42,43,29]. Monte Carlo calculations have been performed for the Surface Array of the Auger detector [37], for mountain ranges [29] and in general for atmospheric τ production rates in general [39], and there are analytical calculations for the Fluorescence part of the Auger experiment [43] and for HiRes and the Telescope Array [38,40]. In this article we calculate both down-going and upcoming rates of neutrino induced air showers relevant for arrays of particle detectors using analytical methods. Comparative studies made by simulation have also been discussed by in the context of shower detection from satellites [39]. Although the detailed rates in an air shower array are expected to be very dependent on detector technicalities we examine the energy spectra of the induced showers to compare the down-going neutrino event rate to the upcoming rate induced by earth skimming neutrinos.

EeV Neutrino Detection: Generalities
Above 1 PeV the Earth becomes opaque to neutrinos and only down-going or earth skimming EeV neutrinos can be detected. The challenge lies in the identification of these showers in the background of down-going cosmic rays and atmospheric muons. Inclined showers in the atmosphere are expected to play a crucial role for the detection of EeV neutrinos.
The main background for the detection of inclined showers produced by neutrinos is due to showers induced by protons and nuclei. These are expected to develop high in the atmosphere so that when the shower front reaches ground level it has very different properties from "ordinary" vertical showers that are observed in the vertical direction. Deep inclined showers induced by neutrinos can develop close to ground level so that their shower front resemble that of a typical vertical proton cosmic ray shower. There is a second type of background to inclined showers induced by neutrinos which is due to deep showers induced otherwise, mostly through hard muon bremsstrahlung. Unlike those produced cosmic rays these cannot be distinguished from down-going neutrino induced showers on the basis of the muon interaction point because they can also be deep. It is difficult to believe this background could be separately identified, alternatively it is thought that it sets the detectability limits for neutrino interactions.
The success in the search for inclined showers induced by neutrinos depends crucially on obtaining reliable information on shower properties related to the depth of the first interaction. Fluorescence detectors actually register the depth development of the air shower and the depth of shower maximum, which is closely related to it, can be directly obtained. For these detectors upcoming showers, such as typical neutrino earth skimming events, can be in principle easily identified. In these respects the fluorescence technique is most reliable and only good reconstruction accuracy is required for neutrino identification. Unfortunately background light introduces a duty cycle which limits its acceptance both for cosmic ray and neutrino detection. Air shower arrays measure the particle densities as the shower front reaches ground level and thus the inference of the depth of the first interaction is a far more indirect measurement. For close to horizontal events it is actually not possible to distinguish upcoming from down-going events from the arrival times of the signals, so that other shower characteristics must be used to identify potential neutrino induced events. The study of the difference between deep and shallow inclined showers has been the subject of much activity in the recent past. It is now believed that there are a number of properties of the shower front that can be easily used to distinguish neutrino induced showers, mostly stemming from the time distribution of the shower particles. We discuss briefly the basis of these differences. However the detailed calculation of the efficiency of a given detector to separate deep showers from those induced by protons and nuclei is quite a technical issue which lies beyond the scope of this article.
The neutrino cross sections have a direct impact on the expected shower rates. Unfortunately the cross section has to be calculated using quite significant extrapolations of data obtained in accelerators and there is uncertainty involved. At the end of this section we briefly review the basics of the standard model neutrino interactions and give the details of the cross sections used for completeness. The effects of non-standard model behaviors of the cross section can be quite dramatic in the shower rates. These have been discussed several times in the literature [44] and will not be further addressed here.

Inclined Showers from Cosmic Rays
Atmospheric air showers are produced by high energy cosmic rays that interact soon after entering the upper part of the atmosphere and they are regularly registered by air shower detectors [45,46,47]. The atmosphere has just the adequate matter depth (∼ 1000 g cm −2 ) so that vertical showers fully develop not too far above ground level, where the shower front contains large numbers of gamma rays, electrons and positrons (the electromagnetic component) as well as muons in rough proportions 100:6:4:1. As the arrival direction of the cosmic ray particles increases in zenith angle, θ, the slant depth to ground level rises, in proportion to sec θ for low zenith angles and with a modified behavior for θ > 60 • because of the Earth's curvature. The maximum slant depth for a completely horizontal (60 • ) shower at sea level is ∼ 36000 (2000) g cm −2 , see Fig. 1, about 36 (2) times larger than for vertical showers. As a result most of the electromagnetic component of showers with θ > 60 • is much absorbed before reaching ground level.
The analysis of the Haverah Park inclined data set established the cosmic ray EeV background for neutrino detection [49]. A fairly general model to describe the inclined showers induced by protons, nuclei and photons [48] was developed which allowed the first analysis of EeV data at zenith angles between 60 • and 85 • on on an event by event basis. The observed rate was shown to be consistent with a baryonic origin of the cosmic rays and a limit on photon abundance above 10 19 eV was set [49,50,21]. The electromagnetic part of cosmic ray showers at high zenith develops and gets practically absorbed in the first 2000 g cm −2 and to a very good approximation only muons reach the ground. The lower energy muons actually decay in flight and can contribute a small electromagnetic component that follows closely that of the muons. As a result the average energy of the muons reaching ground increases rapidly as the zenith angle rises. The magnetic field of the Earth acts as a spectrometer for the surviving muons separating the negative from the positive in inverse proportion to their energy and producing complex density patterns on the ground [48]. The shower front that reaches ground level for inclined showers is very different from vertical showers. Inclined shower fronts practically only contain energetic muons and their density patterns have lost the , Figure 1. Left: Slant depth as a function of zenith angle for three different altitudes: from top to bottom sea level, 1300 and 3000 m. Right: Depth development of electromagnetic (red) and muonic (black) components of a typical shower induced by a proton (qualitatively identical to that induced by a nucleus).
cylindrical symmetry because of the Earth's magnetic field. The degree to which this symmetry is broken depends on the relative strength and orientation of the magnetic field with respect to shower axis as well as the characteristic distance traveled by the muons which is governed by the zenith angle [48].
Studies have been made of the time structure of the shower front induced by inclined cosmic rays and deep showers [52]. While the former are very close to a plane front deep showers display a front with significant curvature, in similar fashion to ordinary vertical showers. This can be characterized by an effective radius of curvature that is directly related to the distance at which shower maximum is achieved. The time structure of the arriving particles has been shown to be very different because it is due to muons that are produced very far away. The muons do not have much multiple elastic scattering and, as a result, their arrival time distribution can be related to a very good approximation only to their geometrical path and the kinematic delays associated with the muon sub-luminal velocity [53]. At large zenith ordinary cosmic rays produce the shower muons very far from ground level and both the front curvature and its width are very small. On the contrary vertical and deep inclined showers, with large numbers of electrons, positrons and photons, display both a small radius of curvature because shower maximum is relatively close as well as a a large signal dispersion in time which increases as the distance to shower axis increases. The time structure of the showers detected with the Auger observatory have shown clearly these effects [54,55].
The inclined shower technique must exploit these differences to search for neutrino induced inclined air showers. A simple signature can be used: inclined showers that look like ordinary vertical showers with a significant electromagnetic component. More sophisticated analysis methods can enhance the performance of the detectors to search for deep showers for instance relating the time structure of the muons to the depth development of the shower [56]. It has been recently shown that it is possible to establish the depth of production of the muons in a shower from the time distributions of the muons which are far from shower axis [57]. If this technique can be effectively exploited it could allow the selection of deep showers also for showers that have had their electromagnetic component absorbed.

Inclined Showers Produced by Hard Muons
High energy muons can produce deep air showers by bremsstrahlung, pair production and nuclear interactions. Although these processes typically contribute similar amounts to muon energy loss, bremsstrahlung is the hardest process and when folded with a steeply falling atmospheric muon spectrum, is the most important for producing high energy showers. This is fact the dominant mechanism for explaining the observed inclined shower rates in the energy range between 100 TeV and 10 PeV [16]. The high energy muons that are mostly traveling along the shower axis must undergo the bremsstrahlung interaction in which a large fraction of the muon energy is transferred to the photon. This must be produced at about 400 g cm −2 of slant depth before the muon reaches ground level, in a relatively narrow range of depth so that the shower reaches its maximum close to ground level. Both bremsstrahlung and pair production will give rise to pure electromagnetic showers that would hardly contain any muons, but the hadronic interactions of the muons will produce showers that are of hadronic type. Since atmospheric muons are produced in air showers these energetic muons must come from extensive air showers which exceed the muon energy by some factor, which is, incidentally, dependent on the nature of the primary particle. The hard process itself generates a sub-shower of a fraction of the muon energy which develops along shower axis. It is even possible that several of these hard muon processes combine in a single larger shower. The softer muons, which are spread laterally in patterns which are governed by the geomagnetic field, could in principle also be detected if the primary shower energy is sufficiently high.
The deep showers induced by these hard processes is thus "embedded" in a larger shower which may or may not be detectable. If they are detected it could in principle be possible to distinguish the deep sub-shower within the cosmic ray shower which would clearly signal this kind of processes but this seems a rather difficult task. If not, it would be impossible to distinguish from a charged current electron neutrino shower for instance. Fortunately the atmospheric muon flux is very soft and (at sufficiently high energy) the expected rate of showers from hard muon processes is expected to be very small. However the muon bremsstrahlung rate is subject to uncertainties, mainly due to the production of prompt muons through charmed mesons, a subject that has not been resolved [58]. Inclined showers of shower size between 10 5 and 10 7 are expected to give valuable information concerning prompt charm production [59]. Backgrounds of hard muon processes for neutrino detection in the context of the Auger detector have been discussed in Ref. [60].

The Neutrino Cross Sections
The dominant interaction for neutrinos in matter is Deep Inelastic Scattering (DIS) on nucleons (ν + N → l + X, where N stands for a nucleon and l for a lepton). The cross sections are given in terms of the standard structure functions. Within the standard model propagator effects associated with the mass of the weak exchange bosons, M B , become important at high energies: Here y is the fraction of energy transferred to the nucleon in the laboratory frame, Q 2 is minus the square of the 4-momentum transfer, m p is the proton mass and x is defined by the relation Q 2 = 2m p Exy. For charged (neutral) current interactions the boson mass is taken to be M B = M W (M Z ), the W ± (Z 0 ) boson mass. The ± sign in Eq. 1 implies a + sign for ν and a − sign forν. F 1,2,3 are the structure functions which can be expressed in terms of universal parton distributions of the nucleons, obtained from fits to accelerator data and QCD evolution equations. For charged current interactions on isoscalar targets at leading order (LO) they are given by: where x is to be interpreted as the momentum fraction carried by the quark distributions (u, d, s, c, b, t, their dependence on x and Q 2 has been omitted for clarity). For neutral currents the structure functions involve different combinations of the parton distributions with neutral coupling factors that depend on z = sin 2 θ W , where θ W is the electroweak mixing angle: which together with the Callan-Gross relation (Eq.6) apply to both ν and ν.
Unfortunately the neutrino cross section at the energies of interest must be deduced from extrapolations of the parton distributions to high Q 2 and low x values, where there is no accelerator data. As accelerator experiments provide more accurate data on parton distributions, new sets of parameterizations are developed which lead to somewhat different cross section calculations [61,62,63,64]. For energies above 10 6 GeV charged current interactions are about 2.5 times more likely that neutral current interactions. Uncertainties can represent up to a factor of two for the total cross section at energies around 10 20 eV as estimated in [65]. Next to leading order [64] and the atomic number of the target [66] have been shown to imply changes in the cross sections which are smaller than the uncertainties due to the extrapolations of the structure functions. In the work that follows here we will use one of the most recently updated sets, CTEQ6 [67] for reference.
This adds an additional interest to neutrino detection because neutrino experiments may provide new insights into particle physics. Besides the obvious role of the neutrino cross section in the actual interaction that leads to the possible neutrino detection there are other more subtle effects. The y distribution of the cross section also has an impact on the detection rate. The average value obtained with standard parameterizations is < y >≃ 0.2 for energies above 10 6 GeV. Since the average energy transfer to the nucleon is very dependent on the low x behavior of the parton distributions in a region which is not yet accessible to accelerator experiments [68], it has been pointed out that simultaneous measurements of different neutrino channels could provide insight on the low-x behavior of the parton distributions [69]. For down-going neutrino induced air showers the interaction rate is proportional to the cross section, for Earth skimming upcoming showers the cross section acts in the opposite way increasing the neutrino absorption in the Earth. A very interesting idea has been discussed that a simultaneous measurement of the two rates could lead to an indirect measurement of the neutrino cross section at energies well beyond accelerator range [70]. The relative rates of down-going and upcoming showers are thus of great relevance from this perspective too.
Finally we also give an approximate cross section for the Glashow resonance, e ν e → W − → l l ′ considering only the first order s-channel diagram and neglecting mass corrections: where here y is to be taken an the fraction of energy carried by the most negative charged lepton, and Γ W is the W boson width. The average fraction of energy transferred to the most negative lepton is < y >= 0.25. In practice electron cross sections are suppressed by the the ratio of the electron to proton mass and this interaction is only relevant for ν e energy of around 6.4 10 6 GeV, where it dominates over DIS interactions. This expression is only adequate for the resonance region. The reactions are explicitely shown in Table 1.

Down-going Neutrino Induced Air Showers
There are multiple channels for producing down-going neutrino induced showers. In order of importance first are charged current interactions with atmospheric nuclei and then neutral current interactions followed by ν GeV and if it decays into a τν τ the shower is produced by τ decay (64% of the times into hadrons and 18% into eν e ). The calculations of the shower rates induced by a differential neutrino flux (per unit energy E ν , area, A, solid angle, Ω and time, t): are based on Ref. [16]. We will used a simple neutrino flux for with a constant spectral index of 2 and compatible in the 10 6 − 10 8 GeV region with the Gamma Ray Burst (GRB) flux given in Ref. [71] and one from a characteristic top down scenario (TD) . This last model has been chosen because it is similar to that used in Ref. [37] and it is consistent with experimental data [72,73]. They are illustrated in Fig. 2. The choices are driven by their different energy dependences in the most interesting region between 10 16 − 10 18 eV and for comparison with other calculations. The total ν flux is obtained multiplying the ν µ andν µ flux by a factor 1.5 to approximately account for ν e production according to naive counting. The basic expression for the calculation of the number of showers, dN sh , induced by this flux interacting in the interval of atmospheric matter depth dx, is simply given by: where N A is Avogadro's number. Depending on the interaction the differential rates in shower energy have slightly different expressions involving the cross sections discussed in Section 2.3. In charged current ν e interactions all the interaction energy is transferred (mass of X particle M X = 3 10 14 GeV, evolution parameter p = 2 and injection normalization Q 0 = 10 −34 erg cm −3 s −1 ) to the shower: E ν = E sh and the differential shower rate produced is: where x max (θ) is the slant depth available in the atmosphere for the neutrino to interact.
Unfortunately not all the showers induced by neutrinos can be detected and distinguished from showers produced by cosmic rays. The shower rate that can be identified to be produced by a neutrino depends crucially on detector performance and capabilities. We can establish a probability for detection P det , a function that necessarily depends on shower energy, arrival direction, dΩ, depth at which the shower is produced, x, fraction of energy transfer, y, and impact parameter, l ⊥ . When this probability is included the shower rate becomes: The above expression can be integrated both in solid angle and in impact parameter over the surface detector of area A d to give: Here P mix det refers to the subsequent shower of mixed electromagnetic and hadronic components respectively carrying fractions (1 − y) and y of the incident neutrino energy. If edge effects are ignored the d 2 l ⊥ integral simply factorizes as A d and introduces a "cos θ" factor in the remaining integral. The acceptance (as expressed between the large brackets) for the surface detector of the Pierre Auger Observatory was calculated in Ref. [24]. Expressions for both purely hadronic and purely electromagnetic showers were obtained integrating in solid angle for θ ≥ 75 • and demanding that the shower intercepts ground level with a significant electromagnetic component.
An important contribution to the inclined shower rate induced by neutrinos is due to neutral current interactions. In these processes the shower energy is simply E sh = yE ν . The result above is modified in the argument of the neutrino flux which takes a dependence on y, in the y integral that must include a y −1 factor from the change of variables and on P had det which now is only for hadronic type showers: The same expression accounts for charged current ν µ interactions (replacing σ nc by σ cc ) because the produced µ is unlikely to decay in the whole atmospheric depth unless its energy is below 1 TeV. In the case of ν τ charged current interactions, the emerging τ is likely to decay provided that its energy is below ∼ 10 9 GeV [33]. Since a radiation length in the atmosphere is of order 300 m, for τ energies below about 10 7 GeV the decay will merge with that from the hadronic debris and a result similar to Eq. 14 can be expected. For the τ decay shower the fraction of τ energy that goes into the shower after its decay (f ) must also be accounted for. Assuming a constant fraction, E sh = yf E ν , the expression becomes: with P τ det modified according to the shower details. For larger energies the double bang nature of the process will be manifest and each of the two showers could be separately detected.
The calculation proceeds in analogous ways for resonant W − production with σ cc replaced by σ res in the corresponding expressions. When the W boson decays into qq pairs it can be assumed that all the energy is channeled after fragmentation into a hadronic shower and Eq. 14 applies with a hadronic detection probability function P had det that does nor depend in y. As a result the integral in y can be trivially performed, what amounts to substituting the integral in y for the corresponding total cross section, σ res . When a eν e pair is produced Eq. 15 must be used provided that y (for the resonant differential cross section) is taken to be the fraction of the neutrino energy carried by the charged lepton in the laboratory frame. When the W decays into τ ν τ pair and the shower is produced by τ decay Eq. 16 applies.

Relative Contributions of different Channels
The precise form of the detection probability P det is crucial in the calculation of the shower rates. This calculation is complex and dependent on technical details of the detector response that must include the capabilities to identify showers produced by Table 1. Contributing channels to down-going shower rates indicating the neutrino reactions, the corresponding neutrino flavors, the normalization factors, the shower type, its fraction of the neutrino energy and the relevant expression applied.

Reaction
Type ν flavor Norm Fact Shower Type Energy Frac Eq.
ν + N → l ± X Chargedν e ν e 1. mixed To establish some level of comparison with upcoming showering events from Earth skimming neutrinos, we make a very simple assumption to evaluate it: provided the zenith angle is greater than 60 • all showers that start developing in the second half of the slant depth of the atmosphere can in principle be distinguished from those produced by cosmic rays. This is not necessarily a conservative assumption. For air shower arrays the particle detectors must be capable of resolving the arrival time and time spread of the signals and this is necessarily challenging near detector threshold. However by demanding that the neutrino travels through at least the first ∼ 1000 g cm −2 without interacting, because of geometrical considerations it can be assumed that the subsequent shower must have different properties (such as radius of curvature) from a cosmic shower of the same zenith. ¶ We then require that showers at ground level have more than a fixed number of electrons and positrons (N e ) to be detected. This requirement will determine the depths at which the neutrino can interact to give a significant signal. We make several choices of N e = 10 5 , 10 6 and 10 7 to simulate different detector responses naively. These are arbitrary (a shower of 10 15 eV has about 10 6 particles at shower maximum) and set a fairly abrupt low energy cutoff for the calculations that follow. By choosing deliberately low values the limits of the technique are explored by examining its potential low energy behavior. In a realistic calculation this is very dependent on the detector, on the geometry, on the lateral distribution of the different particles in the shower front and on many other aspects which have been ignored. However as the energy of the shower rises, the range of depths in which the neutrino can interact to be detected increases to eventually reach half the atmospheric depth. For high energies, when any detector is most likely to detect deep showers, we do not expect this simple calculation to be very ¶ We remark here that the atmosphere is split in halves in depth. This implies that for a detector 1.3 km above sea level and a zenith angle of 60 • the midpoint is 10.5 km away while the typical first interaction point for a proton or nucleus is ∼ 40 km away. A ν shower that develops in the second half is closer to the detector by more than a factor of two. far off.
The calculations have been made using the above described methods together with parameterizations for the depth development of hadronic [74] and of electromagnetic showers [75] to calculate x max . All neutrino flavors are assumed to have the same flux when reaching the Earth after mixing. In Table 1 different channels are summarized, together with the resulting normalization factors (relative to each neutrino flavor), the type of shower generated, the energy fraction transferred to the shower and the appropriate expression in which P has to be chosen according to the shower type. The grouping results from the combinations of the neutrino flavor(s) involved, the shower type, the fraction of energy transferred to the shower and the cross section. The normalization factors are associated to the number of neutrino flavors, and in the case of qq production a factor of six is included to account for the sum over all final state possibilities, (two weak doublets and three colors). , Figure 3. Shower rates produced by the main Deep Inelastic Scattering processes (full) and resonantν e e interactions (dashed) for the GRB model (left) and the TD model (right). DIS interactions are separated in charged current ν e +ν e , charged current ν µ +ν µ + ν τ +ν τ (neglecting τ decay) and neutral current interactions for all neutrino flavors. The qq and eν e channels for the Glashow resonance are shown.
The results of the main contributions are compared in Fig. 3. We have included all the neutral interactions, the ν e ,ν e charged current interaction and the ν µ ,ν µ , the ν τ ,ν τ charged current interactions neglecting τ decay, (a detailed calculation of this effect will be presented elsewhere [76]) as well as the resonant qq and eν e channels. Basically for DIS interactions the calculated shower rate follows closely that of the incident neutrino flux. The relative contributions of each channel depend on its spectral index. This is not surprising since the main difference between channels is the fraction of energy that goes to the shower. If the ν spectrum is steep the charged current electron neutrino interactions dominate because all the neutrino energy goes to the shower, but for hard fluxes, like the topological defect models discussed here, the shift in shower energy for neutral current channels is compensated by the increase in flavors and hence this channel dominates.
The total shower rates produced by the discussed DIS and resonant channels are plotted for the two fluxes chosen in Fig. 4, together with the dominant DIS channel for each case. In this figure we also compare the effects of different thresholds as discussed above. The figure illustrates that the overall shower rate flux is close to a factor of 3 times above that calculated using the dominant neutrino channel only (the precise factor depends on the spectral shape of the flux in question). The figure illustrates the dramatic effect that the low energy behavior of an air shower array can have in the total event rate calculation. Naturally this is particularly significant for the GRB model which is much steeper. , Figure 4. Down-going shower rates produced by all channels and flavors added up, compared to charged current ν e +ν e for the GRB model (left) and to charged current ν µ +ν µ + ν τ +ν τ for the TD model (right). Detector efficiency is roughly estimated requiring shower size (ρ e ) at ground level to be above nominal threshold values.

Earth Skimming Neutrinos
As neutrinos go through the Earth they can undergo charged current interactions close to their exit point below the surface of the Earth. Under certain circumstances these interactions can produce extensive air showers that develop upwards in the atmosphere. Only those produced by ν τ 's are of significance from the quantitative point of view, because they lead to large detectable shower rates in the high energy range of interest. This is due to the combined effect of the lepton lifetime, its energy loss and the neutrino cross section. For ν τ the mechanism can be described in three stages, namely the ν τ enters the Earth and propagates through it, it then has a charged current interaction just below the surface producing a τ that continues the propagation in matter with considerable energy loss and finally the τ exits the surface and decays in the atmosphere inducing an air shower.

Propagation in the Earth
Neutrinos have a cross section which increases as the neutrino energy rises. As a flux of neutrinos goes through a given amount of matter a characteristic energy scale is selected, the energy at which the neutrinos are likely to interact. Above this energy the flux is attenuated while below it is practically undisturbed. Both charged and neutral current interactions contribute to neutrino attenuation. For a neutrino crossing the Earth the amount of matter depth is a rapidly varying function depending on impact parameter. If they travel through the center, that is they enter with a relatively low zenith angle, the matter depth is not much below 6.6 10 9 g cm −2 (the depth of a diameter). Using the cross section discussed in subsection 2.3 this corresponds to the attenuation length of a 100 TeV neutrino. As the impact point gets further away from the Earth's center, and both the incidence (zenith) angle and exit (nadir) angle (which are the same) approach 90 • , the matter depth decreases to zero. At an incidence zenith of 89 • the depth has reduced by over a factor 100 to a value (which depends on local density) typically between ∼ 2 − 5 10 7 g cm −2 . Since the cross section rises with energy approximately as E 1 3 in this range, the characteristic energy scale, at which the neutrinos become attenuated changes from 1 PeV to about 10 19 − 10 20 eV at 89 • . The Earth thus filters the neutrino spectrum in a complex way suppressing the high energy part in a rapidly changing way as the zenith angle approaches the horizontal.
Besides attenuation there is a regeneration effect. Neutral current interactions shift the neutrino energy, transferring part of the absorbed flux into lower energies. In addition, the τ (µ) lepton produced in charged current interactions can in turn decay producing ν τ (ν µ ). This double sequence also regenerates the neutrino flux [34]. Typically the absorbed high energy ν τ flux gives a τ which only decays when it is below a characteristic energy scale ∼ 10 8 GeV (this is described in subsection 4.2). Regeneration is important below the smaller of this scale and the energy scale fixed by the absorption. For low zenith the regeneration is important in the 1 PeV region [34] but for Earth skimming neutrinos regeneration is most important at about 2 10 7 GeV as shown by simulations [77]. These neutrinos have to interact produce a τ , exit the Earth and decay in the atmosphere so that the shower energy is further reduced. Since the electron is stable there is no such effect for ν e while for ν µ the long µ lifetime results in regeneration at a low energy which is not relevant for this work. We will here consider neutral current regeneration but will not further discuss the ν τ regeneration. Earth skimming rates obtained particularly for the topological defect flux will be somewhat higher at around 10 7 GeV.

Lepton Propagation: Energy loss and Decay
Once a ν µ or a ν τ undergoes a charged current interaction, the produced lepton has a probability to decay which increases as it loses energy while propagating in matter. The propagation of the lepton in the Earth determines the effective volume that is available for detection in this channel. Following procedures similar to those described in [43,53] we calculate the survival probabilities for leptons with continuous energy losses. The survival probability has the following differential equation: where l is the distance traveled, βc the lepton velocity, m its mass, τ its decay lifetime and E(l) its energy, which decreases as l increases. The energy loss can be expressed by an approximate expression: in which the energy loss per unit depth (x) has a constant energy loss term associated with ionization (a) and a linear term in energy which is due to hard processes, bremsstrahlung, pair production and hadronic interactions. This is usually expressed in terms of a characteristic grammage, ξ = 1/b, and a critical energy, ǫ = aξ above which radiative processes dominate. These values are slightly varying with energy and somewhat uncertain but they have implications for the τ production rate [37]. We here use the values ξ = 1.25 10 6 g cm −2 and ǫ = 3000 GeV consistent with results for the 10 8 GeV region in Ref. [78]. This equation can be integrated to give [79]: which relates E ′ , the initial lepton of energy, to E, the lepton energy after propagation over depth x. This expression can be used to express E in Eq. 17 in terms of l, assuming the lepton propagates at the speed of light in a constant density medium. With this substitution Eq.17 can be integrated in x to obtain: where the second approximate expression holds for E >> ǫ and the exponent, η, is a constant that depends on the medium density ρ, on the lepton mass, m l , its decay constant, τ l , and the loss parameters ξ and ǫ: The value of η for muons in rock of density 2.65 g cm −3 is ∼ 3.5 10 −4 which implies that decay can be ignored for all cases of interest here (the muon lifetime is very long compared to energy loss time). In the atmosphere η ∼ 0.66 and the muon can decay in flight. However for the energies of interest here the muon exits the atmosphere well before decaying.
For the τ lepton propagating in standard rock η ≃ 3.2 10 4 is very large. This implies that unless the ratio of energies in Eq. 20 is very close to one the τ does not survive. We can then approximate Eq. 20 as: For the τ not to decay the absolute value of the exponent must be small and the fraction of energy lost by the τ is limited: For E is below ǫ ionization losses dominate and the fraction of energy loss must be below a very small value. When E >> ǫ this bound increases linearly with E, the energy of the τ after propagation, until it becomes 1 when E ≃ ǫη. This introduces a new energy scale for the emerging τ : ǫη ≃ 9 10 7 GeV.
, Figure 5. Comparison of the τ range and the ν τ absorption length, using two different neutrino cross sections (CTEQ4 and CTEQ6) and and two different sets of parameters for τ energy loss. The right hand side panel is a blow up of the intersection region.
The fraction of energy loss can be related to the τ range in rock through Eq. 19. Expanding the range for x << ξ we obtain: For E < ǫη ≃ 9 10 7 GeV the fraction of energy loss is below 1 and the above approximation holds. The range increases approximately linearly with E and so does the effective volume for neutrino detection through Earth skimming neutrinos. This is an important energy scale because above it, the effective volume for neutrinos ceases to grow linearly with E. A second energy scale arises because of absorption. In Fig. 5 the range at which the survival probability is 0.5 is shown, plotted as a function of E ′ . It is compared to the matter depth at which the flux has been reduced by a factor of 2 because of absorption. Two different cross sections and two different values of the energy loss parameters have been used to illustrate the uncertainties involved. The crossover energy sets the limiting scale to a value between 10 9 -10 10 GeV above which neutrino absorption suppresses the τ rate.

The Emerging τ flux
The calculation of the emerging τ flux proceeds in a relatively simple manner. Neutrinos arriving at a given zenith angle propagate through the Earth to a point along the corresponding Earth chord at which they interact. The total matter depth along the chord is x T (θ) and the depth of the interaction point, x, is measured along the chord from the exit point. We assume that the neutrino flux at x can be expressed as . The τ rate produced in an interval of distance dx is simply obtained convoluting the differential neutrino flux at the interaction point with the differential cross section dσ cc /dy multiplied by the τ survival probability: where E ′ τ and E τ are respectively the τ energies at the interaction and exit points and ρ is the matter density at the interaction point.
Using Eq. 18 we can express dx in terms of dE τ to get the differential τ rate as the following integral: The interaction point enters the expression through the function g(x T − x, E ν ). Once E ν and y are established, the τ lepton is produced with fixed energy E ′ τ = E ν (1 − y). If E τ is also fixed the interaction point x can be expressed in terms of the three variables x(E ν , E τ , y) as: The limits of the y integral in 26 are: The corresponding fluxes at different zenith close to the horizontal are illustrated in Fig 6 for the GRB flux used. A strong zenith angle dependence is observed which arises in Eq. 26 from x T , both in the function g(x, E ν ) (with x = x T − x) and in the limit y max . To a good approximation the function g(x, E ν ) can be thought of as an attenuation due to the total cross section σ tot which is effectively reduced by Figure 6. Differential τ spectra produced by earth-skimming neutrinos. All τ 's are included and the spectra are differential in τ energy, E τ . Left panel: τ spectra for the GRB flux at different zenith angles. Right panel: τ spectra after solid angle integration for both the TD and GRB fluxes.
regeneration. If we ignore regeneration due to the τ decays in matter we can simply express it as g(x, E ν ) = (x T −x)σ ef f where the effective absorption cross section is given by: and dσ nc /dy is the differential neutral current interaction cross section. It thus depends strongly on the zenith angle. We can finally integrate Eq. 26 over solid angle to obtain the energy spectrum of the total emerging τ rate. The results for both fluxes are also shown in Fig. 6. We note that the bulk of the τ lepton spectrum is contained within about 3 10 5 − 10 7 GeV for the prediction from GRB's, while for topological defects, which is a harder flux, the range is ∼ 10 6 − 3 10 8 GeV.
The earth skimming τ flux as a function of E τ is a complex result from the competition of an increasing effective volume and cross section for the neutrinos and a decreasing solid angle because of neutrino absorption. This is convoluted with a neutrino flux spectrum which often falls with energy like ∼ E −γ ν . The effective volume increases linearly with E τ until the scale of 9 10 7 GeV is reached above which the rise is slow. The cross section increases with energy as E ∼1/3 for E τ > 100 T eV while the solid angle decreases with approximately the inverse energy dependence. If γ > 2 the dominant part of the τ flux will be for energies at the 100 TeV scale which include all directions. But if γ < 2 the rate will be dominated by the turnover scale of 9 10 7 GeV. For hard spectral indices of order one, the dominant part of the flux shifts to the absorption scale ∼ 10 10 eV.

The τ Decay in the Atmosphere
Finally the emerging τ must decay to produce a detectable shower. The decay process is mediated by a W boson and always involves a ν τ which is irrelevant from the point of view of the developed shower. The charge boson couples to both electron and muon lepton pairs with probabilities close to 18%. When a µν µ pair is produced there is no shower, while if a eν e pair is produced about one third of the τ energy is expected to go into an electromagnetic shower. The remaining 64% of the times it couples to a quark antiquark pair that fragments into hadrons, mostly between one and three hadrons. These showers can be considered hadronic. Assuming that the shower carries all the quark energy we can expect about two thirds of the total τ energy to go into the air shower. When all the probabilities are weighted together the average fraction of energy carried by the shower is just below 50%.
The τ travels a distance l in the atmosphere to decay and produce a secondary shower that develops to reach shower maximum at a further distance l sh , which is a slowly varying function of shower energy (of order 10 radiation lengths or ∼ 3 km in air for both electromagnetic and hadronic showers). The altitude at which the shower reaches shower maximum is crucial for its detection and affects different types of detectors in different ways. This altitude is a combination of both the decay length and the distance to shower maximum (l + l sh ), the nadir angle with which the τ exits and the curvature of the Earth. The differential decay probability is again given by Eq. 20, where we can now neglect τ energy loss which is very small. The survival probability after traveling distance l is simply: For a given τ flux exiting the Earth, we can estimate the detection rate, integrating in length the flux multiplied by both the differential decay probability and the probability of detecting the shower: We note that provided that P τ det is known, the l integral can be in principle performed independently of the τ flux and expressed as an overall probability factor for detection, to multiply the rate obtained in Eq. 26. Unfortunately this requires detailed knowledge of the detector response and it is extremely dependent on the threshold behavior of the detector, which depends critically on shower development, geometry and the trigger conditions. The calculation of this is quite technical and difficult to estimate analytically. Precise calculations of P τ det must be preformed for each experiment. These functions are now being studied for detectors such as the surface array and fluorescence detector of the Pierre Auger Observatory, EUSO and others.
To understand the effect of this factor for an air shower array we consider that the shower will be detected provided that the particle density exceeds a given threshold at the point when the plane of shower maximum intercepts ground level. The center of the shower in this plane lies at a distance h = r g sin θ to the ground level where h is the altitude at which it lies above ground. The expression of h is simple in terms of l and l sh : where R ⊕ is the Earth's radius. For a given shower energy we calculate r max , the lateral distance below which the electron density of the shower is above threshold ρ e . We make h = r max sin θ and invert Eq. 32 to obtain the corresponding l. This value of l which we will denote l max is the maximum distance the τ can travel before decaying to be successfully detected according to our simplifying assumptions. The function chosen for P τ det allows analytic integration of the probability factor thus simplifying the calculation: We have expressed the detection probability integrated over decay length as an overall factor to multiply the τ flux rate. The corresponding suppression factor has been plotted in Fig. 7 for different zenith angles and different conditions to fix the minimum particle density to be detected at ground level. The probability is reduced for both low and high energy showers. Low energy showers do not meet the threshold requirement while high energy showers have maximum too high above ground. At energies around 10 8 GeV the detection restriction becomes important even for showers with 87 • nadir angle.

The Resulting Shower Rate for an Air Shower Array
The final shower rate recorded with an air shower array is severely constrained by the probability that the shower has sufficient particle density at ground level to be detected and identified as deep shower. This rate can be estimated taking the τ flux rate obtained in Eq. 26 and multiplying it by the integrated probability (as in Eq. 33).
We have calculated the resulting rates approximately as a function of the τ energy, taking into account that on average the τ decay showers carry about half the energy of the τ . Different detection sensitivities are crudely implemented using ρ e = 0.1, 0.4 and 1 m −2 . The total shower rates for the two models chosen are compared in Figs. 8 to the actual rate when the shower is required to be detected by an air shower array. The total shower rates are just those shown in the right hand panel of Fig. 6 rescaled to be differential in shower energy. When the decay probability and curvature of the Earth are accounted for as described in the previous subsection, the shower rate decreases quite dramatically at low energies when compared to the production rate, as illustrated by the figure. At high energies the detection is suppressed simply because the showers decay too far away from the surface and this is not too dependent on ρ e . At low energies detection details become crucial and constrain the detection in an intricate way. , Figure 8. Comparison of the total τ shower rate produced by earth-skimming ν τ +ν τ flux assuming that all τ ′ s make showers and the detected shower rate in an air shower array accounting for τ decay in the atmosphere and requiring electron number density at ground to be above different thresholds (ρ e > 0.1,0.4,1 cm −2 from top to bottom). The left panel is for the GRB flux and the right panel for the TD model.

Discussion
We have finally plotted the rates of down-going showers induced by neutrinos compared to the up-coming shower rate induced by Earth-skimming neutrinos in Fig. 9 for the two fluxes discussed in this article. The figure compares the results of Figs.4 and 8. In this figure the maximum earth-skimming τ shower rate (assuming all τ 's induce a shower) is clearly shown to be larger than the down-going neutrino induced shower rate in the energy region between 10 5 -10 9 GeV. A fully efficient 1000 km 2 active area operating for a year corresponds to 3 10 20 s cm 2 of aperture. This is the order of magnitude required for detection of the τ produced by earth-skimming τ neutrinos for the lowest of the two discussed models (TD) assuming all τ 's decay and are detected.
Unfortunately it is not possible to detect all these τ decays with most existing detectors. Fluorescence detectors would be limited by duty cycle while detection of these showers by arrays of particle detectors on Earth surface is made complicated by an awkward geometry as indicated by the curves corresponding to the different thresholds. Moreover the energy region where the earth-skimming rates dominate is somewhat below the typical energies for detectors in construction such as the Auger detector. For a 1000 km 2 area the expected peak rate is of order 4 (14) events per year for the TD (GRB) model between 10 6 and 4 10 7 GeV (3 10 5 and 3 10 6 GeV). This implies that if neutrinos are to be detected they would induce showers of energy very close to the detector threshold. As a result the rate calculations are quite uncertain unless very detailed simulations are performed. Some of these have been made and results are in agreement with our calculations [37,43].
The restriction introduced by air shower arrays somewhat compensates the duty factor for fluorescence detection and this explains why roughly similar rates can be obtained for the fluorescence and surface detector calculations of the earth skimming event rates. , Figure 9. Down-going inclined shower rates compared to earth-skimming up-coming shower rates from τ decays. The plot is the same as Fig. 8 but it also includes the shower rates obtained for down-going neutrinos as shown in Fig. 4.
When decay probability, earth curvature and detector sensitivity are included, the obtained rates for down-going and earth skimming events are not very different for an air shower array. The differential shower rate from τ decay (induced by earth-skimming neutrinos) dominates mostly in the energy region around 10 8 GeV and is limited to an interval of about one and a half decades around this value. The relative rates will be very dependent on the detector sensitivity in this energy range. The potential for neutrino detection relies on the detector capability to identify very inclined showers which are almost horizontal in the case of ν τ earth skimming events. Once the zenith angle has been established the detector must also identify them as late development or deep showers. Although τ decay showers induced by earth-skimming neutrinos have higher potential for neutrino detection than down-going neutrino events, the actual relative rates in the air shower arrays is very dependent on detector details and should be quite comparable for detectors having energy thresholds in the 10 8 GeV range.