R-parity Violating Supersymmetry at IceCube

The presence of $R$-parity violating (RPV) supersymmetric interactions involving high-energy neutrinos can lead to resonant production of TeV-scale squarks inside large-volume neutrino detectors. Using the ultra-high energy neutrino events observed recently at the 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 2-body final states involving leptons and quarks through the RPV couplings. With more statistics, we expect these limits to be comparable/complementary to the existing limits from direct collider searches and other low-energy processes.


INTRODUCTION
The detection of ultra-high energy (UHE) neutrinos at the IceCube [1,2] has opened a new window on the Universe. In the 4-year dataset, 54 UHE neutrino events have been observed, which constitute a 6.5σ excess over the expected atmospheric background [2]. It is imperative for both astrophysics and particle physics communities to understand all possible aspects of these UHE neutrino events reported by the IceCube Collaboration. From the astrophysics side, one needs to identify the possible extraterrestrial source(s) [3], and the underlying spectral shape [4] and flavor composition [5], of the neutrino flux. From a particle physics point of view, one can use this as a unique opportunity to test the Standard Model (SM) at energy scales that are otherwise not achievable on Earth [6]. So far, no statistically significant deviations from the SM prediction have been found in the IceCube data [2,[7][8][9][10], although many new physics scenarios have been envisaged (see Ref. [6] for an overview) to explain some peculiar features. With more statistics, if the data remains consistent with the SM predictions, one can put useful, complementary constraints on various new physics scenarios. Anticipating this, we examine the current and future prospects of a well-motivated new physics scenario, namely, R-parity violating (RPV) supersymmetry (SUSY) [11], at the IceCube and beyond.
The most general RPV superpotential in the MSSM is 2) L -singlet chiral superfields, respectively (with i, j, k = 1, 2, 3 being the generation indices) and H u is the up-type Higgs superfield.
Here we have suppressed all gauge indices for brevity. SU (2) L and SU (3) c gauge invariance enforce antisymmetry of the λ ijk -and λ ijk -couplings with respect to their first and last two indices, respectively. Since we are interested in the UHE neutrino interactions with nucleons, we will only focus on the λ ijk -couplings. Any of these 27 new dimensionless complex parameters can lead to resonant production of TeV-scale squarks at IceCube energies [68], thereby making a potentially significant contri-bution to the UHE neutrino events. Note that even without SUSY, similar resonance features in neutrino-nucleon interactions can also occur in models with TeV-scale leptoquarks [69][70][71][72][73][74][75]. 1 The λ-couplings in Eq. (1) can give rise to resonant production of selectrons from neutrino interactions with electrons, reminiscent of the Glashow resonance [76] 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. Note that with non-zero λ -couplings, we need to explicitly forbid the λ -terms, e.g. by imposing baryon triality [77], to avoid rapid proton decay [78,79]. We also ignore the bilinear terms in Eq. (1), since they do not give rise to the resonance feature exploited here. 2 Using the fact that for a given power-law astrophysical neutrino flux, there is no statistically significant resonance-like feature in the current IceCube high-energy starting event (HESE) data, we derive robust upper limits on the RPV couplings |λ ijk | as a function of the resonantly-produced down-type squark mass md k , independent of the other SUSY parameters, provided the squarks decay dominantly through their RPV couplings. With the currently available low statistics of the Ice-Cube HESE data, our bounds turn out to be weaker than the existing indirect constraints [11,[81][82][83] from precision measurements in various low-energy processes. However, with more data pouring in from IceCube, and with the possibility of a second km 3 detector, such as KM3Net [84], and even a multi-km 3 extension, such as IceCube-Gen2 [85], our projected future limits could be comparable to the existing ones and complementary to the direct probes of the λ -type RPV SUSY at the LHC [86,87], as well as other indirect searches at future low-energy neutrino experiments [88]. This should provide yet another science motivation for the nextgeneration neutrino telescopes, or at least should allow for an independent test of a possible finding in the LHC or other experiments.

NEUTRINO-NUCLEON INTERACTIONS
We start with the λ -part of the RPV Lagrangian, after expanding the superpotential (1) in terms of the super-1 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. Moreover, the recent papers analyzing the IceCube data in the context of leptoquark models have not taken into account the LHC and low-energy constraints. 2 However, they could lead to other distinct signatures (e.g. triple bang) relevant for future multi-km 3 neutrino telescopes [80]. field components: 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 down-type squarks are produced resonantly. Similarly, the s-channel process in Fig. 1(e) mediated by a left-handed down-type squark can also give a resonant enhancement to the (anti)neutrino-nucleon cross section. Here we have implicitly assumed that the R-parity conserving (RPC) squark decays to a quark and a gluino, neutralino or chargino are suppressed [89], compared to the RPV decays induced by Eq. (2). 3 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 will be sizable only for the first generation quarks, which are the predominant constituents of the nucleon, and to some extent, for the second-generation quarks. Therefore, we will ignore the contributions from the third-generation quarks. For the SM CC and NC interactions [90,91], 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 neutrinonucleon interactions, written in terms of the Bjorken scal- 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 ν is the incoming neutrino energy, E ν = E ν − E is the energy loss in the laboratory frame, E is the energy of the outgoing lepton, and where s = 2m N E ν is the square of the center-of-mass energy, g is the SU (2) L gauge coupling, and m W is the W -boson mass. It is obvious to see that in Eqs. (4) and (5), 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 involvingd and u-type quarks, which do not have interference with the RPV processes, the coefficients in Eq. (3) are simply obtained by putting λ ijk = 0 in Eqs. (4) and (5).
Similarly, for the NC processes shown in Figs. 1(c)-(f), the only non-trivial coefficients are [68] a NC where L d = −(1/2) + (1/3)x w and R d = (1/3)x w are the chiral couplings, x w ≡ sin 2 θ 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. (3) are simply obtained by putting λ ijk = 0 in Eqs. (6) and (7) and replacing L d → For neutrinoantiquark interactions, the coefficients for the NC processes can be obtained simply by crossing symmetry, Similarly, for antineutrinonucleon interactions, we can just replace the PDFs f ↔f in Eq. (3). Note the Breit-Wigner resonance form of Eqs. (4) and (6), which is regulated by the right-handed down-type squark width assuming that the only dominant decay modes ared kR → ν iL d jL (NC) andd kR → 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 left-handed down-type squark, Γd kL = Γd kR /2, sincẽ d kL → ν iL d jR is the only available decay mode. The resonance condition is satisfied for the incoming energy E ν = m 2 d kR /2m N x, but due to the spread in the initial quark momentum fraction x ∈ [0, 1], the resonance peak will be broadened and shifted above the threshold energy Fig. 2). Nevertheless, one can immediately infer that for md kR ∈ [100 GeV, 2 TeV], E th ν is in the multi TeV-PeV range, and hence, can be probed by the available IceCube HESE data.

