Signatures of Supersymmetry in Neutrino Telescopes

TL;DR

This paper explores how neutrino telescopes can detect signs of R-parity violating supersymmetry through resonant production of SUSY particles by ultra-high energy neutrinos, setting limits on RPV couplings and discussing recent anomalies.

Contribution

It provides the first robust limits on RPV couplings from IceCube data and discusses potential SUSY explanations for ANITA's anomalous events.

Findings

No significant deviation from Standard Model in IceCube data.

Established upper limits on RPV couplings based on neutrino observations.

Proposed SUSY interpretations for ANITA's upward-going air showers.

Abstract

We review the prospects of probing $R$-parity violating Supersymmetry (RPV SUSY) at neutrino telescopes using some of the highest energy particles given to us by Nature. The presence of RPV interactions involving ultra-high energy neutrinos with Earth-matter can lead to resonant production of TeV-scale SUSY partners of the SM quarks and leptons (squarks and sleptons), thereby giving rise to potentially anomalous behavior in the event spectrum observed by large-volume neutrino detectors, such as IceCube, as well as balloon-borne cosmic ray experiments, such as ANITA. Using the ultra-high energy neutrino events observed recently at IceCube, with the fact that for a given power-law flux of astrophysical neutrinos, there is no statistically significant deviation in the current data from the Standard Model expectations, we derive robust upper limits on the RPV couplings as a function of the resonantly-produced squark mass, independent of the other unknown model parameters, as long as the squarks decay dominantly to two-body final states involving leptons and quarks through the RPV couplings. Also, we discuss RPV SUSY interpretations of the recent anomalous, upward-going EeV air showers observed at ANITA, in terms of long-lived charged or neutral next-to-lightest SUSY particles.

Equations

$$(-1{)}^R=(-1{)}^{3(B-L)+2S}.$$

$$W_{\mathrm{MSSM}}=\mu H_u{H}_d+h_{u}Q{H}_u{U}^c+h_{d}Q{H}_d{D}^c+h_{e}L{H}_d{E}^c,$$

$$\displaystyle W_{\rm RPV}\ =\$$

$$\displaystyle{\cal L}_{LQD}\ =\ \lambda^{\prime}_{ijk}\bigg{[}\tilde{\nu}_{iL}\bar{d}_{kR}d_{jL}+\tilde{d}_{jL}\bar{d}_{kR}\nu_{iL}+\tilde{d}_{kR}^{*}\bar{\nu}_{iL}^{c}d_{jL}$$

$$\displaystyle\qquad\qquad-\tilde{e}_{iL}\bar{d}_{kR}u_{jL}-\tilde{u}_{jL}\bar{d}_{kR}e_{iL}-\tilde{d}_{kR}^{*}\bar{e}_{iL}^{c}u_{jL}\bigg{]}+{\rm H.c.}$$

$$\frac{d^2\sigma }{dxdy}=$$

$$\displaystyle a_{d_{j}}^{\rm CC}\$$

$$\displaystyle b_{\bar{u}_{j}}^{\rm CC}\$$

$$\displaystyle a_{d_{j}}^{\rm NC}\$$

$$\displaystyle b_{d_{j}}^{\rm NC}\$$

$${\Gamma }_{{\tilde{d}}_{kR}}\simeq \frac{m_{{\tilde{d}}_{kR}}}{8\pi }\sum _{ij}|{\lambda }_{ijk}^{\prime }{|}^2,$$

$${\mathcal{L}}_{LLE}=$$

$$N_{\mathrm{bin}}=T{N}_A{\int }_{E_{\min }^{\mathrm{bin}}}^{E_{\max }^{\mathrm{bin}}}d{E}_{\mathrm{dep}}{\int }_0^{1}dy{V}_{\mathrm{eff}}\Phi \Omega \frac{d\sigma }{dy},$$

$$L=\prod _{\mathrm{bins}\text{leavevmode}i}\frac{e^{-{\lambda }_i}{\lambda }_i^{n_i}}{n_i!},$$

$$-2\Delta \ln L=-2(\ln L-\ln L_{\max }),$$

$$\displaystyle\frac{d\sigma_{1}}{dt}\$$

$$\displaystyle\frac{d\sigma_{2}}{dt}\$$

$$c\tau \sim {\left(\frac{\sqrt{F}}{10\text{leavevmode}\mathrm{PeV}}\right)}^4{\left(\frac{100\text{leavevmode}\mathrm{GeV}}{m_{{\tilde{\tau }}_R}}\right)}^{5}10\text{leavevmode}\mathrm{km},$$

$${\sigma }_{LLE}=\frac{8\pi s}{m_{\tilde{\tau }}^2}\frac{\Gamma (\tilde{\tau }\to e\nu )\Gamma (\tilde{\tau }\to \tau {\tilde{\chi }}_1^0)}{(s-m_{\tilde{\tau }}^2{)}^2+m_{\tilde{\tau }}^2{\Gamma }_{\tilde{\tau }}^2}{\left(\frac{s-m_{{\tilde{\chi }}_1^0}^2}{m_{\tilde{\tau }}^2-m_{{\tilde{\chi }}_1^0}^2}\right)}^2,$$

$${\sigma }_{LLE}$$

$$\Gamma ({\tilde{\chi }}_1^0\to {\tau }^{-}{e}^{+}\bar{\nu })\simeq \frac{g^{\prime 2}|{\lambda }_{ijk}{|}^2}{512{\pi }^3}\frac{m_{{\tilde{\chi }}_1^0}^5}{m_{\tilde{\tau }}^4}.$$

$$\displaystyle\frac{d\sigma^{\nu}_{LQD}}{dxdy}\ =\$$

$$\displaystyle\frac{d\sigma^{\bar{\nu}}_{LQD}}{dxdy}\ =\$$

$${\Gamma }_{{\tilde{d}}_{kR}}=\frac{m_{{\tilde{d}}_{kR}}}{8\pi }\left[\sum _{ij}|{\lambda }_{ijk}^{\prime }{|}^2+\frac{2}{9}{g}^{\prime 2}\right],$$

$${\Gamma }_{{\tilde{d}}_{kL}}=\frac{m_{{\tilde{d}}_{kL}}}{16\pi }\left[\sum _{ij}|{\lambda }_{ijk}^{\prime }{|}^2+\frac{1}{9}{g}^{\prime 2}\right].$$

Full text

I Introduction

The Standard Model (SM) of particle physics has worked beautifully to predict what experiments have shown so far about the basic building blocks of matter Tanabashi:2018oca . However, physicists have long recognized that it cannot be the ultimate theory of Nature. Experimental evidence in support of this include the observation of nonzero neutrino mass, dark matter (DM), dark energy, and matter-antimatter asymmetry, none of which the SM can explain. Among the theoretical problems the SM fails to address are the quadratic divergence of the Higgs mass (the so-called hierarchy problem), the adhoc choice of the gauge groups and particle representations, unification of gauge couplings, inclusion of gravity, and the origin of electroweak symmetry breaking. While these arguments all point to some beyond the SM (BSM) physics, without further assumptions, they do not point decisively to the energy scale at which it should occur.

