Observation of Seven Astrophysical Tau Neutrino Candidates with IceCube

In 2013 IceCube discovered a flux of neutrinos of astrophysical origin [1, 2, 3]. The astrophysical neutrino ($\nu^{\mathrm{astro}}$) flux normalization and index $\gamma$ carry information about neutrino sources and their environments [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Different $\nu^{\mathrm{astro}}$ production mechanisms lead to different $\nu_e:\nu_\mu:\nu_\tau$ ratios at the sources but, after standard neutrino oscillations over astrophysical distances, detectable numbers of all three neutrino flavors are expected at Earth [16, 17, 18, 19, 20, 21, 22, 23, 24]. Previous measurements at lower energies, using neutrinos produced at accelerators and in the atmosphere ($\nu^{\mathrm{atm}}$), have detected $\nu_\tau$ produced directly [25] and through neutrino oscillations [26, 27, 28]. At the much higher energies accessible to this analysis, $\nu_\tau^{\mathrm{atm}}$ are strongly suppressed relative to $\nu_\tau^{\mathrm{astro}}$ [29], while an unexpected level of presence of $\nu_\tau^{\mathrm{astro}}$ in the $\nu^{\mathrm{astro}}$ flux could be an indication of new physics [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43].

Previous analyses [44, 45, 46, 47] by IceCube to detect $\nu_\tau^{\mathrm{astro}}$ included searches for double-cascade signatures, such as the distinctive “double bang” [16] in the full detector or “double pulse” (DP) waveforms in one or two individual photosensors. The DP signature is produced by the distinct arrival times of light signals at one or more photosensors from the $\nu_\tau$ interaction and $\tau$ decay vertices. IceCube previously observed two candidate $\nu_\tau^{\mathrm{astro}}$, ruling out the null hypothesis of no $\nu_\tau^{\mathrm{astro}}$ at 2.8$\sigma$ [46]. The analysis presented in this Letter reports on the low-background, high-significance detection of seven $\nu_\tau^{\mathrm{astro}}$ candidate events through the use of convolutional neural networks (CNNs).

IceCube [48] is a neutrino observatory with 5160 Digital Optical Modules (DOMs) on 86 strings [48, 49] in a cubic kilometer of ice at the South Pole. Charged particles produced in neutrino interactions emit Cherenkov light [50] while propagating through the ice; photomultiplier tubes in the DOMs convert this light into electrical pulses that are digitized in situ. Light is deposited in the detector in several distinct patterns: long tracks, single cascades, and double cascades. Tracks are produced by muons from, e.g., $\nu_\mu$ charged-current (CC) interactions, and can start or end inside, or pass through, the detector. Single cascades arise from electromagnetic and/or hadronic particle showers produced by deep inelastic neutrino-nucleon interactions in or near the detector, or by $\bar{\nu}_e$ via the Glashow Resonance [51, 52]. Double cascades are formed by high-energy $\nu_\tau$, CC interactions in or near the detector that produce a hadronic shower and a $\tau$ lepton at the interaction vertex, followed by a second electromagnetic or hadronic shower at the $\tau$ decay vertex (BR[$\tau\to(e,h)$] $\simeq83\%$). With a decay length of $\sim50$ m/PeV, the $\tau$ can travel a macroscopic distance in the ice. For $\nu_\tau$ energies satisfying $E_{\nu_\tau} \gtrsim1$ PeV and favorable geometric containment, the double-bang signature can be created, with two energetic and well-separated cascades. Such events are intrinsically rare. In contrast, for lower $E_{\nu_\tau}$ between roughly 50 TeV – 1 PeV, the $\nu_\tau$ flux is expected to be higher but in CC interactions the two cascades are closer together. Two cascades as close as about 10 m can produce distinctive patterns correlated across multiple DOMs and strings as the light from each cascade passes by, as well as DP waveforms in one or more DOMs.

We analyzed 9.7 years of IceCube data from 2011–2020, triggering on approximately $10^{12}$ downward-going cosmic-ray muons, $10^{6}$ $\nu_{\mathrm{atm}}$, and $10^{4}$ $\nu_{\mathrm{astro}}$ [53, 54, 55, 56]. We required that the DOMs on the most illuminated string collected at least 2000 photoelectrons (p.e.) (see Figure 1) and at least 10 p.e. in the two next-highest-charge, nearest-neighbor strings. Signal events will appear more like cascades than tracks in IceCube, so we also selected events whose morphology was better described by the cascade hypothesis. Aside from 0.6% (22 live-days) of the data sample used to confirm agreement between data and simulation (data that were subsequently excluded from our analysis and which contained no signal-like events), we performed a “blind” analysis that only used simulated data to devise all selection criteria and analysis methodologies. After application of these initial selection criteria, there was roughly 300 times more background than signal. The expected number of p.e. on the most illuminated string for $\nu_{\tau}^{\mathrm{astro}}$ CC events after application of these criteria, and additional CNN criteria described below, is shown in Figure 1.

We then created 2-d images of DOM number (corresponding to depth) vs. time in 3.3 ns bins, with each pixel’s brightness proportional to the digitized waveform amplitude in that time bin. Images were created for the 180 DOMs on the most illuminated string and its two nearest and highest-illuminated neighbors, providing three images per event. The image for the highest-charge string on a candidate signal event is shown in Figure 2 (left). The three images were then processed by CNNs, trained to distinguish images produced by simulated signal and background events and based on VGG16 [57], with a total of $O(100\ \mathrm{M})$ trainable parameters for the high-dimensional signal parameter space.

Three separate CNNs were used to distinguish the $\nu_{\tau}$ signal from remaining backgrounds produced by 1) single cascade neutrino interactions such as $\nu_{e,\mu,\tau}$ neutral current (NC) and $\nu_e$ CC, 2) downward-going muons ($\mu_\downarrow$), and 3) both $\nu_\mu$ interactions producing muon tracks and $\mu_\downarrow$; the associated CNN scores are denoted $C_1$, $C_2$ and $C_3$, respectively, with ranges $[0,1]$. Figure 2 (right) shows $S(C_1)$, the saliency [58] for $C_1$, here defined as the magnitude of the gradient of the CNN score (scaled to $[0,1]$) of $C_1$ with respect to the signal amplitude at each pixel. For reference, the contour (solid line) shows where the detected light falls to zero, and is essentially an outline of the plot on the left. (Points outside the contour are variously causally early, very late, or at distances that are many absorption lengths from the event vertices.) Large $S(C_1)$ values indicate where and when changes in light level most effectively change $C_1$. Small $S(C_1)$ values appear in highly-illuminated regions and in regions with no light. Bright regions contribute to $C_1$, but $C_1$ is not as sensitive there to changes in light level as at the leading and trailing edge envelopes of the light from the event, which are roughly coincident with the contour. The saliency thus shows that $C_1$ is sensitive to the overall shape of emitted light in the detector.

