Probing new physics with longlived charged particles produced by atmospheric and astrophysical neutrinos
Abstract
As suggested by some extensions of the Standard Model of particle physics, dark matter may be a superweakly interacting lightest stable particle, while the nexttolightest particle (NLP) is charged and metastable. One could test such a possibility with neutrino telescopes, by detecting the charged NLPs produced in highenergy neutrino collisions with Earth matter. We study the production of charged NLPs by both atmospheric and astrophysical neutrinos; only the latter, which is largely uncertain and has not been detected yet, was the focus of previous studies. We compute the resulting fluxes of the charged NLPs, compare those of different origins, and analyze the dependence on the underlying particle physics setup. We point out that even if the astrophysical neutrino flux is very small, atmospheric neutrinos, especially those from the prompt decay of charmed mesons, may provide a detectable flux of NLP pairs at neutrino telescopes such as IceCube. We also comment on the flux of charged NLPs expected from proton–nucleon collisions, and show that, for theoretically motivated and phenomenologically viable models, it is typically subdominant and below detectable rates.
, , ,
1 Introduction
The fundamental nature of dark matter poses a profound challenge to contemporary theoretical particle physics. Observations constrain the neutrino—the only electrically and color neutral nonbaryonic elementary particle within the Standard Model of particle physics—to have a negligible contribution to the overall dark matter budget [1, 2]. Dark matter is regarded as one of the most compelling hints towards new physics beyond the Standard Model. The question of its elementary essence has triggered enormous theoretical and phenomenological efforts [3].
The existence of suitable dark matter particle candidates in several theoretically cogent extensions of the Standard Model, like lowenergy supersymmetry [4] or extradimensional scenarios [5, 6], focused a strong interest on a class of dark matter candidates known as weakly interacting massive particles (WIMPs). Similarly to other Standard Model particles, WIMPs would fall out of thermal equilibrium and freeze out in the early Universe, leaving a relic abundance compatible with the inferred amount of dark matter [7]. These WIMPs can be directly detected by experiments looking for the minuscule energy deposition caused by dark matter particles scattering nuclei [8]. The pair annihilation of WIMPs into energetic gamma rays, neutrinos, and antimatter, is a second, yet indirect, handle on the presence and potential imprint of galactic particle dark matter [9].
The connection of the aforementioned scenarios to the electroweak scale, soon to be probed at the Large Hadron Collider (LHC), motivated the exploration of complementarity between collider searches for new physics and the question of the elementary nature of dark matter [10, 11, 12, 13, 14]. In most cases, if dark matter is a WIMP, the anticipated experimental signature at LHC would be the production of strongly interacting massive particles which promptly decay to the lightest and stable WIMP, plus a number of energetic jets and leptons. The neutral dark matter particle would escape the detector unobserved, leading to large missing transverse energy as well [15]. Conclusive identification of escaping neutral particles at LHC with dark matter permeating our Galaxy and other cosmic structures would, however, require some evidence from the direct and/or indirect WIMP searches listed before [9].
WIMPs are indeed attractive dark matter candidates, but are not the only theoretically envisioned possibility. The dark matter particle could exhibit even feebler interactions with ordinary Standard Model particles than a WIMP, making direct and indirect searches completely hopeless. For instance, the supersymmetric gravitino [16] or the KaluzaKlein graviton of universal extra dimensions [17, 6] are perfectly plausible ‘superweakly interacting’ [18] dark matter candidates (superWIMPs). If Nature chose an option like this, collider signatures of new physics, if any, would strongly depend upon the nature of the nexttolightest supersymmetric particle (NLSP). Since a superWIMP is also very weakly coupled to the other newphysics heavier states, the NLSP would likely be quasistable. If the NLSP is neutral, the qualitative experimental landscape would look like that of a standard WIMP scenario. However, we would lack the needed proof of a connection between the weakly interacting longlived particles produced at colliders and galactic dark matter.
If the NLSP is instead charged (constituting a charged massive particle, or CHAMP), the LHC would potentially observe the extremely distinct signature of a ‘heavy muon’: charged tracks and penetration of the outer muon subdetector, possibly at very low relativistic beta. CHAMPs are constrained by direct collider searches at LEP2 [19], as well as at the Tevatron [20].^{1}^{1}1These and future searches at LHC are not trivial, but advanced work is being done to ensure that such potential signals would not be missed [21]. Detection for masses TeV is essentially guaranteed.
If a CHAMP were stable on collider scales, it could still decay on cosmological scales, and thus impact precision astrophysical measurements [22], including the chemical potential associated with the cosmic microwave background blackbody spectrum [23], the extragalactic gammaray background [24], the reionization history of the universe [25], the formation of small scale structures [26, 27] and the synthesis of light elements in the early Universe [28, 29] (see, for implications of neutral particle, Ref. [30]). Anomalies in the above mentioned quantities, however, could hardly be considered smoking gun evidence that the collider CHAMPs were indeed related to lighter, superweakly interacting dark matter.
If CHAMPs were stable on collider timescales, but featured a short lifetime on cosmological timescales, say on the order of a year or less, CHAMPs produced in colliders might be trapped in large water tanks surrounding the detectors [31]. The tanks would then be periodically drained to underground reservoirs where CHAMP decays might be observed in lowbackground conditions [31]. While certainly not straightforward experimentally, such a technique might provide important information first on the actual metastability of the charged species, and, secondly on the nature of the superweaklyinteracting particle the CHAMP would decay into. However, even if CHAMP decays were actually observed, this would still not suffice as conclusive evidence that the elusive particle CHAMPs decay into is indeed the dark matter constituent.
To our knowledge, beyond highenergy collider experiments, the only direct experimental handle on a superweakly interacting dark matter particle featuring a heavier, metastable charged partner is CHAMP pair production via neutrino–nucleon collisions, followed by direct observation at neutrino telescopes. This idea, originally proposed in Ref. [32], relies on the fact that the energy losses of CHAMPs in Earth are significantly smaller than those of muons, therefore CHAMP pairs (unlike muon pairs) can reach the detector even if they were produced far away. This makes the relevant target volume for neutrino–nucleon interactions much larger. The CHAMP pairs can be efficiently separated from muon pairs, due to large track separations in the detector. The original proposal was subsequently followed up by related studies [33, 34, 35, 36, 37, 38], which focused on the specific case of a gravitino lightest supersymmetric particle (LSP), and a stau NLSP playing the role of the CHAMP. Among other aspects, these studies investigated in detail stau energy losses in Earth and in the detector, computed expected event rates for a few sample models and the relevant background, and addressed the possibility of discriminating singlestau from singlemuon events.
Given the steady progress in the deployment of nextgeneration kmsize neutrino telescopes—particularly IceCube at the South Pole, already under construction and taking data—we consider it timely to address in detail a few points relevant for improving our understanding of the prospects for implementing the above outlined technique. In particular, as background rejection is not a substantial issue, the crucial point appears to be the evaluation of the CHAMP pair event rate at IceCube. To this end, we focus on the following four aspects:

Incoming neutrino flux: So far, all longlived CHAMP analyses for neutrino telescopes have considered a relatively optimistic flux of astrophysical neutrinos—as large as the WaxmanBahcall (WB) bound [39]—as the primary source for CHAMP pair production. However, the WB bound is only a theoretical upper limit on the flux of astrophysical neutrinos expected from optically thin sources; therefore, the absolute normalization as well as the spectral shape of the true astrophysical neutrino flux remain largely unknown (of course the biggest reason for this is that these neutrinos are as yet undetected!). On the other hand, atmospheric neutrinos have been detected and their flux is rather accurately known (below 100 TeV [40]). At larger energies ( TeV), although there is no detection so far, one can rather reliably extrapolate the flux of conventional atmospheric neutrinos from lower energies. In addition, one also expects a significant flux of socalled promptdecay atmospheric neutrinos, which originate from the decay of shortlived charmed mesons and feature a harder spectral index than the conventional component. While we again have only an upper limit on this prompt decay component, we know from particle physics that it necessarily exists, at some level, in the highenergy regime, and will be measured accurately by IceCube. In any event, atmospheric conventional and promptdecay neutrinos evidently contribute as well to CHAMP pair production in neutrino–nucleon collisions. In this paper, we study the role of these standard, guaranteed neutrino sources, and compare it to the contribution from astrophysical neutrino flux models including the WB bound.

Underlying particle physics model: The event rate depends not only on the incoming neutrino flux, but also on the nature of the assumed particle physics model. Here, we consider generic supersymmetric models featuring a gravitino LSP and stau NLSP, and study how the stau pair production cross section and event rates at neutrino telescopes depend on the given mass spectrum. Neutrino–nucleon interactions produce slepton–squark pair final states, the amplitudes mediated by supersymmetric fermion exchange (neutral and charged gauginos and higgsinos).

Other CHAMP sources: Unlike the highenergy neutrino flux, the flux of very energetic protons is accurately measured and wellknown up to extremely high energies; proton–nucleon interactions feature a supersymmetric pair production cross section significantly larger than that of neutrino–nucleon processes. It therefore seems reasonable to quantitatively assess the flux of stau pairs produced in the interaction of incident primary protons with nuclei in the atmosphere. The tradeoff is the enormous cross section for proton–nucleus scattering into Standard Model particles, which depletes the incoming proton flux, highly suppressing any stau pair production rate. Since quark–antiquark processes can directly produce stau pairs, however, the kinematic threshold for stau pair production, as a function of the incoming primary particle energy, is lower for proton–nucleus than for neutrino–nucleus processes. Yet we find that, for reasonable and phenomenologically acceptable particle models, the expected stau pair event rates from proton–nucleus collisions are experimentally negligible at all energies, leaving neutrinos as the only relevant primary source particles for CHAMP pair production.

Simplified analytic approach: We present the computation of the flux of staus from neutrino–nucleon interactions from first principles, and argue that, to an acceptable degree of accuracy, the total number of expected staus can be computed as one simple integral of three factors. Specifically, we show that the quantities of physical relevance are (1) the incident flux of primary neutrinos, (2) the ratio of the cross sections of neutrino plus nucleon into supersymmetric particle pairs to the total neutrino–nucleon cross section, and (3) a geometric efficiency factor.
Hereafter, we specifically use supersymmetric staus as charged metastable NLPs, but note that the following arguments are applicable to any other possible candidates of longlived CHAMPs. The outline of the remainder of the paper is as follows. We discuss the various components of the high energy neutrino flux and their connected uncertainties in Section 2. We introduce the new physics scenarios and compute the stau pair production cross sections in Section 3. Section 4 outlines the computation of the stau event rate at IceCube, including the abovementioned simplified analytic treatment. The stau rate dependence on the particle physics framework is addressed in Section 5, and we draw conclusions in Section 6.
2 Highenergy neutrino flux
In this Section, we summarize the highenergy neutrino fluxes we consider in the present study. As discussed in Section 1, past works considered only a flux of neutrinos close to, or saturating, the WB upper limit. However, other neutrino sources potentially contribute as well: these include conventional atmospheric neutrinos, promptdecay atmospheric neutrinos and possibly astrophysical neutrinos other than those considered in the WB setup. Since the stau production rate does not depend on neutrino flavor,^{2}^{2}2Flavor is conserved in the underlying supersymmetric particle pair production event, but all supersymmetry particles then decay promptly to a stau plus Standard Model particles. Thus, an electron neutrino eventually produces a stau pair just as a tau neutrino does. all we care is the total neutrino plus antineutrino flux, i.e., the flux of , where each here indicates neutrino plus antineutrino of flavor . From this point on, we mean this combined quantity when we use the term ‘flux’, unless otherwise stated. In Figure 1, we summarize and collect the various neutrino sources we discuss below, and the ranges of normalizations we consider.
2.1 Highenergy astrophysical neutrinos from extragalactic sources
Very powerful astrophysical objects such as active galactic nuclei (AGNs) and gammaray bursts (GRBs) are candidate highenergy neutrino sources. This is because strong gammaray emission detected from these objects can be attributed to the particle acceleration and successive interaction with the surrounding medium, magnetic and photon fields, which might also be a source of neutrinos via charged meson production.
If these neutrino sources are optically thin, then the upper bound on neutrino flux is obtained from the wellmeasured cosmicray flux, because each proton that arrives at Earth should produce no more than a few neutrinos at the source. Based on this argument and assuming that the cosmicray spectrum above GeV is of extragalactic origin, Waxman and Bahcall [39] derived the upper bound for the flux before neutrino oscillation to be – GeV cm s sr, where the range reflects cosmological evolution of source density. As we expect that flavor ratio at production (i.e., before oscillation) is :::2:0 due to meson decays, the WB bound summed over flavors after neutrino oscillation is – GeV cm s sr. Here we adopt GeV cm s sr for our reference value and show this bound in Figure 1.
We stress that although all previous studies adopted the WB upper limit on astrophysical neutrino flux as their incident source, that flux is not a model prediction, but an upper bound. One can never regard the output of a WB maximal neutrino flux as a solid prediction for the stau event rate at IceCube; the resulting stau flux is merely an upper limit. In addition, a few comments are in order on the WB bound’s robustness, especially at our energies of interest. First, the WB bound is valid for GeV, the threshold corresponding to proton energies of GeV (a daughter neutrino carries of its parent proton’s energy [39]). Below this, the WB bound is only an extrapolation, since the cosmicray flux below GeV is totally dominated by the galactic component. Second, by its definition the WB bound is applicable only to optically thin sources. If neutrinoemitting opaque objects existed, their contributions might sum to produce a neutrino flux exceeding the WB bound. Potential sources include baryonrich GRBs [42] and starburst galaxies [43].
2.2 Atmospheric conventional neutrinos
While the astrophysical neutrino flux is totally unknown, there is a guaranteed and wellmeasured neutrino component—atmospheric neutrinos. These arise from the decays of mesons produced by cosmic rays striking the upper atmosphere. Neutrinos coming from pion and kaon decays form the ‘conventional’ component, which is wellstudied both theoretically [44] and experimentally [40]. Although there is no detection of any neutrinos for GeV, the energy range we are mainly interested in, this component should quite easily be extrapolated using measured data at lower energies, thus providing guaranteed seeds for CHAMP production. We use the model of Ref. [45] for the conventional atmospheric flux. As shown in Figure 1, the wellknown spectrum of these neutrinos falls rather steeply with increasing energy.
2.3 Atmospheric promptdecay neutrinos
Hadronic cosmic ray interactions with Earth’s atmosphere also produce shortlived charmed mesons [46, 47, 48, 49, 50, 51, 52, 53, 54]. Even though branching ratios into these final states are not large, the neutrino spectrum from the subsequent decays of charmed mesons is quite hard as they immediately decay, before losing energy. As a consequence, the contribution from this ‘promptdecay’ component to the total flux of the atmospheric neutrinos falls less rapidly with energy than the conventional component. The absolute normalization of this prompt flux is however still unknown, and there is a large range in model predictions. The main source of uncertainty is the proper treatment of nexttoleading order charm meson production cross sections, which strongly depend on the behavior of the nucleon parton distribution functions (PDFs). In our study, we consider a range between a smaller prompt flux from Ref. [52] and a larger prompt flux obtained by using the shape presented in Ref. [53] and normalizing it to IceCube’s experimental upper limit [54]. Such a large flux, just allowed by the data, is in fact characteristic of the largest model predictions among Refs. [46, 47, 48, 49, 50, 51, 52, 53, 54]. We show these fluxes in Figure 1.
3 Interaction cross section
Several TeVscale extensions of the Standard Model feature a metastable massive charged particle. Perhaps the bestmotivated scenario from a theoretical standpoint is supersymmetry, which provides several examples. If the LSP is very weakly interacting—e.g., a gravitino or a right handed sneutrino, where the interactions with the rest of the supersymmetric partners are suppressed by gravitational couplings or a gauge symmetry—the NLSP is generally metastable. Specifically, a charged NLSP can occur in supergravity theories with a gravitino LSP [16, 55, 56], gauge mediated supersymmetry breaking setups [57, 58], scenarios featuring a stau–neutralino neardegeneracy (particularly in the socalled coannihilation region [59, 27]; here one can have a neutralino LSP), or supergravity scenarios with a righthanded sneutrino LSP [60].
Another TeVscale new physics setup that naturally encompasses a metastable NLP is that of universal extradimensions (UED) with a KaluzaKlein (KK) graviton as the lowest mass eigenstate in the KK tower [6]. In the minimal UED setup, the nexttolightest KK state is usually neutral, and corresponds to the KK first excitation of the U(1) gauge boson, , if the Higgs mass is below 200 GeV, for any value of the compactification inverse radius . However, as pointed out in Ref. [61], even in the minimal UED setup the nexttolightest KK particle (NLKP) can be charged, and the LKP can be the KK graviton if GeV and GeV. In this case, the NLKP corresponds to the KK first charged Higgs mode. Alternatively, the boundary conditions at the orbifold fixed points can alter the spectrum, and give rise to scenarios with a KK lepton as the NLKP, and again, a KK graviton LKP [6]. In general, electroweak precision observables constrain the scale where the first excitation of a 5dimensional UED scenario might be expected to a few hundred GeV, depending upon the value of the Higgs mass [62]. The analysis outlined below applies, with the appropriate production cross sections and energy losses, to such UED models and to any other similar framework featuring a metastable charged particle.
We choose to work with two of the wellknown and wellmotivated supersymmetric frameworks mentioned above: gauge mediated supersymmetry breaking (GMSB) [57, 58], and minimal supergravity (mSUGRA) with a gravitino LSP (see, e.g., Refs. [16, 55, 56]). For each framework we examine two models. One is a ‘supersymmetric benchmark’ (the SPS7 point of Ref. [63] for GMSB, and the model of Ref. [64] for mSUGRA with gravitino LSP), and the other is a variant with a lighter spectrum (models I and II). These are essentially rescalings of the first two. By adopting benchmark models, not only do we consider phenomenologically viable and theoretically soundlymotivated setups, but we also make it easier to compare the detection technique discussed here with a wealth of existing phenomenological analyses of the same models (see, e.g., Refs. [63, 64]). As will become apparent below, for the present analysis the details of the spectrum of the heavy supersymmetric particle pair produced is crucial. Assuming a degenerate sfermion spectrum, while potentially useful to get an understanding of the role of the mass scale in the expected size of the signal, can potentially be a misleading oversimplification.
mSUGRA Models  
I  280 GeV  10 GeV  11  0  
[64]  440 GeV  20 GeV  15  GeV  
GMSB Models  
II  70 TeV  35 TeV  15  3  
SPS7 [63]  80 TeV  40 TeV  15  3 
Model  

