R-parity Violating Supersymmetric Explanation of the Anomalous Events at ANITA

The ANITA balloon experiment has observed two EeV-energy, upgoing events originating from well below the horizon. This is puzzling, because (i) no Standard Model (SM) particle is expected to survive passage through Earth at such energies and incident angles, and (ii) no such events were reported by IceCube. In this paper, we address both these issues by invoking a beyond SM interpretation of the EeV events as due to the decay of a long-lived bino in the R-parity violating (RPV) supersymmetry. In particular, a TeV-scale slepton/squark can be resonantly produced through the interaction of the EeV neutrino with electrons/nucleons inside Earth, that decays to a light, long-lived bino, which survives the propagation through Earth matter before decaying back to neutrinos, leptons and/or quarks, thus producing upgoing air showers in the atmosphere near ANITA. We find that the ANITA events can be explained with a GeV-scale bino and ${\cal O}(0.1)$ RPV couplings, which are consistent with all existing high and low-energy constraints. We also find that an isotropic neutrino flux is inadequate for a beyond the SM explanation of this kind, and an anisotropic flux must be invoked. Finally, we also address the apparent tension of these observations with IceCube. Various aspects of our interpretation are testable in the near future at different frontiers, such as by the LHC, Belle II, ANITA-IV and IceCube.


