---
abstract: 'The Pierre Auger (cosmic ray) Observatory provides a laboratory for studying fundamental physics at energies far beyond those available at colliders. The Observatory is sensitive not only to hadrons and photons, but can in principle detect ultrahigh energy neutrinos in the cosmic radiation. Interestingly, it may be possible to uncover new physics by analyzing characteristics of the neutrino flux at the Earth. By comparing the rate for quasi-horizontal, deeply penetrating air showers triggered by all types of neutrinos, with the rate for slightly upgoing showers generated by Earth-skimming tau neutrinos, we determine the ratio of events which would need to be detected in order to signal the existence of new non-perturbative interactions beyond the TeV-scale in which the final state energy is dominated by the hadronic component. We use detailed Monte Carlo simulations to calculate the effects of interactions in the Earth and in the atmosphere. We find that observation of 1 Earth-skimming and 10 quasi-horizontal events would exclude the standard model at the 99% confidence level. If new non-perturbative physics exists, a decade or so would be required to find it in the most optimistic case of a neutrino flux at the Waxman-Bahcall level and a neutrino-nucleon cross-section an order of magnitude above the standard model prediction.'
author:
- 'Luis A. Anchordoqui'
- Haim Goldberg
- Dariusz Góra
- Thomas Paul
- Markus Roth
- Subir Sarkar
- Lisa Lee Winders
title: |
    Using cosmic neutrinos to search for non-perturbative physics\
    at the Pierre Auger Observatory
---

Introduction
============

Ultrahigh energy cosmic neutrinos (UHEC$\nu$) are expected to be produced in association with the observed ultrahigh energy (charged) cosmic rays (UHECR), either at the same sites responsible for UHECR acceleration, or via interaction of the UHECR during propagation, particularly with the cosmic microwave background (CMB). These neutrinos are unique probes of new physics as their interactions are uncluttered by the strong and electromagnetic forces and, upon arrival at the Earth, they may experience collisions with center-of-mass energies up to $\sqrt{s} \alt 250$ TeV. However, rates for new physics processes are difficult to test since the flux of cosmic neutrinos is virtually unknown. Interestingly, it is possible in principle to disentangle the unknown flux and new physics processes by using multiple observables [@Kusenko:2001gj; @Anchordoqui:2001cg; @Anchordoqui:2005ey; @PalomaresRuiz:2005xw; @Anchordoqui:2006ta].

The Pierre Auger Observatory provides a promising way to detect UHEC$\nu$ by looking for deeply–developing, large zenith angle ($\agt
75^\circ$) or “quasi-horizontal” air showers [@Capelle:1998zz]. At these large angles, hadron-induced showers traverse the equivalent of several times the depth of the vertical atmosphere and consequently their electromagnetic component is extinguished before reaching the detector. Only very high energy muons survive past about 2 equivalent vertical atmospheres. Therefore, the shape of a hadron-induced shower front is very flat and prompt in time. In contrast, a neutrino shower exhibits the roughly same morphology as a vertical shower. It is therefore possible to distinguish neutrino induced events from background hadronic showers. Moreover, because of full flavor mixing, tau neutrinos are expected to be as abundant as other species in the cosmic flux. $\nu_\tau$’s can interact in the Earth’s crust, producing $\tau$ leptons which may decay above to the Auger detectors; such events will be referred to as “Earth–skimming” events [@Feng:2001ue; @Bertou:2001vm].

Possible deviations of the neutrino–nucleon cross-section due to new non-perturbative interactions[^1] can be uncovered at the Auger Observatory by combining information from Earth-skimming and quasi-horizontal showers. In particular, if an anomalously large rate is found for deeply developing quasi-horizontal showers, it may be ascribed either to an enhancement of the incoming neutrino flux, or an enhancement in the neutrino-nucleon cross-section (assuming non-neutrino final states dominate). However, these possibilities can be distinguished by comparing the rates of Earth-skimming and quasi-horizontal events. For instance, an enhanced flux will increase both quasi-horizontal and Earth-skimming event rates, whereas an enhanced interaction cross-section will also increase the former but [*suppress*]{} the latter, because the hadronic decay products cannot escape the Earth’s crust. Essentially this approach constitutes a straightforward counting experiment, as the detailed shower properties are not employed to search for the hypothesized new physics.

In this paper, we compute how many Earth-skimming [*vs.*]{} quasi-horizontal showers one would have to collect at the Auger Observatory to convincingly demonstrate new non-perturbative physics in which the final state energy is dominated by the hadronic component. We show that even a small number of events could be sufficient to rule out the standard model (SM). Thus the expected low neutrino “luminosity” is not at all a show-stopper, and the Observatory has the potential to uncover new physics at scales exceeding those accessible to the LHC.

In order to demonstrate this, we first compute acceptances for Earth-skimming and quasi-horizontal events using detailed models of the terrain in the vicinity of the Observatory as well as detailed simulations of the response of the Surface Array to highly inclined air showers. We then perform a likelihood analysis to determine the event counts needed to exclude the SM at various confidence levels. The analysis includes systematic effects both from theoretical uncertainties in the (perturbative QCD) SM cross-section [@CooperSarkar:2007cv] and from uncertainties associated with the detector response.

The outline of the paper is as follows. In Sec. \[sec:new\] we discuss some possible new physics scenarios which could manifest themselves as non-perturbative interactions at LHC energies and beyond. Next, in Sec. \[sec:sims\] we describe the detailed Monte Carlo studies of the acceptance for Earth-skimming and quasi-horizontal showers under the assumption of SM interactions, including systematic uncertainties. Finally in Sec. \[sec:results\] we perform the statistical analysis to ascertain the discovery reach of the Observatory. Our conclusions are collected in Sec. \[summary\]