I  101 GeV  620 GeV  110 GeV  200 GeV 
[64]  153 GeV  940 GeV  180 GeV  340 GeV 
II  101 GeV  800 GeV  140 GeV  240 GeV 
SPS7 [63]  120 GeV  900 GeV  160 GeV  270 GeV 
We specify the mSUGRA and GMSB input parameters in Table 1. Notice in particular that Model II lies on the SPS7 slope defined in Ref. [63]. In Table 2 we detail the four models’ relevant particles masses.
The gravitino mass need not be specified as long as the stau decay length, , is larger than or of the same order as the Earth radius, . This implies a lower limit on the gravitino mass [65],
(1) 
Rearranged, the formula implies, for instance, that for a 100 GeV stau the gravitino mass can be as light as 1 MeV, and for a 1 TeV stau has to be larger than 0.3 GeV. Notice that Equation (1) does not take into account the relativistic boost factor, , which is sizable for staus produced in very high energy neutrino–nucleon interactions. It is therefore a conservative constraint on the gravitino mass.
In addition, the stau lifetime should be short enough to be consistent with limits obtained from the effects on light element abundances processed in bigbang nucleosynthesis [66] (see also Refs. [67]) and from excessive distortions to the cosmic microwave background spectrum [68]. Of particular relevance are constraints resulting from overproduction of Li and Li [28, 69, 70], induced by catalytic effects produced by bound states consisting of light nuclei and the metastable CHAMP (here, the lightest stau). Notice, however, that the details of the estimate of the amount of primordially synthesized lithium are still under debate [71]. Also, the constraints depend upon the fraction of electromagnetic energy released in the decay. Conservatively, if one requires the lifetime of the charged metastable species to be shorter than – s, as implied by the analysis of Ref. [72], the gravitino mass is constrained to be approximately below 1 GeV for a 100 GeV stau, and below 100 GeV for a 1 TeV stau. This evidently leaves a very wide window, of almost three orders of magnitude, for the viable gravitino mass range.
To calculate the stau flux,^{3}^{3}3Here the word ‘stau’ denotes staus and antistaus collectively. This is because both could be produced by the same interaction, but cannot be distinguished at neutrino telescopes. we first need to calculate the stau production cross section as a function of incoming neutrino energy. We do this using the susymadevent package [73, 74], which calculates the differential or total cross section for any scattering process in the MSSM given an SLHAconforming (standardized format for spectra) model input file [75]. We generate model input files using the suspect spectrum generator package [76], which also automatically checks the generated model against various known precision data constraints, such as . Our four models do not conflict with any known constraints. The processes contributing to stau production mainly stem from treelevel and exchange diagrams with a slepton ( and ) plus a squark ( and ) in the final state, through neutralino () or chargino () exchange, as discussed in Ref. [33]. More specifically, they are: , , where and indicate neutralinos and charginos, respectively. The squarks and sleptons produced in the final state promptly cascade decay into metastable staus.
We calculate total cross sections, summing over neutrino and antineutrino inelastic scattering on protons using exact matrix elements to supersymmetric pair final states for both charged current (CC) and neutral current (NC), employing CTEQ6L1 PDFs [77]. The possible final states are (anti)sneutrino or (anti)slepton, plus (anti)squark. Our calculations are leading order, as there are no available NLO QCD corrections for neutrino–proton () scattering. Judging from the known results for proton–proton () scattering, however, we probably make an underestimate of the rate on the order of . In this regard our calculation is conservative. For our purposes at the relevant energies, the neutrino–neutron cross section is sufficiently close to the neutrino–proton cross section that we can ignore calculating it separately.
An interesting side observation is that the and cross sections are not equal except at very large , where low quarks dominate the PDFs and are approximately egalitarian [78]; we reproduce this observation here. As gets within a couple orders of magnitude above threshold, however, dominates because the larger CC process picks out valence quarks, where dominates slightly over . Closer to threshold this remains true for the CC process, but for NC dominates because of chiralityselection for the finalstate squarks. Near threshold these two diverging components accidentally roughly cancel, resulting in once again. For the present purposes, however, we are concerned only with total rates, and ignore chargeseparated subsamples, which would vary somewhat as a function of stau energy.
The left panel of Figure 2 illustrates our results for the neutrino–proton scattering cross section into any supersymmetric particle pairs, for the four models listed in Table 1, as a function of the incident neutrino energy. The general trends in are consistent with those found in other analyses, see e.g., Ref. [33]: the steep rise in the lowenergy end reflects the strong kinematic suppression associated with the squark–slepton pair production threshold. The subsequent rise of the cross section with the incoming neutrino energy depends upon the small behavior of the PDFs. We give in Section 5, where we discuss the role of the specific supersymmetric particle spectrum in the determination of , an analytical interpretation of the specific powerlaw behavior that emerges from the numerical computation. Notice that compared to the optimistic toy model used in Ref. [33], theoreticallymotivated (optimistic) supersymmetric setups appear to give a maximal that is roughly one order of magnitude smaller in the asymptotic highenergy regime.
In the right panel of Figure 2, we show the ratio of the neutrino–nucleon cross section into supersymmetric particle pairs over the total neutrino–nucleon cross section (that is, as apparent from the figure, always close to the purely Standard Model cross section). As we explain in Section 4.2, this is the physical quantity of interest, in the limit of Earth as a thick target for high energy neutrinos, for the computation of the stau flux. Beyond threshold effects, we point out that in the energy range of interest the branching ratio of neutrino–nucleon interactions into supersymmetric particle pairs lies between and .
4 Flux of longlived staus
We devote this section to a detailed analytical treatment of the computation of the flux of staus produced by neutrino–nucleon collisions that might be detected at a neutrino telescope such as IceCube. We start, in Section 4.1, with a derivation of the differential flux of staus from first principles, leading to the result presented in Equation (6). In the following Section 4.2 we assume that Earth is opaque to neutrinos at energies relevant here. In this thick target approximation, we analytically show that the flux of staus at the detector can be computed as a simple integral, over incident neutrino energies, shown in Equation (10), of the product of three factors:

the differential flux of incident neutrinos (shown in Figure 1),

the ratio of the neutrino–nucleon cross section into supersymmetric particles over the total neutrino–nucleon cross section (shown in Figure 2, right), and

a ‘geometric efficiency’ factor, to be defined below and explicitly shown in the right panel of Figure 3.
When a neutrino interaction occurs, the branching ratio for stau production among the final states is given by the above cross section ratio. In the thicktarget approximation, all incoming neutrinos below the horizon will interact in Earth. Since the staus are collinear with the incoming neutrino direction, and are produced with a sizable fraction of the neutrino energy, this leads to a simple but important result about the stau flux. Neglecting stau energy losses in matter for a moment, we see that Earth acts as a neutrinotostau converter, with a probability that is independent of direction (below the horizon). This gives an upper bound on the stau flux through IceCube. The relevant range of energies will be between threshold ( GeV; see Figure 2) and the point at which the stau flux becomes too small ( GeV; note the product of Figures 1 and 2). Taking stau energy losses into account will only reduce the stau flux at the detector. Since staus cannot reach the detector from too far away, this means that only a limited range of nadir angles will be relevant, and this defines our geometric efficiency.
With our approximations, we give a preliminary assessment of the total expected stau flux at the detector (Table 3). Finally, we compute the actual accurate stau flux resulting from the fullglory integration of Equation (6) in Section 4.3. We provide numerical estimates of the integrated flux in Table 4 (that fall within a factor 2 of the approximate results anticipated in Table 3), as well as the actual differential flux of staus from the different primary incident neutrino sources. Finally, we comment in Section 4.4 on the flux of staus predicted from nucleon–nucleon reactions.
4.1 Formulation
In the framework we consider here, all supersymmetric particles produced in neutrino–nucleon interactions promptly decay into the NLSPs—here, the stau, which is metastable: longlived enough to propagate through Earth (with energy loss) and reach the detector. Our objective is to calculate the spectrum of staus after this energy loss, following similar principles for the spectra of neutrinoinduced muons [79, 80].
The stau electromagnetic energy loss rate is given by [32, 35]
(2) 
where is the column depth of matter in units of gcm (i.e., density times distance), and the and terms represent ionization and radiation losses, respectively. We neglect discrete scattering by weak interactions, as the effect is small at the energies we focus on near threshold [36]. Hereafter, we use the subscript to indicate quantities referring to staus. Our coordinate system locates the detector at and particles are produced at , so that the energy is a growing function of ( is positive). For the density profile of Earth, we use the model given in Ref. [81].
The ionization coefficient for staus, , is approximately the same as that of muons [35]; specifically, we use GeV cm g. Radiative losses, on the other hand, depend on particle mass, and we take the corresponding coefficient to be given by cm g [35]. By integrating Equation (2), the ‘distance’ traversed by the stau while its energy decreases from to reads
(3) 
We assume that each stau produced in interactions carries a large fraction of the parent neutrino energy. We denote this as and assume , independent of neutrino energy. (We discuss in Section 4.2 below the dependence of the stau flux on the parameter .) At these high laboratory energies, the staus may be taken to be collinear with the original neutrino direction. The differential flux of produced staus per energy and distance reads
(4) 
where . The exponential factor takes into account the neutrino attenuation, mainly due to the Standard Model interactions (; see Figure 2, right); represents the value of corresponding to the surface of Earth, which depends on the direction. The overall factor comes from the change of variables from the spectrum of neutrinos to that of produced staus, i.e., .
The spectrum of staus at the detector is given by a double integral over all production positions and energies, subject to the constraint of having the stau energy be between and at :
(5) 
where the function is defined by the energy loss, Equation (3). We use the energy constraint to perform the integration, and as a result we obtain the following expression:
(6)  
where and inside the integral have to be evaluated according to the chosen on the lefthand side and the at that step inside the integral. From the energyloss equation, Equation (3), we have
(7) 
where, again from the kinematic definition, .
It is convenient to change the integration variable by dividing the differential by and multiplying the integrand by . Then the integration steps are in , and in Figure 4 we show this new integrand for different nadir angles. In the lefthand panel, the stau energy losses in matter are neglected, so that the neutrino interaction and geometric effects are shown clearly. For each nadir angle ( is at the horizon, and is along the diameter of Earth), the sharp edges at large occur due to the boundary of Earth. The peaks arise because the neutrino interactions occur logarithmically near . Substantial attenuation of the neutrino flux is only seen at small nadir angles; in units of the axis, the exponential scale height is . In the righthand panel, the stau energy losses are now included. As expected, this prevents staus from arriving at the detector from too large of distances (beyond –0.05 in units of the axis, depending on nadir angle due to the radial variation of the density profile). In this panel, the visual area under each curve shows its relative importance to the total stau flux through the detector.
4.2 Thick target approximation
For neutrino energies relevant for our purposes, GeV, Earth is opaque to neutrinos at the most important nadir angles. Thus, to a good approximation, we can assume that staus are produced from neutrino interactions logarithmically near Earth’s surface. In this case, the exponential factor in equation (6) is sharply peaked at , i.e., . For our change of variables, the integrand is proportional to .
We require the staus to be more energetic than a given detector threshold . If the staus are relativistic enough to emit Čerenkov light, they will be detected. We adopt GeV, since the typical stau mass we consider is 100 GeV, but we note that our results change negligibly for higher or lower thresholds. For instance, a shift by one order of magnitude, TeV, affects the final result by only .
To satisfy the detection requirement , the initial stau and corresponding neutrino energies must be larger than some minima, and , given by
(9) 
as evident from Equation (7). We plot this minimal neutrino energy as a function of the nadir angle in the left panel of Figure 3 for the various supersymmetry models under consideration. Note that the crucial quantity here is the stau mass, hence we obtain the same result for models I and II. In the thick target approximation, the flux of staus reaching the detector from below is thus given by
(10)  
where is the step function, and in the second equality, we have simply used definitions given above. In the last equality, we are defining a ‘geometric efficiency’ factor , assuming that the incident neutrino intensity is isotropic. We show this geometric efficiency factor in the right panel of Figure 3. Thus, under the thick target approximation, the detection flux of staus can be divided into three independent factors:

incident neutrino spectrum

cross section ratio

geometric efficiency
These three factors are illustrated in figures 1–3. After integrating over the neutrino energy in Equation (10), we obtain approximate stau fluxes at the detector, summarized in Table 3 for the ranges of incident neutrino fluxes given in Section 2.
Model  WB bound  WB GRB  Atm. Prompt  Atm. Conv. 

I  0.20  0.012–0.73  0.0023–0.0038  
II  0.066  0.0028–0.16  0.00028–0.00048  
SPS7  0.045  0.0018–0.099  0.00014–0.00024  
0.034  0.0011–0.062  0.000069–0.00012 
So far we worked with the assumption that ; i.e., each stau carries half of the incident neutrino energy. However, especially when including the complex pattern of chain decay of squarks and sleptons into staus, a smaller fraction of the maximally available energy is expected to be carried by the staus, implying a larger value for . The precise value for is modeldependent, and its detailed evaluation is beyond the scope of the present analysis. We thus investigate the dependence of the stau flux on the parameter, for simplicity under the thick target approximation. Our calculations show that the flux of staus, with an incoming flux saturating the WB bound, and with our benchmark supersymmetric model I, is 3.2, 2.7, and 1.5 km yr for , 0.7, and 0.95, respectively. This shows that the stau flux changes at most a factor of 2 for a wide range of , which is well within other model uncertainties. The weak dependence we find stems from the fact that a larger value of requires larger neutrino energies to produce staus with a certain energy. In turn, at larger energies the incident flux is smaller, but the cross section ratio is larger. The cancellation of these two effects results in the mild dependence we find. The same argument also applies when we evaluate the flux more accurately in the next subsection, where we thus again use the assumption , for simplicity.
4.3 Results of numerical integration: stau spectrum
We now solve Equation (6) numerically to obtain a more precise spectrum and rate estimate of staus at the detector, as well as to check that the approximation made in the previous subsection is reasonable. Before giving the final flux estimates, we start by investigating the generic structure of the integrand of Equation (6).
Figure 4 shows the integrand of Equation (6) as a function of the column depth for various values of the nadir angle , assuming GeV. For definiteness, we show here the result only for a neutrino injection model saturating the WB upper limit. From the right panel of this figure, one can see that the stau events are dominated by directions with large enough nadir angle, specifically with , so that the staus can reach the detector. On the other hand, for very large nadir angles, the contributions to the total event rate are modest, because Earth is not completely opaque in those directions, even at these high energies.
Figure 5 shows the differential stau fluxes as a function of final stau energy. First, we note that nearly all of the detectable staus are well above the threshold required to be relativistic, which is a stau energy comparable to the stau mass. Only relativistic charged particles produce the Čerenkov light that IceCube can measure.
The energy loss of relativistic particles may be dominated by ionization or radiation, depending on whether the term or the term dominates in Equation (2), respectively. This transition for muons occurs at an energy GeV. For staus, it occurs at an energy a factor higher, i.e., at least GeV. The energy loss associated with Čerenkov radiation is always negligible; on the other hand, the Čerenkov radiation per unit length is the same for all relativistic particles. Thus, for particles at any energy in the relativistic ionizationdominated regime, all tracks will look the same in IceCube.
At higher energies, in the radiation dominated regime, there is additional Čerenkov radiation arising from relativistic electrons and positrons created in hard radiative processes. In this regime, one can indeed tell the energy of the primary particle by the intensity of the total Čerenkov radiation. For most of the relevant final stau energies shown in Figure 1, the staus will at most be only slightly in the radiative regime, and so all stau tracks going through IceCube will be indistinguishable from each other (and from lowenergy muons). While the total energy deposited in the detector is much smaller for staus than it is for lowenergy muons, this is irrelevant for IceCube, which detects only the Čerenkov light.
We argue that lowenergy but relativistic stau pairs could also be detected (albeit without energy measurement), giving a sizable event rate. Recall that while it might be difficult to distinguish between staus and muons on the basis of a singleparticle detection, it would still be possible if we use dualtrack events: since staus propagate over much longer distances in Earth than muons, tracks entering the detector simultaneously are expected to be wellseparated [32, 33, 34]. The careful analysis of Ref. [34] shows that the separation distribution of stau pairs ranges widely from 50 m to 1 km, and that it peaks around 500 m. On the other hand, the separation distribution for dimuons—the main background for the stau pair track—peaks at around 10 m, and essentially no dimuon events with m separation are expected. Therefore, this criterion rejects almost all background dimuon events, but would capture a large fraction of stau events (typically ). Thus, we are interested in the stau flux integrated over energies larger than the relativistic threshold—which is very small, GeV—and where our results are least sensitive. This is clear from Figure 5.
Table 4 shows the expected stau flux at the detector, obtained by integrating the spectrum above 300 GeV. Comparing with the results of Table 3, where we used the thick target approximation, we find the two results are consistent with each other within a factor of 2, justifying the reliability of the thick target approximation.
Model  WB bound  WB GRB  Atm. Prompt  Atm. Conv. 