I. INTRODUCTION
The Antarctic Impulsive Transient Antenna (ANITA) collaboration has recently reported two anomalous upward-going ultra-high energy cosmic ray (UHECR) air shower events with deposited shower energies of 0.6 ± 0.4 EeV and 0.56 +0. 3 −0.2 EeV, respectively [1,2]. Both events, one from ANITA-I [1] and another from ANITA-III [2], originate from well below the horizon, with elevation angles of (−27.4 ± 0.3) • and (−35.0 ± 0.3) • , respectively. They do not exhibit phase inversion due to Earth's geomagnetic effects -a primary characteristic of conventional downgoing UHECR air showers which produce downgoing radio impulses that are reflected off the Antarctic ice surface. Potential background events from anthropogenic radio signals that might mimic the UHECR characteristics, or unknown processes that might lead to a non-inverted polarity on reflection from the ice cap are estimated to be ≤ 0.015. This leads to 3σ evidence for the interpretation of the two anomalous events as due to direct upward-moving Earth-emergent UHECR-like air showers above the ice surface [2].
However, such an interpretation faces severe challenges within the known Standard Model (SM) framework, because no SM particle is expected to survive passage through Earth a chord distance of ∼ 7000 km (corresponding to the observed zenith angles of the two events) at EeV energies. In particular, the interpretation of these events as τ -lepton decay-induced air showers at or near the ice surface arising from a diffuse UHE flux of cosmic ν τ is strongly disfavored due to their mean interaction length of only ∼ 300 km. Even including the effect of ν τ regeneration [3][4][5][6], the resulting survival probability over the chord length of the ANITA events with energy greater than 0.1 EeV is < 10 −6 [7], largely due to τ -lepton energy loss inside Earth because of ionization, e + e − pair production, bremsstrahlung, and photonuclear interactions [8], thereby excluding the SM interpretation at 5.8σ confidence. A possible way out is by invoking significant suppression of the deep-inelastic neutrino-nucleon cross section above EeV [9][10][11][12][13] due to gluon saturation at small Bjorken-x < 10 −6 [14]. This will likely decrease the exponential attenuation of the Earth-crossing neutrino flux by at most a factor of 2-3 [15][16][17], whereas an order of magnitude or more suppression is needed to explain the two ANITA events.
Another explanation of the anomalous events within the SM framework was proposed in terms of the transition radiation from a particle shower crossing the Earthair interface and induced by an Earth-skimming neutrino [18]. In this model, the plane-of-polarization correlation to geomagnetic angles would be coincidental. Since both ANITA events are well-correlated to the local geomagnetic angle, and are consistent within 3 • -5 • of measurement error, coincidental alignment for both is possible only at the few % level [2]. Moreover, the diffuse neutrino flux necessary for this explanation to work is in tension with the current best limits from the Pierre Auger [19], IceCube [20] and ANITA [21] data.
Several beyond SM (BSM) interpretations of the ANITA anomalous events, such as sterile neutrino mixing [22,23], heavy dark matter [24,25] and stau [7,26] decays, have also been discussed. All of these explanations assume that the showers observed by ANITA were initiated by the hadronic decays of a τ -lepton. However, a major challenge for any BSM interpretation in which the ANITA events are initiated by a decaying τ lepton is arXiv:1810.08479v1 [hep-ph] 19 Oct 2018 to explain the apparent discrepancy with the null observation of any comparably energetic and steeply inclined throughgoing track events at IceCube [27] 1 , which has been operating at its design sensitivity for more than nine years, as compared to ANITA's total runtime of merely 67 days during the three flights. With an effective area of 1 km 2 (as compared to ANITA's 4 km 2 ) at EeV energies, the IceCube exposure is almost 12 times that of ANITA. Based on this argument, it was pointed out [23] that the sterile neutrino explanation [22] is in strong tension with IceCube. The same conclusion holds for the hypothesis of quasi-stable dark mater decay inside the Earth [24], and for the stau-based proposals [7,26]. Moreover, as pointed out in Ref. [25], given the local dark matter density of 0.3 GeV.cm −3 , the capture rate of an EeV-scale decaying dark matter is ridiculously low, corresponding to one dark matter decay every 137 years in the entire volume of the Earth.
The assumption that the upward going showers observed by ANITA were initiated by τ -leptons may be premature, since it is not clear how the ANITA experiment would distinguish between showers initiated by different kinds of particle decays on an event-by-event basis. The decays of a highly boosted BSM particle directly into hadrons, electrons, or photons would also result in the production of an impulsive radio cone, and this might give rise to miscalibrated energy measurement or effective area prediction when interpreted in terms of a τ -hypothesis. Moreover, all the BSM scenarios discussed above assume an isotropic flux of incident neutrinos, which has serious problems producing the observed arrival directions at ANITA without overproducing at shallower angles. As we discuss in this paper, it is difficult to account for the ANITA anomalous events using an isotropic Greisen-Zatsepin-Kuzmin (GZK) neutrino flux. We therefore consider an anisotropic flux to fit the ANITA events and show that this is so far consistent with the existing searches for potential candidate transient sources in the northern sky.
We then propose a new BSM solution to the ANITA puzzle in terms of a long-lived neutral particle. In particular, we advocate a GeV-scale bino in R-parity violating supersymmetry (RPV-SUSY) as a natural candidate for this purpose. Our solution has several advantages over the other BSM explanations entertained earlier: (i) The RPV couplings allow the on-shell, resonant production of a TeV-scale squark/slepton from neutrino-nucleon/electron scattering inside Earth, thereby naturally enhancing the signal cross section.
(ii) The squark/slepton decay inside the Earth can produce a pure bino which interacts with the SM 1 Though it is worth noting that there are three IceCube events which could be interpreted as throughgoing τ tracks with energy 10 to 100 PeV and inclined at angles 10 to 30 degrees below the horizontal [7,28].
fermions only via the U (1) Y gauge interactions and heavier supersymmetric particles, and therefore, can easily travel through thousands of km inside the Earth without significant energy loss.
(iii) For a suitable, yet realistic sparticle mass spectrum and couplings consistent with all existing low and high-energy constraints, we find some parameter space where the bino is sufficiently long-lived (with proper lifetime of order of ns) and decays to SM fermions at or near the exit-surface of Earth to induce the air shower observed by ANITA.
(iv) For LLE-type RPV couplings, the bino decays to a τ -lepton (if kinematically allowed) and electron, either of which could induce the air shower seen by ANITA, whereas for LQD-type couplings, the bino directly decays to quarks and neutrino, which mimic the SM τ -decay. In the latter case, there exist some parameter space for which no throughgoing track events are predicted at IceCube.

II. THE MODEL SETUP
We start with the general trilinear RPV superpotential in the Minimal Supersymmetric SM [29]: (2) L -singlet chiral superfields, and i, j, k = 1, 2, 3 are the generation indices.
Here we have suppressed all gauge indices for brevity. Note that SU (2) L and SU (3) c gauge invariance enforce antisymmetry of the λ ijk and λ ijk couplings with respect to their first and last indices, respectively. Since we are interested in the UHE neutrino interactions with matter, we will only consider the λ and λ -terms, one at a time. 2