The scores were calculated for each event, and a signal-to-noise ratio of $\sim14$ was obtained by requiring events to have high scores ($C_1 \ge0.99$, $C_2 \ge0.98$ and $C_3 \ge0.85$). The dominant background came from other $\nu_{\mathrm{astro}}$ flavors and $\nu_{\mathrm{atm}}$. The expected energy spectra for signal and the dominant backgrounds, after application of initial and then final selection criteria (including the high CNN scores), are shown in Figure 3.

A sub-dominant “edge event” background was observed from simulated cosmic-ray muons that deposited most of their Cherenkov light on a single string on the outer edge of the detector. We required $C_{3}>0.95$ for edge events, reducing this background by about an order of magnitude at an estimated 15% signal loss. Table 1 lists the expected number of events, after application of the initial and final sets of selection criteria, assuming the best-fit parameters from two IceCube flux measurements.

	$\nu_{\tau,\mathrm{CC}}^{\mathrm{astro}}$ [59]	$\nu_{\mathrm{other}}^{\mathrm{astro}}$ [59]	$\nu_{\mathrm{conv.}}^{\mathrm{atm}}$ [60, 61, 62, 63]	$\nu_{\mathrm{prompt}}^{\mathrm{atm}}$ [56, 64, 65, 66]	$\mu_{\mathrm{conv.}}^{\mathrm{atm}}$ [67, 68, 69, 70]	all background
initial	$160 \pm0.2\ (190 \pm0.3)$	$400 \pm0.7\ (490 \pm0.8)$	$580 \pm7$	$72 \pm0.1$	$8400 \pm110$	$9450 \pm110\ (9540 \pm110)$
final	$6.4 \pm0.02\ (4.0 \pm0.02)$	$0.3 \pm0.02\ (0.2 \pm0.01)$	$0.1 \pm0.008$	$0.1 \pm0.001$	$0.01 \pm0.008$	$0.5 \pm0.02\ (0.4 \pm0.02)$

