PeV scale Supersymmetry breaking and the IceCube neutrino flux

The observation of very high energy neutrino events at IceCube has grasped a lot of attention in the fields of both astrophysics and particle physics. It has been speculated that these high energy neutrinos might originate either from purely conventional astrophysical sources or from the late decay of a super heavy (PeV scale) dark matter (DM) particle. In order for decaying DM to be a dominant source of the IceCube high-energy neutrinos, it would require an unusually suppressed value of the coupling of DM to neutrinos. We attempt to explain this small coupling in the context of an $R$-parity conserving minimal supergravity model which has right-handed neutrino superfields. With the main assumptions of super-partner masses at the PeV scale and also a reheating temperature not much larger than the PeV scale, we find in our model several natural order-of-magnitude"miracles", (i) the gravitino is produced via freeze-in as a DM candidate with the correct relic density (ii) the right-handed (RH) sneutrino makes up only a tiny fraction ($10^{-6})$, of the present day energy density of the universe, yet its decay lifetime to the gravitino and neutrinos is such that it naturally predicts the right order-of-magnitude for the IceCube neutrino flux. The long lifetime of the RH sneutrino is explained by the existence of a global $R$-symmetry which is only broken due to supersymmetry breaking effects. Our model also predicts a flux of 100 TeV gamma rays from the decaying RH sneutrino which are within the current observational constraints.


Introduction
For the past few decades, supersymmetry (SUSY) has been a dominant paradigm to address almost all the shortcomings of the Standard Model (SM) [1]. However, so far direct searches at the Large Hadron Collider (LHC) [2, 3] for weak scale supersymmetry have produced only null results, which has forced us to reconsider whether naturalness is the prime motivation for predicting the SUSY breaking scale. Indirect searches from flavor physics observables for generic SUSY breaking scenarios with large flavor mixings have long constrained the SUSY breaking scale to be greater than 100 TeV [4]. These strong constraint have forced theorists to work in scenarios with minimal flavor violation that are arguably less generic, such as gauge mediated SUSY breaking (see for example [5][6][7] and references therein).
However, SUSY is still a theoretically attractive and predictive paradigm. If one takes the message from LHC and flavor physics searches seriously, then the SUSY breaking masses should lie beyond the 100 TeV scale. Even if SUSY does not solve the hierarchy problem, it has many other theoretically attractive features, such as the prediction of improved gauge coupling unification compared to the SM [8][9][10] among other features.
If DM decay is the correct explanation for the high energy neutrino events, it forces us to consider a PeV scale decaying DM candidate in the theory. However, a PeV scale DM candidate cannot be produced via thermal freeze-out, as the cross-section would be limited by the unitarity bound [44]. Thus a PeV-scale neutralino would not be an ideal candidate to explain the observed abundance of DM and simultaneously predict the IceCube neutrino flux. This observation compels us to consider other potential DM candidates e.g. gravitinos, axinos etc., which could be produced using a non-thermal mechanism to obtain the correct DM relic abundance. In addition to this, if the decay of DM is considered to be the dominant source of high-energy neutrinos, one requires an extremely high lifetime of the DM (dominantly decaying into one or more neutrinos) in order to predict a neutrino flux that is in good agreement with IceCube data. The requirement of the high lifetime ∼ O(10 28 − 10 30 ) seconds of DM translates into an unusually suppressed value of the DM-neutrino interaction strength which is very challenging to explain theoretically. Refs. [18,[22][23][24][25][26]28] attempted to explain the observed neutrino flux by allowing for a fine-tuned value of the DM-neutrino interaction coupling in their models. There have however been a few attempts which try to explain the origin of a small value of the DMneutrino coupling. In refs. [18,31], the required coupling is obtained by considering tiny Rparity violating interactions which are "technically natural". One can also consider explicit continuous global symmetry breaking parameters that can be chosen to be small (e.g. ref. [25] introduces a small chiral symmetry breaking mass for an electron partner and the right chiral electron field). Ref. [21,27] generate the small coupling through higher dimensional operators which are suppressed by some heavy mass scale, such as the Grand Unified Theory (GUT) scale or Planck scale.
In this work, we propose a simple N = 1 supergravity model which includes the particles of the Minimal Supersymmetric Standard Model (MSSM) [45] and an additional righthanded (RH) neutrino superfield. Our main assumption is that the superpartner masses are at the PeV scale. The gravitino is assumed to be the lightest supersymmetric particle 1 The idea that decay of superheavy DM may lead to high energy cosmic ray signatures predates the IceCube data and has been studied in the literature in both non-supersymmetric [32,33] and supersymmetric models [34].
(LSP) and the RH sneutrino 2 is assumed to be the next-to-lightest superparticle (NLSP). Our model assumes the existence of a global U (1) R -symmetry with specific R-charge assignments. A well known motivation for imposing such a global R-symmetry is that the breaking of this symmetry due to SUSY breaking effects can explain the small value of the µ parameter, hence addressing the infamous "µ-problem" via the Giudice-Masiero mechanism [46]. We will argue that in the presence of R-parity and the additional R-symmetry, the RH sneutrino becomes a quasi-stable candidate with decays forbidden by the approximately preserved R-symmetry. However, after the breaking of R-symmetry at the PeV scale, the non-renormalizable operators involving the sneutrino allow for it to decay into a gravitino and other SM particles with an extremely long lifetime. Since the dominant decay mode of the sneutrino happens to produce a neutrino as one of the final decay products, it is an attractive candidate to explain the PeV scale neutrino events.
We will see several "miracles" that occur naturally under this set-up, (i) the gravitino is predicted to make up most of the dark matter and naturally has the correct order-ofmagnitude relic abundance from freeze-in production, (ii) we will see that the sneutrino is predicted to contribute to less than 10 −6 of the present day energy density of the universe, however the value of the decay lifetime will naturally result in a flux spectrum normalization that can explain the PeV scale neutrino events observed at IceCube. Furthermore, we will show that the same set-up can solve the µ-problem and can also explain the small value of the neutrino masses. This paper is organized as follows: In §2, we first discuss the details of our N = 1 supergravity model. We write down the leading R-symmetry conserving non-renormalizable operators appearing in the superpotential and Kähler potential for this model. In §2.1, we present the spectrum of superpartners along with some benchmark values and we work out the physical mass eigenstates of the sneutrino by diagonalizing the mass matrix involving off-diagonal mixing terms coming from soft SUSY breaking trilinear couplings. In §2.2, we briefly discuss the conditions under which one can obtain a SM-like Higgs with PeV-scale SUSY. In §2.3, we calculate the mass of the SM-like left-handed neutrino in this model, which turns out to be purely Majorana in nature. In §3, we calculate the decay width of the right-handed sneutrino and show that is has a lifetime longer than the age of the universe, indicating that it can act as a viable decaying DM candidate. In §4, we calculate the relic abundance of both the gravitino and the sneutrino produced non-thermally from the scattering and decay of heavier supersymmetric particles. In §5, we compute the expected neutrino flux from the decay of sneutrino DM along with a simple unbroken power-law astrophysical source, and we show that the combination explains the IceCube data well. We also show that our results are consistent with current constraints from high-energy gamma ray observations. In §5, we summarize our results. In appendix §A, we review the mechanism that leads to both R-symmetry and SUSY breaking in our model.