A. LLE Contribution
Let us consider the λ-terms first. Expanding them in Eq. (1), we obtain the Lagrangian With these interactions, we can have new contributions to the (anti)neutrino-electron scattering inside Earth, as shown in Fig. 1 (top panel). In particular, given enough energy of the incoming (anti)neutrino, this will lead to the resonant production of a left-handed (LH) slepton through the second term in Eq. (2), and similarly a right-handed (RH) slepton through the third term in Eq. (2). For an incoming neutrino energy E ν , the slepton mass at which the resonance occurs is simply [33,34], where s is the centerof-mass energy. This is reminiscent of the Glashow resonance in the SM, where an on-shell W boson is produced from theν e − e scattering with an initial neutrino energy of E ν = m 2 W /2m e = 6.3 PeV [35]. Once produced, the slepton can decay back to an electron and neutrino through the same RPV interaction in Eq. (2) or to the corresponding lepton and neutralino through gauge interactions. Since we eventually want a τ -lepton in the final state to induce the upward air shower seen by ANITA, we will assume the slepton to be a stau ( τ ), and also assume the lightest neutralino (χ 0 1 ) to be much lighter than the stau, so that the decay τ → τ χ 0 1 is kinematically allowed. All other sparticles are assumed to be heavier than the stau and do not play any role in our analysis, except for the gravitino ( G), which could be the lightest supersymmetric particle (LSP) and plays the role of dark matter in this scenario. 3 The cross-section for νe → τ → τ χ 0 1 production is given by [33] 3 Gravitino LSP and bino next-to-LSP (NLSP) can be realized, e.g. in natural gauge mediation without gaugino unification [36]. which can be approximated by a Breit-Wigner formula close to the s-pole with s → m 2 τ : where g ≡ e/ cos θ w is the U (1) Y gauge coupling (e being the electromagnetic coupling and θ w the weak mixing angle), and j = 1, k = 3 or vice versa, depending on whether it is the RH or LH-slepton resonance, respectively. The index i for the incoming neutrino is free and we will assume a democratic flux ratio 1 : 1 : 1 for ν e : ν µ : ν τ and similarly for antineutrinos, as expected for a typical astrophysical neutrino flux with 1 : 2 : 0 flavor composition at the source [37]. We assume the bino is long-lived enough to survive its passage through Earth, before decaying close to or at the surface of exit. It can decay back to the τ -lepton and an off-shell stau, leading to a 3-body final state: In principle, the upgoing shower may be initiated either directly by the electron, or by the subsequent decay of the τ , as shown schematically in Fig. 2. In the limit m χ 0 1 m τ , the 3-body decay rate can be estimated as According to the geometry shown in Fig. 2, the incident neutrino travels a short distance l 1 inside Earth, and the remaining distance l 2 is traveled by the bino. In the limit l 1 l 2 , we can approximate l 2 as the chord length of ∼ 7000 km required for the two ANITA events, which translates into the mean lifetime of bino in its rest frame as τ χ 0 1 ∼ 0.22 ns. From the decay kinematics of the event, we estimate that the incoming neutrino energy should be roughly four times the detected energy at ANITA. Given that the two ANITA events had an average energy of 0.5 EeV, the initial neutrino energy should be E ν ∼ 2 EeV. Then the resonance condition fixes the stau mass: m τ = √ 2E ν m e 2 TeV. Substituting this in Eq. (5), we find that for a typical value of |λ| ∼ 0.1, as allowed by current experimental constraints [38], one needs a light bino of mass m χ 0 1 ∼ 2 GeV. A more accurate calculation of the allowed range in the (m χ 0 1 , |λ|) plane, taking into account all statistical interaction/decay probabilities for the neutrino, bino and tau, will be presented in a later section.
We should mention here that for gravitino LSP and bino NLSP, the bino can also have a 2-body decay into a photon and a gravitino via its photino component and the decay rate is given by [39] 2. A sketch of our model setup. The incoming UHE neutrino interacts with electron or quark inside the Earth within a distance l1 and produces a bino through the diagrams shown in Fig. 1. The bino travels a distance l2, before decaying to νe − τ + or νqq close to the surface, which induces the air shower seen by ANITA. Here, ltot is the total distance between the point where the UHE neutrino enters Earth to the ANITA detector, located at a height h above Earth's surface, and θ is the incoming angle of the neutrino with respect to the vertical direction.
where x 3/2 ≡ m G /m χ 0 1 . The photon can also initiate the air shower, but as mentioned above, we will only consider the τ final state. For the parameter space we work with, the 3-body decay rate given by Eq. (6) is larger than the 2-body decay rate given by Eq. (5) for very light gravitino with mass m G 0.1 eV, which is actually preferred by cosmology [40]. In particular, our scenario is naturally free from the cosmological gravitino problem [41] and consistent with cosmological constraints, such as from Lyman-α forest [42], cosmic microwave background (CMB) lensing and large-scale structure [43].