However, given the observed Higgs boson mass of 125 GeV, the quadratically divergent corrections in the SM imply that unless the new physics scale $\Lambda$ is close to the TeV scale, we need unnatural fine-tuning of parameters to keep the Higgs mass at its observed value. For instance, if $\Lambda\sim 10^{16}$ GeV, i.e. at the Grand Unified Theory (GUT) scale, then the mass-squared parameter in the SM Lagrangian will have to be fine-tuned to 1 part in $10^{26}$ to maintain a physical Higgs mass at 125 GeV. Although not logically impossible, such extreme fine-tuning is usually thought to be symptomatic of a deeper issue, often dubbed as the naturalness problem. If taken seriously, this leads to the inevitable conclusion that there must be new degrees of freedom around TeV scale.

Supersymmetry (SUSY), which predicts the existence of a partner differing in spin by $\frac{1}{2}$ to every known particle is widely considered as the most attractive solution to the naturalness and hierarchy problems. Specifically, due to the new symmetry between bosons and fermions, the quadratically divergent loop contributions of the SM particles to the Higgs mass squared are exactly canceled by the corresponding loop contributions of their SUSY partners. There are numerous other theoretical and phenomenological motivations for SUSY, such as (i) the most general symmetries of the $S$-matrix in quantum field theory are a direct product of the super-Poincaré group, which includes SUSY, with the internal symmetry group; (ii) gravity is naturally incorporated if SUSY is made local, as in supergravity models; (iii) SUSY is an essential ingredient of superstring theories, which are by far the best candidates for a consistent quantum theory of gravity; (iv) the SUSY particle spectrum leads to the remarkable unification of the strong and electroweak gauge couplings at an energy scale of $M_{\rm GUT}\sim 2\times 10^{16}$ GeV, which is consistent with proton decay constraints; (v) the lightest SUSY particle (LSP), if electrically and color neutral, makes an excellent candidate for DM; and (vi) the observed Higgs boson mass is consistent with the upper limit of $\lesssim 135$ GeV on the lightest Higgs boson in the minimal SUSY SM (MSSM). All these arguments in favor of SUSY make it arguably the strongest contender for BSM physics Haber:1984rc ; Baer:2006rs .

However, the lack of evidence for superpartners in the Large Hadron Collider (LHC) data so far atlas-susy ; cms-susy has forced the simplest SUSY scenarios toward regions of parameter space unnatural for the Higgs sector Papucci:2011wy ; Hall:2011aa ; Brust:2011tb ; Buckley:2016kvr ; Ahmed:2017jxl . A simple way to preserve the Higgs naturalness by evading the current experimental constraints is by allowing $R$-parity Violation (RPV) Barbier:2004ez in the production and decays of superpartners at the LHC. Apart from significantly lowering the collider bounds on the SUSY spectrum Carpenter:2006hs ; Dreiner:2012np ; Allanach:2012vj ; Dreiner:2012wm ; Asano:2012gj ; Graham:2014vya ; Monteux:2016gag ; Dercks:2017lfq ; Guo:2018hbv ; Bansal:2018dge , RPV SUSY implies the violation of baryon and/or lepton numbers, which has important phenomenological consequences Barbier:2004ez . For instance, one can automatically generate non-zero neutrino masses and mixing Hall:1983id ; Ellis:1984gi ; Dawson:1985vr ; Joshipura:1994ib ; Hempfling:1995wj ; Roy:1996bua ; Mukhopadhyaya:1998xj ; Bhattacharyya:1999tv ; Davidson:2000ne ; Diaz:2003as ; Dreiner:2011ft ; Peinado:2012tp at either tree or one-loop level without introducing any extra particles beyond the MSSM field content. Similarly, the presence of $\Delta L\neq 0$ RPV vertices can lead to observable neutrinoless double beta decay ($0\nu\beta\beta$) Mohapatra:1986su ; Vergados:1986td ; Hirsch:1995zi ; Babu:1995vh ; Pas:1998nn ; Faessler:2007nz ; Allanach:2009xx ; Rodejohann:2011mu , as well as successful baryogenesis Cline:1990bw ; Masiero:1992bv ; Sarkar:1996sn ; Dolgov:2006ay ; Kohri:2009ka ; Cui:2012jh ; Sorbello:2013xwa ; Arcadi:2015ffa ; Farina:2016ndq ; Pierce:2019ozl . Moreover, RPV scenarios also provide a compelling solution to some recent anomalies, such as the muon anomalous magnetic moment Kim:2001se ; Bhattacharyya:2009hb ; Hundi:2011si ; Chakraborty:2015bsk , lepton flavor universality violation in semileptonic $B$-meson decays Deshpande:2012rr ; Biswas:2014gga ; Huang:2015vpt ; Zhu:2016xdg ; Deshpand:2016cpw ; Altmannshofer:2017poe ; Das:2017kfo ; Earl:2018snx ; Sheng:2018ylg ; Trifinopoulos:2018rna ; Hu:2018lmk and anomalous upgoing events at ANITA Collins:2018jpg . If any of these anomalies becomes statistically significant, one should consider RPV SUSY as a strong candidate for the underlying new physics. Hence, it is of paramount importance to find complementary ways at as many different energy scales as possible to test this scenario.

In this review, we discuss one such possibility, namely, testing the RPV SUSY scenario using ultra-high energy (UHE) neutrinos at neutrino telescopes, such as IceCube and KM3NeT. The basic idea is that the presence of RPV SUSY interactions involving UHE neutrinos can lead to resonant production of TeV-scale sfermions inside the Earth Carena:1998gd ; Dev:2016uxj ; Collins:2018jpg , thereby giving rise to distinct features in the detected neutrino spectrum at neutrino telescopes.111Even in the absence of RPV, one could consider pair production of sparticles in neutrino-nucleon interactions Albuquerque:2003mi . In this case, the relatively small cross-section for the production of sparticles, as compared to the SM neutrino-nucleon cross section, can be partially compensated for, provided one of the sparticles is very long lived Albuquerque:2003mi ; Ahlers:2006pf ; Albuquerque:2006am ; Ando:2007ds ; Canadas:2008ey ; Meade:2009mu ; Connolly:2018ewv . We will briefly discuss this possibility in Section V. Note that in the SM, the only observable resonance for UHE neutrinos is the Glashow resonance Glashow:1960zz : $\bar{\nu}_{e}e^{-}\to W^{-}$, which occurs at incoming electron antineutrino energy $E_{\nu}=\frac{m_{W}^{2}}{2m_{e}}=6.3$ PeV. In contrast, the squark/slepton resonance can occur at a different incoming neutrino energy, depending on the mass of the SUSY particles. Since no such resonance feature has been observed in the IceCube data so far Aartsen:2017mau (except for a mild $2\sigma$-level discrepancy around 100 TeV which could very well be an artifact of the fitting procedure or astrophysical flux uncertainty), one can derive robust upper limits on the RPV couplings as a function of the resonantly-produced sfermion mass, largely independent of the other unknown model parameters Dev:2016uxj . With more statistics, we expect these limits to be comparable/complementary to the existing limits from direct collider searches and other low-energy processes. More importantly, if IceCube does observe a statistically significant resonance-like excess feature in the future, it might point to a preferred range of the relevant RPV coupling and sfermion mass that could be looked for more vigorously in the LHC data.

