Exploring a Non-Minimal Sterile Neutrino Model Involving Decay at IceCube and Beyond

We study the phenomenology of neutrino decay together with neutrino oscillations in the context of eV-scale sterile neutrinos. We review the formalism of visible neutrino decay in which one of the decay products is a neutrino that potentially can be observed. We apply the formalism developed for decay to the recent sterile neutrino search performed by IceCube with TeV neutrinos. We show that for $\nu_4$ lifetime $\tau_4/m_4 \lesssim 10^{-16} {\rm eV^{-1}s}$, the interpretation of the high-energy IceCube analysis can be significantly changed.

The outline of this paper is as follows. In section II we introduce the neutrino decay model. In section III we develop our framework for models of neutrino oscillation with neutrino decay for two scenarios. First we consider the case where the neutrino daughters are all Beyond the Standard Model (BSM) particles, which we term "invisible decay". Secondly, we consider the case where one neutrino daughter is a lighter SM neutrino, which we refer to as "visible decay". In section IV we illustrate our model in the IceCube experiment using atmospheric neutrino data. Finally, in section V, we provide concluding remarks.

II. NEUTRINO DECAY MODELS
In this paper, we are interested in neutrino decay via new interactions. The simplest cases are those in which the neutrino decays into two particles. These are expected to be dominant over decays into three or more particles, as these involve high dimensional operators [77]. Fig. 1 shows a diagram of this process. We can organize the two-daughter processes into (I) invisible and (II) visible decays. Note that in (I), if one of the new particles decays into neutrinos, then this can be reduced to an effective three-body decay, which we have chosen to neglect. From these statements, it follows that in (I) ν i → φψ and in (II) ν i → ν j φ, where ψ is a fermion and φ a boson. Due to it being phenomenologically more interesting, we will restrict our examples to the visible decay scenario case, i.e., ψ is one of the SM neutrinos, and consider φ to be either a scalar or pseudoscalar.

A. Lagrangians of simplified models
In this work, we will use simplified models to introduce the aforementioned decay scenarios, as done, for example, in [78,79]. For simplicity, we write this section assuming that the neutrino is a Majorana field. This assumption can be relaxed by introducing right-handed neutrinos [77], which would be related to Dirac mass terms and would give rise to a different phenomenology due to the lack of ν →ν transitions and visible decay. With this assumption, the scalar-neutrino interaction is [80] where g s ij (g p ij ) are the (pseudo)scalar couplings that control the transition strength from parent (i) to daughter (j): ν i → ν j φ. In other words, g carries two indices: one for the parent mass eigenstate and another for the daughter neutrino mass eigenstate. Note that we assume = c = 1 throughout, except where we explicitly restore constants to estimate decay lengths.