I  0.13  0.0070–0.41  0.0010–0.0018  
II  0.049  0.0019–0.11  0.00017–0.00029  
SPS7  0.035  0.0012–0.070  0.000091–0.00015  
0.027  0.00080–0.047  0.000049–0.000082 
Note that the flux of conventional atmospheric neutrinos is not totally isotropic, but peaks in the horizontal direction by almost one order of magnitude compared to other directions [45]. (The prompt flux is isotropic.) This is an important effect because most of the staus reaching the detector arose from neutrinos from horizontal directions (Figure 4). Therefore, our results in Tables 3 and 4 for atmospheric neutrinos would be larger by a factor of 3, as these results were obtained with a directionaveraged incident neutrino flux.
As a consequence, given that the predicted stau flux could be as large as 1 km yr and that one expects essentially no background from muon pair events, the search for stau pair tracks is warranted in the actual data. Our discussion here shows that a significant fraction of stau events could possibly come from atmospheric promptdecay neutrinos (if the flux is close to the current upper bound), regardless of the assumed supersymmetry models.
4.4 Stau flux from nucleon–nucleon collisions
The stau pair flux produced from a differential flux of primary highenergy nucleons^{4}^{4}4Secondary nucleons and other hadrons contribute a small fraction of the incoming flux, and taking them into account does not affect our conclusions. colliding with atmospheric nuclei is given by
(11) 
because the thick target approximation (introduced in Section 4.2) is very good for interaction. The symbols and indicate the total nucleon–nucleon interaction cross section and the cross section into any supersymmetric particle pair, respectively. To reiterate, since direct decays into gravitinos are strongly suppressed by gravitational couplings, all final state parity odd particles decay into the NLSP, i.e., lightest stau pairs. For the total cross section, we assume the parameterization for the hadronair total cross section [82]
(12) 
where is the average number of nucleons in a nucleus of air. We approximate the nucleon–nucleon cross section into supersymmetric particles with the proton–proton cross section, and we compute the latter using Prospino2.0 [83]. For the incoming nucleon flux we use the estimate in figure 1 of Ref. [82].
Quark–antiquark processes can produce stau pairs directly, unlike neutrino–quark processes, where the final state has a larger threshold as the final state must contain a typicallyheavier squark. Hence the kinematic threshold, as a function of the primary particle energy, is lower in nucleon–nucleon collisions than in neutrino–nucleon collisions. Also, the subsequent occurrence of various supersymmetric particle thresholds at larger and larger masses, including particles featuring large degeneracy factors (such as squarks), implies a more rapidly growing behavior for the cross section than that for . As a consequence, we find that for the models under consideration here, is almost constant over several orders of magnitude in .
The final flux of staus is, however, dramatically suppressed by the ratio [84], even taking into account multiplicity effects in stau pair production or proton reinteractions. In particular, for the models we consider here, where the strongly interacting supersymmetric particles are typically much more massive than the NLSP, the combination of threshold effects and of the rapidly decreasing flux of incident primary protons leads to dramatically less optimistic predictions than those recently reported in Ref. [37]. There, the authors considered squarks and gluinos with extremely low masses (150 and 300 GeV), while the theoreticallymotivated benchmark models we use here feature squark and gluino masses between 600 GeV and 1 TeV. While we agree with the numerical results reported in Ref. [37] when making the same assumptions on squark and gluino masses, we obtain much lower figures for the benchmark models we adopt here. Namely, we find that for the two most optimistic models, I and II, we predict a stau flux in IceCube from nucleon–nucleon interactions of and per km per year, respectively—much lower than even the contribution from conventional atmospheric neutrinos. We believe that this relative smallness compared with that from the atmospheric incident neutrinos would be a rather modelindependent feature.
5 Role of the supersymmetric particle spectrum
In the previous sections we focused on specific supersymmetric models. We now wish to address the modelindependent question of how the stau pair rate at neutrino telescopes depends upon the supersymmetric particle masses. As we already pointed out, the stau pair flux depends upon the cross section: the relevant masses entering the cross section are those of the heavy pair produced, and those of the supersymmetric partners of the electroweak gauge bosons which mediate the charged and neutralcurrent interactions responsible for slepton–squark production. In addition, the number of produced stau pairs as a function of incoming neutrino energy crucially depends on the final state kinematic threshold: the larger the latter, the smaller the flux of incoming neutrinos that can lead to stau pair production, hence a smaller expected stau pair rate.
To explore quantitatively the statements above, we employ a convenient phenomenological parameterization of the supersymmetric setup at the lowenergy scale, and no longer rely on a specific supersymmetry breaking framework. The left panel of Figure 6 shows the stau production cross section variation in highenergy neutrino–proton collisions in the plane defined by the slepton (xaxis) and squark (yaxis) masses. For definiteness, we fix the relevant gaugino (respectively, bino and wino) masses to TeV and TeV. We choose TeV, GeV, , set all trilinear scalar couplings to zero, and further assume that the soft supersymmetry breaking scalar masses of sleptons and (separately) of squarks are degenerate, CPconserving and flavordiagonal. In the figure, we show isolevel curves at fixed values of the production cross section for staus at an incident neutrino energy of GeV, in order to avoid production threshold effects. As the figure illustrates, cross section variation is very mild well above threshold. Quantitatively, the effect varies within little more than one order of magnitude for scalar masses varied between 100 and 1000 GeV.
In the inset, we show how the cross section scales the masses of channel supersymmetric particles (neutralinos and charginos) exchanged in squark–slepton pair production from neutrino–proton collisions. Neutralino and chargino masses are entirely determined at tree level by the gaugino soft supersymmetry breaking masses and and the higgsino mass term , and by . For simplicity we assume a common ‘gaugino/higgsino’ mass scale, , defined as the common value of . We employ a common scalar mass for both sleptons and squarks of 300 GeV, and set all other lowscale supersymmetric parameters as in the rest of the figure ( GeV, , all trilinear scalar couplings zero). As we illustrate in the inset, beyond the scalar mass scale () where kinematic effects play a nontrivial role, cross section scaling goes like the gaugino/higgsino mass scale to the power . This can be analytically understood, since
(13)  
where, in Equation (13), we made use of the fact that in the large regime, [85]. One thus gets:
(14) 
which explains both the scaling in the inset of Figure 6 and of in Figure 2.
As the supersymmetric particle spectrum gets heavier, not only does the neutrino–proton cross section become smaller (left panel of Figure 6), but more importantly the shift in incident neutrino energy threshold for stau production strongly suppresses the final stau flux. This depends on the dramatic dependence of the flux of incident neutrinos on energy, as illustrated in our Figure 1. The right panel of Figure 6 quantifies this trend. There we show the relative contribution to the total stau flux from various incident neutrino fluxes as a function of the common squark and slepton masses . The fluxes are normalized to that resulting from the incident WB upper limit of extragalactic neutrino flux and GeV. We set all the supersymmetric parameters to the same values as in the left panel. Comparing the relative flux for various origins summarized in Table 4 with that shown in the right panel of Figure 6, we can roughly estimate that the four benchmark models correspond to GeV in the context of this phenomenological approach. In addition, recalling that the benchmark models predict stau flux of 1 km yr for WB neutrino bound, this suggests that one can expect a stau flux well above 1 km yr, provided slepton and squark masses are smaller than 500 GeV.
The different neutrino flux scaling with energy dictates that the relative importance of the various neutrino sources depends on the mass scale of the particles produced in the neutrino–proton collision. In particular, while with a very light spectrum the contribution from conventional atmospheric neutrinos can be comparable to (or even larger than) the extragalactic neutrino component, for TeV the contribution from atmospheric neutrinos becomes negligible. Note that in that mass range the supersymmetric particles are so heavy that the overall stau flux is extremely suppressed, and likely undetectable. On the other hand, the figure illustrates that, in principle, prompt decay neutrinos can be the dominant source of staus for almost any value of if the extragalactic neutrino flux is close to the GRBderived range (dashed red line in the figure) rather than the WB upper limit. In addition, even conventional neutrinos contribute a stau flux of the same order of magnitude as that expected from the WB upper limit on astrophysical neutrinos, as long as the supersymmetric scalars mass scale is below 200 GeV. If both the promptdecay neutrino flux and the extragalactic neutrino flux are maximal, prompt neutrinos contribute at the same level as extragalactic neutrinos for GeV, and dominate for below GeV. Since a detectable signal is expected only for a light supersymmetric spectrum, this leads us to the following prediction: if the signal discussed here is indeed detected, a very sizable fraction of it will originate from conventional and promptdecay neutrinos.
6 Conclusions
We reassessed the flux of metastable staus produced by neutrino–nucleon and nucleon–nucleon interactions that might be detectable at km neutrino telescopes. We derived the flux of staus from first principles, and showed that, under the approximation that Earth is opaque to very high energy neutrinos, the number of staus at the detector is given by a simple integral over the neutrino energy of the product of three factors: the incident neutrino flux, the ratio of the neutrino–nucleon cross section into supersymmetric particle pairs over the total neutrino–nucleon cross section, and a geometric efficiency factor. We showed that this approximation reproduces an exact numerical computation within a factor 2, which in turn is much better than the level of our knowledge of the first two factors entering the stau flux computation—namely, the incident neutrino flux and the features of the supersymmetric particle setup.
We focused on each of the factors relevant to the computation of the final flux of staus. We concentrated on four wellmotivated supersymmetric benchmark models, and independently evaluated the relevant production cross sections. We pointed out that previouslyneglected atmospheric neutrinos from prompt charmed meson decays could give a potentially large stau flux, even in the absence of a (yet to be discovered) astrophysical highenergy neutrino flux. This will depend on the prompt atmospheric neutrino flux being at the upper end of the theoretically expected range. Nucleon–nucleon processes, even for the most optimistic benchmark models, would not contribute sizably to the final stau flux. Finally, we numerically and analytically studied how the relevant cross sections depend on the supersymmetric model mass spectrum, and how the relative importance of primary neutrino sources depends on the mass scale of supersymmetric scalars. In particular, we predict that if the signal discussed here is indeed detected, a very sizable fraction of it would originate from conventional and promptdecay neutrinos.
Acknowledgments
This work was supported by the Sherman Fairchild Foundation (SA); CCAPP, The Ohio State University, and NSF CAREER grant PHY0547102 (JFB); DoE grants DEFG0392ER40701, DEFG0205ER41361, and NASA grant NNG05GF69G (SP); and DoE grant DEFG0291ER40685 (DR). We would like to thank Francis Halzen, Chris Quigg, and Xerxes Tata for enlightening discussions and Markus Ahlers for useful correspondence.
References
References
 [1] M. Fukugita, K. Ichikawa, M. Kawasaki and O. Lahav, Phys. Rev. D 74, 027302 (2006).
 [2] D. N. Spergel et al. [WMAP Collaboration], Astrophys. J. Suppl. 170, 377 (2007).