The largest backgrounds are due to other astrophysical neutrino interactions, and conventional and prompt atmospheric neutrinos, followed by muons from $\pi^{\pm}/K^{\pm}$ decays in cosmic-ray air showers. The backgrounds listed in Table 1 were estimated using simulation packages for astrophysical neutrinos [59], muons from cosmic-ray air showers [67, 68] (with cosmic-ray primary flux given by [69] and hadronic interaction model by [70]), conventional atmospheric neutrino flux from $\pi^{\pm}/K^{\pm}$ decays [60] following our published $\nu^{\mathrm{atm}}$ flux measurements above $E_{\nu}\sim50$ TeV [61, 62, 63], and prompt atmospheric neutrino flux [56, 64, 65, 66] postulated to arise from the decays of charm or heavier mesons produced in air-showers and modeled following Ref. [64]. Electromagnetic (EM) and hadronic showers below 1 PeV were simulated based on the parameterizations of the mean longitudinal and lateral profiles in Ref. [72] and included fluctuations in the energy of the hadronic component. Above 1 PeV the LPM effect [73, 74, 75] is used for EM showers. (Our treatment of possible prompt atmospheric muons is described in the Appendix.) The total $\nu N$ deep inelastic scattering cross section is from [76].

Additional potential background from muon deep inelastic scattering ($\mu$ DIS), given by $\mu+X\rightarrow\nu_{\mu}+X'$, where the light from the incoming $\mu$ followed by the light from the hadronic cascade could mimic the $\nu_{\tau}$ signature, is estimated from the predicted atmospheric $\nu_{\mu}$ CC background. At energies above roughly 100 TeV, we expect comparable numbers of atmospheric $\nu_{\mu}$ and $\mu$ [71], but the $\mu$ energies will be diminished as they pass through the ice to the detector, decreasing their ability to mimic the $\nu_{\tau}$ signature. We conservatively doubled the estimated background from atmospheric $\nu_{\mu}$ CC interactions, from 0.005 to 0.01, to account for the potential background.

We also estimated the background expected from charmed hadrons produced in energetic $\nu_{e}$ CC and $\nu$ NC interactions. This background component had not initially been considered in designing the analysis. After unblinding, we became aware of recent results [77] that indicate that the strange sea in the nucleon is not as suppressed as had been previously believed, so that charm production would thereby be somewhat enhanced compared to our original estimate. Using a simulated neutrino dataset based on the HERAPDF1.5 [78] parton distribution functions (PDFs), and applying a modest correction to reflect more modern PDFs [77, 79, 80, 81, 82], the estimated background from $\nu^{\mathrm{astro}}$ increases by 23% relative to the simulations excluding these interactions. The theoretical uncertainty from the PDFs at the 100 TeV scale is roughly 3%, so the increase corresponds to only about $(15\pm0.5)\%$ (0.08$\pm$0.002 events) of the total background estimation. We included this additional background directly to maintain our blindness protocol that disallowed retraining the CNNs to reject charm background. Uncertainties in the cross section for the interaction of charmed mesons and baryons with ordinary matter had a negligible impact. Backgrounds from on-shell $W$ production [83] from high-energy $\nu_e/\nu_\mu$ interactions, top-quark decay and Glashow resonance interactions [84] can produce energetic $\nu_\tau$ or $\tau$, but are collectively estimated to contribute roughly an order of magnitude fewer background events than other sources and were not included in our background estimate.

For the range of astrophysical neutrino fluxes measured by IceCube (denoted $\phi^{\mathrm{IC}}_{\mathrm{astro}}$), and for a 1:1:1 neutrino flavor ratio at the detector, we predicted a final sample of 4–8 $\nu_\tau$ CC signal events. Similarly, the predicted total background varied for each $\phi^{\mathrm{IC}}_{\mathrm{astro}}$. Using IceCube’s previous measurements of the spectral index $\gamma_{\mathrm{astro}}$, this relatively small number of events constrains the $\nu^{\mathrm{astro}}$ flux normalization $\phi_{\mathrm{astro}}$. Data satisfying $C_2 > 0.98$ were placed in $4\times4$ bins in their $C_1$ and $C_3$ scores.

