New constraints on supersymmetry using neutrino telescopes

We demonstrate that megaton-mass neutrino telescopes are able to observe the signal from long-lived particles beyond the Standard Model, in particular the stau, the supersymmetric partner of the tau lepton. Its signature is an excess of charged particle tracks with horizontal arrival directions and energy deposits between 0.1 and 1 TeV inside the detector. We exploit this previously-overlooked signature to search for stau particles in the publicly available IceCube data. The data shows no evidence of physics beyond the Standard Model. We derive a new lower limit on the stau mass of 320 GeV (95% C.L.) and estimate that this new approach, when applied to the full data set available to the IceCube collaboration, will reach word-leading sensitivity to the stau mass ( m ˜ τ = 450 GeV).

GeV at 95% C.L. respectively [2]. For these limits in particular, the stau's mass is the only parameter of interest, due to the assumed Drell-Yan production.
Stau searches have also been proposed in the context of megaton-mass neutrino telescopes. Highly energetic cosmic particles (cosmic rays and neutrinos) colliding with nucleons in the Earth's atmosphere are capable of producing staus. These in turn would then appear as charged particles in the detectors. In particular, they would appear as charged particles propagating through the detector, so-called tracks. For such events, other particles pro-* Corresponding author.
ducing charged tracks would act as a background, mainly atmospheric muons and muons produced by neutrinos. Fig. 1 shows the relative fluxes at the surface. The orders of magnitude difference between the stau flux and background makes disentangling them difficult, even for a low stau mass of 100 GeV. The neutrino flux is divided into its primary contributors. These are the astrophysical flux, the conventional flux from the decay of π and K mesons and the prompt flux from the decay of heavier mesons.
One proposed strategy to disentangle staus from the background is to search for stau pairs. Due to the highly relativistic boosting, the two staus would move in the same direction, cross the detector simultaneously, and thus give rise to two parallel tracks -a smoking-gun signature. [5][6][7][8][9][10]. However, current telescopes can not distinguish these events from the overwhelming background, single tracks from cosmic ray muons (hN → μX) and/or from charged-current muon neutrino interactions (ν N → μX), unless the two stau tracks have a large separation (in Ice-Cube by ∼ 135 meters [11,12]). Thus, the majority (∼ 99.9 %) of the potential stau particles would go undetected.
In this letter we discuss how neutrino telescopes, in particular IceCube, can exploit a different signature to observe, on a statistical basis, a stau signal. At the energies of interest, staus are expected to be significantly more penetrating than muons of the same energy, because they essentially do not loose energy through stochastic processes [9]. Hence, for nearly horizontal trajectories, tens of kilometers of ice shield IceCube from the vast majority of muons but not from staus. We thus search for an excess of track events over the background expected from muons crossing the de-  Predicted cosmic ray flux, φ 0 (red), muon flux (φ μ (green), muon neutrino flux φ ν (dashed) and stau flux φτ (blue) at the surface above the IceCube detector at θ Zenith = 88 • . The neutrino flux is further split into its components, the astrophysical flux (orange), conventional flux (purple) and prompt flux (magenta). The stau mass was set here to 100 GeV and the interaction model Sibyll 2.3c [3] and primary model and H4a [4] were used. The low stau flux, when compared to muon and neutrino fluxes, makes searches using neutrino telescopes challenging.
tector horizontally at a zenith angles near ∼ 90 • . In contrast to previous works, this new analysis does not rely on an identifiable double track signature and thus less stringent event selection criteria. We demonstrate the potential strength of the method by analyzing one year of publicly available IceCube data [13,14]. In order to identify this stau component, we utilize differences in the corresponding energy and angular distributions compared to the contributions from muons.
Considering state-of-the-art cosmic-ray flux and hadronic interaction models, we calculate the stau signal and backgrounds using MCEq [15], a tool that solves the cascade equations that govern the interaction of cosmic rays and propagation of the resulting air shower. We modified MCEq to include the generation of staus with a probability of [16] P h where σ h,nucleon τ is the total stau production cross-section from the collisions of a hadron h with a nucleon in the atmosphere, σ h,air total is the total cross-section of h with air, and A = 14.6 is the average number of nucleons in a nucleus of air. The stau production cross-section has been computed assuming the Drell-Yan process through MadGraph [17,18]. More details about these calculations and corresponding assumptions are given in Appendix B. After production, staus are assumed to propagate straight through matter with an energy loss per column density, X , that can be approximated by: where a τ (E) accounts for ionization losses and b τ (E) for energy loss by pair production, hadronic interactions, and bremsstrahlung. Given that the ionization effects for stau are expected to be similar to those of a muon we assume a τ (E) ≈ a μ (E). In contrast, we expect all other effects to depend on the particle speed and hence to be in good approximation given by b τ (E) ≈ b μ (E)m μ /m τ [16].
Our backgrounds are muons that produce a detector response indistinguishable from that of the staus. The muon backgrounds can be divided into two components according to whether muons are produced by a hadronic interaction in a cosmic ray airshower or by a neutrino interacting in the Earth. We simulate these contributions similarly to what we do for the staus. In addition to the cosmic-ray flux and composition used before, we now also include the flux of astrophysical neutrinos measured by IceCube [22]: To calculate event rates from particle fluxes we rely on publicly available IceCube effective areas, as discussed in Appendix C.
The IceCube detector is capable of reconstructing the direction from which a low-energy muon or stau comes with high accuracy. Given the symmetry of the Earth and of the detector, the information of interest about the direction is typically expressed in terms of the angle with respect to the zenith. For the energy range of interest, E ∈ [100 GeV, 1 TeV], the resolution with which IceCube can reconstruct this angle is 1 • [23], and this is folded into our simulations by smearing the particle arrival directions accordingly. The energy released in the detector's active volume is reconstructed with a resolution of 20% for muon events with energy above 1 TeV. Below this threshold, the energy reconstruction is biased towards 700 GeV, see Fig. 2. The origin of this bias is related to the particles being minimally ionizing at these energies. For the energy estimation, current reconstruction methods rely on stochastic losses, which are negligible at low energies. Both the bias and resolution in the energy reconstruction are folded in our simulation. More details are given in Appendix D. Fig. 3 shows the expected angular distribution for signal and background events in the public IceCube data set. Staus are barely stopped by the Earth or the material surrounding the detector and hence their angular distribution appears flat. However, because the amount of material to be traversed increases rapidly towards the horizon, their distribution is slightly peaked at 86 • . The muon contribution due to hadron interactions is maximal at the zenith (i.e. at 0 • ) and decreases steeply for increasing angles as the muon-flux is attenuated by the ice-overburden over the detector. This contribution drops below the rate expected for staus towards the horizon, for angles above 84-85 • . On the contrary, the muon contribution due to neutrino interactions increases towards the horizon because that's where the flux of atmospheric neutrinos is largest. The distribution grows till about 90 • , after which it saturates. The opposite trend of the two background distributions creates a small range of angles in which the background is minimal and the sensitivity to a stau signal is maximal. This range is independent of the stau mass, as shown in Fig. 3. Higher masses scale the distribution of the staus, while leaving the shape unchanged.  The energy distribution of the staus and muon background due to neutrinos for angles between 85 and 90 • is shown in Fig. 4. As a consequence of the low rate of energy loss, staus deposit within the detector as much energy as the lowest-energy muons. Their distribution is sharply peaked as they cross the entire detector depositing always the same energy regardless of their kinetic energy. The full distribution is in the energy range for which the energy reconstruction is biased and the events are given a random energy estimation that is independent by their original energy distribution. After the energy reconstruction, there is a minimal difference between staus of different masses. The reconstructed energy distribution for the muon background partially overlaps with the stau distribution but has a tail extending towards high energies. The different energy distribution provides another handle to separate our signal from the background.
Based on our modeling, we identified the observable space in which the ratio between the stau signal and the muon-induced background is maximal. This space corresponds to angles between 85 and 90 • and energies between 0.1 and 1 TeV. The most efficient way to extract the signal would be to perform a bivariate analysis in energy and angle. Such an analysis should take into account the systematic uncertainties related to the detector response, in particular related to the bias in the energy reconstruction. Given that these uncertainties are not completely available in literature, we opted not to use in the analysis the full shape of the energy distribution but only to apply a loose cut to remove the highenergy part of the background events (E > 1 TeV). This approach reduces the sensitivity of the analysis but makes it more robust.
We applied these simple selection criteria (100 GeV < E < 1 TeV and 85 • < θ Zenith < 90 • ) to the data of IceCube [13,14]. It follows well the background-only distribution. The angular distribution for events with energy between 0.1 and 1 TeV is shown in Fig. 5. Also, in this case the data points follow well the background-only distribution. The p-value of our data given the background-only hypothesis is 0.1, supporting our intuition that there is no evidence for a signal. Assuming the currently leading experimental limit on the stau mass (mτ = 430 GeV), we expect to retain 8 stau events. To extract a limit on the mass of the staus we perform a binned likelihood fit of the data shown in Fig. 5. The fit has one free parameter, the mass of the stau particle. The rates of the background components are fixed by our modeling. The best fit results in 0 stau events. Inverting a standard frequentist hypothesis test [24], we set a lower bound on the stau mass of mτ > 320 GeV at 95% C.L. This lies approximately 10% above the expected limit from simulations.
This result does not include systematic uncertainties. Because of the way we constructed the analysis -i.e. by using only a loose cut on the energy and a bin size equal to the angular resolution of the experiment -we expect systematic uncertainties to be negligible, to first approximation. Future analysis based on larger data samples that exploit fully the energy information will need to refine the treatment of the systematic uncertainties.
Our limit of mτ > 320 GeV is the most stringent constraint on stau masses ever set by a non collider experiment. This proves the potential of our analysis approach. The analysis can be improved significantly by increasing the exposure (see Fig. 6). With the ten years of IceCube data that have been already recorded, and assuming no improvements to the reconstruction, simulation or energy resolution, our analysis approach would provide a sensitivity to stau masses of up to 450 GeV. Thus, neutrino telescopes have the chance to probe an unexplored parameter space, together with future LHC experiments but using a very different and complementary technique.
Furthermore, several neutrino telescopes, such as P-ONE [25], KM3NET [26], GVD [27] and IceeCube-Gen2 [28], are currently taking data or are in preparation. While their designs are different, they are all based on the detection of the Cherenkov light emitted Fig. 6. Projection of the sensitivity to long-lived staus at IceCube as a function of time, assuming Drell-Yan production. The red line corresponds to the current ATLAS limit [2]. The orange cross shows the here obtained limit, while the black line plots the expected one. The purple shaded region shows the 15% uncertainty region. by charged particles traveling through water or ice. The analysis performed in this work using data from IceCube can be extended to these upcoming datasets as well. Finally, this search strategy is not limited to the stau particle but can be utilized to search for other long-lived, charged particles beyond the Standard Model.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Atmospheric shower simulation
The particle interactions are modeled with Sibyll 2.3c [3], EPOS-LHC [29], QGSJET-II [30] and DPMJET-III [31]. For the cosmic ray models we use the Gaisser-Hillas models H3a and H4a [4] and Gaisser-Stanev-Tilav Gen 3 and 4 [32]. The resulting differences are used as an estimate of the uncertainties in our calculation. We found these to be negligible compared to other uncertainties in our analysis that will be discussed later on. We use the NRLMSISE-00 [33,34] model to simulate the atmosphere.

Appendix B. Stau production
We use the built in MSSM model in MadGraph. We furthermore adopt the parton distribution functions (PDF) from LHAPDF6 [35]. Concretely, we use the CT10nlo [36] PDF as well as the NNPDF30_nnlo_nf_5_pdfas from the NNPDF3.0 [37] PDF set.
This allows us to include the uncertainties due to different PDF parametrizations. For the interactions of other hadrons with air, we scale the p − p cross-sections following [16].

Appendix C. Effective areas
To calculate the contribution of neutrino-induced muons, we fold the 2D effective area, as a function of energy and declination, from [14] with the neutrino fluxes. The neutrino energy to muon energy mapping is approximated using the normalized 3D effective areas given in [38].
To make predictions for the stau component, we require a detector response to staus. We use the same approach as for the muons and divide the convolution by the total neutrino crosssection. The resulting efficiency includes effects of muon propagation in the ice. These we compensate by scaling the results, so at 1 TeV the effective area for muons corresponds to the spatial area of the detector, 10 6 m 2 . This results in a signal efficiency of 78%.

Appendix D. Energy reconstruction
We include the effects of energy reconstruction as described in [39,38] by constructing a function, mapping the true particle energy, E true , to the reconstructed energy, E reco , of the form Linear interpolation To get the specific values for σ and μ, we fitted the background prediction to the data for energies above 1 TeV. For such energies we do not expect any stau events. The results are σ 1 = 0.4, μ 1 = 700 GeV, σ 2 = 0.3 and μ = E true . The energy reconstruction distribution is shown in Figure (4). These values agree well with those shown in [39]. To map the stau energies to their reconstructed energies, we map their energy deposit to muon energies with an equivalent loss according to equation (2) and then proceed as with the muons.