[3]
L. Bergstrom,
Rept. Prog. Phys. 63, 793 (2000);
G. Bertone, D. Hooper and J. Silk, Phys. Rept. 405, 279 (2005).  [4] G. Jungman, M. Kamionkowski and K. Griest, Phys. Rept. 267, 195 (1996).
 [5] H. C. Cheng, J. L. Feng and K. T. Matchev, Phys. Rev. Lett. 89, 211301 (2002).
 [6] D. Hooper and S. Profumo, Phys. Rept. 453, 29 (2007).
 [7] S. Wolfram, Phys. Lett. B 82, 65 (1979).
 [8] C. Munoz, Int. J. Mod. Phys. A 19, 3093 (2004).
 [9] L. Bergstrom, New Astron. Rev. 42, 245 (1998).
 [10] M. Drees, Y. G. Kim, M. M. Nojiri, D. Toya, K. Hasuko and T. Kobayashi, Phys. Rev. D 63, 035008 (2001).
 [11] H. Baer, C. Balazs, A. Belyaev, T. Krupovnickas and X. Tata, JHEP 0306, 054 (2003).
 [12] G. Weiglein et al. [LHC/LC Study Group], Phys. Rept. 426, 47 (2006).
 [13] A. Masiero, S. Profumo and P. Ullio, Nucl. Phys. B 712, 86 (2005).
 [14] E. A. Baltz, M. Battaglia, M. E. Peskin and T. Wizansky, Phys. Rev. D 74, 103521 (2006).
 [15] H. Baer and X. Tata, “Weak Scale Supersymmetry: From Superfields to Scattering Events”, Cambridge University Press, 2006.
 [16] G. R. Blumenthal, H. Pagels and J. R. Primack, Nature 299, 37 (1982).
 [17] J. L. Feng, A. Rajaraman and F. Takayama, Phys. Rev. D 68, 085018 (2003).
 [18] J. L. Feng, A. Rajaraman and F. Takayama, Phys. Rev. Lett. 91, 011302 (2003).
 [19] F. Cerutti et al. [LEP2 SUSY Working Group], LEPSUSYWG/0205.1, http://lepsusy.web.cern.ch/lepsusy/www/stablesummer02/stable208.html
 [20] D. Acosta et al. [CDF Collaboration], Phys. Rev. Lett. 90, 131801 (2003).
 [21] M. Fairbairn et al, Phys. Rept. 438, 1 (2007).
 [22] See e.g. X. L. Chen and M. Kamionkowski, Phys. Rev. D 70, 043502 (2004).
 [23] L. Zhang, X. Chen, M. Kamionkowski, Z. g. Si and Z. Zheng, Phys. Rev. D 76, 061301 (2007)
 [24] J. A. R. Cembranos, J. L. Feng and L. E. Strigari, arXiv:0704.1658 [astroph].
 [25] E. Pierpaoli, Phys. Rev. Lett. 92, 031301 (2004).
 [26] K. Sigurdson and M. Kamionkowski, Phys. Rev. Lett. 92, 171302 (2004).
 [27] S. Profumo, K. Sigurdson, P. Ullio and M. Kamionkowski, Phys. Rev. D 71, 023518 (2005).
 [28] See,e.g., K. Kohri and F. Takayama, Phys. Rev. D 76, 063507 (2007); M. Kawasaki, K. Kohri and T. Moroi, Phys. Lett. B 649, 436 (2007); T. Jittoh, K. Kohri, M. Koike, J. Sato, T. Shimomura and M. Yamanaka, arXiv:0704.2914 [hepph].
 [29] K. Jedamzik, arXiv:0707.2070 [astroph].
 [30] D. Cumberbatch, K. Ichikawa, M. Kawasaki, K. Kohri, J. Silk and G. D. Starkman, arXiv:0708.0095 [astroph]; and references therein.
 [31] J. L. Feng and B. T. Smith, Phys. Rev. D 71, 015004 (2005) [Erratumibid. 71, 0109904 (2005)].
 [32] I. Albuquerque, G. Burdman and Z. Chacko, Phys. Rev. Lett. 92, 221802 (2004).
 [33] M. Ahlers, J. Kersten and A. Ringwald, JCAP 0607, 005 (2006).
 [34] I. F. M. Albuquerque, G. Burdman and Z. Chacko, Phys. Rev. D 75, 035006 (2007).
 [35] M. H. Reno, I. Sarcevic and S. Su, Astropart. Phys. 24, 107 (2005).
 [36] Y. Huang, M. H. Reno, I. Sarcevic and J. Uscinski, Phys. Rev. D 74, 115009 (2006).
 [37] M. Ahlers, J. I. Illana, M. Masip and D. Meloni, JCAP 0708, 008 (2007).
 [38] M. H. Reno, I. Sarcevic and J. Uscinski, arXiv:0710.4954 [hepph].
 [39] E. Waxman and J. N. Bahcall, Phys. Rev. D 59, 023002 (1999).

[40]
K. Daum et al. [Frejus Collaboration.],
Z. Phys. C 66, 417 (1995);
J. Ahrens et al. [AMANDA Collaboration], Phys. Rev. D 66, 012005 (2002);
Y. Ashie et al. [SuperKamiokande Collaboration], Phys. Rev. D 71, 112005 (2005).  [41] E. Waxman and J. N. Bahcall, Phys. Rev. Lett. 78, 2292 (1997).
 [42] S. Ando and J. F. Beacom, Phys. Rev. Lett. 95, 061103 (2005).

