Decays of Long-Lived Relics and Their Signatures at IceCube

We consider long-lived relic particles as the source of the PeV-scale neutrinos detected at the IceCube observatory over the last six years. We derive the present day neutrino flux, including primary neutrinos from direct decays, secondary neutrinos from electroweak showering, and tertiary neutrinos from re-scatters off the relic neutrino background. We compare the high-energy neutrino flux prediction to the most recently available datasets and find qualitative differences to expected spectra from other astrophysical processes. We utilize electroweak corrections to constrain heavy decaying relic abundances, using measurements impacted by electromagnetic energy injection, such as light element abundances during Big Bang nucleosynthesis, cosmic microwave background anisotropies, and diffuse $\gamma$-ray spectra. We compare these abundances to those necessary to source the IceCube neutrinos and find two viable regions in parameter space, ultimately testable by future neutrino, $\gamma$-ray, and cosmic microwave background observatories.


Introduction
The IceCube detector, located in the Antarctic ice layer, is sensitive to neutrino energies ranging from 10 -10 10 GeV [1]. Over its six year run, IceCube has detected several neutrinos in the energy range 100 TeV -10 PeV [2][3][4][5]. The measured neutrino flux in this range is significantly larger than that expected from the atmospheric neutrino background [2][3][4][5]. This suggests an alternative source at a significance of 5.6σ [5]. Previously, no statistically significant correlation between the direction of origin of the detected neutrinos and any known high energy γ-ray sources existed, suggesting an isotropic extra-galactic source [6]. Recently however, multi-messenger astrophysics linked one 290 TeV neutrino to a flaring blazar [7]. More data is necessary to determine whether blazars can explain the highest energy events. Other possible astrophysical sources such as Supernova remnants (SNRs), star forming regions, Fermi bubbles, and active galactic nuclei (AGNs), have also been considered in the past [8][9][10][11][12][13][14][15][16]. Beyond the standard model physics (BSM) explanations have been investigated with regard to heavy decaying dark matter (see, for example, [17][18][19][20][21][22][23][24][25][26][27][28]). However, many models of decaying dark matter as a source of the IceCube neutrinos are highly constrained because they are predicted to produce γ-rays in excess of current measurements [23,28,29]. In this paper, we explore the experimental signature of a heavy relic directly decaying to neutrinos, sourcing an isotropic extra-galactic high-energy neutrino flux. We focus on lifetimes that are shorter than the age of the universe. We examine whether this high-energy neutrino flux can fit the excess events seen between 250 TeV -10 PeV. We show that many constraints imposed by γ-ray observations can be avoided under this set of assumptions.
Recently, electroweak corrections at energy scales well above the electroweak (EW) scale have drawn considerable attention [30][31][32][33]. For high-energy scattering and decays, the EW effects significantly impact phenomenology by producing higher multiplicity final states. Different implementation strategies have been explored with regards to heavy decaying DM [30,34,35]. In our analysis we implement a fixed order EW shower. We use the results of the shower to predict a neutrino spectrum and fit it to that detected at IceCube. We also explore how the decaying relic model is constrained by its impact on light element abundances, CMB anisotropies, and diffuse γ-ray spectra, after including the EW shower effects.
A long lived relic has been considered previously as a source for the IceCube neutrinos, and analyzed up to redshifts of z = 1000 [36,37]. We extend this range by including neutrinos arising from re-scatterings off the relic neutrino background in our analysis. Our inclusion of EW corrections further changes the qualitative features of the neutrino flux today, leading us to conclude that EW corrections are a necessary part of an accurate forecast.