Outline of the model
In this section, we start by presenting the details of our supergravity model of PeV scale SUSY breaking. As mentioned in the introduction, the model assumes an additional global U (1) R symmetry. We will first discuss why such a U (1) R symmetry has been considered desirable in previous works.
Usually it is assumed that SUSY is broken in the hidden sector at a very high scale and SUSY breaking effects are transmitted to the visible sector via gravitational interactions. Therefore the coefficient of the SUSY breaking operators, such as the masses of superpartners, are determined by non-renormalizable operators governing the interactions between the visible and hidden sectors. Infamously, there is an ambiguous µ parameter present in the superpotential which does not depend on the dynamics of SUSY breaking. Its natural value can be either zero (if protected by some symmetry) or on the order of the Planck scale. However, the value of the µ parameter needs to be tuned to the order of the SUSY breaking scale (or smaller) in order for electroweak symmetry breaking (EWSB) to occur successfully. This is the well known "µ problem" in SUSY theories. A solution suggested by Giuidice and Maisero [46] in the context of N = 1 supergravity theories imposes an additional global U (1) R symmetry. Their idea was based on the simple observation that if the Higgs superfields are vector-like under the U (1) R symmetry, then the term µH u H d would be forbidden in the superpotential, and EWSB can successfully occur. However, in order to generate higgsino masses, they also included Planck suppressed terms like (X † /M p )H u H d in the Kähler potential, with X being a hidden-sector superfield (SUSY breaking spurion) and M p being the Planck mass scale.
Subsequently, the Giuidice-Maisero mechanism was extended to involve a SM singlet RH neutrino in refs. [47][48][49] with the motivation of realizing a TeV scale see-saw mechanism of neutrino mass generation. This scenario predicted a TeV-scale Majorana neutrino which could be observed at the LHC. In these models, the charges of the right-handed neutrino superfields (N c ) under the global U (1) R symmetry prevented Majorana mass terms from showing up in the superpotential (similar to the way the Higgs bilinear term is forbidden), but allowed for the non-renormalizable operator (X † /M p ) N c N c in the Kähler potential, thus giving rise to a right-handed neutrino mass on the order of the SUSY breaking scale. Interestingly, this mechanism explains the small value of the observed neutrino masses via the see-saw mechanism. In the model discussed in ref. [47], the neutrino mass is generated from non-renormalizable terms in the Kähler potential while the right-handed neutrino Yukawa couplings and corresponding soft SUSY breaking trilinear coupling parameters are generated from other non-renormalizable operators in the superpotential.
We draw motivation from these models above, to build a model of PeV scale SUSY breaking which possesses a U (1) R symmetry and has right-handed neutrino superfields in addition to the superfields of the MSSM. In order to generate a quasi-stable, long-lived righthanded sneutrino, we assign different R-charges to the fields following which all the soft SUSY breaking operators involving the RH neutrino superfield, including Yukawa terms, are generated from higher dimensional operators in the Kähler potential. By imposing the R-charges R( N c ) = R( H u ) = R( H d ) = 1/2, R( X) = 1, R( L) = R( Q) = 0, R( e c ) = R( u c ) = R( d c ) = 3/2, one does not get any renormalizable terms in either the Kähler or superpotential involving direct interactions between RH neutrinos and other MSSM superfields. The non-vanishing interaction terms thus appear from the following R-symmetry conserving non-renormalizable operators in the Kähler potential: where M p is the Planck scale and M * is the scale obtained by integrating out other heavy states (assumed to be present in the theory).
In the 'Planck-scale mediated' SUSY breaking scenario, if SUSY is broken in the hidden sector by an F -term vacuum expectation value (VEV) F X , soft mass terms in the visible sector are of the form The spontaneous breaking of local supersymmetry implies the existence of a massive spin-3/2 gravitino. Its mass is also given by In order for the sparticle masses to be at the PeV-scale, one needs to assume F X ∼ 10 12 PeV 2 . This gives rise to soft SUSY breaking mass terms and a gravitino mass at the PeV scale. The masses of the right-handed neutrino as well as the higgsino appear from the nonrenormalizable operators given by the first two terms of the Kähler potential. Upon considering the F -term VEV, the RH neutrino mass and the higgsino mass parameter (µ) turn out to be of the same order, The Yukawa coupling for the right-handed neutrino appears through the non-renormalizable operator X † M 2 P N c L H u in the Kähler potential. By writing the superfield in terms of component fields, we find that the Yukawa coupling is given by for an assumed coupling α N ∼ O(0.1). Further we can see that the operator ( X † /M 2 p ) N c L H u in the Kähler potential does not produce any trilinear scalar coupling. However, the trilinear term can arise from the sub-leading non-renormalizable operator given by ( For X ∼ m 3/2 (a justification for expecting this particular VEV for the scalar component of the SUSY breaking spurion field is given in appendix A) and F X ∼ 10 12 PeV 2 , A N is given by for an assumed coupling β N ∼ O (1). Similarly the lepton-number violating B-term coefficient B N (which is the coefficient of N N in the Lagrangian) as well the Higgs mixing parameter B µ (which is the coefficient of H u H d in the Lagrangian), will be given by We note that in addition to the soft SUSY breaking terms we have estimated, we also need a mechanism to generate gaugino masses and slepton/squark trilinear interaction terms (A-terms). We will not discuss in detail the mechanism for generation of these terms, but one way these could be generated is by assuming that in addition to hidden-sector field X with R-charge R X = 1 as considered in this paper, there could be another hidden sector superfield Y present in the theory with R Y = 0. The non-renormalizable interactions of this field (in both the Kähler potential and superpotential) would generate the gaugino masses and trilinear terms to be of the order of PeV. We note that this mechanism would not however, generate a trilinear term involving the RH sneutrino (A N ) because the R-charge of the sneutrino superfield forbids the corresponding Kähler terms.
Effect of one-loop corrections: It is evident that values of all the parameters we have estimated are at the hidden sector SUSY breaking scale ( √ F X ∼ 10 6 PeV). In order to obtain their physical values at the low energy (EW/PeV) scale, one has to consider renormalization group (RG) evolution of the mass parameters and couplings. Since all the mass parameters in our model (except B µ and B N ) are around the PeV-scale, it is reasonable to assume that there will only be an O(1) change in their values after RG evolution down to the low energy scale. Similarly, we expect that all the couplings except the suppressed ones involving right-handed neutrinos (y N and A N ) will exhibit an O(1) change after their RG evolution down to the EW scale.
Next we check whether the suppressed parameters y N , A N , B µ , B N get enhanced at the one-loop level due to other dominant couplings and PeV-scale mass parameters present in the model. Using the generic expressions for one-loop RG equations as given in [45], we have checked that the value of Yukawa coupling y N will remain mostly unchanged after its RG evolution down to the low energy scale. However, the radiatively corrected trilinear scalar coupling A N will receive a dominant correction proportional to the Yukawa coupling y N given by, The numerical estimate above holds if we take y N ∼ 10 −13 and the electroweakino mass M a ∼ PeV, with g a being the corresponding gauge coupling. Similarly, by using the generic expressions for one-loop RG equations as given in [45], it can be checked that the value of B N will not get any significant contribution even at one-loop. However, the radiatively corrected soft SUSY breaking mixing parameter B µ receives a dominant contribution given by (2.9)

Parameters
Dependence Numerical values on X and F X Scalar/squark mass (m i ) The values of the higgsino mass parameter (µ) and the electroweak gaugino masses (M a ) can be chosen in such a way that the one-loop corrected B µ turns out to be on the order of the PeV scale. Thus in summary, this model gives PeV scale mass to all supersymmetric partners. The imposition of R-symmetry gives rise to suppressed values of both y N and A N , while the trilinear couplings related to all other supersymmetric partners will still remain on the order of the PeV scale. Table 1 shows the typical scales of the SUSY breaking parameters and the RH neutrino/sneutrino couplings in our model. There are three questions that naturally arise given the parameters of our PeV scale SUSY breaking model in the table.
• What is the mass spectrum of the superpartners? In particular, which are the lightest and next-to-lightest superpartners?
• Can we get electroweak symmetry breaking and the observed Higgs mass of 125 GeV with PeV scale SUSY breaking?
• Given that the Yukawa couplings of the RH neutrinos are extremely suppressed, is it possible to obtain the right order-of-magnitude for the neutrino masses via the see-saw mechanism?
We will address these three issues in the remainder of this section.

Mass hierarchy, sneutrino mass eigenstates and benchmark spectrum
We assume in our model that the gravitino is the LSP after considering RG evolution of the SUSY breaking parameters. Since R-parity is conserved in this model, it is natural to then expect that the gravitino might exist as a viable DM component. We will assume that the lightest RH sneutrino is the NLSP. We will show in §3 and §4 that the lightest RH sneutrino can also co-exist as another quasi-stable, long-lived DM component (albeit with a much smaller contribution to the relic density) because of its almost vanishing Yukawa and trilinear scalar couplings. In §5, we will show that the decays of the RH sneutrino can explain the high energy neutrino flux observed at IceCube. We also assume that the mass hierarchy of the remaining superpartners follows the trend M a > m i > mh u > m N > mg > mÑ > m 3/2 . Here, M a is an electroweakino mass, m i is a slepton/squark mass parameter (other than RH sneutrino), mh u is the uptype higgsino mass, m N is the RH neutrino mass, mg is the gluino mass, mÑ is the RH sneutrino mass 3 and m 3/2 is the gravitino mass. While the exact hierarchy of particles other than the gravitino and RH sneutrino is not absolutely critical for our main results, working with a specific hierarchy allows us to streamline our discussion in subsequent sections.
The masses of all the superpartners including the sneutrino are naturally on the order of the PeV scale. However, in the case of the sneutrino, the small trilinear scalar parameter A N leads to a tiny mixing between left-handed and right-handed sneutrinos. In order to calculate the physical eigenstates and their eigenvalues, we need to diagonalize the sneutrino mass matrix. The full set of SUSY breaking parameters corresponding to the neutrino and other visible superfields (including all three matter generations) is given by where the couplings are 3 × 3 matrices in generation space. The sneutrino mass matrix can be written in the following basis [50], (2.11) In this basis, we find the mass matrix to be of the form Here, ∆ 2 ν = (m 2 Z /2) cos 2β is the D-term contribution. We note that the term v u A N in the mass matrix results inν L ↔Ñ R mixing. Ignoring the phases present in Eq. 2.12 which can lead to CP violating effects, we writeν L andÑ R in terms of real fields as (2.14) With this, the mass matrix reduces to block diagonal form as Diagonalizing the above matrix by performing a unitary transformation, we get the physical mass eigenstates, By considering a single generation of neutrino superfields and neglecting flavor mixings for simplicity, the mixing angle is given by , (2. 16) with the top (bottom) sign for i = 1 (i = 2). Thus, the physical mass eigenstates of the sneutrino are given byÑ 1 ,Ñ 2 ,ν 1 andν 2 . We note that the mixing between left and righthanded sneutrinos is very small, since the value of the mixing angle based on the generic parameter values in Table 1 can be estimated as θν i ∼ 10 −19 . In the next few sections we will show that the lightest eigenstate (presumed to beÑ 1 ) can exist as the appropriate decaying, sub-dominant DM component in this model. In later sections, we will consider the benchmark spectrum: M a = 10 PeV, mh u = 9.5 PeV, m N = 8 PeV, mg = 7.5 PeV, mÑ 1 = 6.5 PeV, and m 3/2 = 1.5 PeV when evaluating numerical results.

EWSB, Higgs mass and PeV scale SUSY breaking
The Higgs potential for the two Higgs doublets takes on the form of the standard MSSM Higgs potential. Given that all the parameters in the Higgs potential are on the order of the PeV scale, it is possible to obtain the SM VEV of v = 246 GeV with fine-tuning [10,51]. In terms of the SUSY breaking parameters, this fine-tuning shows up as a restriction on the allowed region of parameter space. Thus, as one would expect, we can obtain successful EWSB at the cost of a little-hierarchy problem.
In addition, if we require the light Higgs to have a mass of ∼ 125 GeV, this would lead to further restrictions on the space of allowed parameters. Below, we demonstrate the existence of a SM-like light Higgs boson for a suitable choice of SUSY breaking parameters.
The 2 × 2 mass-squared matrix for the neutral CP-even Higgs sector is given by [52], where the entries M 2 ij correspond to the tree level masses of the neutral CP-even Higgs doublets and the terms ∆M 2 ij take into account the loop corrections. The dominant corrections arise from one-loop corrections which are controlled by the top/stop couplings, subleading one-loop corrections which are controlled by bottom/sbottom couplings and two-loop corrections due to the top/bottom Yukawa couplings and the strong coupling constant. The tree level mass-squared matrix for the neutral CP-even Higgs sector is given by where m Z is the Z-boson mass, tan β is the ratio of up-type to down-type Higgs VEVs, m A is the mass of the CP-odd Higgs and λ t,b are the top and bottom Yukawa couplings respectively. The shorthand notation used in these formulas is to be interpreted as follows: Here, M S is the arithmetic mean of the stop squark masses, i.e. M S = 1 2 (mt 1 + mt 2 ), m t is the top quark mass, A t,b are the stop and sbottom trilinear couplings respectively, and X t = A t −µ cot β. The coefficients c ij correspond to the leading two-loop corrections due to the top/bottom Yukawa couplings and strong coupling constant g 3 ; they are given by, 125.

140.
145.  Fig. 1, we have shown the smaller eigenvalue (which corresponds to the mass of the light Higgs) as a function of X t , for this choice of parameters. We find that for X t = ±6 PeV, we can obtain a light Higgs boson of mass ∼ 125 GeV. For this set of parameters, the heavy Higgs boson has a mass of 8 PeV. Thus, we have demonstrated that a light Higgs boson can easily be obtained in this set up. Note that since the superpartners and the other Higgs bosons are very heavy, the light Higgs boson has very nearly Standard Model like couplings. A full scan over the parameter space to find the phenomenologically acceptable SUSY breaking parameters is beyond the scope of this work, but these can be found by using numerical codes such as SUSPECT [53] and SOFTSUSY [54].

Neutrino masses
We would also like to show that our model explains the small value of the observed neutrino masses despite the suppressed neutrino Yukawa coupling y N . In generic models with RH neutrinos, the masses of the neutrinos are controlled by the following mass matrix which mixes left and right-handed neutrinos: (2.20) Here m ν L and M N are Majorana like masses for the left and right handed neutrinos and v u y N is a Dirac mixing induced by a VEV for the up-type MSSM Higgs. In the conventional see-saw mechanism, m ν L is either zero or suppressed. With this, the lightest eigenvalue of the neutrino mass matrix is given by y 2 N v 2 u /M N , with M N being the mass of the (heavy) RH neutrino. This mechanism for neutrino mass generation does not work in our model because of the extremely suppressed Yukawa coupling y N . However in the context of supersymmetric models, various alternatives to the conventional see-saw mechanism have been conjectured to originate from different non-renormalizable Kähler operators present in the theory [47,49,[55][56][57]. In our model, the dominant source of neutrino masses arises from the following non-renormalizable Kähler operator present in the theory (see Eq. 2.1): Here, we have introduced another scale M * based on the assumption that the model exhibits another symmetry at this intermediate scale. Our implicit assumption is that all the visible supersymmetric fields have some non-zero charges under this new symmetry while the supersymmetry breaking field (the hidden-sector field) is uncharged under this symmetry 4 . Therefore, any non-renormalizable interaction term involving only the visible sector superfields (and in particular the term ( L · H † d )( L · H u )), will be suppressed by M * instead of M p . On the other hand, all the non-renormalizable interactions involving the hidden-sector field (for e.g. the term that generates the neutrino Yukawa coupling in our model, X † N c L · H u ) are still suppressed by the Planck scale.

Decay lifetime of the RH sneutrino
In this section, we will calculate the decay width of the lightest right-handed sneutrino. Due to R-parity conservation and the assumed mass hierarchy: M a > m i > mh u > m N > mg > mÑ > m 3/2 , the possible decay modes of the sneutrino include: Figure 2: Feynman diagrams corresponding to the decay modeÑ 1 → ν L h 0 ψ 3/2 . The diagrams (a) and (b) show interactions which proceed through the Yukawa interaction with coupling y N , whereas diagram (c) proceeds through the trilinear interaction with strength proportional to A N . The Yukawa interaction dominates the decay process, and in fact leads to the three-body decay being dominant over the two-body decayÑ 1 → ν L ψ 3/2 which is suppressed by the small mixing angle.
We note that both these decay modes are suppressed. The first decay mode is parameterically suppressed through the small mixing between left-handed and right-handed sneutrinos (which is given by sin θν i ) while the second decay mode (depending on the relevant Feynman diagram) is suppressed either by the small Yukawa coupling (y N ) or by the trilinear interaction coupling (A N ). Of all these suppressed couplings, the Yukawa coupling 5 Other suppressed decay modes through mixing are the three body decays:Ñ1 → νLB * ,Ñ1 → e ± LW * ∓ , where the off-shell gauginosB * /W * ∓ convert into ψ 3/2 and the corresponding SM gauge bosons. These decay widths are extremely suppressed by both three body phase space as well as the small mixing angle. turns out to be the least suppressed, and thus the decay of the RH sneutrino is dominantly to the three body final stateÑ 1 → ψ 3/2 ν L h 0 . We will justify this below.
We will begin by estimating the partial decay widths. Since all the decay channels involve a gravitino as one of the final decay products, we first write the relevant gravitino (ψ µ ) interaction terms in the Lagrangian which are given by [58], where χ i stands for any chiral left-handed Weyl spinor, φ i stands for the scalar partner of χ i in the same chiral supermultiplet, λ α(a) and F α(a) νρ denote the gaugino and gauge field strength in the same vector supermultiplet. At high energies [59,60], the gravitino field can be replaced by the goldstino by a generalization of the Goldstone boson equivalence principle, ψ µ → The non-zero mixing between the left-handed and right-handed sneutrino allows for two-body decays of the typeÑ 1 → ν L ψ 3/2 . Using the interaction terms given in Eq. 3.1, we can estimate the partial decay width for this mode as, where we have used m PeV as a generic representative SUSY breaking PeV mass and in the numerical expression we have used the estimated value of sin θν i ∼ 10 −19 . The three body decay modeÑ 1 → ψ 3/2 ν L h 0 is generally expected to be phase space suppressed relative to the two-body decay. However, looking at the Feynman diagrams for this decay mode, we can see from Fig. 2 (a) and (b), that the decay proceeds through the Yukawa interaction y N ∼ 10 −13 , which is parameterically much larger than the mixing factor sin θν i 6 . The decay width for this three body mode is therefore given by: Comparing Eqs. (3.4) and (3.6), we therefore see that the three-body decay is much more dominant than the two-body decay decay mode. Thus, the lifetime of the right-handed sneutrino turns out to be, Since the lifetime is longer than the age of the universe, the RH sneutrino can act as a quasi-stable DM component co-existing with the stable gravitino DM. As it is evident that the sneutrino decays gives rise to neutrinos, the late decay of a PeV-scale sneutrino could potentially explain the flux of PeV scale neutrino events observed at IceCube. We will discuss this issue in §5.

Relic abundance
In the previous section, we found that in addition to the gravitino LSP, the real component of the lighter right-handed sneutrino (Ñ 1 ) can co-exist as a viable DM component in our model. Before we explore whether sneutrino decays can explain the IceCube PeV events, we compute the relic abundance of both the proposed DM components. Since both components have extremely suppressed interactions with the SM particles, they would not be in thermal equilibrium in the early universe. However, a non-zero relic density for both particle species can be produced from the decay or inelastic scattering of heavier superpartners existing in thermal equilibrium with the SM plasma. This process is known as freeze-in production of DM. In order to keep the discussion self-contained, we briefly review the details of the freezein mechanism from decay and scattering of heavier particles [61,62]. We will consider below the two generic processes: φ → ψ +χ and φ 1 +φ 2 → ψ +χ, with χ being the weakly-coupled DM component. We first discuss the production of χ from the decays of a heavier species φ. The evolution of the number density of χ is described by the following Boltzmann equation [61,62]:ṅ where the factors of dΠ φ j = d 3 p j /((2π) 3 2E j ) correspond to phase space measure and the factors of f j corresponds to phase space density of particles of type j. The sum over the initial and final spins in the squared matrix elements is assumed implicitly in this equation.
Assuming that the initial abundance of χ is zero, we will have f χ = 0. Thus we can neglect Pauli blocking and Bose enhancement effects, and set the factors of (1 ± f χ ) ≈ 1.
Performing the two body phase space integration over the kinematic distributions of ψ and χ and writing the result in terms of the decay width (Γ φ ) and mass (m φ ) of φ, the equation above simply reduces tȯ where g φ is the spin degeneracy of φ. Since the φ particles are assumed to be dilute and in thermal equilibrium, we can approximate f φ e −E φ /T , where T is the temperature of the kinetically-coupled φ particles. Changing the integral over momentum space into an integral over energy, one obtainṡ where K 1 corresponds to the first Bessel function of the second kind. Writing n χ in term of the yield Y χ = n χ /S (where S is the entropy density of the universe at a particular temperature T ) and utilizingṪ ∼ −HT , we obtain where the entropy density, S(T ) and the Hubble rate, H(T ) are given by S(T ) = 2π 2 g S * T 3 /45 and H(T ) = 1.66 g ρ * T 2 /M p in the radiation dominated era, and g S * and g ρ * are the effective entropy and energy degrees-of-freedom present in the radiation bath. After defining x = m φ /T , the integral turns out to be Approximating the limits of integration as x max = ∞ and x min = 0, we get xmax x min K 1 (x)x 3 dx = 3π/2. Using this, the final expression for Y χ φ→ψχ turns out to be Similarly we can calculate Y χ for the scattering process φ 1 + ψ → φ 2 + χ. In particular, we discuss the case where the scattering process involves non-renormalizable operators i.e. of the form L ⊃ 1 Λ φ 1 ψχφ 2 . At temperatures much higher than the masses of all particles involved, the squared-matrix amplitude for the given scattering process is proportional to, where s is the interaction center of mass (c.o.m.) energy squared at a temperature T . The value of κ is determined by other low energy parameters such as masses of the particles involved in the scattering process. The change in number density of χ due to this process is described by the following Boltzmann equation: Assuming that the initial abundance of χ is negligible, we set f χ = 0. Following [62], the collision term can be written as an integral with respect to the c.o.m. energy: For |M| 2 ∼ κ s Λ 2 , the equation above reduces tȯ Expressing n χ in terms of the yield Y χ = n χ /S, we can rewrite the Boltzmann equation as: Integrating this with respect to temperature by using limit of integrations to be T min = 0 and T max = T R (the reheating temperature), the final expression for Y χ φ 1 ψ→φ 2 χ is just given by indicating that the number density produced through the scattering process is linearly proportional to the reheating temperature (T R ). Next we will use the above results to calculate the relic abundance of both the gravitino as well as the RH sneutrino.
Relic abundance of the gravitinos: Since the heavy gravitinos decouple from the thermal plasma in the very early universe because of their Planck suppressed interactions, they do not have a significant population in the early universe from thermal production. The dominant number density of heavy gravitinos will therefore be produced from the decay and/or inelastic scattering of heavier particles present in the early universe. Utilizing the gravitino interaction terms given in Eq. 3.1, the possible processes involving the decay of heavier super-partners (i) into the gravitino include i → ψ 3/2 + SM. (4.14) In the high energy limit, when ψ 3/2(µ) → ∂ µ χ/m 3/2 , the decay width of any heavier particle decaying into a gravitino is given by where m i corresponds to the mass of the heavier superpartner. Using Eq. 4.6, the gravitino yield Y 3/2(decay) will be given by Similarly there are many inelastic scattering processes that contribute to the production of the gravitino. The dominant scattering processes proceed through non-renormalizable operators e.g. production of gravitinos from the process: q +g → ψ 3/2 + q, where q is a SM quark andg is the gluino. The complete list of possible scattering processes is given in [63] and the matrix amplitude is given by where mg corresponds to the mass of gluino. Further using Eq. 4.13, and using κ ∼ (1 + m 2 g /3m 2 3/2 ) ∼ m 2 g /(3m 2 3/2 ) and Λ = M p , the gravitino yield from inelastic scattering of the gluino Y 3/2(scattering) is given by The present day relic abundance of any species χ is related to the yield by [61]: where the present day entropy density, S(T 0 ) = 2π 2 g S * T 3 0 /45 ∼ (2.8 × 10 −4 eV) 3 for the current CMB temperature of T 0 = 2 × 10 −4 eV and the critical density is ρ c /h 2 = (2.95 × 10 −3 eV) 4 .
As is evident from Eq. 4.17 and Eq. 4.20, for a reheating temperature not much larger than the PeV scale, both scattering and decay processes contribute comparably to the gravitino relic abundance via freeze-in. For our benchmark mass parameters m 3/2 = 1.5 PeV, mg = 7.5 PeV, M a = 10 PeV and assuming T R 1000 PeV, electroweakino decay processes dominate over scattering processes and therefore determine the final yield of gravitinos.
The gravitino relic abundance can be estimated from Eq. 4.17 and Eq. 4.21 as This is the first order-of-magnitude miracle that we promised. For superpartners masses near the PeV scale and a reheating temperature less than ∼ 1000 PeV, we automatically get a relic abundance of gravitinos consistent with the observed amount of dark matter in the universe.
Relic abundance of sneutrinos: As we have shown in §2, the right-handed sneutrino (Ñ 1 ) has very suppressed interactions with other MSSM particles. Therefore, as in the case of the gravitino, it would also never attain equilibrium with the thermal plasma. Similar to the gravitino abundance, the relic density of the sneutrino can be produced via decay and scattering of thermally produced heavy MSSM superpartners. However, since the production via scatterings proceeds through renormalizable operators, there is no reheating temperature enhancement and this contribution is always subdominant as compared to decay processes [64]. We will therefore neglect the contribution to RH sneutrino production from the scattering processes. Below we list the possible two-body decay processes that might give a non-zero relic density for the sneutrino (Ñ 1 ): 1.h u (Higgsino) →Ñ 1 ν L (mediated by Yukawa interactions), 2.ν 1 (Heavier sneutrino) →Ñ 1 h 0 (mediated by trilinear scalar interactions), 3.ẽ L (Slepton) →Ñ 1 W ± (mediated by the off-diagonal mixing angle and gauge interactions).
As we have seen before, the Yukawa coupling is the dominant interaction coupling of the sneutrino and it determines the production rate of sneutrinos from the decay of a heavier species. The right-handed sneutrinos are therefore dominantly produced through the decay ofh u →Ñ 1 ν L which is mediated by the Yukawa interaction, and whose rate is given by Incorporating this in Eq. 4.6, the sneutrino yield YÑ 1 (decay) is given by Thus, for PeV masses of the superpartners, the RH sneutrino only makes up a tiny fraction O(10 −6 ) of the energy density of the universe today. We will assume the benchmark value of y N = 0.5 × 10 −13 in the next section.

IceCube neutrino flux
As mentioned in the introduction, the IceCube collaboration has reported 82 ultra-high energy neutrino events with deposited energies in the range from 60 TeV to 10 PeV in their six year data sets [15]. Of these ultra-high energy events, 3 events have been observed with deposited energies greater than 1 PeV. The atmospheric neutrino flux is expected to give a negligible contribution to the event rate above 60 TeV [11,14]. Therefore, the ultrahigh energy events can be interpreted as strong evidence for either hadronic astrophysical processes (pp or pγ) or the decay of superheavy DM, or both. In this paper, we postulate that the PeV scale events observed by the IceCube collaboration arise from the decay of the RH sneutrino, whereas the sub-PeV events arise dominantly from a simple unbroken power law type astrophysical flux. In this section, we calculate the combined neutrino flux from both astrophysical sources and the decay of the RH sneutrino. We will assume the benchmark spectrum for the superpartner masses from the last paragraph of §2.1. We will also see that in addition to explaining the PeV neutrinos seen at IceCube, the decays of the RH sneutrino could give rise to a flux of very high energy gamma-rays.
. Different models for a purely astrophysical origin of the high energy neutrino flux have been tested by the IceCube collaboration in [13], while considering different assumptions about the isotropy of the neutrino flux, as well as its flavor-ratio. The simplest astrophysical assumptions take an isotropic neutrino flux and a flavor-ratio at Earth for ν e : ν µ : ν τ of 1 : 1 : 1, and a simple unbroken power law spectrum (UPL) given by (5.1) Here, J UPL 0 is a normalization of the neutrino flux at 100 TeV and γ is the spectral index. The IceCube collaboration has found a best fit spectral index γ = 0.5 ± 0.09 for neutrinos with energies between 25 TeV to 2.8 PeV [14]. They find that this simple power law flux is consistent with their observations. However, the observations of three high energy events (above 1 PeV) have generated a lot of interest in the particle physics community. Such high energy events, could arise from a new source of neutrinos, such as from dark matter decay. While the data is still insufficient to discriminate between these alternative hypotheses for the origin of the high energy neutrinos, it is interesting to speculate on a possible dark matter origin that future evidence could either confirm or disprove.
We will assume that the (10 − 100) TeV events seen at IceCube arise from an UPL type astrophysical flux with parameters J UPL 0 = 1.8 × 10 −8 GeV cm −2 s −1 sr −1 and γ = 0.9. Next we will try to explain the origin of the PeV scale neutrino events as arising from RH sneutrino decays.  Figure 3: Neutrino source spectra from RH sneutrino decay through the processÑ 1 → ψ 3/2 ν L h 0 , with mÑ 1 = 6.5 PeV, m N = 8 PeV and m 3/2 = 1.5 PeV. The dashed blue line shows the neutrino spectrum obtained from the primary neutrinos produced directly from the decay of the RH sneutrino. The red dashed line shows the neutrino spectrum obtained from the decay of the Higgs into neutrinos. Finally, the dashed black line represents the total neutrino spectrum obtained from the three-body decay of the sneutrino.

Neutrino flux from sneutrino DM decay
To calculate the expected neutrino flux from the decay of the PeV scale RH sneutrino, we first need to obtain the neutrino spectrum at the source of decay. The neutrino spectrum for the processÑ 1 → ν L h 0 ψ 3/2 is determined by both direct neutrino production in sneutrino decay, as well as from the decay of the Higgs particle in the final state. The neutrino spectrum at the source is therefore given by, where the first term corresponds to the direct neutrino spectrum (with neutrino energies E ν ) and the second term gives the neutrino spectrum from Higgs decay (where the Higgs has an energy E h 0 ). The factor of gives the neutrino spectrum from Higgs decays and can be taken from the tabulated values in the PPPC4DMID code [66] which is based on Pythia [67]. In order to evaluate the full spectrum, we need to compute the differential decay widths of the sneutrino dΓ dEν and dΓ dE h 0 , as well as the total decay width Γ.
The decay width can be evaluated in terms of the decay amplitude as, Here, the spin-summed squared matrix amplitude is given by, where E 3/2 = mÑ 1 − E h 0 − E ν is the energy of the gravitino. Using this, the differential decay width with respect to E ν is where the limits of integration for fixed E ν are .
One can now perform the above integration to obtain the differential decay width, however we will not reproduce the full and rather lengthy expression here. Similarly, the differential decay width with respect to E h 0 is given by, where the limits of integration for fixed E h 0 are . We note that the endpoint for the E ν or E h 0 spectra in the above differential decay widths is given by The total decay width can be computed by integrating either of the above differential decay widths. We can numerically compute the partial and total decay widths for our benchmark sparticle spectrum (mÑ 1 = 6.5 PeV, m N = 8 PeV and m 3/2 = 1.5 PeV). The lifetime of N 1 for our benchmark point is found to be τÑ 1 ∼ 10 24 sec, consistent with our estimate in Eq. 3.7. Plugging in the numerical partial decay widths in to Eq. 5.2, we obtain the final form of the neutrino source spectrum of the RH sneutrino which is shown in Fig. 3. Using this source spectrum, we will next compute the terrestrial neutrino flux obtained from the decay of sneutrino DM accounting for both galactic and extragalactic contributions.
Extragalactic contribution: DM decays outside of our galaxy will generate neutrinos at a cosmological distance χ(z) with energy E ν and these neutrinos will then be incident on the Earth with redshifted energy E ν /(1 + z). As we have shown in the previous section, most of the DM relic density is dominated by the gravitino which does not contribute to the neutrino flux. We wish to test the feasibility of the sneutrino as a possible candidate to explain the IceCube PeV-scale neutrino events even though it contributes to only a small fraction of the relic density. The differential neutrino flux observed at Earth from extragalactic RH sneutrino decay is given by [68], where E ν = (1 + z)E ν is the energy of the neutrinos at source points which are at a redshift z, which gives rise to neutrinos at Earth with an energy E ν . Here, mÑ 1 , τÑ 1 and ΩÑ 1 denote the mass, lifetime and present day energy density fraction respectively of the RH sneutrino and ρ c is the critical density of the universe. The Hubble parameter H(z) is given by H(z) = H 0 Ω Λ + Ω m (1 + z) 3 + Ω r (1 + z) 4 and Ω Λ , Ω m and Ω r correspond to the present day energy fractions of dark energy, matter and radiation respectively.
Galactic contribution: Neutrinos from the MilkyWay galactic halo will not experience any cosmological redshift. For the galactic halo, the flux of neutrinos at earth from sneutrino decay per unit energy and time in a volume element located at some point in the halo is given by [68], where s is the distance along the line-of-sight from Earth to the DM decay point and r is distance between galactic center and DM decay point which can be written in the form r = s 2 + r 2 − 2sr cos b cos l, where (s, b, l) represent standard galactic coordinates (distance along line-of-sight, latitude and longitude respectively) and r = 8.5 kpc is the distance between the earth and the center of the MilkyWay. In the formula above, the factor of ΩÑ 1 /Ω total DM accounts for the fact that only a fraction of the DM halo density ρ DM is in the form of the RH sneutrino. We assume the density distribution of DM in the  Figure 4: The dashed blue and red lines show the predicted differential neutrino flux at Earth from RH sneutrino decay accounting for galactic and extragalactic contributions, respectively. The dashed black line indicates the total differential neutrino flux which is a sum of these two. We have used mÑ 1 = 6.5 PeV, m N = 8 PeV, m 3/2 = 1.5 PeV as our benchmark point, which results in τÑ 1 ∼ 10 24 sec.
galactic halo follows the standard Navarro-Frenk-White (NFW) profile [69], which is given by ρ DM (r) = ρ 0 (r/r c ) with r c = 20 kpc and ρ 0 = 0.33 GeV/cm 3 . Notice that this value of ρ 0 accounts for the total DM present in the galactic halo.
Thus, the total neutrino flux predicted from astrophysics and RH sneutrino decay (including both galactic and extragalactic contributions) is given by, For the parameters ΩÑ 1 = 1.5 × 10 −6 , mÑ 1 = 6.5 PeV and τÑ 1 ∼ 10 24 sec, and taking values of all the cosmological parameters from the latest results given by Planck [70], we can compute the final differential flux obtained from sneutrino decay as a sum of galactic and extragalactic contributions. The results are shown in Fig. 4.
In Fig. 5, we compare the inferred neutrino flux from six years of IceCube data with our predictions for the total neutrino flux accounting for both astrophysical UPL (given in Eq. 5.1) and neutrinos from RH sneutrino decay. From the figure, it is apparent that the  Figure 5: This figure shows the total predicted terrestrial neutrino flux coming from both sneutrino decays as well as the nominal UPL astrophysical background. The black data points represent the observed neutrino flux which is inferred from IceCube data [15]. We can see that for our benchmark values mÑ 1 = 6.5 PeV, m N = 8 PeV and m 3/2 = 1.5 PeV and τÑ 1 ∼ 10 24 sec, the predicted total neutrino flux closely matches that seen by IceCube, and in particular the PeV events are well accounted for by the sneutrino decays.
sub-PeV flux is best fitted by astrophysical background while the shape of the high energy PeV flux points fit well with the neutrino spectra coming from PeV scale RH sneutrino decay.
Here, we see the second order of magnitude miracle. The lifetime of the RH sneutrino for PeV scale masses of the superpartners is in the range of 10 24 − 10 26 seconds (see Eq. 3.7). This combined with the suppressed relic abundance of the RH sneutrino from freeze-in processes (Eq. 4.26), automatically gives us a flux of PeV neutrinos consistent with the flux inferred from IceCube data.
Additional constraints: Gamma rays signals provide a complementary test of any interpretation of the IceCube data. Previous work has studied the compatibility of gamma  Figure 6: The predicted terrestrial gamma-ray flux from the RH sneutrino decayÑ 1 → ψ 3/2 ν L h 0 , where the decays of the Higgs give rise to gamma-rays. We have not accounted for secondary gamma-rays here. We have used the benchmark value mÑ 1 = 6.5 PeV and m 3/2 = 1.5 PeV. The current constraints on the diffuse gamma-ray flux from Fermi-LAT, H.E.S.S., CASA-MIA and KASCADE data are shown, and our predicted flux is well below limits set by these experiments. We have also shown the expected diffuse flux sensitivity of CTA with 500 hours of observation, and we can see that our predicted flux is unfortunately below CTA's sensitivity.
ray observations with astrophysical interpretations [71] as well as decaying DM interpretations of the neutrino flux [72][73][74]. In our model the Higgs produced from RH sneutrino decay can further decay to give prompt as well as secondary high energy gamma-rays. To check the compatibility of expected gamma-ray signals with the observed gamma-ray flux measurements, we calculate the prompt gamma-ray flux spectrum from the decay of sneutrinos by using both galactic and extragalactic contributions. Then we check for consistency of these predictions with gamma-ray bounds from the Fermi-LAT measurement of the isotropic diffuse gamma-ray background [75], CASA-MIA [76,77], the KASCADE experiment [78,79] and H.E.S.S. measurements [80]. The comparison between the predicted prompt gamma-ray flux and the constraints from these experiments is shown in Fig. 6. We can see from the figure that the predicted gamma-ray flux is several orders of magnitude smaller than the observational constraints. We also show in Fig. 6 the predicted diffuse gamma ray flux sensitivity of the Cerenkov Telescope Array (CTA) to a 10 • × 10 • box around the galactic center with 500 hours of observation. We have made a naive rescaling of the CTA point source sensitivity as given in [81,82] by assuming background domination in the box, and calibration of the background level by making observations just outside the box. Unfortunately, as we can see from the figure, our predicted gamma-ray flux is below the sensitivity of CTA as well. We note that ref. [74] showed that the Fermi-LAT bounds for a decaying DM interpretation of the IceCube neutrino flux can be improved by up to one-and-a-half orders of magnitude by using spatial information of the gamma ray flux as well as by using more detailed modelling of the astrophysical sources of gamma rays. However, our low energy gamma ray flux for the benchmark model parameters would still be below this improved sensitivity reach of Fermi.

Summary, conclusions and future prospects
The search for supersymmetry at the LHC has so far only given us null results. The absence of flavor changing neutral current signals in other experiments have long hinted that the scale of generic SUSY breaking must lie beyond the 100 TeV scale. Moreover the intriguing observation of unexplained PeV scale neutrino events seen at IceCube may be our first direct hint that the SUSY breaking scale is PeV. While this would leave us with a little-hierarchy problem for the Higgs mass, many other attractive features of SUSY such as improved gauge coupling unification would remain.
In this work, we tried to explore the implications of a generic supergravity model with SUSY breaking at the PeV scale. We have built a model with a U (1) R symmetry with specific R-charge assignments, where all the superpartners are at the PeV scale. We have assumed that the gravitino is the LSP and is therefore a good DM candidate. In addition we have assumed that the lightest RH sneutrinoÑ 1 is the NLSP and in our model, it is also quasi-stable and therefore a candidate for an auxiliary (albeit sub-dominant) DM component. We showed that with our specific R-charge assignments, the RH sneutrino has naturally suppressed couplings and the only interactions that it has are Planck-scale suppressed interactions in the Kähler potential involving the SUSY breaking fields.
We demonstrated that once SUSY, and hence R-symmetry is broken, the RH sneutrino is left with suppressed effective interactions to the other particles through trilinear scalar interactions, small mixing with LH sneutrinos and a tiny Yukawa coupling (y ν ∼ PeV/M p ∼ 10 −13 , see table 1). Of all of these interactions, the Yukawa coupling is the least suppressed and is responsible for both the production of RH sneutrinos in the early universe as well as their eventual decay through the channelÑ 1 → ψ 3/2 ν L h 0 , with a very long lifetime much larger than the age of the universe, O(10 24 ) seconds.
We then calculated the relic abundance of the gravitino and the RH sneutrino in our model via freeze-in processes. This is where we encountered our first surprise. Assuming a reheating tempreature not much larger than the PeV scale, the gravitino abundance au-tomatically had the right order-of-magnitude to make up the bulk of the dark matter of the universe. The RH sneutrino was found to have a suppressed abundance in comparison, making up only 1 part in 10 6 of the present day density of the universe. This led to the assertion that it makes up a sub-dominant part of the DM density.
The second surprise we encountered was when computing the expected flux of neutrinos observed at Earth from the decay of the RH sneutrino. We found that even though the RH sneutrino is expected to form a sub-dominant component of the DM relic density, its lifetime is naturally of the right order that it predicts a PeV scale neutrino flux at IceCube consistent with the level that has been observed. In the literature it has been pointed out that if decaying DM has to explain the IceCube neutrino events, one typically requires the DM to have a lifetime of O(10 29 ) − O(10 31 ) seconds to be consistent with the resulting flux inferred from the IceCube data. This is based on the assumption that the decaying DM accounts for the entire DM relic abundance (i.e. Ω DM h 2 ∼ 0.1) of the universe while in our case, the lifetime and relic density fraction of the sneutrino turn out to be O(10 24 ) seconds and O(10 −6 ) respectively. Thus, the suppression in the magnitude of neutrino flux due to suppressed relic abundance is balanced by the larger decay width of the sneutrino, which gives us the correct magnitude of the neutrino flux required to match IceCube observations.
We also showed that in our model, one naturally solves the µ problem, via the Giudice-Masiero mechanism, and we can also obtain the observed Higgs mass of 125 GeV at the cost of a relative tuning between the weak and PeV scales. Another issue that arose, was whether the observed neutrino masses could arise via the see-saw mechanism, given the suppressed Yukawa coupling of RH sneutrinos. We showed that non-renormalizable Kähler terms with a suppression by an intermediate scale of M * = 10 10 GeV, could explain the observed neutrino masses.
Thus in summary, with the simple assumptions of PeV scale superpartner masses and reheating temperatures not much larger than the PeV scale, we naturally get the right relic abundance for dark matter in the form of gravitinos, and we also naturally predict the right flux for the PeV neutrino events observed at IceCube from RH sneutrino decays.
Prospects for verifying this scenario: A precise measurement of the neutrino flux spectrum from an improved IceCube detector could potentially resolve the differences between simple astrophysical power law backgrounds and the RH sneutrino decay signal that we have postulated in this work. In addition, if the spatial distribution of the high energy neutrino flux is consistent with a DM origin (such as from clustering of events pointing towards the galactic center or dwarf galaxies), this could bolster the evidence for the RH sneutrino decay as the origin of the PeV neutrino flux seen by IceCube.
Since the three-body decay of the sneutrino gives a Higgs boson in the final state, it is possible that DM will also contribute to the prompt gamma ray flux from the decay of the Higgs into gamma-rays. Therefore we have also computed the possible gamma-ray flux expected from the sneutrino decay. Our results indicate that the prompt gamma-ray flux coming from sneutrino decay is well below the bounds set by high energy gamma-ray observatories. We have performed a naive analysis on the prospect of CTA being able to detect this flux, however it seems unlikely that CTA will be sensitive to the predicted flux of 100 TeV gamma-rays.
PeV scale superpartners would be unobservable at the LHC or any foreseeable high energy collider, but future flavor physics measurements could also be sensitive to generic squark/slepton mixings that are expected in gravity mediated SUSY breaking scenarios. We do not expect any dark matter direct detection signals, since most of the dark matter is in the form of gravitinos in our scenario and the sub-dominant RH sneutrino has extremely tiny couplings to SM particles.
Finally, we note that due to the suppressed value of the Yukawa couplings, the righthanded neutrino in our model cannot be responsible for vanilla leptogenesis [83]. Therefore one must consider other origins of the matter-antimatter asymmetry in this scenario such as soft leptogenesis [84].

Acknowledgments
MD would like to acknowledge Stephan West, Gaurav Tomar and Ujjal Dey for useful clarifications. MD would also like to thank Gaurav Goswami for several useful discussions. VR acknowledges useful discussions with Varun Bhalerao, Ranjan Laha, Pratik Majumdar, Danny Marfatia, Tuhin Roy and Christoph Weniger. Both authors would like to thank the organizers and participants of the "Blueprints Beyond the Standard Model" workshop at TIFR for valuable comments and suggestions. VR would like to thank ICTP, Trieste where a part of this work was completed. VR is supported by a DST-SERB Early Career Research Award (ECR/2017/000040) and an IITB-IRCC seed grant.

A Supersymmetry breaking mechanism
In this appendix we briefly discuss a mechanism of supersymmetry and R-symmetry breaking. We will see that in addition to the usual F -term VEV, this mechanism gives rise to a non-zero VEV for the scalar component of the hidden-sector field. Another interesting feature of the SUSY breaking mechanism presented here is that it also leads to a vanishing cosmological constant.
Typically, supersymmetry is broken in the hidden sector and then communicated to the visible sector via Planck scale suppressed operators (sequestering). Assuming that the theory has a supersymmetric hidden SU (N c ) gauge theory with N f pairs of hidden quark superfields Q andQ in the fundamental representations N c andN c , respectively, there will be a superpotential interaction of the hidden-sector field with hidden quarks superfields given by [85], This is consistent with both global R-symmetry as well as axial symmetry. The global axial symmetry and R-symmetry have a hidden QCD anomaly which generates a nonperturbative superpotential for X given by [85]: where Λ s corresponds to the scale of strong coupling of the hidden QCD while w is a constant term added to the superpotential. This constant term is necessary in any realistic model to obtain a vanishing cosmological constant after SUSY breaking.