We calculated confidence intervals following Ref. [85] and using the test statistic defined as $\mathrm{TS}(\lambda_\tau)=\ln L(\hat{\lambda}_\tau)-\ln L(\lambda_\tau)$, where $\lambda_\tau=\phi(\nu^{\mathrm{astro}}_\tau)/\phi_{\mathrm{nom.}}(\nu^{\mathrm{astro}}_\tau)$, the measured-to-nominal flux ratio. Here the nominal flux is one of the four IceCube measured values, and $\hat{\lambda}_\tau$ the value of $\lambda_\tau$ that maximizes the Poisson likelihood $L$ across all 16 bins. Critical values were extracted at the desired confidence level using the TS distributions from a range of $\lambda_\tau$ values, each of which were simulated with $10^4$ pseudo experiments. This procedure incorporated as nuisance parameters the systematic uncertainties in the estimated fluxes for each background component (prior width of 30% for $\nu^{\mathrm{atm}}$ and $\nu^{\mathrm{astro}}$; 50% for cosmic-ray muons), the detection efficiency of the DOMs (10%), and the optical scattering properties of the ice (5%). Since many of these parameters are degenerate in their effect on the analysis observables, and we expected fewer signal events than nuisance parameters, we estimated their impact by incorporating randomized versions of the parameters for each of the pseudo experiments used to calculate the critical value of our TS. This procedure increased the critical values relative to their values in the absence of the systematic uncertainties, widening the extracted confidence intervals.

Seven events remained after applying the final set of selection criteria to the data, consistent with expectation. Figure 4 shows the final expected signal and background, assuming IceCube’s GlobalFit flux, as a function of $C_3$ vs. $C_1$. Five of the candidate $\nu_\tau$ events are in the upper right bin and two are in the bin just below it. Three of the seven events were seen in previous IceCube analyses [1, 46, 47, 86], and one of these three had previously been identified [46, 47] as a candidate $\nu_\tau$. For each candidate event we evaluated the “tauness” as $P_\tau(i)=n_s(i)/(n_s(i)+n_b(i))$, where $n_s$ and $n_b$ are the expected signal and background in bin $i$ (see Figure 4). $P_\tau$ ranges from 0.90–0.92 for two of the candidate $\nu_\tau$, and 0.94–0.95 for the other five, depending on the assumed $\phi^{\mathrm{IC}}_{\mathrm{astro}}$. For IceCube’s GlobalFit flux we predict a total background of $0.5\pm0.02$ events (see Table I); using the distribution of the seven observed events and expected backgrounds in the 16 bins in Figure 4, we exclude the null hypothesis of no $\nu_\tau^{\mathrm{astro}}$ at a (single-sided) significance of $5.1\sigma$. Under the other three flux assumptions [54, 55, 56], the significances are $5.2\sigma$, $5.2\sigma$ and $5.5\sigma$, respectively. The best-fit $\nu_\tau$ flux normalizations are all within the 68% frequentist confidence intervals of the four IceCube fluxes.

We performed multiple checks on the candidate events to ensure they were consistent with expectation. For simplicity and to avoid introducing additional systematic uncertainties, the analysis did not employ a tailored $\nu_\tau^{\mathrm{astro}}$ reconstruction. However, as a post-unblinding check we used a reconstruction for single-cascade events [87] to estimate the energies and directions (Fig. 5) and vertex positions (see the Appendix). The median expected $E_{\nu_\tau}$ was roughly 200 TeV (for the flux in Ref. [56]). The dominant up-down asymmetry is due to Earth absorption and consistent with expectation. Other polar angle effects such as higher vertical vs. horizontal DOM density and $\nu_\tau$ regeneration [88] are also included. (Simulations predict that for $E_\tau\lesssim0.5$ PeV, the selected events are biased toward higher average $\tau$ decay lengths $\langle L_\tau\rangle$; e.g., for $E_\tau\sim100$ TeV, $\langle L_\tau\rangle\sim10$ m.) The events were more clustered in depth than expected but were consistent with a statistical fluctuation, as discussed in the Appendix. We observed no significant coincident activity in the IceTop cosmic-ray air-shower surface array for any of the events. We tested the robustness of the CNNs to hypothetical improperly modeled uncertainties by evaluating their susceptibility to correlated and uncorrelated variations of the raw data underlying the images. We found that the probability of background-to-signal migration was $< (2 \pm0.2)\cdot10^{-5}$ and of signal-to-background migration $< (3 \pm0.8)\cdot10^{-3}$. We also employed targeted tests to estimate the CNNs’ robustness against less likely changes in the underlying raw data, including adversarial attacks [89] against candidate signal and simulated background events. We found that events only migrated in response to changes outside our uncertainty envelope. These tests are described in the Appendix. We conclude that the CNNs are robust against detector systematic effects that could present as either uncorrelated or correlated changes in light levels in one or more DOMs, or entire strings, in the detector.