The rest of this article is organized as follows: in Section II, we briefly discuss the RPV interactions in the MSSM. In Section III, we focus on the UHE neutrino interactions involving trilinear RPV. In Section IV, we calculate the event rate at IceCube and derive constraints on the $\lambda^{\prime}$-type RPV couplings as a function of the squark mass. In Sections V and VI, we consider the possibility of long-lived charged and neutral SUSY particles respectively, and implications for neutrino telescopes, as well as balloon experiments. We conclude in Section VII.

II RPV Interactions

In the SUSY extension of the SM, all the fields, both fermionic and bosonic, get promoted to become chiral superfields which have both a bosonic and a fermionic part. The bosonic partners of known SM fermions are called sfermions (i.e. squarks and sleptons), whereas the fermionic partners of the Higgs fields are called Higgsinos, and those of the SM gauge fields are called gauginos (i.e., gluino, bino, wino, photino).

$R$-parity is associated with a $Z_{2}$ subgroup of the group of continuous $U(1)$ $R$-symmetry transformations acting on the gauge superfields and the two chiral doublet Higgs superfields $H_{u}$ and $H_{d}$ responsible for electroweak symmetry breaking Fayet:1974pd . $R$-parity appears as the discrete remnant of this continuous $U(1)_{R}$ symmetry that must be broken to allow for the gravitino (the superpartner of graviton in supergravity models) and gluinos to acquire masses. There is an intimate connection between $R$-parity and baryon ($B$) and lepton ($L$) number conservation laws Farrar:1978xj , which is made explicit by the following expression for the $R$-parity of a particle with spin $S$:

[TABLE]

Under this symmetry, all the SM fields are even, while all the superpartner fields are odd. Thus, the conservation of $R$-parity implies that all superpartner fields must be produced in pairs and must decay to states containing at least one superpartner field. This also implies that the LSP is absolutely stable, which has two cosmological implications: (i) it cannot be electrically charged, as it is cosmologically disfavored DeRujula:1989fe ; Dimopoulos:1989hk ; Pospelov:2006sc , and (ii) if electrically and color neutral, it is a good candidate for DM Jungman:1995df .

The MSSM superpotential invariant under the SM gauge group $SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$ that preserves $R$-parity is written as

[TABLE]

where $L\ni(\nu,e)_{L}$ and $Q\ni(u,d)_{L}$ are the $SU(2)_{L}$-doublet and $U^{c},D^{c},E^{c}$ are the $SU(2)_{L}$-singlet chiral superfields, respectively. Since SUSY is not an exact symmetry of nature, one must also add the SUSY breaking terms, which consist of soft mass terms for superpartners and trilinear $A$-terms of the type in Eq. (2) but with the scalar components of the superfields and mass terms for gauginos, all of which still respect the gauge symmetry.

As mentioned in Section I, in light of the current experimental results, the case for $R$-parity conservation (RPC) in effective low-energy SUSY is less compelling than RPV. We have a few additional points in favor of RPV:

• •

$R$-parity by itself is not enough to protect against dangerous proton decay operators Nath:2006ut , as RPC but non-renormalizable superpotential interactions of the general form $U^{c}U^{c}D^{c}E^{c}$ can give rise to proton decay. So the standard proton-stability motivation for $R$-parity is actually gone.

• •

Although RPV interactions spoil the stability of the DM, it is still possible to have the LSP long-lived enough to play the role of DM, e.g. gravitino because of its Planck-suppressed gravitational interactions Takayama:2000uz ; LopezFogliani:2005yw ; Buchmuller:2007ui ; Monteux:2014tia . Another possibility is to have a different symmetry than the $R$-parity, under which the SM is inert.

• •

Introducing RPV at low scale still preserves the gauge coupling unification Altmannshofer:2017poe – a hallmark prediction of SUSY. In addition, RPV can be made compatible with the unification of quarks and leptons in GUT models DiLuzio:2013ysa ; Bajc:2015zja . In fact, large RPV couplings are favorable to realize the $b-\tau$ Yukawa unification Hu:2018lmk .

• •

RPV is mandatory in minimal theories that extend MSSM to include local $B-L$ symmetry Mohapatra:1986su ; Martin:1992mq ; Mohapatra:2015fua . Since the nonzero neutrino masses seem to suggest that the underlying BSM physics must include a broken local $B-L$ symmetry, the fact that its supersymmetrization leads to RPV makes it even more compelling.

There are two classes of RPV models: (i) spontaneous breaking by giving vacuum expectation values to neutral superpartners, such as sneutrino Aulakh:1982yn , or (ii) explicit breaking by adding RPV terms in the superpotential (2) Hall:1983id . In this review, we will focus on the latter case, where the explicit RPV terms in MSSM can be written as

[TABLE]

where $i,j,k=1,2,3$ are the family indices and we have suppressed all gauge indices for brevity. Note that $SU(2)_{L}$ and $SU(3)_{c}$ gauge invariance enforce antisymmetry of the $\lambda_{ijk}$- and $\lambda^{\prime\prime}_{ijk}$-couplings with respect to their first and last two indices, respectively.

Since we are mostly interested in the UHE neutrino interactions with nucleons, we will only focus on the $\lambda^{\prime}_{ijk}$-couplings. Any of these 27 new dimensionless complex parameters can lead to resonant production of TeV-scale squarks at IceCube energies Carena:1998gd ; Dev:2016uxj , thereby making a potentially significant contribution to the UHE neutrino events.222Even without SUSY, similar resonance features in neutrino-nucleon interactions can also occur in models with TeV-scale leptoquarks Robinett:1987ym ; Doncheski:1997it ; Anchordoqui:2006wc ; Alikhanov:2013fda ; Barger:2013pla ; Dutta:2015dka ; Dey:2015eaa ; Dey:2017ede ; Becirevic:2018uab . However, there are subtle differences between scalar leptoquark and RPV SUSY models, e.g. due to the presence of additional decay channels and chiral mixing between squarks in the RPV case. The $\lambda$-couplings in Eq. (3) can give rise to resonant production of selectrons from neutrino interactions with electrons, reminiscent of the Glashow resonance Glashow:1960zz in the SM; however, for TeV-scale selectrons, their contribution to the total number of IceCube events is negligible in the energy range considered here and we will comment on this possibility later in the context of ANITA events. Note that with non-zero $\lambda^{\prime}$-couplings, we need to explicitly forbid the $\lambda^{\prime\prime}$-terms, e.g. by imposing baryon triality Ibanez:1991hv , to avoid rapid proton decay, as the product $\lambda^{\prime}\lambda^{\prime\prime}$ is constrained to be $\lesssim 10^{-24}$ from proton lifetime bounds Smirnov:1996bg . We also ignore the bilinear terms in Eq. (3), since they do not give rise to the resonance feature exploited here, although they could lead to other distinct signatures (e.g. triple bang) relevant for future multi-km3 neutrino telescopes Hirsch:2007kx .

III Neutrino Interactions involving Trilinear RPV

We start with the $\lambda^{\prime}$-part of the RPV Lagrangian, after expanding the superpotential (3) in terms of the superfield components:

[TABLE]

At the IceCube, these interactions will contribute to both charged-current (CC) and neutral-current (NC) processes mediated by either $s$-channel or $u$-channel exchange of a down-type squark, as shown in Fig. 1.