Models
In this paper, we consider two models in which our relic, X, directly decays to neutrinos. In our analysis, the PeV-scale neutrinos observed at IceCube are assumed to come from these direct decays. Naively, one may wish to consider a toy-model decay: X → νν [38]. However, implementing an EW shower highlights the inconsistency of this treatment. At ultra-high energies, the final state radiation includes many soft W 's, which turn charged leptons into neutral ones and vice versa. This leads to the production of roughly the same amount of neutral and charged leptons for center-of-mass (COM) energies far beyond the EW scale. This is a side effect of unbroken isospin in the high-energy limit. Model-independently, this implies that any high-energy neutrino spectrum sourced directly from a heavy relic decay will be accompanied by a spectrum of electromagnetically interacting particles, which will carry roughly the same amount of energy as the neutrino spectrum.
At energy scales much above the EW scale, Sudakov logarithms contribute to highermultiplicity final states. These corrections grow logarithmically as the mass increases. Effectively this leads to the production of EW jets. To quantify the neutrino spectrum arising from these jets, we implement a fixed order EW shower. The qualitative features of the EW jets are model independent, as any heavy particle that decays to neutrinos will also radiate gauge bosons. To zeroth order, this effect takes a delta function centered around M X 2 , and smears it towards lower energies. The energy lost by the neutrinos is carried away by gauge bosons, which themselves can decay into neutrinos, and contribute to the neutrino spectrum at lower energies. We describe the implementation of the EW shower in detail in Appendix A.
We consider two benchmark models that produce neutrinos through direct decays while remaining consistent with the isospin structure dictated by the Standard Model. We do not study a specific production mechanism for the heavy relic abundance, and assume it is cold. We note that inflationary dynamics can trivially produce such a particle during the reheating period [39]. Model-dependent constraints on these production mechanisms exist based on measurements such as isocurvature; however, these are not stringent enough to rule out the small abundance of decaying relics necessary to source the IceCube neutrinos [40,41].

Model I: Heavy Scalar X 1
We consider a heavy scalar X 1 , that couples to the standard model lepton doublets L i . Here i = 1, 2, 3 indexes the generation. For simplicity, we assume flavor universality: The zeroth order decays are given by: The ratio of branching ratios is essentially 1 : 1 at tree level. We will refer to the above decay model I as X → νν.

Model II: Heavy Fermion X 2
In our second model we consider a heavy Dirac fermion, that couples to the standard model lepton doublets (L i ) and Higgs doublet (φ).
We assume relic and its anti-particle have the same number density. The zeroth-order decays of X 2 are given by: Again, the decays to W ± , Z, h have equal branching ratios at tree-level in the high mass limit. We will refer to the above decay model II as X → V .

Derivation of the Present-Day Neutrino Flux
To derive the shape of the differential flux today we extend the analysis performed in [36]. We consider a number density of cold heavy relic Xs that decay with a given lifetime τ X : where n X,0 (t) is the number density in the limit τ X → ∞. For any given decay, high-energy neutrinos are injected into the thermal bath. The maximum possible energy is set by the mass of the heavy relic: The fractional energy distribution f Emax (x) of these neutrinos is determined by the EW shower, where x = E Emax and E is the injection energy of the neutrino. This decay gives rise to the following source term: Depending on when they were produced, the neutrinos may free-stream or scatter off the relic neutrino background. The cross sections for all relevant (anti-)neutrino-(anti-)neutrino scattering processes are listed in [36]. The total scattering rate is determined by the thermally averaged cross section: In the massless neutrino limit, the relative velocity simplifies to: v rel = s 2Ek , where s is the squared COM energy and k is the energy of the relic background neutrino. The scattering rate can then be written as [36]: where the second line is achieved via integration by parts. Neutrinos that scatter off the relic neutrino background, at the COM energies we consider, may produce two energetic neutrinos, two charged leptons, or two quarks. We define Γ ν and σ ν as the scattering rate and cross-section for 2 → 2 neutrino scattering. We account for this re-injection of neutrinos by adding an additional source term. This is sometimes referred to as a tertiary source term [42].
where Φ(t, E ) is the differential neutrino flux defined in terms of the neutrino number density n ν (t) = dE Φ(t, E ), and E is the scattered neutrino energy. We simplify equation (3.6) by rewriting the differential cross section in terms of the injection energy, E , and the fractional scattered energy y = E/E .
We can make the approximation in equation (3.7) because, for large boosts (γ > 100), g(y) becomes independent of E . We derive g(y) by boosting the relevant differential cross sections from the COM-frame to the laboratory frame: Defining separate functions g(y) for each scattering independently -for neutrino-neutrino scattering (and its conjugate scattering), g νν (y) and for the θ-dependent anti-neutrinoneutrino scattering (and its conjugate), g νν (y) -we write: where the ratios of scattering rates of νν and νν are 0.6 and 0.4, respectively. We now can rewrite equation (3.6): We can now set up the Boltzman equation which describes the thermal evolution of the differential neutrino flux: This partial differential equation can be solved numerically to obtain the present-day differential flux. In our analysis, we implement a propagation code to track the cosmological evolution of individual neutrinos, which is equivalent to solving equation (3.12) in small time steps. For a given lifetime τ X we generate events over the appropriate distributions of redshifts, a decaying exponential. The energy spectrum of the injected neutrino is determined by the EW shower. Based on the decay redshift, z, we divide the total traveling time of the neutrino into intervals such that the average number of scatterings within the interval is much smaller than one. If a scattering event occurs within a time step, the probability of re-injecting two neutrinos with energy g(y) and g(1−y) is weighted by Γν Γtot . If two neutrinos are re-injected, they undergo the same treatment as the primary injection, starting at redshift z , where the scattering has occurred. This process iteratively continues until the neutrinos either arrive today or scatter into charged leptons or quarks.
The output of the simulation is a histogram Φ h (t 0 , E), shown in Figure 1, which is related to the differential neutrino flux described in equation (3.12) by dividing by X's number density: where Φ(t 0 , E) is the solution to (3.12) at t = t 0 and thus accounts for tertiary neutrinos and EW effects. For short lifetimes, including the tertiary neutrinos significantly enhances for a heavy decaying relic X for two different lifetimes and two different decay models. The mass of X is set to M X = 2.4 (1 + z τ X ) PeV, the best fit mass of the neutrino spectrum measured at IceCube in the energy range 0.25 − 4 PeV [4]. z τ X is the redshift z at the decay lifetime τ X .
the flux of lower energy neutrinos, whereas for long lifetimes, these have negligible impact, since almost no scattering occurs. In the limit of negligible tertiary neutrinos, equation (3.12) can be solved analytically [36].
In our analysis, we only account for neutrino fluxes emerging from extragalactic relic decays. Extragalactic decays are the only relevant neutrino source for relics with lifetimes τ X ≤ 8 * 10 16 s, while galactic decays become important when considering longer lived relics [23,36,43]. We leave a detailed investigation of that region of parameter space to future work.