Energetic astrophysical sources, in conjunction with neutrino oscillations over cosmic baselines, provide the only known way to produce large numbers of $\nu_\tau$ energetic enough to create the observed event morphologies. The result presented in this Letter demonstrates that astrophysical $\nu_\tau$ consistent with this hypothesis are present in the IceCube data and provides powerful confirmation of the earlier IceCube discovery of astrophysical neutrinos [2, 3, 90].

Acknowledgments

The IceCube Collaboration acknowledges the significant contributions to this manuscript from the Pennsylvania State University. We dedicate this paper to the memory of Lovisa Arnesson-Cronhamre, a young graduate student whose tragic and untimely passing stole from us all a promising neutrino physicist. USA – U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, U.S. National Science Foundation-EPSCoR, U.S. National Science Foundation-Major Research Instrumentation Program, Deep Learning for Statistics, Astrophysics, Geoscience, Engineering, Meteorology and Atmospheric Science, Physical Sciences and Psychology (DL-SAGEMAPP) at the Institute for Computational and Data Sciences (ICDS) at the Pennsylvania State University, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin–Madison, Open Science Grid (OSG), Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS), Frontera computing project at the Texas Advanced Computing Center, U.S. Department of Energy-National Energy Research Scientific Computing Center, Particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, and Astroparticle physics computational facility at Marquette University; Belgium – Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany – Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen-Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden – Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; European Union – EGI Advanced Computing for research; Australia – Australian Research Council; Canada – Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark – Villum Fonden, Carlsberg Foundation, and European Commission; New Zealand – Marsden Fund; Japan – Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea – National Research Foundation of Korea (NRF); Switzerland – Swiss National Science Foundation (SNSF); United Kingdom – Department of Physics, University of Oxford.

Appendix

The following appendix includes event displays and saliency maps not shown in the main text, and a more detailed discussion of backgrounds and data-driven studies of the CNN performance.

Event Displays

Images and saliency maps for all seven candidate $\nu_{\tau}^{\mathrm{astro}}$ events are shown in Figure A1.

Candidate Event Spatial Distribution

The seven $\nu_\tau^{\mathrm{astro}}$ candidate events were reconstructed a posteriori using the algorithm described in Ref. [87]. We also applied the same reconstruction to the $\nu_{\tau,\mathrm{CC}}^{\mathrm{astro}}$ signal simulation. Figure A2 shows a top view of the reconstructed vertices for the seven events and Figure A3 shows a side view of the events, superimposed on the expected distribution from simulated signal. The numbers 1–7 next to each data point correspond to the images in Figure A1 (moving left to right and top to bottom).

From the depth distribution, it becomes apparent that the events seem to cluster close to a prominent, 100 m thick, dust layer centered around $z=-80$ m [91, 92]. This region has high optical absorption, reducing detectable light. Photon scattering and absorption through all ice layers and their subsequent effects on the detection efficiency were included in our simulation. A Poisson goodness-of-fit test on the entire 2d histogram in Figure A3 gives a moderate $p$-value of 0.38, based on a suite of pseudo-trials. On the other hand, the projected $z$-distribution indicates a clustering that is inconsistent with expectation at the $3\sigma$ level, according to a Kuiper test [93].

To investigate further, we explored the possibility that the discrepancy might be due to mismodeling of the dust layer or due to atmospheric muons from prompt decays of charm or unflavored vector mesons.

Impact of the Dust Layer

At the energies to which this analysis is sensitive, the selected events have $\tau$ leptons that travel 10–20 m. At roughly 100 m thickness, the dust layer is too thick for it to enable background events to mimic signal events. At our energies, the effect of the dust layer would be to obscure the light from one or both of the $\nu_\tau$ cascades, making the event look instead like a single cascade, or simply too dim, respectively.