[43]
A. Loeb and E. Waxman,
JCAP 0605, 003 (2006);
F. W. Stecker, Astropart. Phys. 26, 398 (2007);
T. A. Thompson, E. Quataert, E. Waxman and A. Loeb, astroph/0608699. 
[44]
G. D. Barr, T. K. Gaisser, P. Lipari, S. Robbins and T. Stanev,
Phys. Rev. D 70, 023006 (2004);
G. D. Barr, T. K. Gaisser, S. Robbins and T. Stanev, Phys. Rev. D 74, 094009 (2006);
M. Honda, T. Kajita, K. Kasahara and S. Midorikawa, Phys. Rev. D 70, 043008 (2004);
M. Honda, T. Kajita, K. Kasahara, S. Midorikawa and T. Sanuki, Phys. Rev. D 75, 043006 (2007).  [45] J. Candia and E. Roulet, JCAP 0309, 005 (2003).
 [46] E. Zas, F. Halzen and R. A. Vazquez, Astropart. Phys. 1, 297 (1993).
 [47] L. Pasquali, M. H. Reno and I. Sarcevic, Phys. Rev. D 59, 034020 (1999).
 [48] G. Gelmini, P. Gondolo and G. Varieschi, Phys. Rev. D 61, 056011 (2000).
 [49] C. G. S. Costa, Astropart. Phys. 16, 193 (2001).
 [50] L. V. Volkova and G. T. Zatsepin, Phys. Atom. Nucl. 64, 266 (2001) [Yad. Fiz. 64, 313 (2001)].
 [51] G. Fiorentini, V. A. Naumov and F. L. Villante, Phys. Lett. B 510, 173 (2001).
 [52] J. F. Beacom and J. Candia, JCAP 0411, 009 (2004).
 [53] A. D. Martin, M. G. Ryskin and A. M. Stasto, Acta Phys. Polon. B 34, 3273 (2003).
 [54] A. Achterberg [IceCube Collaboration], astroph/0611597.
 [55] J. L. Feng, S. Su and F. Takayama, Phys. Rev. D 70, 075019 (2004).

[56]
J. R. Ellis, K. A. Olive, Y. Santoso and V. C. Spanos,
Phys. Lett. B 588, 7 (2004);
J. R. Ellis, A. R. Raklev and O. K. Oye, JHEP 0610, 061 (2006);
A. Ibarra and S. Roy, JHEP 0705, 059 (2007). 
[57]
D. A. Dicus, B. Dutta and S. Nandi,
Phys. Rev. D 56, 5748 (1997);
K. Cheung, D. A. Dicus, B. Dutta and S. Nandi, Phys. Rev. D 58, 015008 (1998);
J. L. Feng and T. Moroi, Phys. Rev. D 58, 035001 (1998);
P. G. Mercadante, J. K. Mizukoshi and H. Yamamoto, Phys. Rev. D 64, 015005 (2001).  [58] G. F. Giudice and R. Rattazzi, Phys. Rept. 322, 419 (1999).

[59]
S. Ambrosanio and B. Mele,
Phys. Rev. D 55, 1399 (1997)
[Erratumibid. 56, 3157 (1997)];
A. V. Gladyshev, D. I. Kazakov and M. G. Paucar, Mod. Phys. Lett. A 20, 3085 (2005);
T. Jittoh, J. Sato, T. Shimomura and M. Yamanaka, Phys. Rev. D 73, 055009 (2006). 
[60]
T. Asaka, K. Ishiwata and T. Moroi,
Phys. Rev. D 73, 051301 (2006);
S. K. Gupta, B. Mukhopadhyaya and S. K. Rai, Phys. Rev. D 75, 075007 (2007).  [61] J. A. R. Cembranos, J. L. Feng and L. E. Strigari, Phys. Rev. D 75, 036004 (2007).

[62]
T. Flacke, D. Hooper and J. MarchRussell,
Phys. Rev. D 73, 095002 (2006)
[Erratumibid. 74, 019902 (2006)];
I. Gogoladze and C. Macesanu, Phys. Rev. D 74, 093012 (2006).  [63] B. C. Allanach et al., in Proc. of the APS/DPF/DPB Summer Study on the Future of Particle Physics (Snowmass 2001) ed. N. Graf, In the Proceedings of APS / DPF / DPB Summer Study on the Future of Particle Physics (Snowmass 2001), Snowmass, Colorado, 30 Jun  21 Jul 2001, pp P125 [hepph/0202233].
 [64] A. De Roeck, J. R. Ellis, F. Gianotti, F. Moortgat, K. A. Olive and L. Pape, Eur. Phys. J. C 49, 1041 (2007) [arXiv:hepph/0508198].
 [65] J. R. Ellis, K. A. Olive, Y. Santoso and V. C. Spanos, Phys. Lett. B 588, 7 (2004).

[66]
R. H. Cyburt, J. R. Ellis, B. D. Fields and K. A. Olive,
Phys. Rev. D 67, 103521 (2003);
K. Jedamzik, Phys. Rev. Lett. 84, 3248 (2000); J. L. Feng, S. f. Su and F. Takayama, Phys. Rev. D 70, 063514 (2004); F. D. Steffen, JCAP 0609, 001 (2006).  [67] M. Kawasaki, K. Kohri and T. Moroi, Phys. Rev. D 71, 083502 (2005); M. Kawasaki, K. Kohri and T. Moroi, Phys. Lett. B 625, 7 (2005).
 [68] W. Hu and J. Silk, Phys. Rev. Lett. 70, 2661 (1993).
 [69] M. Pospelov, Phys. Rev. Lett. 98, 231301 (2007).
 [70] M. Kaplinghat and A. Rajaraman, Phys. Rev. D 74, 103004 (2006)
 [71] K. Jedamzik, arXiv:0707.2070 [astroph].
 [72] R. H. Cyburt, J. R. Ellis, B. D. Fields, K. A. Olive and V. C. Spanos, JCAP 0611, 014 (2006).
 [73] G. C. Cho, K. Hagiwara, J. Kanzaki, T. Plehn, D. Rainwater and T. Stelzer, Phys. Rev. D 73, 054002 (2006).
 [74] F. Maltoni and T. Stelzer, JHEP 0302, 027 (2003).

[75]
B. Allanach, P. Skands et al,
JHEP 0407, 036 (2003);
T. Hahn, hepph/0408283;
J. A. AguilarSaavedra et al., Eur. Phys. J. C 46, 43 (2006).  [76] A. Djouadi, J. L. Kneur and G. Moultaka, hepph/0211331.
 [77] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. Nadolsky and W. K. Tung, JHEP 0207, 012 (2002).
 [78] M. Carena, D. Choudhury, S. Lola and C. Quigg, Phys. Rev. D 58, 095003 (1998).
 [79] T. K. Gaisser and T. Stanev, Phys. Rev. D 31, 2770 (1985).
 [80] M. D. Kistler and J. F. Beacom, Phys. Rev. D 74, 063007 (2006).
 [81] R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Astropart. Phys. 5, 81 (1996).
 [82] J. I. Illana, M. Masip and D. Meloni, Phys. Rev. D 75, 055002 (2007).
 [83] W. Beenakker, R. Hopker and M. Spira, hepph/9611232.
 [84] M. Byrne, C. F. Kolda and P. Regan, Phys. Rev. D 66, 075007 (2002).
 [85] R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Phys. Rev. D 58, 093009 (1998).