The $s$-channel processes in Figs. 1(a) and (c) involve valence quarks, thus giving the dominant contributions to the (anti)neutrino–nucleon cross sections, provided the right-handed (RH) down-type squarks are produced resonantly. Similarly, the $s$-channel process in Fig. 1(e) mediated by a left-handed (LH) down-type squark can also give a resonant enhancement to the (anti)neutrino–nucleon cross section. Here we have implicitly assumed that the RPC squark decays to a quark and a gluino, neutralino or chargino are suppressed Butterworth:1992tc , compared to the RPV decays induced by Eq. (4). This is the case, for instance, in the region of RPV MSSM parameter space, where the gaugino masses $M_{1}$, $M_{2}$, as well as the $\mu$-term, are larger than the squark masses, thus kinematically forbidding the two body RPC decays of squarks. The 3-body RPC decays via virtual gauginos will in general be smaller compared to the 2-body decays through RPV couplings, as considered here.

On the other hand, the contributions from the $u$-channel processes in Fig. 1(b), (d) and (f) are much smaller, since they do not have a resonant enhancement, and in addition, for (b) and (d), due to the sea quark involvement. Moreover, the RPV contributions are the largest for the first generation quarks, which are the predominant constituents of the nucleon, and to some extent, for the second and third generation quarks, which become increasingly important at higher energies. At the same time, the existing indirect constraints from precision measurements in various low-energy processes are the strongest for the couplings involving the first generation and weaker for the second and third generations Barbier:2004ez ; Allanach:1999ic ; Kao:2009fg ; Dreiner:2012mx ; Domingo:2018qfg . For the SM CC and NC interactions Gandhi:1995tf ; CooperSarkar:2011pa , which must be included in the total neutrino-nucleon cross section giving rise to the IceCube events, we take into account all valence and sea quark contributions.

The total differential cross section for the neutrino-nucleon interactions, written in terms of the Bjorken scaling variables $x=Q^{2}/2m_{N}E^{\prime}_{\nu}$ and $y=E^{\prime}_{\nu}/E_{\nu}$, is given by

[TABLE]

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 lepton, $E_{\nu}$ is the incoming neutrino energy, $E^{\prime}_{\nu}=E_{\nu}-E_{\ell}$ is the energy loss in the laboratory frame, $E_{\ell}$ is the energy of the outgoing lepton, and $f(x,Q^{2}),\bar{f}(x,Q^{2})$ are the parton distribution functions (PDFs) within the proton for $f$-quark and anti $f$-quark, respectively. For the CC processes shown in Figs. 1(a) and (b), induced by an incoming neutrino of flavor $i$, the only non-trivial coefficients in Eq. (5) are respectively Carena:1998gd

[TABLE]

where $s=2m_{N}E_{\nu}$ is the square of the center-of-mass energy and $g$ is the $SU(2)_{L}$ gauge coupling, and $m_{W}$ is the $W$-boson mass. It is obvious to see that in Eqs. (6) and (7), the first term on the right-hand side is the SM CC contribution, whereas the second term is the RPV contribution. So for all SM CC processes involving $\bar{d}$ and $u$-type quarks, which do not have interference with the RPV processes, the coefficients in Eq. (5) are simply obtained by putting $\lambda^{\prime}_{ijk}=0$ in Eqs. (6) and (7).

Similarly, for the NC processes shown in Figs. 1(c)–(f), the only non-trivial coefficients are Carena:1998gd

[TABLE]

where $L_{d}=-(1/2)+(1/3)x_{w}$ and $R_{d}=(1/3)x_{w}$ are the chiral couplings, $x_{w}\equiv\sin^{2}\theta_{w}$ is the weak mixing angle parameter and $m_{Z}$ is the $Z$-boson mass. For all SM NC processes involving $u$-type quarks, which do not have interference with the RPV processes, the coefficients in Eq. (5) are simply obtained by putting $\lambda^{\prime}_{ijk}=0$ in Eqs. (8) and (9) and replacing $L_{d}\to L_{u}=(1/2)-(2/3)x_{w}$, $R_{d}\to R_{u}=-(2/3)x_{w}$. For neutrino-antiquark interactions, the coefficients for the NC processes can be obtained simply by crossing symmetry, i.e. $a_{f}\leftrightarrow b_{f},\ xs\leftrightarrow Q^{2}-xs$. Similarly, for antineutrino-nucleon interactions, we can just replace the PDFs $f\leftrightarrow\bar{f}$ in Eq. (5).

Note the Breit-Wigner resonance form of Eqs. (6) and (8), which is regulated by the RH down-type squark width

[TABLE]

assuming that the only dominant decay modes are $\tilde{d}_{kR}\to\nu_{iL}d_{jL}$ (NC) and $\tilde{d}_{kR}\to e_{iL}u_{jL}$ (CC), and the masses of the final state fermions in these 2-body decays are negligible compared to the parent squark mass. For the LH down-type squark, $\Gamma_{\tilde{d}_{kL}}=\Gamma_{\tilde{d}_{kR}}/2$, since $\tilde{d}_{kL}\to\nu_{iL}d_{jR}$ is the only available decay mode. The resonance condition is satisfied for the incoming energy $E_{\nu}=m^{2}_{\tilde{d}_{kR}}/2m_{N}x$, but due to the spread in the initial quark momentum fraction $x\in[0,1]$, the resonance peak will be broadened and shifted above the threshold energy $E_{\nu}^{\rm th}=m^{2}_{\tilde{d}_{kR}}/2m_{N}$. Nevertheless, given that the maximum contribution to the cross section happens around $x\sim 10^{-3}-10^{-4}$, one can immediately infer that for $m_{\tilde{d}_{kR}}\in$ [100 GeV, 2 TeV], $E_{\nu}^{\rm th}$ is in the multi TeV–PeV range, and hence, can be probed by the available IceCube data. Interestingly, this is exactly the mass range currently being probed at the LHC, and therefore, the neutrino telescope probes allow for an independent, complementary test of the SUSY parameter space.

Now let us consider the $\lambda$-couplings in Eq. (3), which lead to the $LLE$-type RPV Lagrangian

[TABLE]

With these interactions, we can have new contributions to the (anti)neutrino-electron scattering inside Earth. In particular, given enough energy of the incoming (anti)neutrino, this will lead to the resonant production of an LH slepton through the second term in Eq. (11), and similarly an RH slepton through the third term in Eq. (11). For an incoming neutrino energy $E_{\nu}$, the slepton mass at which the resonance occurs is simply $m_{\widetilde{e}_{i}}=\sqrt{s}=\sqrt{2E_{\nu}m_{e}}$, where $s$ is the center-of-mass energy Carena:1998gd ; Dev:2016uxj . This is reminiscent of the Glashow resonance in the SM, where an on-shell $W$ boson is produced from the $\bar{\nu}_{e}-e$ scattering with an initial neutrino energy of $E_{\nu}=m_{W}^{2}/2m_{e}=6.3$ PeV Glashow:1960zz . However, the corresponding threshold energy $E_{\nu}^{\rm th}=m^{2}_{\tilde{e_{k}}}/2m_{e}$ is beyond 10 PeV for selectron masses above 100 GeV or so. Since smaller selectron masses are excluded from the LEP data Abbiendi:2003rn ; Abdallah:2003xc , we cannot probe the $\lambda_{ijk}$-couplings with the current IceCube data. Nevertheless, if future data reports any events beyond 10 PeV, the $LLE$-type RPV scenario could in principle provide a viable explanation, given the fact that it would be difficult to explain those events within the SM and with an unbroken power-law flux, without having a significantly larger number of events in all the preceding lower-energy bins. We will come back to Eq. (11) in Section VI in a different context.