Nevertheless, we performed several tests to verify that the CNNs were not unduly sensitive to light signals in the dust layer. For candidate events near the dust layer, shifting the waveform times for all DOMs in the dust layer by $\pm300$ ns, or even removing those DOMs entirely from the event, did not change the CNN response. In another test, the times of individual pixels in the candidate events were shifted, and migration out of the signal region occurred only for shifts exceeding about 100 ns, well in excess of uncertainties expected from mismodeling of ice or DOM properties. (Each pixel holds the PMT charge in a 3.3 ns bin in the DOM’s waveform. The timing uncertainty due to ice properties is estimated to be about 20 ns for distances of about 100 m.) We also altered simulated background events near the dust layer and found that the CNN scores were similarly robust (see next section for more details).

Impact of Prompt Cosmic-Ray Muons

The flux of muons from prompt decays of heavy mesons differs from the conventional flux as a function of both energy and arrival direction. To check whether the unanticipated $z$-distribution might be due to the (un-simulated) prompt muon component of cosmic-ray showers, a dedicated simulation was performed sampling muons from an $E^{-2}$ power-law spectrum and using a parameterized prompt model based on DPMJet-2.55 [94]. The targeted simulation was performed at depth, efficiently sampling muons from the parameterized model. The simulated muon flux was then subjected to an analysis where we loosened the cut on the CNN $C_3$ score, designed to distinguish $\nu_\tau^{\mathrm{astro}}$ signal from muon tracks produced by $\nu_\mu$ interactions and downward-going muons, from 0.85 (0.95 with charge asymmetry requirements for outer strings) to 0.75. If the candidate data events were contaminated by muons, loosening $C_3$ as described would result in more simulated muons entering the signal region. In this context, we performed a comparison of the spatial distribution across the detector of the generated muon events and that expected for data. Since all of the simulated events that survive this CNN selection reconstruct as not fully contained, either near the outer edge or top of the detector, for this comparison we excluded the three candidate $\nu_\tau^{\mathrm{astro}}$ events that were well-contained within the fiducial region of IceCube. (We note that one of the remaining four candidate $\nu_\tau^{\mathrm{astro}}$ events in this subsample reconstructs as having passed through hundreds of meters of active detector volume above its interaction vertex, and the absence of appreciable light in this region further supports its neutrino hypothesis, but in what follows we do not incorporate this information.)

Figure A4 shows a scatter plot of $z$ vs. $\rho\equiv\sqrt{x^2+y^2}$ (as measured from the center of IceCube at a depth of 1950 m) for data and unweighted simulated events classified as not fully contained. It is apparent that the additional data events allowed in by loosening the $C_3$ score mainly appear at the top of the detector. Comparing these data to the simulation of cosmic-ray muons, a Kolmogorov-Smirnov consistency test [95] in $z$ gives a $p$-value of 0.1 for the hypothesis that all of the data events are muons. Moreover, reverting to the original definition of the $C_3$ selection criterion leaves only eight unweighted simulated cosmic-ray muons, along with the four candidate edge events. Performing again a Kolmogorov-Smirnov consistency test in $z$ gives a $p$-value of 0.004, indicating that it is unlikely that all of the four candidate $\nu_\tau^{\mathrm{astro}}$ edge events are muons.

Muons from cosmic-ray air showers enter primarily at the top of the detector. Muon rates were estimated from measurements of the cosmic-ray muon flux by IceCube and other instruments [53, 97, 98, 99]. These measurements are compatible with muons solely from $\pi^{\pm}/K^{\pm}$ decays, although upper bounds on the prompt contribution exist. We therefore did not include muons from charm decays in the background simulations used here. For the sake of completeness, however, we calculated the effect of arbitrarily increasing the atmospheric muon background by an order of magnitude with respect to the originally expected value of $0.005 \pm0.004$ (post-unblinding, using the sampling flux from Ref. [69]). We find that the significance remains above $5\sigma$, even with such an increase.

Impact of Relaxing CNN scores C1 and C2