Estimating X's Number Density
We use Φ h (t 0 , E) to estimate the number density n X,0 (t 0 ) needed to roughly produce the excess number of events seen in the high energy bins at IceCube [4]. The number of predicted events in this range at the IceCube detector is obtained by integrating over the differential flux times the effective area A eff (E), which is provided by the IceCube collaboration [4], and multiplying by the detection time T (2078 days), and solid angle 4π, as well as the flux velocity v = c to restore SI units.
Based of the total number of events (N t = 5) in the range 0.25 − 4 PeV in [4] we estimate the number density n X,0 (t 0 ) that is needed to produce the observed number of events:

Constraints
In the following sections we consider different observables that can be used to constrain heavy decaying relics, and how these constraints affect the relic models best suited to generate the PeV neutrinos observed at IceCube. The summary of our findings appear in Figure 2. The shortest lived relics, those with τ X ≤ 10 12 s, are most strongly constrained by their impact on the abundance of light elements generated during big bang nucleosynthesis (BBN). Relics with intermediate lifetimes, 10 12 s < τ X ≤ 5 * 10 15 s, are most strongly constrained by their impact on the CMB anisotropy power spectrum. Relics with slightly longer lifetimes, 5 * 10 15 s < τ X ≤ 8 * 10 16 s, are most strongly constrained by the γ-ray spectrum they generate. These constraints all depend on the amount of energy injected into the thermal bath in the form of electromagnetically interacting (EM) particles. In order to explore constraints on our relic models we define Ξ, the EM energy density produced by relic decays divided by the energy density of cold dark matter ρ CDM : Here f int is the fraction of the relic energy density that becomes EM energy and should in principle be redshift-dependent due to rescattering. However, for the parameter range we are considering, the dominant source of EM energy is from the decay shower where this fraction is largely M X -independent. We take M X = 2.4 (1 + z τ X ) PeV, which gives the best fit mass for the two particular lifetimes shown in Figure 3, where z τ X is the redshift z at the decay lifetime τ X . We use this mass as a benchmark for evaluating the constraints for all lifetimes shown in Figure 2.
Based on the results of the EW shower we estimate a conservative lower bound of f int = 0.25 for both decay models. This estimate assumes that about one third of all hadronic energy is electromagnetically interacting, as well as one third of the energy coming from muon and tau decays. This is the number we use for all constraints below.