IV Event Rate

To give an illustration of the power of neutrino telescope in testing the RPV SUSY parameter space, let us calculate the event rate for IceCube including the RPV effect. In general, the expected number of high-energy starting events (HESE) in a given deposited energy bin at IceCube due to the modified cross section (5) can be estimated as

[TABLE]

where $T$ is the exposure time, $N_{A}$ is the Avogadro number, $E_{\rm dep}(E_{\nu})$ is the electromagnetic-equivalent deposited energy which is always smaller than the incoming neutrino energy $E_{\nu}$ in the laboratory frame by a factor depending on $E_{\nu}$ and the type of interaction (CC or NC) and lepton flavor, $V_{\rm eff}(E_{\nu})$ is the effective target volume of the detector, $\Phi(E_{\nu})$ is the incident neutrino flux, $\Omega(E_{\nu})$ is the effective solid angle of coverage, and we have integrated the differential cross section in Eq. (5) over $x\in[0,1]$, including both neutrino and antineutrino initial states with all flavors. For more details of the event rate calculation, see e.g. Refs. Laha:2013eev ; Chen:2013dza ; Aartsen:2013vja ; Aartsen:2013jdh ; Vincent:2016nut ; Palladino:2018evm ; Sui:2018bbh .

We calculate the predicted number of events with and without RPV interactions for the publicly available 6-year IceCube HESE data corresponding to $T=2078$ days. Here we have assumed the IceCube best-fit Aartsen:2017mau for a single-component unbroken astrophysical power-law flux $\Phi(E_{\nu})=2.5\times 10^{-18}(E_{\nu}/100\leavevmode\nobreak\ {\rm TeV})^{-2.9}\leavevmode\nobreak\ {\rm GeV}^{-1}\>{\rm cm}^{-2}{\rm s}^{-1}{\rm sr}^{-1}$ with a standard (1:1:1) flavor composition ratio on Earth, and have used the NNPDF2.3 leading order PDF sets Ball:2012cx for the cross section calculations.333The PDF uncertainties on the total cross-section are at most at 5% level Chen:2013dza for the energy range of interest. The flux uncertainties, on the other hand, are currently at the level of 15% Aartsen:2017mau . In addition, there is a tension between the HESE and throughgoing muon samples, which cannot be resolved by the given single-component unbroken power-law HESE fit and might require some BSM explanation by itself Sui:2018bbh . The expected background due to atmospheric neutrinos and muons are also taken from the IceCube analysis Aartsen:2017mau . Note that atmospheric $\nu_{e}$ events will also get modified due to nonzero RPV couplings $|\lambda^{\prime}_{11k}|$. However, since the atmospheric background is dominated by the $\nu_{\mu}$-induced events and the event rate for atmospheric $\nu_{e}$ is much smaller, we can safely ignore the $|\lambda^{\prime}_{11k}|$ effects on the background and just assume it to be basically the same as in the SM case. Also one might wonder whether the source flavor composition and flux of the neutrinos could be modified due to the new RPV interactions. However, this effect is expected to be small for the values of squark masses and RPV couplings considered here, since the SM weak interactions with strength $G_{F}$ (the Fermi coupling constant) will be dominant over the RPV interactions with relative strength of $|\lambda^{\prime}_{ijk}|^{2}/m_{\tilde{d}_{k}}^{2}$. The RPV effect at the IceCube detector could get enhanced only due to the resonant production of squarks in a conducive range of the incoming neutrino energy. Thus, adding the flux uncertainty in our analysis will equally affect the events due to both SM and RPV interactions, without changing the relative enhancement of the events in presence of the RPV interactions with respect to the SM prediction. This justifies our use of the IceCube best-fit value for the flux.

Since no statistically significant excess over the SM prediction is seen in the current IceCube data, we use this information to put an upper bound on the $|\lambda^{\prime}_{ijk}|$ couplings. To this effect, we perform a binned likelihood analysis with the Poisson likelihood function

[TABLE]

where the observed count $n_{i}$ in each bin $i$ is compared to the theory prediction $\lambda_{i}$, including the RPV contribution induced by $\lambda^{\prime}_{ijk}$. We then construct a test statistic

[TABLE]

from which a $1\sigma$ upper limit on $|\lambda^{\prime}_{ijk}|$ corresponding to the value of $-2\Delta\ln L=1$ can be derived. Here $L_{\rm max}$ represents the likelihood value for $\lambda^{\prime}_{ijk}=0$, i.e. with only the SM contribution to the interaction. Our results are shown in Fig. 2 (blue solid curves) for $|\lambda^{\prime}_{11k}|$ and $|\lambda^{\prime}_{12k}|$, i.e. for electron-type neutrinos interacting with the 1st and 2nd generation quarks, respectively.

For comparison, we also show the direct limits on $|\lambda^{\prime}_{11k}|$ from direct searches in $e^{\pm}p$ collisions at HERA with $\sqrt{s}=319$ GeV Collaboration:2010ez ; South:2016cmx , as shown by the magenta-shaded region in Fig. 2 (left). Squark masses below 100 GeV or so are disfavored from direct searches for RPV SUSY at LEP Abbiendi:2001aj ; Heister:2002jc ; Abbiendi:2003rn , Tevatron Abbott:1999nh ; Abe:1998gu and LHC atlas-susy ; cms-susy , and therefore, are not considered here.444In the absence of the possibility of a resonant production (as e.g. in the sneutrino case) in $e^{+}e^{-}$ and hadron-hadron collisions, it is difficult to cast most of the collider limits onto the $m_{\tilde{d}}-|\lambda^{\prime}_{ijk}|$ plane in a model-independent manner, and therefore, we do not attempt to show them in Fig. 2. In addition, the recent search for scalar leptoquarks at the 13 TeV LHC Aaboud:2019jcc is relevant for our RPV scenario, since $\lambda^{\prime}_{ijk}$-couplings also give rise to the same $e_{i}e_{i}jj$ final states via pair-production of down-type squarks (from gluon fusion), followed by $\tilde{d}_{kR}\to e_{iL}u_{jL}$ which has a branching ratio of 0.5. The corresponding 95% confidence level (CL) ATLAS limit of $\sim 1.2$ TeV on the first-generation scalar leptoquark mass can be directly translated into a lower bound on the down-type squark mass, as shown by the vertical dot-dashed line in Fig. 2. There also exist indirect constraints on $|\lambda^{\prime}_{11k}|$ from lepton universality in pion decay, measured by the ratio $R_{\pi}=\frac{{\rm BR}(\pi^{-}\to e^{-}\bar{\nu}_{e})}{{\rm BR}(\pi^{-}\to\mu^{-}\bar{\nu}_{\mu})}$, unitarity of the CKM element $V_{ud}$ and atomic parity violation Kao:2009fg , the most stringent of which is shown in Fig. 2 (left) by the red dotted curve. Other low-energy constraints, such as neutrino mass Bhattacharyya:1999tv , electric dipole moment Yamanaka:2014nba and flavor-changing $B$-decays Ghosh:2001mr ; Dreiner:2013jta , always involve the product of two independent RPV couplings, and hence, are not applicable in our case. For $k=1$, we have an additional constraint from $0\nu\beta\beta$, which however depends on the masses of other SUSY particles Rodejohann:2011mu , unlike the other limits discussed above, which are independent of the rest of the SUSY spectrum, as long as the 2-body RPV decay modes of the squark are dominant. Just for the sake of comparison, we translate the $0\nu\beta\beta$ half-life limit on 76Ge from GERDA phase-II Agostini:2018tnm onto our RPV parameter space in Fig. 2 (left) for two benchmark points with gluino and neutralino masses of 10 TeV (gray, dashed) and 10 PeV (gray, dotted), while keeping all sfermion masses equal to $m_{\tilde{d}}$. In the former case, the $0\nu\beta\beta$ limit is the most stringent one, whereas for either heavier gaugino masses or $k\neq 1$, the pion decay constraint is stronger than the limit obtained from IceCube in the entire mass range considered. Similarly, as shown in Fig. 2 (right), for $|\lambda^{\prime}_{12k}|$, the indirect constraints from lepton universality in neutral and charged $D$-meson decays, measured by the ratio $R_{D}=\frac{{\rm BR}(D\to Ke\nu_{e})}{{\rm BR}(D\to K\mu\nu_{\mu})}$, are stronger than the IceCube limit derived here. This rules out the possibility of any $|\lambda^{\prime}_{1jk}|$-induced observable excess in the 8-year IceCube data.