Finally, we investigated the effects of less strict requirements on CNN scores $C_2$ and $C_1$. Loosening just the $C_2$ requirement only lets in events near the top of the detector, and therefore cannot provide an explanation for the observed $z$-clustering. Figure A5 shows a scatter plot of $z$ vs. $\rho\equiv\sqrt{x^2+y^2}$, demonstrating the effect of loosening just the $C_1$ score criterion to $C_1 > 0.90$ to determine if a larger population of events is concentrated near the dust layer, the top of the detector, or at the detector perimeter, as would be expected for an enhanced background. The additional events (red squares) instead broaden the spatial distribution, consistent with a statistical explanation for the originally observed clustering. From simulations, we expected that loosening $C_1$ would yield 9.4 signal and 2.9 background events, for a total of 12.3 events (assuming the IceCube GlobalFit [53] flux), consistent with the 12 events observed. The additional five events have an average “tauness” $\langle P_{\tau} \rangle= 0.49$. The 12 events also exclude the null hypothesis at approximately $5\sigma$.

CNN Robustness

Pre-Unblinding Data vs. Simulation Agreement

Prior to unblinding we investigated the agreement in the CNN scores between data and simulation in the regions populated by background events. Figure A6 shows a cumulative plot of the number of data and expected signal and background events vs. the CNN score $C_1$; CNN scores $C_{2,3}$ show similar levels of agreement.

Post-Unblinding Tests

Here we describe in more detail the various data-driven and simulation-based tests of CNN robustness that we performed. For the first suite of tests, we define the background region as $C_{1,2,3} < (0.9, 0.9, 0.75)$, comprising 8,175 of the original 8,188 events passing the preliminary selection criteria. We applied randomized scale factors to DOM waveforms that artificially increased or decreased the magnitude of the detected light level within expected systematic uncertainties, in five distinct patterns:

• each of 180 DOMs were randomly scaled independent of one another,

• dividing the detector into regions in depth, the group of DOMs in each region was randomly scaled by the same factor, as follows:

  • 20 groups of 9 DOMs each in regions in depth of about 50 m,

  • 15 groups of 12 DOMs each in regions in depth of about 68 m,

  • 12 groups of 15 DOMs each in regions in depth of about 85 m, and

• on each of the two less-illuminated strings, all 60 DOMs were randomly scaled by the same factor (a total of two distinct factors were used).

The first pattern addresses relative DOM detection efficiencies [48], the next three address ice optical properties as a function of depth, and the fifth addresses ice birefringence [100] as a function of azimuthal angle.

For each pattern, we performed 750 trials per event, for a total of $(5 \times8175 \times750)$ or about $3 \cdot10^7$ trials. We found that the probability of background-to-signal migration did not exceed $(2 \pm0.2) \cdot10^{-5}$ in any of the tests, corresponding to $0.16 \pm0.02$ events, and that migration occurred only when events were already close to the signal region. The final significance remains above $5\sigma$ whether we use simulated data sets (described earlier) or this data-driven approach to handle these detector systematics. The same tests performed on the seven signal events showed a signal-to-background migration probability of $(3 \pm0.8) \cdot10^{-3}$.

We also employed targeted tests to estimate the CNNs’ robustness against less likely changes in the underlying raw data. In candidate events with prominent double pulse waveforms, we interpolated between the two peaks to merge together the first and second pulses, and found that this did not cause any candidate events to migrate out of the signal region. As mentioned earlier, shifting pixel arrival times in individual DOM waveforms in candidate events by up to 100 ns did not appreciably change the CNN response.

We applied adversarial attacks [89] against the candidate events, using an optimization algorithm to find the pixel(s) whose physically reasonable alterations resulted in the largest changes to the CNN scores. Just one of the seven candidate events could be forced to migrate, and only when the average change over all pixels was at least 2.5%, a situation that is well outside our estimated uncertainties. Similarly attacking simulated background events, we found that in no particular region of the detector did the CNNs exhibit heightened susceptibility, and generally the changes required to induce migration were much larger than allowed by our uncertainties. We also attacked 634 simulated astrophysical $\nu_e$, allowing the individual pixel uncertainties to be as high as 10%, finding in this harsh test that only one simulated $\nu_e$ was misclassified as a $\nu_\tau$. Finally, we attacked the candidate events after randomly varying their pixel values with 10% uncertainty. Using $10^4$ trials per event, only one event was found to have a (2.1 $\pm$ 0.14)% migration probability. These targeted tests, and other studies described earlier in this Supplemental Material, indicate that even under quite harsh conditions, the CNNs remained capable of rejecting the background events while retaining the candidate signal events.