Visible Decay of Astrophysical Neutrinos at IceCube

Neutrino decay modifies neutrino propagation in a unique way; not only is there flavor changing as there is in neutrino oscillations, there is also energy transport from initial to final neutrinos. The most sensitive direct probe of neutrino decay is currently IceCube which can measure the energy and flavor of neutrinos traveling over extragalactic distances. For the first time we calculate the flavor transition probability for the cases of visible and invisible neutrino decay, including the effects of the expansion of the universe, and consider the implications for IceCube. As an example, we demonstrate how neutrino decay addresses a tension in the IceCube data. We also provide a publicly available code to calculate the effect of visible decay.

Neutrino decay modifies neutrino propagation in a unique way; not only is there flavor changing as there is in neutrino oscillations, there is also energy transport from initial to final neutrinos. The most sensitive direct probe of neutrino decay is currently IceCube which can measure the energy and flavor of neutrinos traveling over extragalactic distances. For the first time we calculate the flavor transition probability for the cases of visible and invisible neutrino decay, including the effects of the expansion of the universe, and consider the implications for IceCube. As an example, we demonstrate how neutrino decay addresses a tension in the IceCube data. We also provide a publicly available code to calculate the effect of visible decay.

I. INTRODUCTION
IceCube's discovery of a new flux of astrophysical neutrinos at TeV-PeV energies [1][2][3] of extragalactic origin [4,5] opens up the window to many new physics probes. Due to the extremely long distances involved in neutrino propagation, this flux of neutrinos provides a unique opportunity to probe new physics models, in particular neutrino decay. While neutrinos do decay in the Standard Model (SM) 1 , the decay is highly suppressed and unapproachable by any foreseen experiment [6,7]. In this work, we consider the well-studied Majoron model, which postulates the existence of a singlet Higgs-like scalar with non-zero lepton number [8][9][10][11]. Such a scalar could generate the Majorana mass of a right-handed sterile neutrino through a lepton number violating vacuum expectation value. The spontaneous breaking of the global lepton number symmetry would then provide us with a massless Goldstone boson, φ -the Majoron. Through mixing between active states and the sterile, we could obtain off-diagonal couplings between the neutrino mass eigenstates and φ. Generically, we are interested in Lagrangian terms of the form where the g ( ) ij are the scalar (pseudo-scalar) couplings. While Majorons are potential dark matter candidates [12][13][14][15][16][17][18][19][20][21][22], in those scenarios the Majorons are too heavy for neutrinos to decay.
Such a term in the Lagrangian would then induce the decays of the light neutrino mass eigenstates, ν i → ν j +φ. Both ν i → ν j and ν i →ν j decays can occur. In the case * asli.abdullahi@durham.ac.uk; 0000-0002-6122-4986 † pdenton@bnl.gov; 0000-0002-5209-872X 1 While the SM does not tell us how neutrino masses are generated, once they are neutrinos will decay. The details of neutrino mass generation do not affect this work.
of Dirac neutrinos, the ν i →ν j decays result in neutrinos in the mostly sterile direction (assuming a light sterile) which are undetectable (unless there is a secondarȳ ν j → ν k decay), while in the case of Majorana neutrinos such a channel is detectable leading to the possibility of differentiating Majorana neutrinos from Dirac neutrinos by the measurement of neutrino decay [23][24][25]. In this work we will focus on the Majorana case, although our results can be easily extended to the Dirac case as well.
Phenomenologically there are two kinds of decay, each of which leave distinct imprints on the measured flux. Invisible Decay In these decays, the final state neutrino is not observable in the detector, either because it is sterile or because its energy is too low to produce a signal through scattering. Invisible decay results in a depletion of the number of events below a fiducial energy given by the coupling and distance traveled. Typically this results in fewer ν µ 's and ν τ 's assuming the normal mass ordering. While we include results for invisible neutrino decay for comparison with other works in the literature, we are agnostic about the model details and focus mostly on the effect of visible decay.
Visible Decay Here, the final state neutrino is observable in the detector. The decay of neutrinos during propagation appears as a depletion of the heavier mass eigenstates as with invisible decay. Simultaneously, we get an increase in the number of lightest eigenstates, and hence more ν e 's in the normal mass ordering. Another signature of visible decay is in the energy of the observed neutrinos. As the neutrinos produced in decay have less energy than the source neutrinos, the increase in lighter mass eigenstates is in lower energy bins than the parent neutrinos.
Neutrino decay is constrained in a wide range of experiments. In general, the longer the baseline and the lower the neutrino energy, the stronger the constraint. The weakest constraints come from atmospheric, longbaseline accelerator, and reactor neutrinos which are mainly sensitive to ν 3 decays [23,[26][27][28][29][30][31][32][33][34][35][36][37][38][39][40]. Hints of ν 3 decay has been found in NOvA and T2K data [35,41]. The next strongest constraints come from solar neutri-  [31,56,61,66,75]. Each constraint only applies to certain mass/flavor eigenstates except for the CMB constraint which is essentially independent of which or how many states decay, and the IceCube constraint which only applies when all states decay with the same lifetime. The vertical blue line is a hint of invisible partial neutrino decay from IceCube data [62].
While the above constraints are directly related to the depletion of neutrinos, it is also possible to constrain the process of neutrino decay via the νν → φφ diagram in the early universe leading to very tight bounds [70][71][72][73][74][75]. Early universe constraints do not depend on neutrinos explicitly decaying, thus it is in principle possible to construct a model that evades the early universe constraints and still predicts neutrino decay, for one such realization see ref. [76]. All of these constraints are summarized in fig. 1.
In addition, neutrino decay can be further probed in the future in the case of a galactic supernova [24,[77][78][79], a measurement of the diffuse supernova neutrino background [80,81], improved solar neutrino measurements [82], or a measurement of the cosmic neutrino background [83]; see ref. [84] for a review of new physics searches in upcoming neutrino experiments.
Due to the large mixing in the neutrino sector, neutrino decay may not significantly deplete the flux of neutrinos from SN1987A [62,67] making IceCube the most sensitive direct probe of neutrino decay. In addition, since Ice-Cube is sensitive to high energy neutrinos well above the τ production threshold, IceCube is sensitive to six of the nine flavor changing channels P ee , P µe , P eµ , P µµ , P eτ , P µτ (it is unlikely that ν τ 's are produced in high energy astrophysical sources) unlike supernova neutrino probes which can only detect ν e and thus are sensitive to only three 2 The track to cascade spectral difference was also investigated in [64] which pointed out that a neutron decay source scenario somewhat relaxes the tension, although such astrophysical models are quite unlikely [65].
flavor changing channels: P ee , P µe , P τ e . The sensitivity to different flavors provides a handle to break degeneracies with astrophysical uncertainties. In this article we further explore the phenomenology of neutrino decay including the regeneration component relevant for visible neutrino decay.
This paper is set out as follows: in section II, we describe in more detail the different decay scenarios, focusing on visible decay, and set up the formalism required to study the impact of neutrino decay on the flavor ratios. We also consider the added effect of cosmology on the observed spectrum of the neutrinos. We then analytically evaluate the integrals in some cases and take some interesting limits in section III. In section IV, we present the results of implementing the formalism in section II. We study how the parameters in our model, namely our couplings with the Majoron, the absolute mass scale of the neutrinos, and the energy spectrum at the source, impacts visibility at Earth. We also suggest some benchmark points that IceCube, and future experiments, could probe. We discuss the possibility of mitigating the > 3 σ tension in the track and cascade events observed at IceCube in section V, and finish up with a discussion of our results and our conclusions in sections VI and VII. The code for this article is publicly available and can be found at github.com/PeterDenton/Astro-Nu-Decay [85] and the data files for all the figures can be found at peterdenton.github.io/Data/Visible_ Decay/index.html.

A. Invisible Decay
Invisible decay refers to the absence of decay products in the detector, either because they are sterile or too low energy to be detected (at lower energies the cross section is lower and, in the case of the high energy astrophysical flux, both the astrophysical flux and the atmospheric backgrounds are higher at lower energies). This means that if we start with a flux of ν α , due to the neutrinos decaying along the baseline, the flux of ν β is expected to be less than that produced in oscillations alone.
To begin, we define the total transition probability as the ratio of the observed spectrum at the Earth for a given channel, Φ E αβ , over the initial spectrum at the source, Φ S α . That is, the probability is the observed spectrum at the Earth divided by the spectrum at the source in the case of no oscillations or decay, Alternatively, the spectrum at the Earth can be simply computed as the probability calculated here multiplied by the initial spectrum. We take a single power law spectrum (SPL) in the remaining work, where Φ S 0 is the flux normalization. While a broken power law may also be considered [62,96,97], it was not particularly preferred by IceCube data to an SPL source. We note that the quantityP αβ defined in eq. 2 is not a true probability in that it can be larger than 1 depending on the spectrum. A clear definition of what we take to be the transition probability is necessary when the flux has non-trivial dependencies, as seen in sec. II C where the flux depends on the source redshift.
In the following, the energy E refers to the energy at the source and that detected at Earth, as there is no change in energy in the SM, or in invisible decay.
Starting with the SM transition probability for ν α → ν β , under the assumption of relativistic neutrinos with equal momentum, We can account for the effect of decay by making the substitution where Γ i is the decay width for mass eigenstate i in the lab frame. The partial width in the lab frame for mass eigenstate ν i to decay to ν j of either helicity is, where E i is the energy of the initial neutrino, ν i , x ij ≡ m i /m j , and we have explicitly included both decay channels: ν i → ν j and ν i →ν j , and, where the f and h equations represent the ν i → ν j case, for scalar and pseudo-scalar mediators respectively, and the k equation represents the ν i →ν j case, for both mediators. The full width 4 is then Γ i = j Γ ij . Unless otherwise specified, we will consider the case of g ij = g ij , the impact of which is shown in the middle left panel of fig. 3. As the baselines under consideration are very large, we can average over the oscillations to get, where the bar denotes the oscillation averaging since the states have decohered. The SM component is, and the depletion component is, We note that whileP dep αβ < 0,P inv αβ > 0. The smallest thatP dep αβ can be is −|U α2 | 2 |U β2 | 2 − |U α3 | 2 |U β3 | 2 , thus the smallest thatP inv αβ can be as a function of g ij , m 1 , and L is |U α1 | 2 |U β1 | 2 since ν 1 is stable. For the nonoscillation averaged depletion component, see appendix A 1.

Limiting Cases
As a benchmark, we will assume that the neutrinos detected by IceCube result from full pion decay wherein 4 The lifetime of neutrino ν i its rest frame is where we recall that Γ i is in the lab frame. the relative flux normalized to the ν e flux is (Φ νe : Φ νµ : Φ ντ ) = (1 : 2 : 0), with all three neutrinos from the π ± decay carrying nearly the same energy. In the SM this results 5 in a flavor ratio of (1 : 1.152 : 1.117) after propagation. It is interesting to consider some limits to demonstrate the effect of the decays on the flux observed at Earth. In the case where Γ 2 = 0 (that is, g 21 = 0) and Γ 3 → ∞, so that all ν 3 decay at the source, we havē P inv αβ → |U α1 | 2 |U β1 | 2 + |U α2 | 2 |U β2 | 2 and we end up with (1 : 0.504 : 0.610). This time assuming both ν 3 and ν 2 decay promptly, we getP inv αβ → |U α1 | 2 |U β1 | 2 , and a final flavor ratio of (1 : 0.140 : 0.342). We can clearly see the relative depletion in the flux of ν µ and ν τ due to the invisible decays.

B. Visible Decay
Here we follow and extend upon refs. [23,98]. The case of visible decay can be qualitatively understood as the depletion of mostly ν µ 's and ν τ 's, as in invisible decay, accompanied by regeneration consisting of mostly ν e 's with lower energy in the normal mass ordering. We now need to account for the regeneration of the neutrino flux with a term describing the appearance of decay products, This section is devoted to calculating this regeneration term which is considerably more complicated, even in the oscillation averaged case, than the SM or depletion terms. The complexity arises from the fact that the initial and final state neutrino energies, E i and E f , are now distinct and thus, the observed spectrum depends on the initial spectrum. Additionally, the possibility of consecutive decays must now be accounted for. Although the effect of multiple decays is smaller than that of single decays, it is not completely negligible. We first shift to the mass basis since this is the basis in which the decays happen, The transition probability involves an integral over both the decay location, L 1 , and the initial neutrino energy, 5 Here and throughout we will use the global fit numbers from nu-fit 4.1 [90].
where the mass eigenstate spectrum is related to the flavor eigenstate spectrum by Φ S i = α |U αi | 2 Φ S α . In this article we will be assuming that initial spectrum of each flavor is the same up to the overall flavor ratio normalization of (1 : 2 : 0). The integral limits, , are fixed by kinematics due to the allowed final state energies.
Assuming a neutrino ν i decays once at L 1 ≤ L, for a baseline L, our regeneration term is then constructed from the following five factors: 3. the decay of ν i and appearance of ν j , Γ ij W ij , 4. for unstable ν j , survival until Earth, e − 1 2 Γj (L−L1) , 5. and the phase accumulation of ν j until Earth, is the normalized energy distribution of the decay products, where The superscripts in eq. 18 include both ν i → ν j and ν i →ν j channels. To account for specifically helicity conserving (flipping) channels then eqs. 16-17 should be used as appropriate, although the full expression for the total width in eq. 6 should always be used. The regeneration amplitude is given by, where we have removed overall phases that only contribute to oscillations. The differential decay probability is the amplitude squared, We use the fact that neutrinos promptly (on scales relative to the total propagation) lose coherency, so we do not include any interference terms. Such terms, however, are important in the context of decays over local distances, and must be accounted for when considering the effects of decay on the measured event rates at, for example, long-baseline neutrino experiments. See appendix A 2 for details. Averaging over oscillations and assuming all decays happen incoherently, we analytically perform the integral over L 1 , (22) Up to this point we have only considered single decays. If both g ( ) 32 = 0 and g ( ) 21 = 0, the neutrinos will experience consecutive decays, ν 3 → ν 2 → ν 1 . This additional decay pathway modifies the ν 3 → ν 1 transition probability. The initial neutrino decays at a distance L 1 ≤ L, after which the intermediary neutrino, with energy E int , propagates and decays at L 2 , with L 1 ≤ L 2 ≤ L. We construct the regeneration amplitude as described above, The total probability is given by, We can again perform the L integrals analytically, followed by integration over the initial and intermediate energies. The ν 3 → ν 1 transition probability now reads

C. Cosmology
So far we have assumed that the neutrinos only propagate over local distances. In the context of high energy astrophysical neutrinos, it is known that their propagation distances are long enough that the expansion of the universe must be accounted for. Decay is dependent on the duration of travel, and so light-travel distance is the correct measure for this purpose. Baselines and neutrino energies are now functions of the redshift to the source, z, where The energy E(z) is the energy at production, and E 0 the observed energy at detection. The function L(z a , z b ) computes the distance between two redshift points z a and z b . We assume a flat ΛCDM cosmology, with Ω m = 0.315, H 0 = 67.4 km/s/Mpc [99] and With the introduction of cosmology, both the invisible and visible scenarios are more subtle and require some care. In the following section, we will continue to work in the mass basis, but all expressions can be re-expressed in the flavor basis using eq. 14. Since the neutrino energies are now a function of redshift, the spectrum at the source and at the point of decay are no longer the same, but are modified by a factor, In accordance with the definition given in eq. 2, the SM transition probability is also modified by the same factor To compute the depletion and regeneration terms, we make the substitutions: As we frequently encounter the product Γ i L in the arguments of our exponential functions, we make a redefinition of the distance integral, eq. 28, .
The depletion term is now given bȳ For the regeneration term, we start from eq. 20 and make the relevant substitutions. The modified version of eq. 15 now involves an integral over z 1 , the redshift at which the neutrino decays, where After making the correct substitutions and taking the squared amplitude, we obtain In the case of consecutive decays, as in sec. II B, we account for integration over the redshift at which our second decay occurs by adding a second term for ν 3 → ν 1 , where the first term is from eq. 36 and the second is, where we note that Finally, since not all neutrinos are coming from the same redshift, we must integrate the total (in)visible probability over the redshift distribution of the source population R(z), where we note that we must include depletion and regeneration components as appropriate (the SM part pulls out of the integral since it does not depend on redshift). For simplicity we take R(z) = δ(z − 1) as our benchmark redshift evolution as this roughly reproduces typical redshift evolution functions, see below.

III. ANALYTIC SOLUTIONS
While a closed form solution for the regeneration term including cosmology likely does not exist, it is possible to find useful expressions without cosmology. These provide intuition for what is going on and are relevant, for example, for neutrinos propagating over galactic distances such as from a galactic supernova wherein cosmology does not play a role. For simplicity we consider only one of the g ( ) ij = 0. We carry out the integrals discussed in II B including a power law spectrum and find, We validate this expression in appendix B. While eq. 42 is relatively intractable, it does allow for several interesting limits to be evaluated.

Limiting Cases
First, we confirm that as E f → ∞ (which takes us back to the SM) the regeneration probability goes to zero. To see this we use the fact that which when inserted into eq. 42, gives zero.
Next, as E f → 0 (the full decay scenario, Γ i → ∞) the regeneration probability tends to a non-zero constant. In this limit T → ∞ and thus Within this limit, we consider the heavy mass case, m 1 → ∞ (x → 1) and find Taking only g ( ) 31 = 0, without cosmology and as E f → 0 and m 1 → ∞, and similar for the other g ( ) ij . It is noteworthy that in this limit the probability is no longer dependent on γ. This makes sense as the integral over the spectrum is over a region [E f , x 2 E f ]. As x → 1, the spectrum is asymptotically constant and the final neutrino has the same energy as the initial neutrino.
Finally, taking the small mass case of the low energy limit, Unlike the previous case, we have dependence on the spectral index γ. As expected, a steeper spectrum at the source depresses the low energy regeneration. That is, if only g ( ) 31 = 0, without cosmology and as E f → 0 and m 1 → 0, We also note that in this limit the probability is independent of m 1 when γ = 1. Some closed form solutions exist for arbitrary E f and m 1 that do not depend on incomplete gamma functions. In our benchmark case of γ = 2, we can use the fact that Γ(2, x) = e −x (1 + x) to writē Additional closed form solutions exist for other integer values of γ; e.g. Γ(3, x) = e −x (x 2 + 2x + 2) and so on. While a full expression for the regeneration term with cosmology is extremely complicated, the various limits evaluate quite simply as for the case without cosmology. We find that, including cosmology, in the heavy mass case which if only g ( ) In the light mass case we have and, if again only g ( ) One factor of (1 + z) −γ comes from the definition of the probability in eq. 2 and the other from the contribution due to cosmology. This allows us to easily write down the flavor ratios in the full decay limit. The flavor ratios from pion decay in the full decay case with cosmology and for m 1 = 0 is shown in table I for various spectral indices in the case of only ν 3 → ν 1 decay. The flavor ratios in the SM and for invisible decay do not change subject to cosmology. We see that as γ → ∞ we recover the invisible decay case regardless of m 1 . This is because the regeneration term goes to zero relative to the SM and the depletion terms, as shown in eqs. 54 and 56.
To summarize, in the full decay limit, there are four main cases as discussed above: m 1 → 0 or m 1 → ∞ and with or without cosmology. Only in the case of m 1 → ∞ and without cosmology does the spectral index not affect the final flavor ratio.

IV. RESULTS
We now numerically integrate the expressions in section II. As there are numerous parameters in our model, we consider benchmark parameters shown in table II and then numerically show how deviations from the benchmark parameters affect the transition probability. We take all six couplings to be 10 −6 , including both scalar and pseudo-scalar couplings; this puts the decay features in IceCube's region of interest. We take the lightest neutrino to be massless, the initial spectral index to be 2, and the redshift evolution to be a delta function at z = 1. We All FIG. 2. The oscillation averaged invisible decay probability for benchmark parameters (see table II) including both scalar and pseudo-scalar interactions and decays to both neutrinos and antineutrinos. The orange, green, and red curves are all for one channel at a time with gij = g ij , and the purple curve is for all six couplings non-zero. The (31) and (32) cases only slightly differ. The light blue region is the widest energy range that IceCube is likely to be sensitive to.
consider both ν → ν and ν →ν channels as is relevant for Majorana neutrinos. To better illustrate the effect of our parameters on the probability, we select a single channel, ν µ → ν e , to focus on and show the rest in fig. 4.
In fig. 2 we show the oscillation averaged probability for ν µ → ν e in the case of invisible decay. This shows that each specific flavors will receive a unique energy dependent modification that also depends on the structure of the coupling matrices, g ij and g ij . In fig. 3, we show the visible decay case for the benchmark point, and then show the effect of varying which coupling is turned on, the absolute neutrino mass scale, the initial spectrum, and the redshift evolution of the sources. Varying each parameter yields a unique modification of the probability and, in principle, can be identified up to experimental precision and astrophysical uncertainties. For the redshift evolution, in addition to R(z) = δ(z − 1), we also consider two from various fits to astrophysical objects which could potentially be sources. We refer to the first one as HERMES [100], with m 3 and z c 1.5, and the second as YKBH [101], with p 1 = 3.4, p 2 = −0.3, p 3 = −3.5, and k = 10 up to z max = 5. We see that R(z) = δ(z − 1) fits inbetween m1 γ S,PS R(z) ν,ν 10 −6 10 −6 10 −6 0 eV 2 S+PS δ(z − 1) ν → ν, ν →ν these two redshift evolution functions in the bottom right panel of fig. 3, justifying its use as a benchmark parameter; it is considerably simpler computationally.
The ν µ → ν e channel is not the only relevant channel for astrophysical neutrino decay; there are nine channels governing neutrino decay in total, but six are relevant in our context. Since no ν τ 's are produced in sources, the ν τ → ν α channels are not relevant. In addition, we note that even though this is an oscillation averaged calculation and we have assumed a loss of coherency throughout (hence no CP violating interference terms), P αβ = P βα for α = β except at high energies (that is, the SM). This is, of course, because while ν 3 → ν 1 , ν 3 → ν 2 and ν 2 → ν 1 decays may proceed, the reverse are not kinematically allowed. There are thus six main channels of interest. Fig. 4 shows the probabilities for these six channels in the SM and for our benchmark visible decay scenario. In this figure we can see that P µe > P eµ for all energies as expected in the normal mass ordering where a regeneration of electron neutrinos is expected to be larger than a regeneration of muon neutrinos since electron neutrinos contain more ν 1 than muon neutrinos do.

V. ICECUBE
IceCube is sensitive to the high energy astrophysical neutrino flux over a range of energies spanning from ∼ 10 TeV to ∼ 10 PeV, although most of their sensitivity is between 100 TeV and 1 PeV. Below ∼ 100 TeV the backgrounds from atmospheric neutrinos and muons become dominant and above 1 PeV the flux has fallen off to just a few events thus far.
IceCube has some sensitivity to the flavor of the neutrino, primarily through the identification of track and cascade topologies. If a neutrino interacts via a charged current (CC) interaction as a ν µ , a muon will be produced leaving a long track. On the other hand, a ν e or ν τ CC interaction will result in an electromagnetic or hadronic shower, making them very hard to differentiate, although in principle possible [102]. However, neutral current (NC) neutrino interactions also result in hadronic showers, although with ∼ 1/3 of the total neutrino energy. It is also possible for a ν τ CC interaction to result in a track if the tau decays to a muon, although this branching ratio is 17.4% [103] and the resulting muon carries about ∼ 1/3 the energy as it would from a ν µ CC interaction. Finally, tracks and cascade topologies are sometimes misidentified [104]. The corrections due NC    interactions and misidentification are very small and so, to an excellent approximation, the track flux is just the ν µ flux and the cascade flux the sum of the ν e and ν τ fluxes [62]. In addition, IceCube has some, albeit very limited, sensitivity in identifying ν τ CC interactions by the spatial or temporal separation between the initial hadronic shower and the subsequent tau decay, although only 1-2 events have been detected using these methods so far [105][106][107][108].
An experimentally motivated useful quantity to consider is the track to cascade ratio which we approximate by, There are some corrections to this as mentioned above, but their effects are very small primarily due to the steeply falling spectrum [62].
In the SM, this ratio R tc should be constant in energy and independent of the source spectrum or the number of source classes 6 . The track to cascade ratio, R tc , is shown in fig. 5 for the benchmark parameters assuming full pion decay (1 : 2 : 0) at the source which oscillates/decoheres to ∼ 1 2 without neutrino decay. While both the visible and invisible cases are quite similar, specific structure in g ij can change this such as if only g ( ) 21 is non-zero. 6 The only exception to this is in the case of damped muon decay where the source is dense enough that high energy pions can be produced and decay, but the muons produced from pion decay lose a significant amount of energy before decaying. We do not consider this case here because a) it will only happen within Ice-Cube's region of interest for a fairly narrow range of astrophysical parameters and b) it yields a fairly small effect (∆γ f 0.2) at the Earth after oscillations even for optimal parameters [62].  While neutrino decay leads to a non-trivial spectrum in each flavor, since IceCube fits a power law to their spectrum, we follow the same procedure to determine the effect of neutrino decay at IceCube. We fit an SPL to the track and cascade spectra in the energy range of 100 TeV to 1 PeV 7 and plot the apparent spectral index as a function of the coupling in fig. 6. The feature at lower couplings comes from the atmospheric g ( ) 31 and g ( ) 32 couplings while the feature at higher couplings come from the solar g ( ) 21 couplings. This is because the solar couplings lead to an effect at lower energies than the atmospheric couplings for a given value of g ( ) ij as seen in fig. 2 and the top left panel of fig. 3. Then, since we are looking at a fixed energy range for IceCube, in order to keep ΓL ∼ 1 fixed, g and E must be changed together. Thus for a fixed energy range the atmospheric features will occur at smaller couplings than the solar feature. Fig. 6 shows a somewhat stronger effect from invisible decay than visible decay. This is because the depletion term is non-positive and the regeneration term is nonnegative. In the invisible case the flux is only lower than in the SM, while the visible case may partially cancel that out, depending on the parameters; compare e.g. only g ( ) 32 > 0 in fig. 2 and the top left panel of fig. 3. We note that this trend is generally true as the regeneration term can be considerably larger than the depletion term in some cases.
To fully investigate this effect we plotted the difference in track and and cascade spectral indices, ∆γ f ≡ γ f,c − γ f,t in fig. 7 for invisible and visible decay as a function of both the coupling and the mass scale. We found that the initial spectral index has almost no effect on ∆γ f for visible decay (it has exactly no effect for invisible decay) as it shifts the spectrum of both tracks and cascades nearly equally. The maximum difference in final track and cascade spectral indices for various coupling configurations are listed in table III. For invisible decay ∆γ f ≥ 0 always holds. For visible decay it nearly only holds, however we find that for the cases when only g ( ) 21 (g ( ) 31 ) is non-zero that the minimum is −0.011 (−0.01). That is, other than a few very small exceptions, it is always the case that the cascade spectrum is softer than the track spectrum in the presence of neutrino decay. This is conditioned upon the normal mass ordering, in the inverted ordering the situation is partially reversed.
IceCube has reported a measurement of the track spectral index of γ f,t = 2.13 ± 0.13 [3] and a measurement of the cascade spectral index of γ f,c = 2.67 ± 0.07 [110] over somewhat different energy ranges. A simple statistical test yields a difference of χ 2 = 11 which can be interpreted as 3.3 σ tension and, in turn, 3.3 σ preference but should be accounted for in any fit to the data.
for invisible partial neutrino decay [62]. To illustrate the preferred parameters the simple χ 2 test statistic we use is, where We take the track and cascade spectral indices over E f ∈ [100 TeV, 1 PeV] as before which is slightly different from the experimental data but does not significantly affect our results. We also considered a constraint on the sum of the neutrino masses from cosmology, χ 2 m1 = ( i m i /0.06 eV) 2 based on the constraint from Planck TT, TE, EE+lowE+lensing+BAO data sets [99] which is that i m i < 0.12 eV at 95% CL. The effect on the preferred regions was negligible, so we fixed m 1 = 0.
The preferred regions for invisible and visible decay are shown in fig. 8 for the benchmark decay parameters and varying the coupling (all six couplings are taken to be equal) and the initial spectral index. Consistent with our previous simple estimate, we find some places where the χ 2 is smaller in the case of decay than in the SM; these regions are indicated in blue with the darker blue regions being the regions that are most preferred. In addition, we note that there is a horizontal line right on the edge of the allowed region for both large and small couplings in the fit at γ = 2.51 which is the best fit spectral index in the SM. The scenario of full decay might seem like an extreme one, but since the spectra of each flavor return to the same one as initially it is difficult to probe this. While the flavor ratio is different in this scenario (see B), since the measured spectral indices are different any comparison with the flavor data only makes sense within a certain energy range. For example, in the full decay scenario the flux at the Earth is pure ν 1 which would predict a deficit of tracks and a dominantly ν e flux. This is compatible with the IceCube flavor data at lower energies which is essentially the reason why the data prefers a neutrino decay scenario over the SM.
Assuming that these spectra are confirmed with future data including e.g. IceCube-Gen2, KM3NeT, and Baikal, we would find ∆γ f = 0.54. At the largest values of ∆γ f from neutrino decay shown in table III of 0.44 (0.40) for invisible (visible) decay which leads to a consistent result at χ 2 = 0.14 (χ 2 = 0.25) which can be interpreted as ∆χ 2 → 3.3 σ model preference in either case assuming Wilks' theorem. A more detailed analysis including more spectral information than the slope for an SPL fit would be helpful in identifying the features of neutrino decay; these analyses require unfolding the final neutrino spectrum from the data under a given assumption and thus are difficult to perform given publicly available data. In addition, a partial deficit in the ν τ flux at lower energies would be clear evidence of neutrino decay [62].

VI. DISCUSSION
While the effects of invisible neutrino decay are fairly straightforward -a depletion of events below a certain energy -and have been previously investigated, the effects of visible decay including a proper treatment of the spectra are much more complicated. Visible decay depends not only on the couplings, neutrino energy, and the distance traveled, it also depends on the absolute neutrino mass scale, the initial neutrino spectrum, the nature of the couplings, and the nature of the masses of the neutrinos. We now discuss the several trends seen in the previous sections. We see that the higher the mass scale, the larger the flux. This is because when m 1 is larger, the final neutrinos have energies closer to the initial neutrinos. Since we are considering steeply falling spectra, this results in a larger effect from regeneration compared with a smaller mass scale wherein the final neutrinos may have much lower energy and become swamped out by the steeply falling flux. In the limit m 1 → ∞ we recover the case where E f = E i even for regeneration since x ij → 1 in this limit and the range of the integral over initial (and intermediate for the case with two decays) energy becomes vanishingly small. However, as shown in fig. 3, even for the largest physically allowed values of m 1 from cosmology deviations from this limit already exist.
The source spectral index also has a very large effect on the flux. If the spectrum is very hard (say, γ = 1 within IceCube's region of interest) then the regeneration component is extremely large. Softer and softer spectra result in less contribution from the regeneration term and, as γ → ∞ we would recover the invisible decay case, although this is not a particularly physical limit. If, however, the neutrino mass scale is arbitrarily large, (m 1 1 eV) then the regeneration term is independent of the initial spectrum. These features are shown analytically in section III without cosmology, but the trends still hold with cosmology.
We also investigated the effect of scalar interactions, pseudo-scalar interactions, or both and found only small differences among those cases. In addition, the channel ν → ν is quite similar to the ν →ν channel. The first channel is relevant for Dirac neutrinos as it is the only one that contributes and the second is relevant for Majorana neutrinos in the case where one can differentiate neutrinos and antineutrinos in the source. Since IceCube has quite limited capabilities to differentiate neutrinos from antineutrinos, we consider both simultaneously: ν → ν and ν →ν. This also has the added benefit of ensuring that the double decay case ν 3 → ν 2 → ν 1 is handled correctly.
Depending on the couplings the effect of neutrino decay falls into one of three categories. For small couplings such that ΓL 1 in the energy and distance ranges of interest, the SM is recovered. For large couplings such that ΓL 1 in the energy and distance ranges of interest, we enter the full decay scenario where the spectrum is the same as in the SM, but the relative normalizations of the different flavors will be different depending on which g ( ) ij couplings are large. For cases when the ΓL ∼ 1 for the relevant energies and distances, we have partial neutrino decay which, in addition to modifying the normalizations of each flavor, can also have significant effects on the shape for each flavor. This opens up the possibility for detection under the reasonable assumption that astrophysics and standard oscillations cannot reproduce such an effect. We see in figs. 7-8 that the expected allowed regions in parameter space for invisible and visible decays are fairly similar. It is, in principle, possible to differentiate between the two, but only under certain assumptions.
If the spectra of two (or three) different flavors is measured in detail and they are found to recover the same spectrum at high energies, visible and invisible decay can be distinguished if one assumes that the initial source spectrum continues down throughout the entire detection regime since invisible decay requires a softer initial spectrum than visible decay to fit the same data. Alternatively, a very fine grained measurement could reveal a slight bump like structure in the spectrum (see many of the curves in fig. 3, although this could be misinterpreted as evidence for dark matter decaying or annihilating to neutrinos [62,[111][112][113][114][115][116][117][118][119][120]).
The effect of neutrino decay can be identified at Ice-Cube and is a possible explanation [62] for the 3 σ tension in the track to cascade spectra [3,62,64,95]. If the spectra continue to disagree with more data, this may become compelling evidence for neutrino decay. A careful investigation of the spectrum of each flavor, including new information about ν τ 's, would be necessary to determine the particular parameters, in particular looking at the texture of the g ij and g ij matrices. It is very interesting to note that neutrino decay in the normal mass ordering predicts a cascade spectrum that is softer than that of tracks; if IceCube had measured the reverse then neutrino decay would not explain the data (or we would have to invoke neutrino decay in the inverted mass ordering). If such a signal of neutrino decay is confirmed, we anticipate it to have couplings of g ( ) 3j ∼ 10 −7 and g ( ) 21 ∼ 10 −6 as shown in fig. 8 in order to have the desired feature within IceCube's region of interest. We also see in fig. 8 that the preferred initial spectral indices are softer in the presence of neutrino decay (2.80 and 2.77 for invisible and visible decay respectively) than in the SM (2.51). This is consistent with expectations as neutrino decay tends to harden a spectrum so matching the data requires a softer initial spectrum when decay is included. Moreover, we see that invisible decay prefers a slightly softer initial spectrum than visible which is due to the regeneration term in visible decay partially canceling out the hardening from the depletion term.
While next-generation experiments such as GRAND and POEMMA [121,122] (in addition to current experiments such as ANITA [123]) will have unique sensitivity to tau neutrinos thus providing important flavor information that is generally difficult for any experiment (including IceCube) to detect, their higher energy thresholds makes them sensitive only to larger couplings than Ice-Cube. Due to the facts that they will only be measuring a single flavor and the source flux is relatively unknown, disentangling neutrino decay features from astrophysical features will be extremely difficult.

VII. CONCLUSIONS
Visible neutrino decay has a rich phenomenology that can be probed at IceCube by simultaneously measuring the flavor and energy of the high energy astrophysical neutrino flux. Due to the extremely long distances the neutrinos at IceCube travel en route to the Earth, Ice-Cube provides one of the most sensitive direct probes of neutrino decay itself. In addition, while the initial flux is unknown, IceCube compensates for this by detecting all three flavors and partially differentiating among them. In this article we have calculated the flavor transition probabilities for all channels as a function of the initial neutrino spectrum, the absolute neutrino mass scale, and the various couplings for the first time. This calculation involves an integral over the distance at which the decay happens and the initial neutrino energy which is larger than the final neutrino energy. We have also included the effect of multiple decays involving extra integrals over the second decay point and the energy of the intermediate neutrino.
Finally, we incorporated the effect of cosmology and the expansion of the universe.
Neutrino decay in the partial decay regime modifies the spectrum of neutrinos seen at IceCube differently for different flavors, while in the SM there can be nearly no difference. This makes IceCube a powerful probe for neutrino decay. Since IceCube has measured this flux and has some sensitivity to flavors, it provides a powerful probe of neutrino decay. As IceCube sees some hints at ∼ 3.3 σ that the spectra may be different for different flavors, neutrino decay offers a mechanism for explaining this difference. The fact that it is the cascade spectrum that is softer than the track spectrum is consistent with neutrino decay in the normal mass ordering. The preferred neutrino decay parameters are shown in fig. 8 and more detailed measurements from IceCube, KM3NeT, and Baikal can further test neutrino decay in coming years.

Regeneration Component
Starting with eq. 20, we reinsert the phases and move to the flavor basis. We can construct the full regeneration transition probability as, As the decay ν 3 → ν 2 → ν 1 does not interfere with the single decays, it does not contribute to the oscillations. We can therefore compute the first term, and simply add on the expression given in eq. 25.
To evaluate the integrals, we use that, where φ = arg (z). After some algebra, we obtain  where A, B, C are prefactors, and φ, ξ, ψ are CP violating phases, given as: Note that we have not integrated over the initial, or intermediate, energies here. We would have to integrate the terms separately, as described in sec. II B, as the integration limits would vary according to the decay in question.