B. LQD Contribution
Now we consider the λ -terms in Eq. (1) which, when expanded, lead to the Lagrangian These interactions will contribute to the neutrinonucleon scattering mediated by either s or u-channel exchange of a down-type squark, as shown in Fig. 1 (bottom panel). After being resonantly produced, the bino will have a 3-body decay via off-shell down-type squark: χ 0 1 → ddν and χ 0 1 → ude. The final-state quarks readily hadronize, and for first-generation squark mediation Since the down-quark parton distribution function (PDF) inside the proton is larger than the strange and bottom PDFs, we will only consider the down-quark in the initial state. Similarly, we only consider the firstgeneration squark in the intermediate state, so that the final-state quarks from the 3-body bino decay hadronize to pions, mimicking the hadronic shower induced by the τ . Therefore, we will only consider the λ ijk couplings with j = 1 and k = 1. All other supersymmetric particles (except for the bino NLSP and gravitino LSP) are assumed to be heavier and not to play any role in our analysis.
The total differential cross section for the (anti)neutrino-nucleon interactions can be written in terms of the Bjorken scaling variables x = Q 2 /2m N E ν and y = E ν /E ν , where m N = (m p +m n )/2 is the average mass of the proton and neutron for an isoscalar nucleon, −Q 2 is the invariant momentum transfer between the incident neutrino and outgoing bino, and E ν = E ν − E χ 0 1 is the energy loss in the laboratory frame. Keeping only the dominant s-channel contributions, we obtain [33] dσ where s = 2m N E ν is the squared center-of-mass energy, and f d , fd are the PDFs for down and anti-down quark within the proton, respectively. The Breit-Wigner resonance is regulated by the squark widths 4 Note that the resonance condition is satisfied for the incoming neutrino energy E ν = m 2 d /2m N x, but due to the spread in the initial quark momentum fraction x ∈ [0, 1], the resonance peak is broadened and shifted above the threshold value E th ν = m 2 d /2m N , unlike the LLE case where the resonance was much narrower. This is one of the reasons why the LQD case allows for a larger parameter space than the LLE case in explaining the ANITA events, as we show in the next section.