We note that the IceCube constraints are mostly limited by statistics, which is expected to improve significantly with more exposure time. To illustrate this point, we just scale the current 8-year dataset by a factor of 4 (roughly corresponding to 30 years of actual data taking) in all the bins analyzed here and derive projected limits on $|\lambda^{\prime}_{1jk}|$ (with $j=1,2$), following the same likelihood procedure described above. This conservative estimate of the future limits is shown in Fig. 2 by the blue dashed curves. We find that the limit on $|\lambda^{\prime}_{1jk}|$ can be improved roughly by up to a factor of 3, and it might even surpass the current best limit in the sub-TeV squark mass range for $j=2$, although the indirect limit from lepton universality could improve by an order of magnitude at Belle II Kou:2018nap . In practice, however, we may not have to wait for 30 years, since a number of unforeseen factors could improve the conservative projected IceCube limits shown here, e.g. the future data in all the bins may not scale proportionately to the current data and may turn out to be in better agreement with the SM prediction. Similarly, other large-volume detectors like KM3NeT Adrian-Martinez:2016fdl and IceCube-Gen2 Aartsen:2014njl might go online soon, thus significantly increasing the total statistics.

We also note that a similar analysis could be performed for incident neutrinos of muon and tau flavors at IceCube, though for muon neutrinos, one has to carefully reassess the atmospheric background including the RPV effects and also take into account the large uncertainty in converting the deposited muon energy to actual incoming $\nu_{\mu}$ energy. In any case, we expect the corresponding limits on $|\lambda^{\prime}_{ijk}|$ (with $i=2,3$ and $j=1,2$) to be weaker than the limits on $|\lambda^{\prime}_{1jk}|$ shown in Fig. 2 simply due to the fact that the effective fiducial volume at the IceCube is the largest for $\nu_{e}$ Aartsen:2017mau . Nevertheless, with more statistics, one could in principle consider the $|\lambda^{\prime}_{ijk}|$ couplings for all neutrino flavors. Also, one could improve the analysis presented here by taking into account the showers and tracks individually. Since the $|\lambda^{\prime}_{1jk}|$ couplings preferentially enhance only one type of events, viz. showers for $i=1,3$ and tracks for $i=2$, a binned track-to-shower ratio analysis is expected to improve the limits on the corresponding $|\lambda^{\prime}_{ijk}|$. In fact, by examining the track-to-shower ratio in future data, one might be able to distinguish between different new physics contributions to the IceCube events, provided the source flavor composition of the neutrinos is known more accurately.

V Long-lived Charged Sparticles

In SUSY theories where the LSP is the gravitino and the next-to-LSP (NLSP) is a long-lived charged particle (typically a slepton), e.g. in gauge mediation models Giudice:1998bp , the diffuse flux of UHE neutrinos interacting with the Earth could produce pairs of slepton NLSPs which could be detected in neutrino telescopes Albuquerque:2003mi ; Ahlers:2006pf ; Albuquerque:2006am ; Ando:2007ds ; Canadas:2008ey ; Meade:2009mu . The basic process involves the $t$-channel production of a slepton and a squark via gaugino exchange, see Fig. 3. Note that this does not require any RPV interactions. The neutrino can interact either with an LH down-type quark, or with an RH up-type anti-quark, with the resulting partonic cross sections respectively given by

[TABLE]

where $s$, $t$ and $u$ are the usual Mandelstam variables, $\alpha\equiv e^{2}/4\pi$ is the fine structure constant, and $m_{\widetilde{\chi}}$ is the relevant gaugino (charged/neutral) mass. The heavier LH slepton and the squark promptly decay to the lighter RH slepton plus SM particles, thus giving rise to a pair of slepton NLSPs. The SUSY cross sections are considerably suppressed with respect to the SM neutrino-nucleon interactions, even when well above the threshold, unlike in the resonant case. This is because the SM cross section is dominated by very small values of $x$, whereas for the SUSY processes, we need $x>\frac{m^{2}_{\rm SUSY}}{2m_{N}E_{\nu}}$, where $m_{\rm SUSY}$ is the typical mass of the SUSY particle being produced. However, if the NLSP slepton happens to be long-lived, as is the case in gauge-mediated SUSY breaking with stau NLSP, for instance, where the decay length in the lab frame is typically

[TABLE]

with $\sqrt{F}$ being the SUSY breaking scale, the long range of slepton compared to the background muon range of a few km compensates for the small SUSY production cross section. The high boost and large range of the NLSPs implies that its upward going tracks could in principle be detected in large ice or water Cherenkov detectors like IceCube or KM3NeT, provided one applies specialized cuts in the track separation to distinguish the signal events from the di-muon background.

In analyzing this type of signal, it is very important to include the energy loss during the propagation of sleptons through the Earth. It turns out that for heavy charged particles, the bremsstrahlung and pair-production energy losses become less important, and the photo-nuclear energy loss is mass-suppressed for slepton masses much larger than the tau mass Reno:2005si . Thus, neutrino telescopes provide a complementary probe of long-lived charged slepton NLSP and intermediate-scale SUSY breaking, beyond the current collider limits of about 450 GeV Khachatryan:2016sfv .

The long-lived stau hypothesis has gained some recent interest in view of two, anomalous upward-going ultra-high energy cosmic ray (UHECR) air shower events, with deposited shower energies of $0.6\pm 0.4$ EeV and $0.56^{+0.3}_{-0.2}$ EeV, observed by the Antarctic Impulsive Transient Antenna (ANITA) balloon experiment Gorham:2016zah ; Gorham:2018ydl . This is puzzling, because no SM particle is expected to survive passage through Earth a chord distance of $\sim 7000$ km (corresponding to the observed zenith angles of the two events) at EeV energies. In particular, the interpretation of these events as $\tau$-lepton decay-induced air showers at or near the ice surface arising from a diffuse UHE flux of cosmic $\nu_{\tau}$ is strongly disfavored due to their mean interaction length of only $\sim 300$ km. Including the effect of $\nu_{\tau}$ regeneration, the resulting survival probability over the chord length of the ANITA events with energy greater than $0.1\;\text{EeV}$ is $<10^{-6}$, thus excluding the SM interpretation at $5.8\sigma$ CL Fox:2018syq ; Romero-Wolf:2018zxt .