EVENT RATE AT ICECUBE
The expected number of HESE events in a given deposited energy bin at IceCube due to the modified cross section (3) can be estimated as [9] where T is the exposure time, N A is the Avogadro number, E dep (E ν ) is the electromagnetic (EM)-equivalent deposited energy for a given incoming neutrino energy E ν in the laboratory frame, V eff (E ν ) is the effective target volume of the detector, Φ(E ν ) is the incident neutrino flux, Ω(E ν ) is the effective solid angle of coverage, and we have integrated the differential cross section in Eq. (3) over x ∈ [0, 1], including both neutrino and antineutrino initial states with all flavors; for details, see Refs. [1,2,[7][8][9][10]. As an illustration, we show in Fig. 2 [92] for the cross section calculations. 4 To illustrate the RPV contribution, we have considered md L = md R = 400 GeV, |λ 11k | = 0.4, and all other |λ ijk | = 0 as our benchmark point.
In Fig. 2, the IceCube data points, as well as the background due to atmospheric neutrinos and muons, are taken from Ref. [2]. 5 Note that atmospheric ν e events will also get modified due to nonzero RPV couplings |λ 11k |. However, since the atmospheric background is dominated by the ν µ -induced events and the event rate for atmospheric ν e is much smaller [94], we can safely ignore the |λ 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 |λ ijk | 2 /m 2 d k . 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.
From Fig. 2, one can see the small enhancement in the total number of events over the SM prediction in presence of RPV SUSY. A larger |λ 11k | will result in a more pronounced excess in some of the energy bins. Our benchmark point was partly motivated by the fact that there seems to be a small excess in the data around 100 TeV, which could in principle be explained by our RPV scenario, if it becomes statistically significant. There seems to be another excess just above 1 PeV, which is however difficult to explain in this scenario, since this would 4 The PDF uncertainties on the total cross-section are at most at 5% level [9] for the energy range of interest. Therefore, we only consider the central values of the cross sections. The flux uncertainties, on the other hand, are currently at the level of 15% [2]. 5 We have not considered here the latest through-going track signal with E dep = 2.6 ± 0.3 PeV [93], since it was not included in the analysis of Ref. [2]. require a squark mass above TeV, for which the production cross section is already small. Moreover, since the neutrino flux has a strong power-law dependence of E −2.58 , the resulting number of events in the higherenergy bins will be further suppressed, thus requiring a non-perturbative value of |λ 11k | to fit any PeV-excess.

CORRELATION WITH 0νββ
For k = 1, a larger |λ 11k |-coupling will also enhance the rate of 0νββ in nuclei, thus giving a smaller lifetime and a negative correlation with the event rate at IceCube. Including only the |λ 111 |-diagrams and ignoring all other RPV contributions [42] for simplicity, we can write down the expression for the 0νββ half-life as where G 01 = 5.77 × 10 −15 yr −1 is the phase space factor for 76 Ge [95] (which is taken here as our benchmark nucleus), m e is the electron mass, m ββ is the effective mass corresponding to the light neutrino contribution with the nuclear matrix element (NME) M ν and φ is a relative phase between the light neutrino and RPV contributions. The explicit form of the NME for the RPV contribution is [37,96] where the amplitudes of the different RPV contributions are given by [41] ηg = πα s 6 and 's denote rotations between the mass and gauge eigenbasis in the gaugino-fermion-sfermion vertices. For 76 Ge, the NMEs for the light neutrino, 2-nucleon, 1-pion and 2-pion exchange modes are respectively given by [37,96,97] Using Eqs. (9) and (10), we show in Fig. 3 the correlation between the 0νββ lifetime and the total number of excess events over the SM expectations (summed over all the bins shown in Fig. 2), as predicted with 4-year exposure at the IceCube in our RPV scenario. For illustration, we have chosen here m ββ = 5 meV so that the RPV contribution is the dominant one for the experimentally accessible range and the relative phase φ in Eq. (10) does not play any role. For larger values of m ββ , one could have either a constructive or destructive interference between the light neutrino and RPV contributions in Eq. (10), depending on the phase φ. For the SUSY spectrum, we have fixed the squark and slepton masses entering into Eqs. (12) at a common value of 400 GeV, whereas the gluino and neutralinos are assumed to be much heavier around 10 PeV, in order to avoid the stringent limits on 0νββ half-life in at least part of the parameter space shown in Fig. 3, which could still yield an observable excess at the IceCube. The blue, solid curve is obtained by varying λ 111 (from 0.1 to 2), whereas the band around this curve shows the 1σ uncertainty in the best-fit neutrino flux [2]. 6 We find that for our chosen benchmark values of the squark and gaugino masses, |λ 111 | 0.7 (marked by the intersection of the blue and red solid lines in Fig. 3) is excluded from the combined lower limit of T 0ν 1/2 > 3 × 10 25 yr from GERDA phase-I+Heidelberg-Moscow+IGEX [99], and this could be improved up to |λ 111 | 0.5 (marked by the intersection of the blue solid and red dashed lines in Fig. 3) with the projected sensitivity of GERDA phase-II [100]. Thus, Fig. 3 implies that we cannot expect a large λ 111contribution to the current or future IceCube HESE data due to the 0νββ constraints. Similar conclusions can be derived for other λ ijk -contributions by considering the corresponding limits from other low-energy processes, as demonstrated in the following section.
In principle, one can also consider the RPV-assisted long range contributions to 0νββ, mediated by a light neutrino and a squark, which are independent of the gaugino mass. This contribution is mostly relevant for the coupling product λ 113 λ 131 , which depends on the left-right sbottom mixing matrix. However, it is strongly constrained by B-physics and light neutrino mass observables [41], and therefore, can not give rise to a significant effect at IceCube.

UPPER LIMIT ON |λ ijk |
Since no statistically significant excess over the SM prediction is seen in the current IceCube data (cf. Fig. 2), we use this information to put an upper bound on the |λ ijk | couplings. To this effect, we perform a binned likelihood analysis [101] with the Poisson likelihood function where the observed count n i in each bin i is compared to the theory prediction λ i , including the RPV contribution induced by λ ijk . We then construct a test statistic from which a 1σ (2σ) upper limit on |λ ijk | corresponding to the value of −2∆ ln L = 1 (2.71) can be derived. Here L max represents the likelihood value obtained with λ ijk = 0, i.e. including only the SM contribution in the analysis.
Our results for the conservative 1σ limits are shown in Fig. 4 (blue solid curves) for |λ 11k | and |λ 12k |, i.e. for electron-type neutrinos interacting with the 1st and 2nd generation quarks, respectively. As expected, the limits for the 1st generation are stronger, since the corresponding cross sections are larger due to the large valence-quark content of the nucleon at the high values of Bjorken-x required to resonantly produce squarks of significant mass.
There exist stringent limits on |λ 11k | from direct searches in e ± p collisions at HERA with √ s = 319 GeV [102,103], as shown by the magenta-shaded region in Fig. 4(a). Squark masses below 100 GeV or so are disfavored from direct searches for RPV SUSY at LEP [104][105][106], Tevatron [107,108] and LHC [13,14,87], and therefore, are not considered here. 7 In addition, the recent search for scalar leptoquarks at the 13 TeV LHC with 3.2 fb −1 data [109] is relevant for our RPV scenario, since λ ijk -couplings also give rise to the same e i e i jj final states via pair-production of down-type squarks, followed byd kR → e iL u jL which has a branching ratio of 0.5. The corresponding 95% CL ATLAS limit of 900 GeV on the first-generation scalar leptoquark mass can 7 In 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 md − |λ ijk | plane in a model-independent manner, and therefore, we do not attempt to show them in Fig. 4. be directly translated into a lower bound on the downtype squark mass, as shown by the vertical dot-dashed line in Fig. 4. There also exist indirect constraints on |λ 11k | from lepton universality in pion decay, measured by the ratio R π = BR(π − →e −ν e ) BR(π − →µ −ν µ) , unitarity of the CKM element V ud and atomic parity violation [82], the most stringent of which is shown in Fig. 4(a) by the red dotted curve. Other low-energy constraints, such as neutrino mass [31], electric dipole moment [110] and flavorchanging B-decays [111][112][113], always involve the product of two independent RPV couplings, and hence, are not applicable in our case. Moreover, for k = 1, we have an additional constraint from 0νββ, as discussed in the previous section, which however depends on the masses of other SUSY particles, 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 show the 0νββ limits in Fig. 4(a) for two benchmark points with mg = m χi = 10 TeV (gray, dashed) and 10 PeV (gray, dotted), while keeping all sfermion masses equal to md. In the former case, the 0νββ limit is the most stringent one, whereas for either heavier gaugino masses or k = 1, the pion decay constraint is stronger than the limit obtained from IceCube in the entire mass range considered. Similarly, as shown in Fig. 4(b), for |λ 12k |, the indirect constraints from lepton universality in neutral and charged D-meson decays, measured by the ratio R D = BR(D→Keνe) BR(D→Kµνµ) , are stronger than the IceCube limit derived here. This rules out the possibility of any |λ 1jk |-induced observable excess in the 4-year IceCube data.
However, one should note that these are the first-ever IceCube constraints on RPV couplings, and at present, are mostly limited by statistics, which is expected to improve significantly with more exposure time. To illustrate this point, we just scale the current 4-year dataset by a factor of 4 (roughly corresponding to 15 years of actual data taking) in all the bins analyzed here and derive projected limits on |λ 1jk | (with j = 1, 2), following the same likelihood procedure described above. This conservative estimate of the future limits is shown in Fig. 4 by the blue dashed curves. We find that the limit on |λ 1jk | can be improved roughly by up to a factor of 3 with 15-yr IceCube data, 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 in a future super-tau-charm factory [114]. In practice, however, we may not have to wait for 15 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 and  IceCube-Gen2 might go online at some point, 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. However, we expect the corresponding limits on |λ ijk | (with i = 2, 3 and j = 1, 2) to be weaker than the limits on |λ 1jk | shown in Fig. 4 simply due to the fact that the effective fiducial volume at the IceCube is the largest for ν e [1]. Nevertheless, with more statistics, one could in principle consider the |λ 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 |λ 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 |λ 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.
This will give rise to a selectron resonance from (anti)neutrino-electron interactions at IceCube. However, the corresponding threshold energy E th ν = m 2ẽ 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 [106,115], we cannot probe the λ 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 powerlaw flux, without having a significantly larger number of events in all the preceding lower-energy bins.

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. Using the 4-year IceCube HESE data in the multi TeV-PeV range, we have derived the first IceCube upper limits on the RPV couplings |λ ijk | as a function of the mass of the resonantly-produced downtype squarks (see Fig. 4). Although weaker than the existing limits from low-energy processes, our limits are expected to be significantly improved with more statistics in future, thereby complementing the RPV SUSY searches at the energy and intensity frontiers. ACKNOWLEDGMENTS B.D. is grateful to Chien-Yi Chen and Amarjit Soni for many useful discussions on IceCube. The work of B.D. is supported by the DFG grant RO 2516/5-1 and W.R. is supported by the DFG in the Heisenberg programme with grant RO 2516/6-1. B.D. also acknowledges partial support from the TUM University Foundation Fellowship, the DFG cluster of excellence "Origin and Structure of the Universe", and the Mainz Institute for Theoretical Physics during various stages of this work. D.K.G. would like to thank Texas A & M University for the hospitality, where the final stage of this work was done.