III. EVENT RATE
We estimate the total number of expected events in the following way: where ∆E ≡ E f − E i is the incident neutrino energy range that gives rise to the resonance, Φ ν (E ν ) is the incoming neutrino flux, T = 67 days is the total runtime for the reported three flights of ANITA, and A eff ·∆Ω is the effective area integrated over the relevant solid angle, averaged over the probability for interaction and decay to happen over the specified geometry shown in Fig. 2.
For the LLE case, since the resonance is very narrow, we can evaluate the energy integral dE ν at the resonance energy E ν = m 2 τ /2m e with the energy spread ∆E = 2Γ τ m τ /m e ∼ 0.05 EeV. The integrated effective area is defined as where l 1 , l 2 are the distance traveled by ν and χ 0 1 respectively, and l tot is the total distance from neutrino entering Earth to the detector, θ is the angle between the travel path and the vertical direction as defined before in Fig. 2, D is the distance between the bino's decay point to the detector, θ c ∼ 0.015 is the Cherenkov cone angle, l SM and l BSM stand for the neutrino interaction lengths in the SM and BSM case, respectively. The interaction length can be generically written as l int ∼ 1/(σN A ρ), where N A is the Avogadro number, ρ is the density and σ is the interaction cross section. In our case, l BSM is a function of the new physics parameters λ ijk and m τ . Including the τ -lepton decay probability would modify Eq. (13) to the following: where l decay,τ is the decay length of τ . A quick test could be done to see if the GZK flux would be strong enough to provide the required number of events. Taking a benchmark point m χ 0 1 = 2 GeV and λ ijk = 0.2, the flux needed to generate 2 events at 3σ confidence (with Poisson distribution) is shown in Fig. 3. The blue shaded region corresponds to the isotropic flux needed to match ANITA events within 3σ. Thus, we need an isotropic flux as strong as ∼ 5 × 10 −23 (GeV · cm 2 · s · sr) −1 , at least two orders of magnitude larger than the GZK upper limit ∼ 2.2 × 10 −26 (GeV · cm 2 · s · sr) −1 at EeV level [20], shown as the grey curve in Fig. 3 . Therefore, under this RPV-SUSY bino scenario, GZK flux is disfavored as the source of UHE neutrinos at more than 3σ CL.
The challenge in explaining the ANITA event rate in terms of an isotropic GZK flux is quite general. Given the GZK upper limit of ∼ 50 EeV beyond which the UHECR flux is suppressed due to interactions of UHE-CRs with relic photons [44,45] and noting that the average neutrino energy is roughly 5-10% of the primary CR energy [46], we can integrate the projected differential GZK flux given in Ref. [47] from 0.1 EeV up to 5 EeV, and use Eq. (12), with ∆Ω = 2π and writing A eff = (4 km 2 ) × where 4 km 2 is the inferred area of the radio cone for the observed ANITA events and is interpreted as a conversion efficiency for incoming neutrinos into observed upward going events at ANITA. We find a predicted event rate of N ∼ 600 . Two events at ANITA therefore suggests ∼ 3 × 10 −3 . Under the hypothesis that the event is initiated by a long-lived particle with γcτ ∼ R Earth which decays below an altitude of 10 km after emerging from Antarctica (otherwise the air density drops rapidly and the shower does not fully develop before it reaches ANITA), there is already a suppression factor in which goes like (10 km)/R Earth 2 × 10 −3 . Similarly, if the events are due to a decaying τ which has been produced in the collision between a highly energetic weakly interacting particle with scattering length 1/(σN Av ρ) ∼ R Earth and a nucleus in the Antarctic crust or ice within 10 km of the surface, there is also the same suppression factor. This is before considering the production cross section for the long-lived particle in a ν-N collision in the northern hemisphere, and any additional branching ratio suppression.
We can still consider anisotropic sources for the incoming neutrinos. In this case, the solid angle considered now (∼ 0.0007π, defined using the uncertainty of angles in the ANITA events [2]) is much smaller than the solid angle in isotropic case (∼ 1.3π). Therefore, to get the same amount of events, the angular averaged anisotropic flux needed will become even larger than the isotropic flux, which is shown in Fig. 3 as the red area. Due to the angular average effect, we can see that the required anisotropic flux is at least two orders of magnitude larger than the corresponding isotropic flux, i.e. ∼ 10 −20 (GeV · cm 2 · s · sr) −1 . According to Refs. [48][49][50][51][52], such a large flux could in principle come from a transient point source, such as blazar, supernova burst, gammaray burst, or starburst galaxy, in the northern sky. In Refs. [50,51], upper bounds on the strength of neutrino flux from point sources in the north sky are given as ∼ 3.2 × 10 −20 (GeV · cm 2 · s · sr) −1 around 0.5 EeV. We also need to note that along the line of sight for the two events observed by ANITA, there might be more than one possible point source aligning in the same direction, making our current estimate a conservative one. The ANITA collaboration [2] also considered the transient source possibility for their anomalous events and, in fact, found a supernova candidate well within their expected angular uncertainty in the sky. The current and future ground arrays, such as Pierre Auger, Telescope Array (TA), Auger-Prime and TA×4 should be able to shed more light on these transient sources [53].
Assuming the average strength of the anisotropic sources being Φ ν ∼ 2×10 −20 (GeV·cm 2 ·s·sr) −1 with the mass of the RPV-SUSY mediator stau at m τ = 2 TeV, we use Eqs. (4) and (12) and perform a statistical analysis to find the 2σ and 3σ favored regions of the parameter space in the (λ ijk , m χ 0 1 ) plane, which is shown in Fig. 4 left panel as the green and yellow shaded regions, respectively. The dashed blue line shows the relation between the parameters once we set the bino decay length to be exactly the maximum distance traveled inside Earth, which corresponds to a mean lifetime of τ χ 0 1 ∼ 0.22 ns. The horizontal purple and red shaded regions in Fig. 4 are excluded from the R τ measurements [38]. The vertical shaded region is the kinematically forbidden region for the bino to decay into a τ -lepton. Combining all the constraints, we find a window with λ i31 ∼ 0.1 and m χ 0 1 ∼ 2 GeV for the new physics parameters to explain the events observed by ANITA. Changing the stau mass to lower values would shrink the allowed parameter space, as it decreases the incoming neutrino energy range allowed by the resonance condition. We cannot increase the stau mass to higher values either, as that would require a higher energy for the incoming neutrino, which is already at the edge of the GZK limit. Therefore, we predict that a ∼ 2 TeV stau should be observed at the LHC, if the LLE-type RPV-SUSY interpretation of the ANITA events is correct. The current LHC lower limits on the stau mass in the RPV-SUSY scenario is about 500 GeV, derived from multilepton searches [54].
As for the LQD case, we can do a similar calculation as for the LLE case described above, except that we can no longer replace the energy integral in Eq. (12) by a delta function, due to a much broader resonance [cf. Eqs. (8) and (9)]. our result is shown in the right panel of Fig. 4. Here we have chosen the down-squark mediator mass as m d = 1 TeV. This is consistent with the current LHC limit, which is around 900 GeV [55]. Increasing the squark mass moves the 2σ and 3σ contours to larger λ values, which are excluded by the V ud and Q w measurements [38], shown by the horizontal red and purple shaded regions, respectively. Thus, we predict that if our LQD-type RPV-SUSY interpretation of the ANITA events is correct, then a TeV-scale squark should be soon found at the LHC.
Another independent test of the allowed parameter space shown in Fig. 4 might come soon from the Belle II upgrade [56], which could significantly improve the R τ measurements.
Note that there is no LEP limit on our light bino scenario, because the Z decay to binos is forbidden at the tree level. In fact, there is no lower limit for the bino mass, as long as it is not the dark matter candidate [57]. Similarly, the stringent cosmological bounds on RPV couplings from the requirement that any baryon or lepton number violating interactions in equilibrium down to the electroweak scale could spoil the mechanism of baryogenesis due to fast electroweak sphaleron processes [58] can be avoided in our setup, because the mediator slepton/squark mass is at the TeV-scale and a low-scale baryogenesis could happen after they freeze out.