In the gauge-mediated SUSY breaking models, the long-lived NLSP stau seems like an ideal candidate to explain the ANITA anomalous events Connolly:2018ewv ; Fox:2018syq . Its cross section for nuclear interactions at $\sim$EeV energies is roughly 10 pb, three orders of magnitude less than the total neutrino cross section of 15 nb at 1 EeV, so that its mean free path through the Earth is several thousand km, while allowing for a reasonable branching ratios (BR $\lesssim 10^{-4}$) via UHECR neutrino-nucleon interactions Albuquerque:2003mi ; Albuquerque:2006am . The stau NLSP then propagates across the required chord length of $\sim$ 7000 km and, for appropriate mass and lifetime, decays to a tau lepton plus LSP prior to Earth emergence and subsequent extensive air shower.

The main problem however is the energy loss of a heavy, charged particle like stau along its journey through the Earth due to ionization and radiative losses. At production, due to its high boost factor of ${\cal O}(10^{7})$, the stau will have a lifetime that is heavily dilated in the lab frame, which however gradually decreases as it loses energy along its path. This was used as an advantage in Connolly:2018ewv to effectively reduce the boost factor, thus lowering the lab-frame decay time, causing the particle more likely to decay at the end of the chord length close to the surface of the Earth. At about 7000 km chord length (corresponding to the incident angle of the ANITA events), the rate of decay probability per distance $\frac{d\rm P_{\rm decay}}{dx}\sim 10^{-4}\ {\rm km}^{-1})$ for a TeV-scale stau Connolly:2018ewv . This could be translated into a decay probability around $0.01\%$ assuming that the decay happens within 40 km above the surface (corresponding to the flight height of ANITA) after exiting the Earth. However, this probability rate was calculated under the assumption that the incoming stau has an initial energy of $10^{22}\ \rm eV$, which is roughly two orders of magnitude above the Greisen-Zatsepin-Kuzmin (GZK) cut-off Greisen:1966jv ; Zatsepin:1966jv . Lowering the initial energy of stau would bring down the decay probability even more. Thus, there is a competing effect between the decay probability rate and the energy loss rate, which makes it difficult, if not impossible, for the long-lived stau scenario to explain the anomalous ANITA events. In any case, this is yet another interesting possibility for probing SUSY at neutrino telescopes that should be further looked into.

VI Long-lived Neutral Sparticles

If the NLSP is electrically neutral (e.g. neutralino or sneutrino) and long-lived, then it has some advantages over the charged NLSP, as far as its detection at neutrino telescopes goes. In particular, a light bino in RPV SUSY turns out to be a natural candidate for explaining the anomalous ANITA events Collins:2018jpg . The $LQD$ (LLE)-type RPV couplings allow the on-shell, resonant production of a TeV-scale squark (slepton) from neutrino-nucleon (electron) scattering inside Earth, see Fig. 4, thereby naturally enhancing the signal cross section. The squark/slepton decay inside the Earth can produce a pure bino which interacts with the SM fermions only via the $U(1)_{Y}$ gauge interactions and heavier SUSY particles, and therefore, can easily travel through thousands of km inside the Earth without significant energy loss, unlike tau, stau or any other electrically-charged particle. 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 Collins:2018jpg . For $LLE$-type RPV couplings [cf. Eq. (11)], the bino decays to a $\tau$-lepton (if kinematically allowed) and electron, either of which could induce the air shower seen by ANITA, whereas for $LQD$-type couplings [cf. Eq. (4)], the bino directly decays to quarks and neutrino, which mimic the SM $\tau$-decay. In the latter case, there exist some parameter space for which no throughgoing track events are predicted at IceCube. This might explain why IceCube has not seen either of the EeV events seen by ANITA, despite of its higher exposure time.

Considering the $LLE$-type interactions first [cf. Eq. (11)], the resonantly-produced slepton can decay back to an electron and neutrino through the same RPV interaction in Eq. (11) or to the corresponding lepton and neutralino through gauge interactions. The slepton might in principle be from any generation, though here we will make the assumption that the slepton is a stau, and also assume the lightest neutralino ($\chi_{1}^{0}$) to be much lighter than the stau, so that the decay $\widetilde{\tau}\to\tau\widetilde{\chi}_{1}^{0}$ is kinematically allowed. All other sparticles are assumed to be heavier than the stau and do not play any role, except for the gravitino ($\widetilde{G}$), which could be the LSP and plays the role of DM in this scenario. Gravitino LSP and bino NLSP can be realized, e.g. in natural gauge mediation without gaugino unification Barnard:2012au .

The cross-section for $\nu e\to\widetilde{\tau}\to\tau\widetilde{\chi}_{1}^{0}$ is given by Carena:1998gd

[TABLE]

which can be approximated by a Breit-Wigner formula close to the $s$-pole with $s\to m^{2}_{\widetilde{\tau}}$:

[TABLE]

where $g^{\prime}\equiv e/\cos\theta_{w}$ is the $U(1)_{Y}$ gauge coupling, 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 $\nu_{e}:\nu_{\mu}:\nu_{\tau}$ and similarly for antineutrinos, as expected for a typical astrophysical neutrino flux with $1:2:0$ flavor composition at the source Learned:1994wg .

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 $\tau$-lepton and an off-shell stau, leading to a 3-body final state: $\widetilde{\chi}_{1}^{0}\to\tau^{+}\widetilde{\tau}^{*-}\to\tau^{+}e^{-}\bar{\nu}$. In principle, the upgoing shower may be initiated either directly by the electron, or by the subsequent decay of the $\tau$. In the limit $m_{\widetilde{\chi}_{1}^{0}}\ll m_{\widetilde{\tau}}$, the 3-body decay rate can be estimated as

[TABLE]

The chord length of $\sim 7000$ km required for the two ANITA events translates into the mean lifetime of bino in the lab frame as $\tau_{\widetilde{\chi}_{1}^{0}}^{\rm lab}\sim 0.022$ s. 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_{\nu}\sim$ 2 EeV. Then the resonance condition fixes the stau mass: $m_{\widetilde{\tau}}=\sqrt{2E_{\nu}m_{e}}\simeq 2$ TeV. Substituting this in Eq. (20), we find that for a typical value of $|\lambda|\sim 0.1$, as allowed by current experimental constraints Kao:2009fg , one needs a light bino of mass $m_{\widetilde{\chi}_{1}^{0}}\sim 8$ GeV.

