Search for Unstable Sterile Neutrinos with the IceCube Neutrino Observatory

We present a search for an unstable sterile neutrino by looking for a resonant signal in eight years of atmospheric ν μ data collected from 2011 to 2019 at the IceCube Neutrino Observatory. Both the (stable) three-neutrino and the 3 þ 1 sterile neutrino models are disfavored relative to the unstable sterile neutrino model, though with p values of 2.8% and 0.81%, respectively, we do not observe evidence for 3 þ 1 neutrinos with neutrino decay. The best-fit parameters for the sterile neutrino with decay model from this study are Δ m 241 ¼ 6 . 7 þ 3 . 9 − 2 . 5 eV 2 , sin 2 2 θ 24 ¼ 0 . 33 þ 0 . 20 − 0 . 17 , and g 2 ¼ 2 . 5 π (cid:2) 1 . 5 π , where g is the decay-mediating coupling. The preferred regions of the 3 þ 1 þ decay model from short-baseline oscillation searches are excluded at 90% C.L.

We present a search for an unstable sterile neutrino by looking for a resonant signal in eight years of atmospheric ν μ data collected from 2011 to 2019 at the IceCube Neutrino Observatory. Both the (stable) three-neutrino and the 3 þ 1 sterile neutrino models are disfavored relative to the unstable sterile neutrino model, though with p values of 2.8% and 0.81%, respectively, we do not observe evidence for 3 þ 1 neutrinos with neutrino decay. The best-fit parameters for the sterile neutrino with decay model from this study are Δm 2 41 ¼ 6.7 þ3.9 −2.5 eV 2 , sin 2 2θ 24 ¼ 0.33 þ0.20 −0.17 , and g 2 ¼ 2.5π AE 1.5π, where g is the decaymediating coupling. The preferred regions of the 3 þ 1 þ decay model from short-baseline oscillation searches are excluded at 90% C.L. DOI: 10.1103/PhysRevLett.129.151801 Long-standing anomalies in short-baseline (SBL) neutrino experiments [1,2] have been interpreted in the standard oscillation framework of three known flavors and one or more hypothetical sterile neutrinos, referred to as "3 þ N" models. The "3 þ 1" model, which involves only one sterile neutrino, has been extensively studied through global fits to datasets sensitive to vacuum oscillations involving a dominant mass splitting of ∼1 eV 2 [3][4][5]. These fits find a strong preference for 3 þ 1 over the three neutrino hypothesis [4]. However, the allowed regions from these fits suffer from internal inconsistencies between datasets [4,6]. In particular, no experiment has found evidence of ν μ disappearance, which is expected in a 3 þ 1 model. This is one motivation to consider alternative models to the 3 þ 1; another is to evade cosmological bounds on light sterile neutrinos [7][8][9][10][11] and possibly resolve the Hubble tension [12][13][14][15][16].
The IceCube Neutrino Observatory has the unique capability of performing such a test. IceCube is a cubickilometer neutrino detector buried 1.5-2.5 km beneath the surface of the Antarctic glacier at the South Pole [37]. Muon tracks from charged current (CC) muon (anti)neutrino interactions are reconstructed based on observation of emitted Cherenkov light that is collected by "digital optical modules" (DOMs) [38] arranged in vertical strings on a hexagonal lattice. Specifically, the track fitting [39] utilizes signals from two detector arrays: (i) the main array of 78 strings spaced 125 m apart, each carrying 60 DOMs with a vertical separation of 17 m between them; and (ii) the DeepCore [40] eight-string array, with lateral spacing varying from 42 to 72 m, and vertical DOM separation of 7 m.
The existence of an eV-scale sterile neutrino can manifest itself as a resonant, matter-enhanced flavor transition for either muon antineutrinos or muon neutrinos traversing the core of the Earth [41][42][43][44][45][46]. This causes a deficit of "up-going" muon (anti)neutrinos at TeV-scale energies. IceCube cannot distinguish between neutrinos and antineutrinos. Therefore, the only signature is a deficit in the combined muon neutrino and muon antineutrino (ν μ þν μ ) CC event distribution at TeV energies. A search in the framework of the 3 þ 1 model using eight years of IceCube data has recently been published [47,48]. This dataset offers an excellent platform to test the hypothesis that the 3 þ 1 þ decay model provides a better description of the data than the 3 þ 1 model without relying upon the vacuum oscillation signature.
In the 3 þ 1 þ decay model, the three-neutrino mixing matrix, U PNMS , which is parametrized by three mixing angles and one CP-violating phase, δ CP , is extended by one row and column, adding one sterile flavor state, ν s , and one heavy mass state, ν 4 . This introduces three new mixing angles, θ 14 , θ 24 , and θ 34 , two new CP-violating phases, δ 14 Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP 3 . and δ 24 , and one additional mass splitting, Δm 2 41 . Lastly, instability of the fourth mass state is introduced as in Ref. [4], governed by the strength of a coupling constant g.
For nonzero values of g, ν 4 can decay into invisible particles beyond the standard model, while g ¼ 0 returns the 3 þ 1 model. The relationship between this coupling, g, the ν 4 mass, m 4 , and its lifetime, τ, is [49] τ ¼ 16π Most of the parameters involved in three-neutrino mixing are well known; these include the light, active neutrino mass splittings and the PMNS matrix elements [50]. However, this Letter is insensitive to these parameters as well as to the neutrino mass ordering and δ CP because, for the relevant neutrino energies (E ν > 100 GeV) and baselines (order the diameter of the Earth or smaller), oscillation probabilities between the three active flavors are insignificant. For the present study, the normal mass ordering is assumed and δ CP is assumed to be zero. Furthermore θ 14 , δ 14 , and δ 24 are set to zero since they have subleading effects [51]; θ 34 is set to zero as this yields conservative results [51,52] and m 1 is set to zero since only mass differences are relevant. This leaves three free parameters in the model to be tested: Δm 2 41 , sin 2 2θ 24 , and g 2 . It is assumed that θ 24 < π=4, which causes the resonance to appear in the antineutrino flux; larger values of θ 24 are heavily constrained [4] and since there are more atmospheric neutrinos than antineutrinos [53], this choice is also conservative.
At IceCube, the ν μ andν μ disappearance probabilities vary as a function of energy and zenith angle (θ z ), where cos θ z ¼ 0 corresponds to neutrinos arriving from the horizon and cos θ z ¼ −1 corresponds to neutrinos arriving from the direction of the North Pole. This study uses a dataset collected over a live time of 2786 days and an event selection that has been described in detail in Refs. [48,54]. The predicted resonance occurs at TeV scales, hence the analysis focuses on muons from CC neutrino interactions with energies between 500 GeV and 10 TeV. Relevant neutrino interactions occur below or within IceCube. Because the signature of the analysis relies on matter effects in the Earth, this analysis requires muon tracks to have upgoing zenith angle (−1.0 < cos θ z < 0.0). The angular resolution of the tracks, σ cos θ z , lies between 0.005 and 0.015, and the track energy resolution, σ log 10 ðE μ =GeVÞ , is ∼0.5 [39]. The selected sample comprises 305 735 CC ν μ andν μ events.
The physics under study affects the flavors of the neutrinos as they propagate through the Earth. This is described using the nuSQuIDS neutrino evolution code [71,72] which accounts for both coherent and noncoherent interactions [73][74][75][76][77][78], as well as tau neutrino regeneration [79,80]. This analysis uses nuSQUIDSDecay, which incorporates the effect of ν 4 decay [34]. The Earth density profile is parametrized by the spherically symmetric PREM model [81]. The CSMS [82] neutrino-nucleon cross section is used to describe the CC interactions below and within the detector.
This analysis builds on the 3 þ 1 analysis in Ref. [48]. The data are binned in reconstructed muon energy and cos θ z , and a modified Poisson likelihood that accounts for finite simulation statistics is used to evaluate the data given sterile neutrino parameters [83]. Eighteen systematic effects related to the atmospheric and astrophysical flux, detector, and cross section uncertainties are incorporated into the likelihood function as nuisance parameters; these are described further in Ref. [48]. The treatment of most systematic uncertainties is unchanged. The dominant category of uncertainties had been identified as those associated with the atmospheric neutrino flux.
Two improvements were made over the 3 þ 1 analysis: (i) the uncertainty in the atmospheric neutrino flux corresponding to the uncertainty in the production of charged mesons in atmospheric showers is calculated using atmospheric data from the NASA Atmospheric InfraRed Sounder satellite [84], rather than the atmospheric model from Ref. [85]; and (ii) the astrophysical and prompt neutrino fluxes are calculated using a corrected depth setting of the glacial ice, compared to Ref. [48], which had little impact on the current or previous results. Combined, these changes increase the likelihood of the data for the three-neutrino model and best-fit 3 þ 1 model by, respectively, 0.09 and 0.18 log-likelihood (LLH) units.
Both a frequentist parameter estimation and a pointwise Bayesian model comparison [86] are performed, following the same procedure as in Ref. [48]. The likelihood function and Bayes factor [87] are evaluated over a grid scan of the three physics parameters-Δm 2 41 , sin 2 2θ 24 , and g 2 -where Δm 2 41 and sin 2 2θ 24 are sampled log uniformly with ten samples per decade in the ranges 0.01-47 eV 2 and 0.01-1.0, respectively, and the parameter g 2 is sampled in steps of ðπ=2Þ in the range 0 − 4π.
The best-fit parameters are found to be Δm 2 41 ¼ 6.7 þ3.9 −2.5 eV 2 , sin 2 2θ 24 ¼ 0.33 þ0.20 −0.17 , and g 2 ¼ 2.5π AE 1.5π. Theν μ disappearance probabilities are given in Fig. 1 for the parameters Δm 2 41 ¼ 6.7 eV 2 and sin 2 2θ 24 ¼ 0.33, and for two values of g 2 ; the top panel shows the situation for g 2 ¼ 0, which corresponds to the 3 þ 1 model, while the bottom panel is for the case g 2 ¼ 2.5π. The bottom panel PHYSICAL REVIEW LETTERS 129, 151801 (2022) represents the best-fit point of the frequentist analysis. Muon neutrino disappearance probabilities do not feature the resonant deficit and make subleading contributions to the sterile signature, so they are not shown.
The best-fit signal expectation and data are both compared to the three-neutrino model expectation in Fig. 2. In these plots, both the signal expectation and the threeneutrino model expectation include systematic uncertainties estimated adopting the respective physics parameters. Both the data and the best-fit signal shapes have a deficit of events for through-going neutrinos at the highest energies and a relative excess for horizon-skimming events at the highest energies. The fit values of all systematic uncertainties are within 1σ of their prior centers, with the exception of the cosmic-ray spectral index. The fit value of this systematic uncertainty deviates by 2.4σ, which is similar to both the result from the 3 þ 1 search [47,48], as well as the fit value assuming no sterile neutrino.
The frequentist confidence regions sliced at the best-fit value of g 2 are shown in Fig. 3. The contours are drawn assuming Wilks' theorem and three degrees of freedom (DOF). Fits to simulated datasets for several points in the parameter space showed the effective DOF was consistent with three or fewer. The slices of the confidence regions for the other values of g 2 are approximately the same in the 2D space of [Δm 2 41 , sin 2 2θ 24 ], with two deviations: the 90% C.L. (confidence level) region for g 2 ¼ 0 excludes any point with sin 2 2θ 24 ≳ 0.2, and above Δm 2 41 ∼ 7 eV 2 , the confidence regions extend to higher Δm 2 41 values for larger values of g 2 . This is shown in the Supplemental Material [88]. The effective DOF at the null hypothesis (only three neutrinos) was determined to be 2.86 AE 0.14, obtained by fitting 300 simulated datasets generated assuming this hypothesis. The null hypothesis is disfavored in favor of the 3 þ 1 þ decay model with −2ΔLLH ¼ 9.1 and a p value of 2.8%. This p value was obtained using Wilks' theorem and three DOF, which is conservative and consistent with the DOF assumed for the contours.
The Bayesian analysis finds the best model to have the parameters Δm 2 41 ¼ 6.7 eV 2 , sin 2 θ 24 ¼ 0.33, and g 2 ¼ 1.5π; this model has a Bayes factor (BF) with respect to the three-neutrino model of 0.025. The Bayes factor of FIG. 1. Muon antineutrino disappearance probabilities for the sterile parameters Δm 2 41 ¼ 6.7 eV 2 , sin 2 2θ 24 ¼ 0.33, and two values of g 2 . Top: g 2 ¼ 0, which corresponds to infinite ν 4 lifetime, i.e., the 3 þ 1 model. Bottom: g 2 ¼ 2.5π; this is the best-fit point. the frequentist best-fit point is 0.027. The Bayesian result for g 2 ¼ 2.5π is shown in Fig. 4. As with the frequentist confidence regions, Bayesian preferred regions sliced at the varying values of g 2 are very similar, with a few exceptions. For g 2 ¼ 0, the region log 10 ðBFÞ ≤ −0.5 excludes points with sin 2 2θ 24 ≳ 0.2. The regions with log 10 ðBFÞ ¼ −1.5 only occur for 1.5π ≤ g 2 ≤ 2.5π. This is shown in the Supplemental Material [88].
The frequentist and Bayesian results profiled over the parameters Δm 2 41 and sin 2 2θ 24 are shown in Fig. 5. Both analyses find some preference for nonzero g 2 . In the frequentist analysis, g 2 ¼ 0 is disfavored in favor of nonzero g 2 with −2ΔLLH ¼ 3.9 and a p value of 0.81%. This p value was obtained using Wilks' theorem and 0.26 effective DOF. The effective DOF was determined by fitting 500 simulated datasets generated assuming the best-fit 3 þ 1 parameters, i.e., fixing g 2 ¼ 0. The uncertainty of the value of the effective DOF is 0.02.
The 95% C.L. allowed region found in this Letter overlaps that of the SBL fits, as seen in Fig. 3. This overlap occurs to some extent for all nonzero values of g 2 , but is larger for g 2 values above π. This overlap remains fixed in Δm 2 41 and sin 2 2θ 24 for varying g 2 . At and above g 2 ¼ π, there is some overlap between the 95% C.L. region of this Letter and the 90% C.L. allowed region from the SBL fits.
In conclusion, we have found no substantive evidence for the 3 þ 1 þ decay model. The null hypothesis of only three neutrinos is disfavored with a p value of 2.8%, and the 3 þ 1 model disfavored with a p value of 0.81%. The FIG. 3. The result of the frequentist analysis for g 2 ¼ 2.5π. The 90%, 95%, and 99% C.L. contours are shown as blue dotted, dashed, and solid curves, respectively. The best-fit point is marked with a blue star. The median sensitivity at 99% C.L., determined from 300 simulated datasets, is shown as a red curve. The medium and light pink bands indicate the 1σ and 2σ regions for the sensitivity. The 2D projection of the SBL fit results from [4] for the range 2.25π ≤ g 2 ≤ 2.75π at 90% C.L., 95% C.L., and 99% C.L. are shown as the solid yellow, green, and purple islands around best-fit parameters are Δm 2 41 ¼ 6.7 þ3.9 −2.5 eV 2 , sin 2 2θ 24 ¼ 0.33 þ0.20 −0.17 , and g 2 ¼ 2.5π AE 1.5π. While we have reported valuable new input to global studies, further work, both within and beyond IceCube, is needed to clarify the picture [89]. In particular, in IceCube, the track energy reconstruction can be improved by using machine learning algorithms [90] and the dataset can be expanded to use a new event morphology [91,92].