IV. ICECUBE SIGNATURES
Explanations for the ANITA anomalous events which proceed via the decay of a τ lepton predict the presence of throughgoing track events at IceCube. While a few events exist which may be interpreted as being of this kind [7,28], their energies are one to two orders of magnitude smaller than those inferred for the events at ANITA. It is therefore worth exploring models which predict fewer or no throughgoing track events at IceCube. For both of the models presented in Section II, only a fraction of events will proceed via a τ decay. A variation on the LQD model may produce no ice-penetrating charged lep- 1 decay length to be exactly the same as its maximum travel length. The horizontal shaded areas are the excluded regions for single RPV couplings from low-energy experiments [38]. The vertical shaded regions are the kinematically forbidden regions for the bino decay considered here.
ton signature at all. For example, a model making use of a L i Q 3 D k vertex would mean no χ 0 1 → td k i decay for m χ 0 1 < m t , while a L 1 Q j D k would lead to an electron which does not penetrate ice. In this case, the leading IceCube signature is χ 0 1 decay within the volume of the IceCube detector, which is suppressed in rate compared to ANITA by an additional factor of h IC /h ANITA ∼ 1/10, in comparison with the throughgoing track signature.

V. CONCLUSION
We have explored a RPV-SUSY interpretation of the two anomalous upgoing air showers seen by ANITA. In our framework, the UHE neutrino interacts with Earth matter to resonantly produce a squark/slepton, which then decays to a long-lived bino, whose decay products are responsible for the upgoing air shower. We considered both LLE and LQD-type interactions and our main results are given in Fig. 4. We find that a light bino of a few GeV mass and the RPV couplings of order 0.1 provide the best-fit solution to the ANITA events. In the LLE case, a stau of mass around 2 TeV, and in the LQD case, a down-type squark of mass around 1 TeV are predicted, which should be accessible by the next run of the LHC. The Belle II upgrade will provide a complementary low-energy probe of the allowed parameter space. Our hypothesis could be completely tested with more events at ANITA-IV (and beyond), as well as by IceCube in the future. It would be remarkable if weak-scale supersymmetry was discovered in such an unexpected way!