Light Element Abundances
Helium-3 (He 3 ) and Deuterium (D) are produced during Big Bang Nucleosynthesis (BBN) and their measured abundances are in general agreement with the predictions of BBN (see review in [44]). Decays of heavy relics can initiate EM cascades that interact with the light elements and alter their abundances. Injected EM particles with energies above 27 MeV can participate in all of the photodisintegration processes pertinent to producing excess He 3 and D by destroying larger nuclei, primarily Helium-4 (He 4 ), as well as those that break He 3 and D down into protons [45][46][47]. Constraints arise from numerically following the evolution of the abundances of all light elements involved in the creation or destruction of He 3 and D, and comparing the end predicted abundances to the measured He 3 and D abundances [45,46,[48][49][50]. This process, and the resultant constraints on a decaying particle injecting EM energy into the thermal bath, have already been worked out in detail by [45,46,[48][49][50]. We utilize those constraints on the allowed energy density and lifetime of a heavy decaying particle [50,51].

CMB Anisotropies
EM energy injection by heavy decaying relics with lifetimes in the range 10 12 s τ X 5 * 10 15 s can increase the free electron fraction around recombination, thereby distorting the CMB anisotropy power spectrum. Detailed constraints have been worked out in [51,52] and we rely heavily on their results, which utilize Monte Carlo Markov chains to calculate the effect of EM energy injection on the CMB anisotropy power spectrum. This study [51] rules out relics that inject enough EM energy at specific redshifts to produce power spectra inconsistent with current measurements.
The injection of EM energy increases the free electron fraction via ionization and collisional excitation. For relics with lifetimes of 10 14 s τ X 10 18 s, the decays enhance the optical depth of the universe after recombination, leading to an additional suppression of the CMB temperature angular power spectra (TT) and polarization power spectra (EE) at small angular scales [51]. Additionally, the increase in the free electron fraction at times between recombination and reionization increases the probability that photons scatter before reionization. This leads to extra polarization, which creates a bump in the EE spectrum at smaller angular scales than the usual reionization bump [51].
Relics with lifetimes of ∼ 10 13 s are the most strongly constrained by the CMB anisotropy. The EM particles released by relics with lifetimes 10 13 s delay recombination. This widens the last scattering surface, damping the temperature power spectra at small angular scales. Like the longer lived particles mentioned above, particles with lifetimes τ X 10 13 s also generate a bump in the EE spectra to smaller angular scales than the usual reionization bump, though the effect is weaker than that generated by longer lived particles [51]. At times much earlier than 10 13 s , the universe is fully ionized, and injection of electromagnetic particles, which increase the ionization fraction, have little impact. As a result, distortions to the CMB anisotropy spectum are exponentially suppressed for relics with τ X much less than 10 13 s. At lifetimes ∼ 10 12 s, only a fraction of the relics decay late enough to alter the CMB and the constraints from CMB anisotropies become weaker than those that arising from BBN.
The analysis done by [51] only considers the effects of particles with kinetic energies in the range [10 keV, 1 TeV], well below the energies of EM particles relevant to our models. We argue that the bounds also apply to injected EM particles with E ≥ 1 PeV because, around recombination, EM particles at these energies scatter off the CMB quickly enough to redistribute their energy to many particles with energies below 1 TeV, well within one Hubble time -energetic photons scatter off CMB photons via pair production extremely efficiently at these energies and redshifts. Electrons and positrons scatter off of the CMB through inverse Compton scattering, which while less efficient than pair production at these energies, is still much faster than the Hubble expansion rate for electrons of all energies considered in this paper, as can be verified.
Different injection energies in the range between [10 keV, 1 TeV] have different efficiency factors determined by their interactions with the thermal bath, which sets the width of the constraints in [51]. To know exactly where within this band our injection energies lie one would have to do a dedicated study. Here, we conservatively apply the least stringent bounds, which correspond to the lowest efficiency of dumping the electromagnetic energy into the thermal bath, noting that a dedicated study for our particular injection energies may improve these bounds by up to a factor of five.  Figure 2: Constraints on a wide array of different lifetimes for a heavy decaying relic X, releasing EM energy into the thermal bath. All constraints are at 95% confidence level. The light red shaded area is excluded by measurements of light element abundances and their agreement with BBN predictions. The blue shaded area is excluded by bounds from CMB anisotropies. The gray shaded region is excluded by diffuse γ-ray observations. All of these constraints are for injections of EM energy above some threshold value unique to the constraint and described in their respective sections of this paper. The cyan line is the forecast from the proposed PIXIE experiment, which could place more stringent bounds from y-distortion [53]. The black and brown lines indicate the abundance necessary to produce the excess IceCube neutrinos for models I and II, based off equation (3.15) and (4.1), assuming M X = 2.4 (1 + z τ X ) PeV. The black and brown (red) markers indicate the data points corresponding to the IceCube spectrum shown in Figure 3 (Figure 4). The dotted lines indicate M X = 2.4 (1 + z τ X ) PeV transitioning to an approximation rather than a best fit, as rescattering effects can change the electromagnetic fraction by O(1).

γ-Ray Constraints
When the heavy relic decays, the EW shower and decays of the showering products produce energetic photons, electrons, and positrons. These are reprocessed, producing a lower energy γ-ray distribution, primarily by inverse Compton scattering and pair production [54]. The γ-rays in this reprocessed spectrum lie in the energy range visible to the Fermi telescope, between 0.1 GeV and 820 GeV [55]. We derive constraints by requiring that the reprocessed spectra of heavy relic decays produce a γ-ray flux that is, in any bin, no more than 2σ above the flux presented in the Fermi Pass 7 Isotropic Extragalactic Gamma Ray (IGRB) spectrum [55].
In order to derive the reprocessed γ-ray spectrum resulting from a heavy relic decay, we follow [54,56]. Processes by which γ-rays can lose energy include photoioization, Compton scattering, photon matter pair production, and scattering off of the CMB. In this analysis, we approximate the γ-ray spectra as if EM particles are only reprocessed by the dominant scattering mechanisms for a particular redshift and energy. We also assume that a photon does not scatter if it has an optical depth d τ < 1. In this context, the optical depth can roughly be thought of the average number of times a photon scatters as it travels toward the Earth.
For redshifts 0 < z ≤ 700, EM particles are reprocessed by initiating cascades with CMB photons through pair production and photon-photon scattering [54]. Pair production is generally more efficient at reprocessing EM particles, except in a small range of energies for 300 ≤ z ≤ 700, in which photon-photon scattering is more efficient. Photon-photon scattering has a negligible effect on the constraints of relics with τ X ≥ 5 * 10 15 s, so we only consider the effect of pair production cascades in this analysis. In pair production cascades, photons pair produce electrons and positrons with CMB photons. The resulting electrons and positrons then upscatter CMB photons by inverse Compton scattering. These two processes continue until the COM energy falls below the pair production threshold. EM particles with energies above the threshold [54,56]: have an optical depth d τ > 1. Particles with energies below E th have optical depths d τ < 1, in which case we assume they free-stream toward the earth. At z > 700, additional scattering processes become important and all EM particles relevant to this analysis thermalize and do not produce any γ-rays observable today [54].
Since the particle cascades occur quickly compared to the expansion rate of the universe [56], we define a universal 'instantaneously' generated differential γ-ray cascade spectrum per unit injection energy 1 L(E, z), such that N γ = E inj L(E, z)dE, where N γ is the number of γ-rays produced when E inj of energy is injected into the thermal bath at redshift z. L(E, z) is built into the source term in the Boltzmann equation describing the evolution of the differential γ-ray flux: where S γ is the source term and H is the Hubble parameter. For a heavy relic, whose decays initiate EM cascades, the source term is: Here we use M X f int to denote the total EM energy injected per relic decay 2 . Solving the above Boltzmann equation (4.3), gives the diffuse γ-ray flux for any given z.
For observational purposes, we are interested in the flux at z = 0. L(E, z) depends on the dominant scattering process for a given redshift. The cascade spectrum for pair production was numerically calculated by [57]. Here, we use an approximate result only taking into account pair production (as the effects where photon-photon scattering is dominant are negligible): We derive constraints by comparing γ-ray spectrum that results from relic decay to the Fermi IGRB spectrum [55] requiring the predicted relic contribution produces less than a 2σ contribution in any one bin. Our results are shown in Figure 2.

Other Constraints
Spectral distortions to the CMB are often used to constrain the release of EM energy in the early universe [58,59]. These constraints can be derived by requiring that the decaying relic not produce µ-and y-distortions larger than the detection limit of COBE-FIRAS [60]. These are weaker than the constraints that arise from the light element abundances for the same redshifts, and thus not relevant for this analysis. However, as shown in figure  2, the proposed Primordial Inflation Explorer (PIXIE) [53], with projected sensitivities to µ-and y-distortions ∼ 1000x better than those of COBE-FIRAS, could detect y-distortions generated by almost all of the heavy relic models considered in the shorter lifetime parameter space window. Other works consider constraints on BSM physics from the 21 cm spin temperature signal [51,61,62]. A heavy decaying relic would heat the intergalactic medium, resulting in a positive change to the differential brightness temperature. We do not consider these constraints in detail in this paper because rough estimates in [51] indicate that they are not currently powerful enough to be relevant. However, more data and improvements in the uncertainty of the differential brightness temperature measurement could eventually provide stronger constraints [61,63]. Figure 2 shows that there are two windows in which a heavy decaying relic could be the source of the PeV neutrinos observed at IceCube, one with longer lifetimes from 5 * 10 14 s to 8 * 10 16 s, and one with shorter lifetimes between 7 * 10 10 s and 10 12 s . Here, we show the full neutrino spectrum predicted by the decay of a heavy relic, including neutrinos that result from EW-showering and re-scattering off of the relic neutrino background, for two sample lifetimes within these two allowed ranges. We compare these spectra to six years of IceCube data and we consider data from two different datasets. The first dataset (DS1) includes neutrinos of all flavors that deposited their energy within the detector [4]. The second dataset (DS2) considers six years of IceCube data on upward going muon neutrinos, where the interaction vertex was also allowed to be outside of the detector, significantly enhancing the effective area [5]. Both datasets are complementary, and predict roughly the same neutrino fluxes for energies above 3 * 10 5 GeV [5]. The main focus on our analysis has been on DS1. We still include DS2 in our analysis because it contains the highest energy neutrino event measured to date. The event, which deposited 4.5 PeV in the detector, has a 88% probability of being caused by a muon-neutrino, in which case IceCube predicts a reconstructed energy of 7.5 PeV [5]. All other anomalous high-energy neutrino events in both data sets have energies below 2.5 PeV. We show the best fit for two allowed sample lifetimes for DS1 [4]. We also show the best fit to all of the data by combining both datasets. We want to stress that the IceCube collaboration has not published a combined measurement, and thus our second set of fits should be taken as purely illustrative. For the following comparison, we choose to fit only to the highest energy events even though there also exists an excess in the lower energy range. We make that choice because a decaying relic cannot comfortably explain both of these excesses at the same time. One should note that systematic uncertainties and atmospheric backgrounds are much higher in the lower energy range than the higher energy range. Additionally, DS1 and DS2 are in tension for bins below 3 * 10 5 GeV [5]. The excess of events in the lower energy range is larger in DS1 than in DS2. A better understanding of the tension between the two datasets in this range may be able to give additional insight into the source of the lower energy excess.

Dataset 1
Here, we compare our forecast to DS1 [4]. Figure 3 shows the neutrino spectrum forecast with the best fit mass to DS1 for two different allowed lifetimes τ Xs = 5 * 10 11 s and τ X l = 5 * 10 15 s. We choose the mass such that the chi-squared is minimized within the range 2.5 * 10 5 GeV -4 * 10 6 GeV [4]. We can see that for both allowed lifetimes, the spectrum resulting from a heavy decaying relic can reproduce the four highest energy non-zero bins reasonably well. Qualitatively, the spectra do not differ much between the different lifetimes. The shorter lifetime τ Xs predicts slightly more events between 2.5 * 10 5 GeV -4 * 10 5 GeV, which is an indicator of tertiary neutrinos contributing to the lower tail of the spectrum. While overall there is some contribution to the lower energy bins between 6 * 10 4 GeV -2.5 * 10 5 GeV, which relieves some of the tension between the expected background and the measurement, it is still an order of magnitude too small to be in agreement with the data. This suggests that different sources or systematic backgrounds would be needed to explain the excess seen between 5 * 10 4 GeV -2.5 * 10 5 GeV.

Combined Datasets 1 and 2
To combine both datasets, we rearrange equation (3.14) to find the average flux per bin Φ a as predicted by the number of events per bin, N b , in DS1: E min and E min correspond to the lower and upper limit in each bin in DS1. We consider all bins between 2.5 * 10 5 GeV -10 7 GeV. We then calculate how many events per bin, N p , the average flux Φ a predicts in DS2: E min and E max correspond to the lower and upper limit in each bin in DS2. A ef f (E) is the effective detection area for DS2 provided in [5]. η f = 1 3 is the flavor efficiency factor, accounting for DS2 only being sensitive to muon-neutrinos.
Combining both datasets shifts the best mass fit from M X = 2.4 (1 + z τ X ) PeV to M X = 8. 1 + z τ Xs PeV and M X = 4.4 1 + z τ X l PeV for the short (τ Xs = 5 * 10 11 s) and long (τ X l = 5 * 10 15 s) lifetimes, respectively. Non-surprisingly, including the higher energy event shifts the mass fits towards higher masses. While for the longer lifetime the spectrum shape still shows the remains of a peak centered around E max , the spectrum for the shorter lifetime does not show this feature. This is due to the spectrum being dominated by tertiary neutrinos, which leads to a power-law shape with a hard cut-off. This effect is more pronounced in the combined dataset fit for τ Xs , because higher M X enhances the scattering rate off of relic neutrinos.

Conclusion
We utilize EW corrections to constrain heavy decaying relic abundances, using measurements impacted by EM energy injection, such as light element abundances during BBN, CMB anisotropies, and diffuse γ-ray spectra. Beyond the scope of our application, in the future EW corrections may be a useful tool to better constrain BSM physics beyond collider reach, using cosmological and astrophysical data.
We derive a precise forecast of neutrino spectra produced by direct decays from heavy relic particles with lifetimes smaller than τ X < 8 * 10 16 s. Due to our analysis including EW showers and tertiary neutrinos, our forecast accurately captures the shape of the possible spectra, and thus can be used as a powerful discriminant against other astrophysical explanations. This will prove useful as IceCube collects more data. IceCube has plans for a large expansion of its detecting abilities referred to as IceCube Gen 2 [64]. These include plans to increase IceCube's detection area by a factor of 10, which is expected to improve IceCube's detection sensitivity for neutrinos with energies in the range 100 TeV -1 EeV by a factor of 10 [64].
Further, we can expect future experiments to shed more insight into the decaying relics proposed here. PIXIE should be able to detect y-distortions from relics with τ X 5 * 10 11 s.
Any isotropic, long lived decaying relic heavy enough to generate 1 PeV neutrinos will also produce a unique γ-ray spectrum. While Figure 2 reveals that none of the decaying relics considered in this paper are ruled out by their γ-ray spectra, this may change as Fermi's detection resolution improves, and as more sensitive γ-ray telescopes come online. The latest Fermi data analysis, Pass 8, is far more sensitive to point sources than Pass 7, and additional analysis seem to indicate that much of the IGRB flux derived from Pass 7 may actually be unresolved point sources [65]. Considering the analysis done by [65], we conservatively estimate that at least half of the IGRB flux measured in Pass 7 is actually unresolved point sources. This would tighten the γ-ray constraints by at least a factor of 2. However, [65] contends that the entire IGRB measured in Pass 7 could in principle be explained by unresolved point sources, suggesting that the γ-ray constraints could become significantly tighter, depending on what fraction of the IGRB is eventually found to be unresolved point sources. These constraints will also improve as γ-ray telescopes with better point source resolution, such as the High Energy Cosmic Radiation Detection facility (HERD) and the Chernekov Telescope Array (CTA), come online in 2020 with expected 10x more sensitivity than current γ-ray detectors [66].
If most point-source contributions to the IGRB are identified, what remains might be a truly isotropic spectrum from a model such as those described here. Thus, more IceCube data, paired with improved γ-ray detection sensitivity, may provide a smoking gun for confirming a heavy decaying relic as source of the IceCube neutrinos and thus physics beyond the standard model.

A.1 Electroweak Splitting Functions
At energies much larger than the electroweak scale, electroweak radiative corrections have a large impact on decay and scattering processes. This has been explored in the literature with regards to a 100 TeV particle collider [31], and indirect dark matter detection spectra [30]. These radiative corrections can be approximated by factorizing the differential cross section (or decay width) into the original 2 → 2 (1 → 2) process, times the differential probability that one of the final states will emit an additional gauge boson (or split into two different particles altogether).
At very high energies there will be more splittings, which requires a summation for a full treatment. However, the majority of the higher order splittings are soft, which means they only carry a small fraction of the total energy. To compare our prediction to the spectrum at IceCube we are only interested in neutrinos within two orders of magnitude of the highest energy neutrinos. This allows us to use a cutoff above which the splitting probability does not exceed 1. In our EW shower we only consider 'hard' first order splittings, in which the gauge boson carries more than 10 −2 of the maximum possible energy. This treatment captures how the resulting spectrum today is affected by the additional particles produced by a decay. However, at these high COM energies, many soft W 's can populate the final state, which can turn a charged particle into a neutral one and vice versa. Therefore at high energies we keep track of all leptons and scalars (for the Higgs and the longitudinal components of the gauge bosons), and do an isospin average in the end.
We use the following splitting functions in the implementation of the EW shower. Here we follow the notation in [30]. D A→B (x) gives the differential probability that a particle A turns into another particle B, with fraction x of the initial energy. Equations (A.3)-(A.9) show the splitting functions for scalars such as the Higgs h, and the longitudinal components of the gauge bosons, W L and Z L . Equations (A.10)-(A.14) show the splitting functions for fermions. All couplings are renormalized. L(x) and l below are the universal kinematical functions [30]: Splitting Functions for h/Z L In the following equations h may be replaced with Z L . In the high-energy limit h and Z L are not distinguishable, which is why they have the same splitting functions.
Notice that the initial particle spin stays the same when emitting a gauge boson. Here for example we start out with a Higgs H, which can emit a W T , which turns the Higgs into a W L , or it can emit a Z L , in which case it remains a Higgs. In either case the mother-particle, which carries the majority of the energy after the splitting, remains a scalar. The Higgs can also split into two top quarks, in which case neither of them have a higher probability of carrying the majority of the energy. This can be seen by A.5 being independent of x.
(Splittings into other quarks and leptons are negligible because their yukawa couplings are small.) Splitting functions for W L

Splitting functions for fermions
Here are the splitting functions for a charged fermion: Here are the splitting functions for a neutral fermion: (A.14)

A.2 Description of Included Processes
In model I, each decay produces two leptons. In model II, each decay of a heavy X-particle produces one scalar and one lepton. We consider one hard splitting off of both daughter particles. We decay all top quarks, keeping track of all gauge bosons. We combine the energy spectrum of the primary lepton with subsequent decays from any gauge bosons (V L and V T ) to secondary leptons. We consider only direct leptonic decays, as neutrinos resulting from hadronic decays are much less likely to be energetic enough to be above our Figure 5: The final neutrino spectrum of decay model I and II considering EW showers. The spectrum includes decays to neutrinos from any gauge bosons, taus and muons produced in the EW shower. We can see that for higher COM energies the peak decreases, which demonstrates how more energy is distributed to EW radiation.
set threshold of x > 0.01. Included gauge bosons come from the primary scalar, radiation off of either leptons or scalars, and subsequent decays from top quarks to W 's.
The total lepton spectrum f tot (x) is the combination of the primary and secondary lepton spectrum. f tot (x) is the probability distribution of producing a lepton with fraction x of M X 2 . Since we have to average over charged and neutral leptons due to the possibility of soft W -emission, the probability distribution of a neutrino with energy fraction x of M X 2 is given by 1 2 f tot (x). The other half of f tot (x) results in charged leptons: electrons, muons, and taus. While electrons are stable, muons and taus decay further before interacting with the thermal bath.
Neutrinos from primary muon and tau decays will also contribute to the measured spectrum today. We assume an isotropic three-body decay, and decay all muons into three particles, two of which contribute to the neutrino spectrum. The tau-decays are more subtle as there is a greater variety of possible final states. We treat the leptonic tau decays in the same manner as the muon decays. We also include other tau-decays with up to three particles in the final state and add the resulting neutrinos to the spectrum, without further decaying any resulting mesons. The final neutrino spectrum, which is shown in Figure 5, is denoted by f Emax (x).
Neglecting hadronic decays may slightly underestimate the low energy tail of our distribution. In the future, it may be worth integrating an EW-shower formalism with a hadronic shower. For our purposes, the accuracy of the high energy tail of the neutrino distribution is most important, to which hadronic decays will not significantly contribute.