A more accurate calculation, taking into account all statistical interaction/decay probabilities for the neutrino, bino and tau, yields the $3\sigma$ allowed range in the $(m_{\widetilde{\chi}_{1}^{0}},|\lambda|)$ plane that explains ANITA, as shown by the yellow shaded region in Fig. 5 (left panel) Collins:2018jpg . Here we have assumed an anisotropic UHE neutrino source with flux $\Phi_{\nu}\sim 2\times 10^{-20}\ (\rm GeV\cdot cm^{2}\cdot s\cdot sr)^{-1}$ and the mass of the RPV-SUSY mediator stau at $\rm m_{\widetilde{\tau}}=2$ TeV. 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 lab-frame lifetime of $\tau_{\widetilde{\chi}_{1}^{0}}^{\rm lab}\sim 0.022$ s. The horizontal purple and red shaded regions are excluded from the $R_{\tau}$ measurements Kao:2009fg . The vertical shaded region is the kinematically forbidden region for the bino to decay into a $\tau$-lepton. The stau mass is roughly fixed to lie in the $1\text{--}2\;\text{TeV}$ region by the requirement that $m_{\tilde{\tau}}=\sqrt{s}=\sqrt{2m_{e}E_{\nu}}$, and that $E_{\nu}$ should be a few times larger than the observed cosmic ray energy of $0.2\text{--}1\;\text{EeV}$. Such a particle may be observed in current or future collider experiments. The current LHC lower limits on the stau mass in the RPV-SUSY scenario is about 500 GeV, derived from multilepton searches Aad:2014iza .

Now consider the $LQD$-type interactions [cf. Eq. (4)]]. For simplicity, we only take the first-generation squark in the intermediate state. As for the initial state quarks, both $d$ and $s$ quark contributions turn out to be comparable. However, due to stringent constraints on the product $\lambda^{\prime}_{i1k}\lambda^{\prime}_{j2k}\lesssim 5\times 10^{-5}$ from $K$-meson decays Barbier:2004ez , we will consider either the down-quark or the strange-quark in the initial state separately, but not both simultaneously. In particular, we will only consider the $\lambda^{\prime}_{ijk}$ couplings with $j=1,2$ and $k=1$ for RH down-squark. After being resonantly produced, the bino will have a 3-body decay via off-shell down-type squark: $\widetilde{\chi}_{1}^{0}\to d\bar{d}\nu$ and $\widetilde{\chi}_{1}^{0}\to u\bar{d}e$. In this case, the final-state quarks from the 3-body bino decay hadronize to either pions or kaons, mimicking the hadronic shower induced by the $\tau$.

The total differential cross section for the neutrino-nucleon and antineutrino-nucleon interactions are respectively given by [cf. Eq. (5)] Collins:2018jpg

[TABLE]

The Breit-Wigner resonance is regulated by the squark widths555The RH down-squark has two RPV decay modes: $\widetilde{d}_{kR}\to\nu_{iL}d_{jL},\,e_{iL}u_{jL}$, whereas the LH down-squark has only one: $\widetilde{d}_{kL}\to\nu_{iL}d_{jR}$. Similarly, for the $R$-parity conserving decays $\widetilde{d}\to d\widetilde{\chi}_{1}^{0}$, the RH squark coupling is twice that of the LH squark (due to different hypercharges).

[TABLE]

Note that the resonance condition is satisfied for the incoming neutrino energy $E_{\nu}=m^{2}_{\widetilde{d}}/2m_{N}x$, but due to the spread in the initial quark momentum fraction $x\in[0,1]$, the resonance peak is broadened and shifted above the threshold value $E_{\nu}^{\rm th}=m^{2}_{\widetilde{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 can be seen from Fig.5 (right panel).

Our results for the $3\sigma$ preferred region in the $LQD$ case are shown in the right panel of Fig. 5 (yellow shaded region) for both $\nu-d$ (solid contours) and $\nu-s$ (overlapping dashed contours) initial states. The vertical shaded regions are the kinematically forbidden regions for the bino to decay into pions or kaons, corresponding to the $\lambda^{\prime}_{i11}$ or $\lambda^{\prime}_{i21}$ scenario, respectively. The horizontal shaded regions bounded by the purple and blue solid lines are excluded from the $V_{ud}$ and $R_{\tau\pi}$ measurements, respectively Kao:2009fg , which constrain the $\nu-d$ scenario. Similarly, the shaded regions bounded by the red and yellow dashed lines are excluded from the $Q_{w}$ and $R_{D_{s}}$ measurements, respectively Kao:2009fg , which constrain the $\nu-s$ scenario. Here we have chosen the RH down-squark mediator mass as $m_{\widetilde{d}}=1$ TeV. We do not include the LH squark contribution, because according to our estimates, the production cross section for a 1-TeV RH down-squark at the $\sqrt{s}=13$ TeV LHC is 6.2 fb, which is safely below the current upper limit of 13.5 fb Aaboud:2017vwy , whereas including the LH squark contribution increases the cross section to 15.5 fb. The black and green shaded regions in Fig. 5 are the exclusion regions based on a recent update of the LHC constraints on the LQD couplings $\lambda^{\prime}_{211}$ and on $\lambda^{\prime}_{221}$, respectively Bansal:2018dge .

Based on these bounds, we find that there is allowed parameter space at both $2\sigma$ and $3\sigma$ for the $\lambda^{\prime}_{i21}$ scenario ($\nu-s$ initial state), whereas for the $\lambda^{\prime}_{i11}$ scenario ($\nu-d$ initial state), there is only a small $3\sigma$ range with $\lambda^{\prime}_{\rm i11}\sim 0.07-0.1$ and $m_{\widetilde{\chi}_{1}^{0}}\sim 7-10$ GeV allowed. Increasing the squark mass moves the $2\sigma$ and $3\sigma$ contours to larger $\lambda^{\prime}$ values, which are excluded by the $V_{ud}$ and $Q_{w}$ measurements Kao:2009fg . 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. 5 might come soon from the Belle II upgrade Kou:2018nap , which could significantly improve the $R_{\tau}$ measurements.

VII Conclusion

RPV SUSY is a well-motivated candidate for TeV-scale new physics beyond the SM, while being consistent with the null results at the LHC so far. Therefore, it is important to test this hypothesis at different energy scales available to us. We have reviewed the prospects of probing RPV SUSY at neutrino telescopes using the highest-energy neutrinos given to us by Nature. We discussed three different aspects of the SUSY models that can be tested in this way: (i) If there exists a TeV-scale squark/slepton, it could be resonantly produced via neutrino interactions with Earth matter, thus leading to a potentially observable excess in the UHE neutrino events at IceCube; (ii) If there exists a long-lived charged NLSP, it could lead to distinct upward-going tracks or air showers that could be detected by IceCube and balloon-borne experiments like ANITA, respectively; and (iii) If there exists a long-lived neutral NLSP, it could also lead to extensive air shower events at ANITA-like experiments, while mostly evading the IceCube-detection (until we accumulate more data). Thus, neutrino telescopes complement the RPV SUSY searches at the energy and intensity frontiers. It would be remarkable if weak-scale SUSY was discovered in such an unexpected way!

Acknowledgments

BD thanks Chien-Yi Chen, Jack Collins, Dilip Ghosh, Werner Rodejohann, Amarjit Soni, and Yicong Sui for collaboration on projects related to the topics discussed here. He also thanks Zackaria Chacko and Lucien Heurtier for useful discussions. Special thanks are due to Yicong Sui for help with Figure 2. This work was supported by the US Department of Energy under Grant No. DE-SC0017987 and also by the Neutrino Theory Network Program under Grant No. DE-AC02-07CH11359.