New Non-Perturbative Physics at the LHC Energy Scale and Beyond {#sec:new}
===============================================================

The analysis techniques described herein constitute an entirely general approach to searching for non-perturbative interactions in which the final state is dominated by hadrons, without any dependence on what hypothetical mechanism might actually cause the ‘hadrophilia’. In order to illustrate some possible new physics signals which may be accessible using these techniques, we consider below two interesting possibilities.

TeV-scale mass black holes
--------------------------

In $D$-dimensional scenarios with large-compact-extra-dimensions (of common linear size $2\pi r_{\rm c}$) the Planck scale is related to the fundamental scale of gravity ($M_D$) according to [@Arkani-Hamed:1998rs] $$M_{\rm Pl}^2 = 8 \pi \,r_{\rm c}^{D-4} M_D^{D-2} \,\, .$$ If $M_D \agt M_W = G_{\rm F}^{-1/2} \simeq 300$ GeV, microscopic black holes (BH) can be produced [*gravitationally*]{} in particle collisions with center-of-mass energies $s \agt 1~{\rm
  TeV}$ [@Banks:1999gd]. Subsequently a TeV-scale BH would promptly decay via thermal Hawking radiation [@Hawking:1974rv] into observable quanta [@Emparan:2000rs]. (For $M_D = 1$ TeV, the lifetime of a BH of mass 10 TeV is less than $10^{-25}$ s.) Since gravitational coupling is flavor blind, a BH emits all the $\approx
120$ SM particle and anti-particle degrees of freedom with roughly equal probability. Accounting for color and spin, we expect $\approx
75\%$ of the particles produced in BH evaporation to be quarks and gluons, $\approx 10\%$ charged leptons, $\approx 5\%$ photons or $W/Z$ bosons, and $\approx 5\%$ neutrinos. Thus, TeV BH production and evaporation constitutes a clear example of beyond-SM non-perturbative physics.

Although such BH production cross-section $\sim M_W^{-1}$ is 5 orders of magnitude smaller than the QCD cross-section $\sim \Lambda_{\rm
  QCD}^{-1}$, it was proposed [@Dimopoulos:2001hw] that such BHs could be produced copiously at the LHC, and that these spectacular events could be easily filtered out of the QCD background. This is possible by triggering on BH events with prompt charged leptons and photons, each carrying hundreds of GeV of energy.

Cosmic ray collisions, with center-of-mass energies ranging up to $10^5$ GeV, certainly produce BHs if the LHC does. The question is, can they be detected? Most cosmic rays are protons or heavier nuclei, which collide with hadrons in the upper atmosphere, producing cascading showers which eventually reach the Earth’s surface. At energies of interest, however, the cosmic ray luminosity ($L \sim 7
\times 10^{-10}~(E/{\rm PeV})^{-2}$ cm$^{-2}$ s$^{-1}$, taking a single nucleon in the atmosphere as a target and integrating over $2
\pi$ sr), is about 50 orders of magnitude smaller than the LHC luminosity, thus making it futile to hunt for BHs in hadronic cosmic ray interactions. On the other hand, neutrino interaction lengths are far longer than the Earth’s atmospheric depth, although they would be greatly reduced by the cross-section for BH production [@Feng:2001ib]. Cosmic neutrinos therefore could produce BHs with roughly equal probability at any point in the atmosphere. As a result, the light descendants of the BH may initiate low-altitude, quasi-horizontal showers at rates significantly higher than SM predictions.

Analytic and numerical studies have revealed that gravitational collapse takes place at high energies and small impact parameters [@Eardley:2002re; @Yoshino:2002tx]. In the course of collapse, a certain amount of energy is radiated in gravitational waves, leaving a fraction $y \equiv {M_{\text{BH}}}/\sh$ available for Hawking evaporation. Here, ${M_{\text{BH}}}$ is a [*lower bound*]{} on the final mass of the BH and $\sh = 2 x m_N E_\nu$ is the center-of-mass energy of the colliding particles, taken to be partons. This ratio depends on the impact parameter of the collision, as well as on the dimensionality of space-time.

The inclusive production of BHs proceeds through different final states for different classical impact parameters $b$ [@Yoshino:2002tx]. These final states are characterized by the fraction $y(z)$ of the initial $\sqrt{\hat s}$ which is trapped within the horizon. Here, $z= b/b_{\rm max},$ and $b_{\rm max}= \sqrt{F} \,
r_s$ is the maximum impact parameter for collapse, where $$\label{schwarz}
r_\mathrm{s} =
\frac{1}{{M_D}}
\left[ \frac{\sqrt{\hat s}}{{M_D}} 
\, \frac{2^{D-4} \pi^{(D-7)/2}\Gamma({D-1\over 2})}{D-2}
\right]^{\frac{1}{D-3}}$$ is the radius of a $D$-dimensional Schwarzschild BH [@Myers:un], and $F$ is the form factor [@Yoshino:2002tx].

The $y$ dependance complicates the parton model calculation, since the production of a BH of mass $M_{\rm BH}$ requires that $\hat s$ be $M_{\rm BH}^2/y^2(z)$, thus requiring the lower cutoff on parton momentum fraction to be a function of impact parameter. Because of the complexity of the final state, we assume that amplitude intereference effects can be ignored and we take the $\nu N$ cross-section as an impact parameter-weighted average over partonic cross-sections, with the lower parton fractional momentum cutoff determined by $x_{\rm min} = M_{\rm BH}^{\rm min}/M_D$. This gives a lower bound ${\cal X} = ({x_{\text{min}}}M_D)^2/[y^2(z)s]$ on the parton momentum fraction $x$. All in all, the $\nu N \to {\rm BH}$ cross-section reads [@Anchordoqui:2003jr] $$\sigma (\nu N \to {\rm BH}) =  \int_0^1 2 z \,dz 
\int_{{\cal X}}^1 dx \, F \pi r_s^2 \,\sum_i f_i(x,Q) \ ,
\label{sigma}$$ where $i$ labels parton species and the $f_i(x,Q)$ are parton distribution functions (PDF).

As an illustration, we consider the $D=10$ string inspired scenario. For, $M_D = 1$ TeV, $x_{\rm min} = 1$, and primary neutrino energy $E_\nu = 10^{10}$ GeV, we obtain $\sigma (\nu N \to {\rm
  BH})\sim 2 \times 10^6~{\rm pb}$ [@Anchordoqui:2003jr]. This is about two orders of magnitude above SM predictions. The BH production cross-section by UHEC$\nu$ scales as $$\sigma(\nu N \to {\rm BH}) \propto 
 \left[\frac{1}{M_D^2}\right]^{\frac{D-2}{D-3}}  \,  .$$ A further suppression arises if $x_{\rm min}$ is increased. For parameters in the semiclassical regime ($x_{\rm min} \agt
3$ [@Anchordoqui:2003ug]) the BH cross-section becomes comparable to the SM cross-section at $M_D \sim 2$ TeV; this determines the multidimensional Planck scale to which Auger may be sensitive. However, the LHC will also be sensitive to extra-dimensional effects at a similar scale. It is interesting to consider whether Auger may have access to new physics [*beyond*]{} the reach of the LHC and we now discuss such a possibility.

Sphalerons
----------

In the electroweak theory, non-trivial fluctuations in $SU(2)$ gauge fields generate an energy barrier interpolating between topologically distinct vacua [@'tHooft:1976up]. An index theorem describing the fermion level crossings in the presence of these fluctuations reveals that neither baryon nor lepton number is conserved during the transition, but only the combination $B-L.$ Inclusion of the Higgs field in the calculation modifies the original instanton configuration [@Klinkhamer:1984di]. An important aspect of this modification (called the “sphaleron”) is that it provides an explicit energy scale of about 10 TeV for the height of the barrier. This barrier can be overcome through thermal transitions at high temperatures [@Kuzmin:1985mm], providing an important input to any calculation of cosmological baryogenesis. More speculatively, it has been suggested [@Aoyama:1986ej] that the topological transition could take place in two particle collisions at very high energy. The anomalous electroweak contribution to the partonic process can be written as $$\hat\sigma_i(\hat s)= 5.3\times 10^3 {\rm
    mb} \cdot e^{-(4\pi/\alpha_W)\ F_W(\epsilon)}\ \
  , \label{sigmahat}$$ where $\alpha_W\simeq 1/30$, the tunneling suppression exponent $F_W(\epsilon)$ is sometimes called the “holy-grail function”, and $\epsilon \equiv \sqrt{\hat s} /(4\pi
m_W/\alpha_W)\simeq \sqrt{\hat s}/30$ TeV. Thus, it is even possible that at or above the sphaleron energy the cross-section could be of ${\cal O}(\rm mb)$ [@Ringwald:2003ns]. Of particular interest to cosmic ray physicists would be enhancement of the neutrino cross-section over the perturbative SM estimates, say by an order of magnitude in the energy range $9.5 < \log_{10}
(E_\nu/{\rm GeV}) < 10.5$. With the methods outlined in this paper, this can be detected as an anomalous ratio of quasi-horizontal Earth-skimming showers to upcoming showers.

It was shown [@Ringwald:2003ns] that for the simple sphaleron configuration $s$-wave unitarity is violated for $\sqrt{\hat s}> 4\pi
M_W /\alpha_W\sim 36$ TeV. For lower parton subenergies, the cross-section is exponentially damped to values well below the perturbative SM electroweak value [@Ringwald:2003ns; @Bezrukov:2003er]. On the other hand, at the higher parton subenergies, the cross-section may well be dominated by non-spherically symmetric classical field configurations [@Gould:1993hb]. If for $\sqrt{\hat s} > 36$ TeV we saturate unitarity in each partial wave then this yields a geometric parton cross-section $\pi R^2$, where $R$ is some average size of the classical configuration. As a fiducial value we take the core size of the Manton-Klinkhamer sphaleron, $R\simeq 4/M_W\simeq
10^{-15} $ cm. In this simplistic model, the $\nu N$ cross-section is $$\sigma^{\nu N}_{\rm black \, disk} (E_\nu) 
 = \pi R^2\ \int_{x_{\rm min}}^1 \sum_{partons} f(x)\ dx \,,$$ where $x_{\rm min} = \hat s_{\rm min}/ s = (36)^2/ 2ME_\nu \simeq
0.065$. In the region $0.065 < x < 3(0.065)$ the PDF for the up and down quarks is well approximated by $f\simeq 0.5/x$, so the expression for the cross-section becomes $$\begin{aligned}
\sigma^{\nu N}_{\rm black \, disk} (E_\nu) 
&\simeq &\pi R^2\,\, (0.5)\,\, (\ln 3)\,\,  (2/2) \nonumber \\
& \simeq & 1.5\times 10^{-30} {\rm cm}^2  \,,\end{aligned}$$ where the last factor of 2/2 takes into account the (mostly) 2 contributing quarks $(u,d)$ in this range of $x,$ and the condition that only the left-handed ones contribute to the scattering. This is about 80 times the SM cross-section [@CooperSarkar:2007cv; @Gandhi:1998ri]. Of course this calculation is very approximate and thecross-section can easily be smaller by a factor of 10 (e.g., if $R$ is 1/3 of the fiducial value used).

Acceptance and Systematic Uncertainties {#sec:sims}
=======================================

To calculate the acceptance we perform detailed Monte Carlo simulations. The incoming neutrinos are propagated through the Earth’s crust, Andes mountains, and the atmosphere using an extended version [@Gora:2007nh] of the code ANIS  [@Gazizov:2004va]. For fixed neutrino energies, $10^6$ events are generated with zenith angles in the range $60^\circ-90^\circ$ (down-going showers) and $90^\circ-95^\circ$ (upgoing showers) and with azimuth angles in the range $0^\circ-360^\circ$. Neutinos are propagated along their trajectories of length $\Delta L$ from the generation point on the top of the atmosphere to the detector in steps of $\Delta L/1000 \, (\geq 6~{\rm km}).$ At each step of propagation, the $\nu N$ interaction probability $$\label{prob}
  P (E_{\nu},E_{l},\theta)  \simeq N_{A} \ \sigma^{\nu N}_{\rm SM}(E_{\nu}) \
  \rho(Z)  \ \Delta L \ ,$$ is calculated using the cross-section ($\sigma_{\rm SM}^{\nu N}
(E_{\nu})$) estimates of Ref. [@CooperSarkar:2007cv], where $\rho(Z)$ is the local medium density, $E_l$ the energy of the outgoing lepton, and $N_A \simeq 6.022 \times 10^{23}~{\text{g}}^{-1}$. The outcoming particle spectrum from $\nu N$ interactions is simulated with PYTHIA [@Sjostrand:2006za] and tau decays are simulated using the package TAUOLA  [@Jadach:1993hs].

The flux of otutgoing leptons as well as their energy and the decay vertex positions are calculated inside a defined detector volume. The geometrical size of the detector volume is set to $3000 \times 10~{\rm
  km}^3$ and it includes the real shape of the Auger Observatory on the ground. A relief map of the Andes mountains was constructed according to a digital elevation data of the Consortium for Spatial Information (CGIAR-CSI) [@CGIAR-CSI]. The map of the area around the Auger site is shown in Fig. \[augermap\].

The detection volume corresponds to the so called [*active volume*]{} in which potentially detectable neutrino interactions are simulated. For a given incoming neutrino with energy $E_{\nu}$ the active volume is defined by a particular plane $A_{\rm gen}$ and distance $\Delta L$. The plane $A_{\rm gen}$ is the cross-sectional area of the detector volume and it is used as a reference plane for the generation of incoming neutrinos. The area depends on the zenith angle $\theta$ of the incoming neutrino. The distance $\Delta L$ is the multiple, $n$, of the average lepton range $\langle
L^l(E_l)\rangle$.

Earth-skimming events occur in the Earth’s crust, and so the relevant neutrinos and taus sample only the Earth’s surface density, $\rho_s
\approx 2.65~{\text{g}}/{\text{cm}}^3$. At these energies, the tau’s propagation length is determined not by its decay length but by its energy loss. The energy loss per unit length of crossed matter is usually approximated by a linear equation (continuous energy loss approach). The $\tau$ lepton loses energy in the Earth according to $$\frac{dE_{\tau}}{dz} = -(\alpha_{\tau} + \beta_{\tau} E_{\tau})
\rho_s \ ,$$ where the factor $\alpha_\tau$ parametrizes the ionisation losses and $\beta_\tau$ the energy losses through bremsstrahlung, pair production, and hadronic interactions. For $E_\nu = 10^7~{\rm GeV},$ $\alpha_\tau$ is negligible and $\beta_{\tau} \approx 0.8\times
10^{-6}~{\text{cm}}^2/{\text{g}}$ [@Dutta:2000hh]. Hadronic interactions (i.e., lepton-nucleus inelastic interactions dominated by small values of the squared momentum transfer $Q^2$) are responsible for the largest and the most uncertain contribution [@Armesto:2007tg]. Such an uncertainty in $\beta_\tau$ dominates the systematic errors in the estimate of the neutrino event rates.

To investigate the response of the Auger detector, we generate the lateral profiles of the shower development using the output of PYTHIA and/or TAUOLA as input for AIRES [@Sciutto:2001dn]. The showers induced by the products of up-going decaying tau leptons, with energies from 0.1 EeV to 100 EeV and decay position at altitudes ranging from 0 to 3500 m above sea leavel, are simulated in steps of 100 m. At each altitude 40 events are generated to cover the tau decay channels implemented in ANIS [@Gora:2007nh]. In the case of down-going showers, the decay altitudes range from ground level up to the upper atmosphere.

The response of the surface detector array is simulated in detail using the simulation package [@Argiro:2007qg]. Besides the standard procedure to simulate the spacial and temporal signal response we have added the simulation of atmospheric background muons to study the impact on the neutrino identification, since such accidental muons might be wrongly classified as shower particles. The background from hadronic showers above $10^8$ GeV is estimated to be ${\cal O} (1)$ in 20 years. At $E_\nu \sim 10^{10}$ GeV the cosmic ray flux is $\approx 10^6$ time smaller, so the expected background for the energy bin considered here ($9.5 < \log_{10} (E_\nu/{\rm GeV})
< 10.5$) is negligible.

The expected neutrino event rate (of flavor $\alpha$) in the detector volume is found to be $$\label{nrate}
   N_{\nu_\alpha}=F^{w}_{\nu}   \, \sum_{i=1}^{N_{\rm acc}}  P_{i},$$ where $N_{\rm acc}$ is the number of events triggering the detector and passing all quality cuts of the cascade analysis. Here, $$\label{flux}
  F^{w}_{\nu}=  N_{\rm gen}^{-1}  \, \Delta T \,\int_{E_{\rm min}}^{E_{\rm max}}
  \Phi_0^{\nu_\alpha}(E_{\nu})  \, dE_{\nu}   \int_{\theta_{\rm min}}^{\theta_{\rm max}}
  A_{\rm gen}(\theta) \, d\Omega \,,$$ $d\Omega$ is the solid angle, $\Delta T$ the observation time, $N_{\rm
  gen}$ is the number of generated events from surface $A_{\rm gen}$, and we take the neutrino flux $\Phi_0^{\nu_\alpha}(E_{\nu})$ to be isotropic. We further assume $\nu_e : \nu_\mu : \nu_\tau \simeq
1:1:1$, which is generally thought to be the case if the neutrinos are produced predominantly thorugh pion decay. In order to ascertain the systematic uncertainties associated with our lack of knowlege of the dependence of the flux on energy, we consider three scenarios which plausibly bracket the range of possibilities:

1.  $  \Phi_0^{\nu_\alpha}(E_\nu) = ( {\cal C}/E_0) \, E_\nu^{-1}$,

2.  $ \Phi_0^{\nu_\alpha}(E_\nu) = {\cal C} \, E_\nu^{-2}$,

3.  $\Phi_0^{\nu_\alpha}(E_\nu) = ( {\cal C}/E_0) \, E_\nu^{-3}$,

4.  $\Phi_0^{\nu_\alpha}(E_\nu) = {\cal C} E_\nu^{-2} \, 
    {\rm exp} [-\log_{10} (E_\nu/E_0)^2/(2 \sigma^2)]$,

where ${\cal C} = 2.33 \times 10^{-8}~{\rm GeV}\, {\rm s}^{-1}\, {\rm
  cm}^{-2}\, {\rm sr}^{-1}$, $E_0 = 10^{10}$ GeV, $\sigma =
0.5$ GeV. This normalization (2) constitutes the common benchmark, the so-called ‘Waxman-Bahcall bound’ [@Waxman:1998yy]. The factor $F^{w}_{\nu}$ of Eq. \[flux\] is chosen to yield the total number of events per year. The expected rates for the entire range over which Auger is sensitive are given in Table \[table1\] and the rates for the high energy bin considered in the following study are given in Table \[table2\].

   flux                                 ratio                                                                           
  ------- ---------- ---------------- ---------- ------------- ---------------- --------------- ----------------------- ----------------------------
           $\theta$   $N_{\nu_\tau}$   $\theta$   $N_{\nu_e}$   $N_{\nu_\tau}$   $N_{\nu_\mu}$   $N_{\nu_{\rm all} } $   $N_\tau/N_{\nu_{\rm all}}$
   $(2)$     90-95        0.68           60-90       0.134          0.109            0.019              0.262                       2.58
   $(2)$     90-95        0.68           75-90       0.075          0.071            0.011              0.157                       4.27

  : Expected events per year ($N_i$) at Auger in the energy range $8 < \log_{10} (E_\nu/{\rm GeV})$, for various incident zenith angle ($\theta$) ranges, assuming the Waxman-Bahcall flux.[]{data-label="table1"}

   flux                                 ratio                                                                          
  ------- ---------- ---------------- ---------- ------------- ---------------- --------------- ---------------------- ----------------------------
           $\theta$   $N_{\nu_\tau}$   $\theta$   $N_{\nu_e}$   $N_{\nu_\tau}$   $N_{\nu_\mu}$   $N_{\nu_{\rm all}} $   $N_\tau/N_{\nu_{\rm all}}$
   $(1)$     90-95        0.14           60-90       0.059          0.049            0.011              0.12                       1.14
   $(2)$     90-95        0.15           60-90       0.059          0.049            0.096              0.11                       1.33
   $(3)$     90-95         0.23          60-90       0.079          0.062           0.0123              0.15                       1.53
   $(4)$     90-95         0.12          60-90       0.046          0.037           0.0080              0.091                      1.33
   $(1)$     90-95         0.14          75-90       0.027          0.031           0.0056              0.064                      2.14
   $(2)$     90-95         0.15          75-90       0.026          0.029           0.0048              0.060                      2.47
   $(3)$     90-95         0.23          75-90       0.036          0.041           0.0062              0.083                      2.75
   $(4)$     90-95         0.12          75-90       0.021          0.024           0.0040              0.049                      2.45

  : Expected events per year ($N_i$) at Auger in the energy range $9.5 < \log_{10} (E_\nu/{\rm GeV}) < 10.5$, for various incident zenith angle ($\theta$) ranges and the 4 flux models considered.[]{data-label="table2"}

Hereafter we consider $\Phi_0^{\nu_\alpha} (E_\nu) \propto E_\nu^{-2}$ as our nominal spectrum. We then estimate systematic uncertainties associated with: different assumptions of the spectrum shape, different parton distribution functions (GRV92NLO [@Gluck:1998js] and CTEQ66c [@Nadolsky:2008zw]), and different estimates on $\beta_\tau$ [@Armesto:2007tg]. The contribution of different systematic errors are listed in Table \[table3\].

       ratio           flux           PDF           $\beta_\tau$                     sum    
  --------------- -------------- ------------- ------------------------ ------------------------------
                    $+ 11\% $        $0\%$             $+24\%$           $\phantom{2.4666666} + 26\%$
       2.47                                                                          2.47
                     $-13\%$        $-21\%$            $-25\%$           $\phantom{2.4666666} - 35\%$

  : Contributions to the systematic uncertainty on the Earth-skimming to quasi-horizontal event ratio. We have considered the energy range $9.5 < \log_{10} (E_\nu/{\rm GeV}) < 10.5$ and the zenith angle range $75^\circ < \theta < 90^\circ$.[]{data-label="table3"}

Auger Discovery Reach {#sec:results}
=====================

Consider a flux of neutrinos with energy in the range $10^{9.5}~{\rm
  GeV}<E_\nu< 10^{10.5}~{\rm GeV}$. In the SM, the interaction path length is $$L_{\rm CC}^{\nu} = \left[ N_A \rho_s \sigma^{\nu}_{\rm CC}
\right] ^{-1} \ ,$$ where $\sigma^{\nu}_{\rm CC}$ is the charged current cross-section for $E_{\nu} = E_0$. (We neglect neutral current interactions, which at these energies serve only to reduce the neutrino energy by approximately 20%, which is within the systematic uncertainty.) For $E_0 \sim 10^{10}~{\rm GeV}$, $L_{\rm CC}^{\nu} \sim {\cal
  O}(100)~{\text{km}}$. Supplemented by the possibility of new non-perturbative physics, the interaction path length is $$L_{\rm tot}^{\nu} = \left[ N_A \rho_s (\sigma^{\nu}_{\rm CC} +
\sigma^{\nu}_{\rm NP}) \right] ^{-1} \ ,$$ where $\sigma^{\nu}_{\rm NP}$ is the new physics contribution to the cross-section for $E_{\nu} = E_0$.

The maximal path length for a detectable $\tau$ is given by $$L^{\tau} = \frac{1}{\beta_{\tau} \rho_s} \ln \left( E_{\rm max} /
E_{\rm min} \right) \ ,
\label{ltau}$$ where $E_{\rm max} \approx E_0$ is the energy at which the tau is created, and $E_{\rm min}$ is the minimal energy at which a $\tau$ can be detected. For $E_{\rm max} / E_{\rm min} = 10$, $L^{\tau} =
11~{\text{km}}$.

Given an isotropic $\nu_{\tau} + \bar{\nu}_{\tau}$ flux, the number of taus that emerge from the Earth with sufficient energy to be detected is proportional to an “effective solid angle” $$\Omega_{\rm eff} \equiv \int d\cos\theta\, d\phi\, \cos\theta\,
P(\theta,\phi) \ ,$$ where $$P(\theta,\phi) = \int_0^{\ell} \frac{dz}{L_{\rm CC}^{\nu}}
e^{-z/L_{\rm tot}^{\nu}} \
\Theta \left[ z - (\ell - L^{\tau} ) \right]
\label{P}$$ is the probability for a neutrino with incident nadir angle $\theta$ and azimuthal angle $\phi$ to emerge as a detectable $\tau$. (In Eq., for the reasons noted above, we have neglected the possibility of detectable signals from new non-perturbative physics by Earth-skimming neutrinos.) Here $\ell = 2 R_{\oplus} \cos\theta$ is the chord length of the intersection of the neutrino’s trajectory with the Earth, with $R_{\oplus} \approx 6371~{\text{km}}$ the Earth’s radius. Evaluating the integrals, we find [@Kusenko:2001gj] $$\begin{aligned}
\Omega_{\rm eff} & = & 2 \pi
\frac{L_{\rm tot}^{\nu}}{L_{\rm CC}^{\nu}}
\left[ e^{L^{\tau} / L_{\rm tot}^{\nu}} - 1 \right] 
\left[ \left( \frac{L_{\rm tot}^{\nu}}{2 R_{\oplus}} \right)^2 \right. \nonumber \\
&- & \left. \left( \frac{L_{\rm tot}^{\nu}}{2 R_{\oplus}} +
\left( \frac{L_{\rm tot}^{\nu}}{2 R_{\oplus}} \right)^2 \right)
e^{-2R_{\oplus} / L_{\rm tot}^{\nu}} \right] \ .
\label{Omegaeff}\end{aligned}$$ At the relevant energies, the neutrino interaction length satisfies $L_{\rm tot}^{\nu} \ll R_{\oplus}$. In addition, for $L_{\rm
  tot}^{\nu} \gg L^{\tau}$, valid when the cross-section enhancement is significant but not so large as typical hadronic cross-section, Eq. simplifies to [@Anchordoqui:2001cg] $$\Omega_{\rm eff} \approx
2\pi \frac{L_{\rm tot}^{\nu\, 2} L^{\tau}}{4 R_{\oplus}^2
L_{\rm CC}^{\nu}} \ .
\label{omegaeffapprox}$$

Equation (\[omegaeffapprox\]) gives the functional dependence of the Earth-skimming event rate on the new physics cross-section. This rate is, of course, also proportional to the source neutrino flux $\Phi^{\nu_{\rm all}}$ at $E_0$. Given these inputs, $$N_{\rm ES} \approx
C_{\rm ES}\, \frac{\Phi^{\nu_{\rm all}}}{\Phi^{\nu_{\rm all}}_0}
\frac{\sigma^{\nu\, 2}_{\rm CC}}
{\left( \sigma^{\nu}_{\rm CC} + \sigma^{\nu}_{\rm NP} \right)^2} \ ,
\label{ES}$$ where $C_{\rm ES} = 0.15$ is the number of Earth-skimming events expected for a fiducial flux $\Phi^{\nu_{\rm all}}_0$ in the abscence of new physics.

In contrast to Eq., the rate for quasi-horizontal showers has the form $$N_{\rm QH} = C_{\rm QH} \frac{\Phi^{\nu_{\rm all}}}{\Phi^{\nu_{\rm all}}_0}
\frac{\sigma^{\nu}_{\rm CC} + \sigma^{\nu}_{\rm NP}}
{\sigma^{\nu}_{\rm CC}} \ ,
\label{QH}$$ where $C_{\rm QH} = 0.06$ for the Auger Surface Array, as determined in Sec. \[sec:sims\].

Given a flux $\Phi^{\nu_{\rm all}}$ and new physics cross-section $\sigma^{\nu}_{\rm NP}$, both $N_{\rm ES}$ and $N_{\rm QH}$ are determined. On the other hand, given just a quasi-horizontal event rate $N_{\rm QH}$, it is impossible to differentiate between an enhancement of the cross-section due to new physics and an increase on the flux. However, in the region where significant event rates are expected the contours of $N_{\rm QH}$ and $N_{\rm ES}$, given by Eqs. (\[ES\]) and (\[QH\]), are more or less orthogonal and provide complementary information. With measurements of $N_{\rm
  QH}^{\rm obs}$ and $N_{\rm ES}^{\rm obs}$, both $\sigma^{\nu}_{\rm
  NP}$ and $\Phi^{\nu_{\rm all}}$ may be determined independently, and neutrino interactions beyond the SM may be unambiguously identified.

We now turn to determining the projected sensitivity of Auger to neutrino fluxes and cross-sections. The quantities $N_{\rm ES}$ and $N_{\rm QH}$ as defined in Eqs. (\[ES\]) and (\[QH\]) can be regarded as the theoretical values of these events, corresponding to different points in the $\Phi^{\nu_{\rm all}}/\Phi_0^{\nu_{\rm
    all}}-\sigma_{\rm NP}/\sigma_{\rm CC}$ parameter space. For a given set of observed rates $N^{\rm obs}_{\rm ES}$ and $N^{\rm
  obs}_{\rm QH}$, two curves are obtained in the two-dimensional parameter space by setting $N^{\rm obs}_{\rm ES} = N_{\rm ES}$ and $N^{\rm obs}_{\rm QH} = N_{\rm QH}$. These curves intersect at a point, yielding the most probable values of flux and cross section for the given observations. Fluctuations about this point define contours of constant $\chi^2$ in an approximation to a multi-Poisson likelihood analysis. The contours are defined by $$\chi^2  =  \sum_i 2 \,\left[ N_i  -  N_i^{\rm obs}\right] +
2\,  N_i^{\rm obs}\,
\ln \left[ N_i^{\rm obs}/ N_i \right] \,\,,
\label{baker}$$ where $i=$ ES, QH [@Baker:1983tu]. In Fig. \[statistics\], we show results for three representative cases. Assuming ($N^{\rm
  obs}_{\rm ES} = 1$, $N^{\rm obs}_{\rm QH} = 10$), ($N^{\rm obs}_{\rm
  ES} = 1$, $N^{\rm obs}_{\rm QH} = 7$), and ($N^{\rm obs}_{\rm ES} =
1$, $N^{\rm obs}_{\rm QH} = 5$) we show the 90%, 95%, 99% and 3$\sigma$ CL contours for 2 d.o.f. ($\chi^2 = 4.61,\ 5.99,\ 9.21,$ and 11.83, respectively). For $N^{\rm obs}_{\rm ES} = 1$ and $N^{\rm
  obs}_{\rm QH} = 10$, the possibility of a SM interpretation along the $\sigma^{\nu}_{\rm NP} = 0$ axis (taking into account systematic uncertainties) would be excluded at greater than 99% CL for [*any*]{} assumed flux. The power of the Earth-skimming information is such that the best fit consistent with the SM would require a flux of about 50 times the Waxman-Bahcall flux, which is already excluded by present limits [@Anchordoqui:2009nf].

Summary
=======

We have re-examined a technique to search for new physics at sub-fermi distances. The strategy involves determining the ratio of quasi-horizontal to Earth-skimming showers initiated by cosmic neutrinos which would need to be detected by the Pierre Auger Observatory in order to signal the existence of exotic non-perturbative interactions beyond the TeV-scale. We perform Monte Carlo simulations of neutrino interactions in the Earth and in the atmosphere, and realistic simulation of the detector acceptance using the Auger software. We find that observation of 1 Earth-skimming and 10 quasi-horizontal events would exclude the standard model at the 99% confidence level. If new non-perturbative physics exists, a decade or so would be required to uncover it in the most optimistic case (cosmic neutrino flux at the Waxman-Bahcall level and $\nu N$ cross-section about an order of magnitude above the standard model prediction). The proposed Northern Auger site [@Abraham:2009en] (which has not been optimized for neutrino studies) would reduce this time by about a factor of 2. Any hint of such an important signal would provide an impetus to infill the array to increase the neutrino acceptance.

Acknowledgements {#acknowledgements .unnumbered}
================

We would like to thank Mandy Cooper-Sarkar for discussions. L.A.A. is supported by the U.S. National Science Foundation (NSF) Grant No PHY-0757598, and the UWM RGI. H.G. is supported by the NSF Grant No PHY-0757959. D.G. is supported by the HHNG-128 grant of the Helmholtz association and the Polish Ministry of Science and Higher Education under Grant 2008 No. NN202 127235. T.P. is supported by the NSF Grant No PHY-0855388. M.R. is supported by the HHNG-128 grant of the Helmholtz association. S.S. acknowledges support by the EU Marie Curie Network “UniverseNet” (HPRN-CT-2006-035863). L.L.W. is supported by the UWM RGI. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

[99]{}

A. Kusenko and T. J. Weiler, Phys. Rev. Lett.  [**88**]{}, 161101 (2002) \[arXiv:hep-ph/0106071\]. L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, *Phys. Rev.  D* **65**, 124027 (2002). L. Anchordoqui, T. Han, D. Hooper and S. Sarkar, Astropart. Phys.  [**25**]{}, 14 (2006) \[arXiv:hep-ph/0508312\]. S. Palomares-Ruiz, A. Irimia and T. J. Weiler, Phys. Rev.  D [**73**]{} (2006) 083003 \[arXiv:astro-ph/0512231\]. L. A. Anchordoqui, A. M. Cooper-Sarkar, D. Hooper and S. Sarkar, Phys. Rev.  D [**74**]{}, 043008 (2006) \[arXiv:hep-ph/0605086\]. K. S. Capelle, J. W. Cronin, G. Parente and E. Zas, Astropart. Phys.  [**8**]{}, 321 (1998) \[arXiv:astro-ph/9801313\]. X. Bertou, P. Billoir, O. Deligny, C. Lachaud and A. Letessier-Selvon, Astropart. Phys.  [**17**]{}, 183 (2002) \[arXiv:astro-ph/0104452\]. J. L. Feng, P. Fisher, F. Wilczek and T. M. Yu, Phys. Rev. Lett.  [**88**]{}, 161102 (2002) \[arXiv:hep-ph/0105067\]. A. Cooper-Sarkar and S. Sarkar, JHEP [**0801**]{}, 075 (2008) \[arXiv:0710.5303 \[hep-ph\]\]. N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, Phys. Lett. B [**429**]{}, 263 (1998) \[arXiv:hep-ph/9803315\]; I. Antoniadis, N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, Phys. Lett. B [**436**]{}, 257 (1998) \[arXiv:hep-ph/9804398\]. T. Banks and W. Fischler, arXiv:hep-th/9906038. S. W. Hawking, Nature [**248**]{}, 30 (1974). R. Emparan, G. T. Horowitz and R. C. Myers, Phys. Rev. Lett.  [**85**]{}, 499 (2000) \[arXiv:hep-th/0003118\]. S. Dimopoulos and G. Landsberg, Phys. Rev. Lett.  [**87**]{}, 161602 (2001) \[arXiv:hep-ph/0106295\]; S. B. Giddings and S. D. Thomas, Phys. Rev. D [**65**]{}, 056010 (2002) \[arXiv:hep-ph/0106219\]. J. L. Feng and A. D. Shapere, Phys. Rev. Lett.  [**88**]{}, 021303 (2002) \[arXiv:hep-ph/0109106\]; L. Anchordoqui and H. Goldberg, Phys. Rev.  D [**65**]{}, 047502 (2002) \[arXiv:hep-ph/0109242\]; A. Ringwald and H. Tu, Phys. Lett. B [**525**]{}, 135 (2002) \[arXiv:hep-ph/0111042\]; L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, Phys. Rev.  D [**66**]{}, 103002 (2002) \[arXiv:hep-ph/0207139\]; E. J. Ahn, M. Ave, M. Cavaglia and A. V. Olinto, Phys. Rev.  D [**68**]{}, 043004 (2003) \[arXiv:hep-ph/0306008\]. D. M. Eardley and S. B. Giddings, Phys. Rev. D [**66**]{}, 044011 (2002) \[arXiv:gr-qc/0201034\]. H. Yoshino and Y. Nambu, Phys. Rev. D [**67**]{}, 024009 (2003) \[arXiv:gr-qc/0209003\]. R. C. Myers and M. J. Perry, Annals Phys.  [**172**]{}, 304 (1986). L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, Phys. Rev. D [**68**]{}, 104025 (2003) \[arXiv:hep-ph/0307228\]. L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, Phys. Lett. B [**594**]{}, 363 (2004) \[arXiv:hep-ph/0311365\]. G. ’t Hooft, Phys. Rev. Lett.  [**37**]{}, 8 (1976); G. ’t Hooft, Phys. Rev.  D [**14**]{}, 3432 (1976) \[Erratum-ibid.  D [**18**]{}, 2199 (1978)\]. F. R. Klinkhamer and N. S. Manton, Phys. Rev.  D [**30**]{}, 2212 (1984). V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett.  B [**155**]{}, 36 (1985); M. Fukugita and T. Yanagida, Phys. Lett.  B [**174**]{}, 45 (1986); P. B. Arnold and L. D. McLerran, Phys. Rev.  D [**36**]{}, 581 (1987). H. Aoyama and H. Goldberg, Phys. Lett.  B [**188**]{}, 506 (1987); A. Ringwald, Nucl. Phys.  B [**330**]{}, 1 (1990); O. Espinosa, Nucl. Phys.  B [**343**]{}, 310 (1990). A. Ringwald, JHEP [**0310**]{}, 008 (2003) \[arXiv:hep-ph/0307034\]. F. L. Bezrukov, D. Levkov, C. Rebbi, V. A. Rubakov and P. Tinyakov, Phys. Rev.  D [**68**]{}, 036005 (2003) \[arXiv:hep-ph/0304180\]; F. L. Bezrukov, D. Levkov, C. Rebbi, V. A. Rubakov and P. Tinyakov, Phys. Lett.  B [**574**]{}, 75 (2003) \[arXiv:hep-ph/0305300\]. T. M. Gould and S. D. H. Hsu, Mod. Phys. Lett.  A [**9**]{}, 1589 (1994) \[arXiv:hep-ph/9311291\]. R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Phys. Rev.  D [**58**]{}, 093009 (1998) \[arXiv:hep-ph/9807264\]. D. Gora, M. Roth and A. Tamburro, Astropart. Phys.  [**26**]{}, 402 (2007). A. Gazizov and M. P. Kowalski, Comput. Phys. Commun.  [**172**]{}, 203 (2005) \[arXiv:astro-ph/0406439\]. T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP [**0605**]{}, 026 (2006) \[arXiv:hep-ph/0603175\]. S. Jadach, Z. Was, R. Decker and J. H. Kuhn, Comput. Phys. Commun.  [**76**]{}, 361 (1993). Consortium for Spatial Information (CGIAR-CSI). [http://srtm.csi.cgiar.org/]{}

S. I. Dutta, M. H. Reno, I. Sarcevic and D. Seckel, Phys. Rev.  D [**63**]{}, 094020 (2001) \[arXiv:hep-ph/0012350\]. N. Armesto, C. Merino, G. Parente and E. Zas, Phys. Rev.  D [**77**]{}, 013001 (2008) \[arXiv:0709.4461 \[hep-ph\]\]. See also, E. V. Bugaev and Yu. V. Shlepin, Phys. Rev.  D [**67**]{}, 034027 (2003) \[arXiv:hep-ph/0203096\]; H. Abramowicz and A. Levy, arXiv:hep-ph/9712415; A. Capella, A. Kaidalov, C. Merino and J. Tran Thanh Van, Phys. Lett.  B [**337**]{}, 358 (1994) \[arXiv:hep-ph/9405338\]. S. J. Sciutto, arXiv:astro-ph/0106044. S. Argiro [*et al.*]{}, Nucl. Instrum. Meth.  A [**580**]{}, 1485 (2007) \[arXiv:0707.1652 \[astro-ph\]\]. E. Waxman and J. N. Bahcall, Phys. Rev.  D [**59**]{}, 023002 (1999) \[arXiv:hep-ph/9807282\]. M. Gluck, S. Kretzer and E. Reya, Astropart. Phys.  [**11**]{}, 327 (1999) \[arXiv:astro-ph/9809273\]. P. M. Nadolsky [*et al.*]{}, Phys. Rev.  D [**78**]{}, 013004 (2008) \[arXiv:0802.0007 \[hep-ph\]\]. S. Baker and R. D. Cousins, Nucl. Instrum. Meth. A [**221**]{}, 437 (1984). A. M. Cooper-Sarkar, private communication.

L. A. Anchordoqui and T. Montaruli, arXiv:0912.1035 \[astro-ph.HE\]. J. Blümer [*et al.*]{} \[The Pierre Auger Collaboration\], New J. Phys. [**12**]{}, 035001 (2010).

[^1]: Throughout this paper we use this term to describe neutrino interactions in which the final state energy is dominated by the hadronic component. We are [*not*]{} considering here new “perturbative” physics e.g. (softly broken) supersymmetry at the TeV scale which would have quite different signatures in cosmic ray showers.