B. Decay rates
In the case where the boson is massless, our scenario is analogous to neutrino-Majoron decays. As in that case, there are two scenarios induced by Eq. (1). The first one is a chirality-preserving process, whose partial decay rate is [78,81] Γ(ν where we have introduced the parent-to-daughter mass ratio, x ij = m i /m j > 1, and have labeled the energy of the parent neutrino, ν i , by E p . The second one is a chirality-violating process (3) In both cases, the auxiliary functions are given by [81] f (x) = where we have dropped the indices for clarity. If the lightest neutrino is massless, then for approximately eV sterile neutrinos, the smallest x ij ∼ 10 2 1, while if the lightest neutrino saturates the current kinematic limits, then the smallest x ij ∼ 10. In the limit of m i m j , the partial decay rates given in Eq. (2) and Eq. (3) are just [81] In any case, the neutrinos that we are considering are relativistic, thus the products of the decay will travel along the direction of the beam [78] and the relevant quantity is the energy distribution of the daughter in the lab frame. In the relativistic limit, the expression for this quantity simplifies if we assume either pure scalar (g p ≡ 0) or pure pseudoscalar (g s ≡ 0) cases. We will make this assumption throughout the rest of the paper, and in particular, we will assume a pure scalar couplings in our analysis. For the chirality-preserving process, the energy distribution takes the form [79] where + corresponds to the pure scalar case and − to the pure pseudoscalar one, and E d labels the energy of the daughter neutrino, ν j . For the chirality-violating process, the distribution takes the form In both cases, due to kinematic constraints, the daughter energy is bounded to: [79]. In the limit with x ij 1, which is the case for approximately eV sterile neutrinos, there is no distinction between the scalar and pseudoscalar scenarios. As can be seen from the energy dependence of these relationships, the chirality-preserving processes produce harder daughters, while the chirality-violating ones produce softer daughters (see Fig. 2). It is important to note that the chirality-preserving processes dominate the visible decay  due to the fact that the atmospheric neutrino flux is a steeply falling power law, φ(E) ∼ E −3.7 [82,83]; the soft daughters in the chirality-violating case are hidden below a large flux which reduces the visible case to the less interesting invisible decay scenario. The chiralityviolating process in the Majorana context induces transitions between neutrinos and antineutrinos, but not if the term arises from Dirac neutrinos. In this latter case, the right-handed neutrinos are sterile and do not interact. Note that x ij diverges as m j → 0, but all decay rates and differential decay rates remain finite in this limit, so decays to massless daughters are well-defined in this framework. Finally, to easily compare with existing constraints on neutrino lifetime, we introduce the partial lifetime τ ij = 1/Γ ij .

C. Existing constraints
Constraints exist on the lifetime of neutrino mass eigenstates ν i for i < 4, i.e., for the active neutrinos. The constraints depend on the neutrino mass ordering as well as the absolute neutrino mass. In the normal ordering (NO), ν 3 and ν 2 are unstable and can decay to ν 1 , which is stable. On the other hand, in the inverted ordering (IO), ν 3 is stable and the others unstable. Cosmology constrains both the sum of neutrino masses, such that i m i 0.12 at 95% C.L. [84], and the radiative neutrino decay lifetime, to 10 19 s [85,86]. However, constraints from cosmology on the presence of a fourth neutrino are model-dependent; assumptions required include the thermal history of the Universe and the influence of dark energy on the expansion history [87]. Moreover, sterile neutrino thermalization can be suppressed by a number of new physics scenarios [88][89][90][91][92][93][94][95][96]; if sterile neutrinos do not thermalize, the bounds do not apply.
On the other hand, given results from neutrino oscillation experiments [97], all neutrino masses cannot be zero. The largest decay rate will be achieved when the stable neutrino is assumed to be massless, i.e., m 1 = 0 or m 3 = 0, for NO or IO, respectively. With this optimistic assumption, in IO the unstable neutrino masses cannot be less than ∼ 0.05 eV and in NO m 2 cannot be smaller than ∼ 0.008 eV and m 3 satisfies the same bound as for IO. The lifetime of a neutrino of mass m i is given by τ −1 i = j τ −1 ij , and the total decay rate in the lab frame is a function of τ i /m i . Constraints on this quantity will depend on the mass ordering: for NO, ν 2 decay is constrained by solar experiments, with τ 2 /m 2 7 · 10 −4 s eV −1 [98], and ν 3 is limited by atmospheric and long-baseline experiments, with τ 3 /m 3 9 · 10 −11 s eV −1 [60]; for IO, τ 1 /m 1 4 · 10 −3 s eV −1 and τ 2 /m 2 7 · 10 −4 s eV −1 are constrained by solar experiments [98]. In contrast, direct constraints on ν 4 have not been set.
Constraints on neutrino decay can be obtained indirectly from measurements of meson decays [99][100][101]. These set strict limits on the neutrino decay process but are flavor-dependent; thus they must be used with care. For example, from kaon decay, the following combination is constrained [99] where α runs over all neutrino flavors. Constraints from supernova 1987A are of similar order [100]. These flavor couplings are related to the ones in Eq. (1) by For simplicity, let's first consider the case where only one g 4j is nonzero. In this case g αβ = g 4j U α4 U * βj , where for short-baseline motivated sterile neutrinos U α4 ∼ O(0.1) [102] and from standard neutrino measurements U βj ∼ O(0.1) for j < 4 [3], which implies that g 4j ≤ O(0.1). For this size of coupling, an eV-scale sterile neutrino lifetime satisfies τ 4 > O(10)/eV. With this lifetime constraint, as we will see in a later example, no interesting interplay exists between oscillations and decay. This is due to the fact that the scales of oscillation and decay, which are given by E/∆m 2 4j and Eτ 4 /m 4 , respectively, are very different. This implies that to have an interesting interplay, more than one of the g 4j needs to be nonzero so that a cancellation can occur in Eq. (9), leading to a decreased bound in τ 4 . We summarize the constraints discussed in this section in Fig. 3 for both active mass states and the mostly sterile state, ν 4 .

III. NEUTRINO OSCILLATIONS AND DECAY
We will now develop a formalism that incorporates both oscillations and decay in a consistent way. We will follow the calculation in [60,103] for invisible decay, and for visible decay we follow [64,78,79] but in the density matrix representation. We will work with a (3+1) model, which adds one sterile state to the three active neutrinos, although the model generalizes readily to (3 + N ) for generic N .

A. Neutrino oscillations
We will begin with the most familiar model, which includes vacuum oscillations exclusively, ignoring decay and matter effects. We have where ∆M 2 is a diagonal matrix with entries ∆M 2 ii = ∆m 2 i1 : the neutrino mass-squared splittings. Conjugation into the flavor basis with the PMNS matrix U ν (appropriately extended to include the sterile state) mixes the mass states into flavor states as follows This mixing gives rise to vacuum oscillations when at least one of the mass splittings is non-zero. We are interested in analyzing TeV-scale neutrinos in IceCube. At these energies, matter effects become important, so we introduce a term modeling neutrino scattering off of electrons and nucleons as an effective matter potential [104][105][106]. The flavor states are eigenstates of the corresponding Hamiltonian term V matter . In the flavor basis, this term takes the form V matter = diag(V e , V µ , V τ , 0), where the fourth eigenvalue is zero because it corresponds to a sterile neutrino state. Thus the total Hamiltonian in the mass basis is given by where the neutrino potentials depend on the electron (N e (l)) and neutron (N n (l)) number density profiles, and l is the position of the neutrino ensemble along the baseline.

B. Invisible neutrino decay
Modeling decay in an evolving neutrino system is not a trivial matter, because it involves the addition of daughter particles to the state space. In the visible decay case, we will treat these new states explicitly. In the case of invisible decay, the daughter states are irrelevant, and one can get around this added complexity by introducing an effective non-Hermitian Hamiltonian, H, which gives rise to a non-unitary evolution when restricted to the state space of the parent neutrinos. Intuitively, the nonunitarity of the associated time evolution operator corresponds to the loss of probability current from the parent state space into the daughter space. Following [103] we can construct H by the explicit addition of an anti-Hermitian term where Γ is a Hermitian operator that is diagonal in the mass basis, where it is just Γ = diag(Γ i (E), ...), where Γ i (E), the decay rate of the i th neutrino, is given by Eq. (5). The factor of 1 2 is necessary for the survival probability of a neutrino born in the ν i mass eigenstate to follow the exponential decay formula In a two-flavor system, the active neutrino vacuum survival probability can written, under the assumption that the lighter of the two neutrinos is stable, as [60] P α→α = cos 4 θ + 1 2 e − Γ 2 l 2 cos ∆m 2 l 2E sin 2 2θ where θ is the two-flavor mixing angle and Γ 2 is the decay rate of the heavier mass state with mass m 2 . Equations (14) and (15) are valid only in the case of invisible decay in vacuum and are included here for completeness.

C. Visible neutrino decay
In this paper, we will concentrate on visible neutrino decay. In this case, Eq. (13) cannot describe the full system evolution because it ignores the evolution of the daughter states. To treat the daughter states explicitly, we first need to extend our formalism from a single, monoenergetic neutrino state to a set of neutrino states indexed by the set of energies relevant to our analysis. We then promote each of these states to an ensemble of states, described by an n × n density matrix ρ(E), where n is the number of neutrino species under consideration. In the density matrix formalism, evolution due to the Hamiltonian given in Eq. (13) can be written as where l is the position of the neutrino ensemble along the baseline, square brackets indicate the commutator, and curly brackets indicate the anticommutator. The factor of 1 2 appears for a similar reason as in Eq. (13). This formalism has already been implemented in an efficient way in the context of neutrino oscillations in a package called nuSQuIDS [107,108]. The decay rate operator is given by where Π i is the projector to the i th mass eigenstate and Γ i is given by the sum of partial rates Γ ij over all daughter states ν j lighter than the parent state ν i . The Γ ij are given by the sum of the chirality-preserving Γ ij from equation (2) and the chirality-violating Γ ij from equation (3) in the case of Majorana neutrinos, or by partial rates Γ ij corresponding to purely chirality-violating processes in the case of Dirac neutrinos. Thus far, we have modeled only the loss of neutrinos using the Γ term in Eq. (16). The advantage of this extended formalism is that it can accommodate terms that generate transitions between neutrinos at different energies. We then complement the Γ term describing the loss of neutrinos from an ensemble at one energy with a "regeneration" term, R, describing the appearance of these neutrinos in ensembles at lower energies. In this way, we keep track of the daughter states instead of neglecting them, and we account for their subsequent evolution.
Let us assume that neutrinos are Majorana. For clarity, we will separate the contributions to regeneration into the chirality-preserving processes (CPP) and the chirality-violating processes (CVP). In the first case, we add the following term to the right-hand side Eq. (16) where E p is the parent energy, E d is the daughter energy, Tr is the trace operation, and the differential decay rate is given in Eq. (6). When this term is added to Eq. (16), we set E d = E. We have an additional contribution from ν →ν in the chirality-violating process, in which case we should add whereρ corresponds to antineutrinos and the differential rate is given by Eq. (7). Because we consider either purely scalar or purely pseudoscalar cases, the differential rates for CVP and CPP must be chosen accordingly when constructing R. This boils down to getting the signs right in equations (6) and (7). Similar equations hold for antineutrinos by replacing ρ byρ and changing the projectors Π toΠ, and vice versa. Note that there is no regeneration in the Dirac case, because the decay is through a chirality-violating process that produces sterile daughters. We have implemented this formalism in nuSQuIDS [107,108]; see Appendix A for details.

IV. EXAMPLE OF APPLICATION TO ICECUBE
The IceCube experiment is an ideal testing ground for the model presented here. The collaboration has released a 400 GeV to 20 TeV single-year data set associated with a sterile neutrino search [9], and an additional six years of data are available to the collaboration. In this section, we use the released data set to explore some of the nuances of including decays as additional phenomena in the IceCube search. We will only consider the decay channel ν 4 → ν 3 φ, with the following ν 4 lifetimes: (A) τ 4 = 10 eV −1 , (B) τ 4 = 1 eV −1 , (C) τ 4 = 0.1 eV −1 , and τ 4 = ∞ eV −1 , corresponding to no decay. For τ = 1 eV −1 , cτ ≈ 0.2 µm.
A. The IceCube experiment and released data set The IceCube detector, located in the Antarctic ice below the South Pole station, observes Cherenkov radiation from interactions of neutrinos that have traversed the Earth. In the 400 GeV to 20 TeV energy range, muon neutrinos are primarily due to decays of kaons produced by cosmic rays impinging on the atmosphere [83]. This analysis makes use of through-going muons that are produced in charged-current ν µ interactions in the ice or bedrock below the detector. The muons produced by neutrinos in this energy range are above critical energy. Hence, the muon energy can be determined from the stochastic light emission profile. The direction of the muon can be determined through reconstruction of the Cherenkov-light time and spatial distribution and is expressed as an angle, θ Z , with respect to the zenith.
The Cherenkov light is observed via digitial optical modules (DOMs), with photomultiplier tubes as the light sensors [109,110]. The detector consists of 5,160 DOMs on 86 vertical strings. The intrastring DOM separation is 17 m and the interstring separation is approximately 125 m. The energy resolution of the detector is σ log10(Eµ/GeV) ∼ 0.5 and the angular resolution, σ cos(θ Z ) , varies from 0.005 to 0.015.
The released data set contains 20,145 wellreconstructed events, described in [9,111,112]. This release is associated with a sterile neutrino search with null results at 90% CL [9]. The power of the sterile neutrino analysis arose from the matter effects that were expected for neutrinos that crossed the core and mantle [113]. This was predicted to lead to an observable deficit of ν µ events in IceCube for parameters that were consistent with short baseline anomalies [16,17]. The deficit was expected to be localized to specific regions of E µ and cos(θ Z ), depending on the sterile neutrino parameter space. As a result of the striking signature, IceCube was able to perform a powerful search. The null result significantly changed the parameter-landscape for 3+1 searches [102].
The released data provides energy and angle information for each data event. It also provides Monte Carlo information on an event-by-event basis. Information on handling of multiple error sources is provided. Using this data set, one can reconstruct the IceCube search results well, as was demonstrated in [102]. The appendix of that paper provides step-by-step instructions for use of the data release, which we follow here.

B. IceCube oscillograms
We refer to an "oscillogram" as a plot of the expected change in neutrino flux from creation in the atmosphere to arrival at IceCube, as a function of true neutrino zenith angle and true neutrino energy. Effects that may change the neutrino flux include oscillation, matter effects, absorption, and decay. The minimal 3+1 ster-ile neutrino model is parameterized by a mass-squared splitting, ∆m 2 41 , and a mixing angle, θ 24 . Fig. 4 shows the shape effects in theν µ spectrum for the 3+1 sterile neutrino model with parameters ∆m 2 41 = 1 eV 2 and sin 2 2θ 24 = 0.1, for the four lifetimes listed previously. In Fig. 4, for no decay scenario as well as for examples points A and B, the depletion ofν µ s in the region E truē νµ ∼ 300 GeV and cos θ true Z ∼ -1.0 is due to matter effects. Decreasing the lifetime of ν 4 decreases the magnitude of this feature and shifts its position. This is due to the fact that, as the ν 4 lifetime becomes smaller, its decay length is smaller than the oscillation scale. In other words, the decay operation breaks the coherence of the system, projecting into the mostly active mass states, preventing the development of oscillation. However, the final flux is still different from the flux in the absence of oscillations. The depletion in the top-left corner of each plot in Fig. 4 is due to absorption of high-energy neutrinos crossing the Earth.

C. Data analysis and systematic uncertainties
The data is binned in reconstructed energy proxy and zenith angle. We use a binned Poisson log-likelihood function, log L, for the data, and incorporate systematic uncertainties by means of nuisance parameters, η, with Gaussian priors. For each lifetime considered and for each point in a fine, logarithmic grid of [sin 2 (2θ 24 ), ∆m 2 41 ], log L(sin 2 (2θ 24 ), ∆m 2 41 , η) is maximized with respect to a set of nuisance parameters.
Along with the "conventional flux" from pion and kaon decay, in principle a "prompt" contribution arises from the decays of heavier mesons, but it has yet to be observed [114]. For this reason, the prompt flux is neglected [115]. The atmospheric neutrino flux is the sum of the neutrino component and the antineutrino component, where the neutrino component is parameterized as where φ π ν and φ K ν are the fluxes of neutrinos originating from decays of pions and kaons, respectively. The overall flux normalization, N 0 ; variations to the spectral index, ∆γ; and the ratio of neutrinos originating from kaons and pions, R K/π , are nuisance parameters. The pivot point for the spectral index change, E 0 , is at a midpoint energy such that changes to the spectral index do not dramatically change the overall flux normalization. The antineutrino flux is parametrized identically to the neutrino flux, up to a relative normalization factor, which is another nuisance parameter. The cosmic ray and hadronic models used are Poly-gonato [116] and QGSJET-II-4 [117], respectively. Uncertainty in the atmospheric density profile is accounted for in a linear parameterization, F(δ) [111,112].
The DOM efficiency is a final nuisance parameter. Monte Carlo data sets corresponding to several discrete values of DOM efficiency are publicly available [118]. A piecewise linear interpolation was fit to them, allowing us to treat the DOM efficiency as a continuous nuisance parameter.

D. Results
We have chosen a few specific parameters to illustrate the effect of this model. To demonstrate that our analysis technique is robust and properly implemented, we first perform the analysis without decay. Fig. 5 shows the best-fit nuisance parameters as functions of the sterile parameters.
The oscillograms shown in Fig. 4 indicate that introducing neutrino decay diminishes the strength of the sterile neutrino effect. In order to compare how much the model power changes when we introduce decay, we compare the profile likelihood with and without decay for the scenarios discussed in the previous section. The difference in profile likelihood as a function of sterile neutrino parameters is shown in Fig. 6. In this figure, red colors indicate that the no-decay scenario is preferred whereas blue colors show preference for the decay solution. As the lifetime decreases, the decay scenario is preferred over the no-decay scenario for ∆m 2 41 ∼ 0.1 − 1 eV 2 . It is worth noting that the saturated blue color does not necessarily indicate a good fit to data.
Our main result is illustrated in Fig. 7, which compares the sterile neutrino hypothesis with decay to the standard three neutrino scenario. In order to quantify the difference between models, we used the approximate Bayes factor, B 01 , as a function of the sterile mixing angle and mass difference for the lifetimes discussed in this paper. The Bayes factor is approximate as we have used a profile instead of a marginal likelihood. In comparing hypotheses, we use the Jeffreys scale, where 0 < log B 01 < 1, 1 < log B 01 < 2.5, and 2.5 < log B 01 < 5 correspond to weak, strong, and decisive evidence against a hypothesis, respectively [119]. We observe that for lifetimes on the order of 0.1 eV −1 , in the regions of parameter space that correspond to the allowed regions from sterile neutrino global fits to short-baseline data assuming no decay, the sterile neutrino hypothesis is not disfavored. For lifetimes greater than 1 eV −1 , the IceCube exclusion is robust under this new physics scenario.

V. CONCLUSIONS
We have reviewed the framework of neutrino decay with oscillations and presented it in a consistent manner using the density matrix formalism. We have further implemented this new physics scenario in the nuSQuIDS software package [120]. We have implemented the highenergy IceCube sterile analysis and then introduced the decay of ν 4 as an additional effect. We show that for small values of the lifetime, τ 0.1 eV −1 , the IceCube  results interpretation can be significantly changed. Thus, neutrino decay can dramatically alter the landscape of eV-sterile neutrinos and will need to be studied in the context of global fits. The dotted, dashed, and solid lines correspond to the Jeffreys scale criterion for weak, strong, and decisive evidence against the alternative hypothesis. In the left-most plot, the grey regions are the 90% confidence level credible intervals from global fits to short-baseline data [17]. In all these plots we assume the decay channel ν4 → ν3φ. For τ = 1 eV −1 , cτ ≈ 0.2 µm.