Abstract High-energy tau neutrinos are rarely produced in atmospheric cosmic-ray showers or at cosmic particle accelerators, but are expected to emerge during neutrino propagation over cosmic distances due to flavor mixing. When high energy tau neutrinos interact inside the IceCube detector, two spatially separated energy depositions may be resolved, the first from the charged current interaction and the second from the tau lepton decay. We report a novel analysis of 7.5 years of IceCube data that identifies two candidate tau neutrinos among the 60 "High-Energy Starting Events" (HESE) collected during that period. The HESE sample offers high purity, all-sky sensitivity, and distinct observational signatures for each neutrino flavor, enabling a new measurement of the flavor composition. The measured astrophysical neutrino flavor composition is consistent with expectations, and an astrophysical tau neutrino flux is indicated at 2.8σ significance.

The discovery of a diffuse flux of astrophysical neutrinos, using High-Energy Starting Events (HESE) observed by Ice-Cube [1] opened the possibility to study the Universe's most powerful cosmic accelerators [2,3]. HESE is an all-flavor, all-sky selection of events of predominantly astrophysical origin, with an analysis region above 60 TeV in deposited electromagnetic-equivalent energy in the detector. Tau neutrinos are expected to be produced only in tiny fractions at neutrino sources, but emerge due to neutrino oscillations over cosmic baselines [4]. For neutrinos from distant sources, the probability of a neutrino created with flavor ν α to reach the detector as ν β is P ν α →ν β = i |U αi | 2 |U βi | 2 [5,6]. Thus, the neutrino flavor composition at Earth depends on the neutrino mixing matrix elements, U αi , and the source flavor composition. For neutrinos from the decay of charged pions produced in hadronic interactions, with a source flavor composition of ν e :ν μ :ν τ = 1/3:2/3:0, we expect ν e :ν μ :ν τ = 0.30:0.36:0.34 at Earth (using the oscillation parameters from [7]), i.e., very close to equipartition (1/3:1/3:1/3). However, the environment at the neutrino production sites may influence the flavor composition, due to cooling or interactions of the charged particles produced in the hadronic interactions [8][9][10][11]. Therefore, the flavor composition of astro-physical neutrinos is a powerful probe of the environments of cosmic accelerators and can help constrain the source populations contributing to the observed neutrino flux. The neutrino flavor composition on Earth is also a sensitive probe of physics beyond the Standard Model (BSM) affecting neutrino propagation and modifying the flavor composition [12][13][14][15][16]; see [17] for BSM-constraints derived using the HESE selection.

Atmospheric neutrinos are a background to astrophysical neutrino searches. As atmospheric neutrinos are accompanied by muons born in the same cosmic-ray-induced shower, their contribution to a sample can be suppressed by muonrejecting event selection criteria, e.g. by using the outer parts of the detector as a vetoing region. This effect, called atmospheric neutrino self-veto [18], is used in HESE [19]. Conventional atmospheric neutrinos are ν e and ν μ from the decay of π ± and K 0,± produced in the atmosphere by cosmic-ray interactions. At energies above ∼ 100 TeV, the atmospheric flux is expected to be increasingly dominated by the prompt component, originating from the decays of charmed hadrons (e.g. [20]). Tau neutrinos, produced from rare decays of D s and D 0,± , contribute only up to 5% to the yet unobserved prompt atmospheric neutrino component [21,22]. This makes the observation of high-energy tau neutrinos a smoking-gun signature of cosmic neutrinos, but so far, none have been identified [23][24][25]. Previous flavor studies only separated the charged-current ν μ contribution from other flavors, leading to a significant degeneracy between the ν e and ν τ flavors [26][27][28]. Here, we present a new flavor composition measurement of astrophysical neutrinos with direct sensitivity to each of the neutrino flavors, performed on the HESE sample. A detailed description of the characteristics of the HESE sample and spectral fits to a diffuse astrophysical neutrino spectrum assuming flavor equipartition, as well as a detailed description of systematic uncertainties and their treatment are provided in [19]. There, the astrophysical neutrino spectrum was fit as a single power law,

where φ ν is the all-flavor ν + ν flux at E 0 = 100 TeV and γ astro is the spectral index. Their best fit values are φ ν = 6.4 +1.5 -1.6 × 10 -18 GeV -1 cm -2 s -1 sr -1 , and γ astro = -2.87 +0.21 -0.19 . The sample and associated results have been made available publicly through a dedicated data release [29].

This manuscript is structured as follows: Sect. 2 describes the signatures of neutrino interactions detected in IceCube and how they map to neutrino flavors; Sect. 3 illustrates the selection and classification of the detected events according to these various signatures; Sect. 4 summarizes the outcome of the classification and the characteristics of the found ν τ candidates; in Sect. 5 the flavor composition constraints from this analysis are derived.

In IceCube, neutrinos are detected by collecting the Cherenkov light emitted by charged secondary particles created in neutrino interactions. All neutral-current (NC) interactions produce showers of hadrons and are indistinguishable between flavors. In a charged-current (CC) interaction, the neutrino flavor can be inferred from the distinct Cherenkov light pattern produced by each flavor of charged lepton. Light depositions from a muon traversing the detector are called tracks and stem from ν μ CC interactions, atmospheric muons, and ν τ CC interactions where the tau decays to a muon (17% branching ratio). Single cascades consist of energy depositions at a single vertex and are produced by ν e CC and NC interactions of all flavors. At PeV energies, both tracks and single cascades can also emerge from the decays of W-bosons produced in resonant neutrino-electron scattering [30]. Double cascades are two energy depositions connected by a track of comparatively low light emission. They are produced by ν τ CC interactions where the first cascade originates from the hadronic interaction of the ν τ producing a tau, and the second cascade stems from the tau decaying to a hadronic or electromagnetic cascade (83% branching ratio) [4]. Due to their short livetime, taus have a short decay length of L τ ∼ 50 m • E τ / PeV, where E τ is the tau energy. This makes the distinction between single and double cascades challenging in IceCube, where the mean horizontal distance between light sensors, called Digital Optical Modules (DOMs), is 125 m. The HESE analysis defines a lower threshold on the deposited electromagnetic-equivalent energy in the detector of events, E tot , of 60 TeV (see below). Above this threshold it is possible to identify some of the ν τ events as double cascades, if L τ 10 m, breaking the degeneracy between ν e and ν τ flavors present at lower energies. 1

Using the HESE selection, we have performed a new analysis of the IceCube data that incorporates major improvements with respect to previous publications [2,3] in our understanding of the detector and the modelling of atmospheric backgrounds. The HESE selection is described in [2]. To pass, an event has to (1) start inside of the outermost layer of DOMs making up the "veto" layer, and (2) deposit more than 6000 photoelectrons in the detector. Muons radiate away Fig. 1 Length resolution for simulated tau neutrinos classified as double cascades for the best-fit spectrum [19]. With increasing length, the resolution improves while the event expectation drops. The inlay shows a schematic of a ν τ CC interaction producing a double cascade with associated reconstruction parameters energy throughout their passage through the ice, with the amount of light deposited increasing with increasing muon energy. It is thus extremely unlikely for atmospheric muons to pass the HESE selection criteria. Due to the atmospheric self-veto ( [18], see also Sect. 1), accompanying muons also greatly reduce the number of downgoing atmospheric neutrinos present in the sample. To further enhance the fraction of astrophysical neutrinos in the sample, the analysis is restricted to events with a reconstructed total deposited energy E tot above 60 TeV. Data collected between 2010 to 2017 using the original HESE selection [2], with a total livetime of 2635 days, have been reprocessed using a new and improved detector calibration. An improved model of the optical properties of the South Pole ice sheet [32], critical to the reconstruction of event properties, has been incorporated into the simulation and reconstruction, and an updated calculation of the atmospheric neutrino self-veto [33] is used. This new HESE sample has 60 events in the analysis region, i.e. with E tot > 60 TeV, and is described in detail in [19].

We use a classification algorithm developed on Monte Carlo (MC) simulated events and first applied to the 6-year HESE sample [25,34] to classify the 60 events as single cascades, double cascades, or tracks (ternary event classification). It was developed with the goal of achieving a high ν τ purity for the events assigned a double-cascade topology, while keeping misclassification fractions low for all topologies [35]. All events are reconstructed using maximum likelihood fits with different hypotheses: single cascade [36], track [37], and double cascade [36,38]. For the fits, the timing and spatial information of the light collected in an event is used. The parameters of the double-cascade fit are (see also A preselection removes events with a failed doublecascade fit from being further considered as double cascades. After preselection, three event properties are used to classify each event: double-cascade length, energy asymmetry, and energy confinement. The double-cascade length is a proxy for the tau lepton's decay length with an average resolution of ∼ 2 m over the entire analysis range at the best-fit spectrum with flavor equipartition [19]. Figure 1 shows the reconstructed double cascade length as a function of the true double cascade length; the length resolution improves with increasing length as the cascades get better separated. The energy asymmetry is defined as

It can take values -1 ≤ A E ≤ 1, with the boundary values corresponding to single cascades. The energy confinement is defined as E C = (E C1 + E C2 )/E tot , where E Ci are the energy depositions within 40 m of the i-th cascade vertex position. For the purpose of this calculation the vertices of the two cascades are taken directly from the doublecascade reconstruction, while the energy depositions are obtained through a track energy unfolding algorithm [36], and thus the confinement can take values 0 < E C < 1. Thus, for double cascades separated by 80 m the relation

Events passing the requirements shown in the second column of Table 1 are classified as double cascades.

True single cascades typically have a small reconstructed double-cascade length and a large, positive energy asymmetry. True tracks typically have energy depositions all along their tracks, leading to low energy confinement values. True double cascades have E C values very close to 1 even for separation lengths in excess of 80 m, due to the low relative brightness of the tau. By choosing a conservative requirement of E C > 0.99 for double cascades, the performed analysis does not lose sensitivity towards higher-energy ν τ producing longer-lived τ leptons. True double cascades show a flat dis-tribution in A E with a resolution of ∼ 0.1 at negative values of A E and worsening towards positive values. Their doublecascade length is correlated to their total deposited energy and follows the exponentially falling distribution seen in the energy spectrum. Events failing the double-cascade requirements are classified according to the procedure shown in the last column of Table 1.

Note that the requirement of L dc ≥ 10 m for double cascades leads to the majority of ν τ induced events to be classified as single cascades. At the best-fit spectrum with flavor equipartition [19], we expect ∼ 15 ν τ events, of which ∼ 12 interact via the double cascade channel. But only ∼ 3 (22.7%) of those are expected to produce a tau that travels at least 10 m before decaying. 42.3% of simulated double cascades with tau decay lengths above 10 m pass the double cascade requirements in Table 1. The total efficiency of the ternary topological classification chain for double cascades is 12.2%. 1.9% of all ν e and ν μ induced events are expected to be misclassified as double cascades.

Glacial ice at South pole flows at a rate of ∼ 10 m per year. It has recently been observed [39] that the optical properties of glacial ice at South Pole vary as a function of the direction with respect to the flow of the glacial ice. This ice anisotropy is one of the limiting factors on the selection of double cascades: directional distortions of Cherenkov light patterns can lead to a misclassification of single cascades as double cascades. See Appendix C for details on the ice anisotropy treatment.

The 60 events are classified into 41 single cascades, 17 tracks, and 2 double cascades. These are the first double cascades in the signal region and the first astrophysical tau neutrino candidate events. The reconstructed properties of the double cascades are shown in Table 2. As the average tau decay length scales with the tau energy L τ ∼ 50 m • E τ / PeV, which depends on the energy of the incoming ν τ as E τ ∼ 0.7• E ν τ , the double cascades length L dc scales with the total deposited energy E tot . Two-dimensional MC probability distribution functions (PDFs) of reconstructed total deposited energy E tot versus reconstructed double cascade length L dc for signal and background contributions to events classified as double cascades are shown in Fig. 2 with the data events overlaid as white circles. For ν τ -induced double cascades (top panel), a clear correlation between E tot and L dc can be seen. Background events (bottom panel) cluster at the thresholds in E tot due to the falling spectrum and in L dc since single cascades typically have very small reconstructed L dc . The regions containing 68%, 90%, and 95% of true single cascades misclassified as double cascades are marked by vertical white lines, i.e. 68% of the true single cascades misclassified as double cascades have L dc < 14.4 m, while 90% have L dc < 20.4 m. The tilted white lines show the region within which 95% of the signal are contained. Few events are expected in the parameter space of event #1, while there are contributions expected from both signal and background in the parameter space of event #2. For single cascades and tracks, the properties total deposited energy, E tot , and cosine of the zenith angle in detector coordinates, cos(θ z ), are used to distinguish atmospheric and astrophysical contributions. The PDFs shown in Fig. 2 and the corresponding PDFs for single cascades and tracks described above are used in the all-flavor analyses presented in [19].

An event view of event #1, observed in 2012 and nicknamed "Big Bird" [3], is shown in Fig. 3. For several DOMs, the photon counts as a function of time are displayed alongside the predicted photon count distributions for single-and double-cascade hypotheses. The double-cascade hypothesis fits the observed data better than the single-cascade hypothesis. However, this event has several saturated and bright DOMs that were excluded from the analysis, a standard procedure for high-energy IceCube analyses [40,41]. A DOM is called saturated if the signal in the PMT exceeds the dynamic range of the readout electronics. A DOM is called bright if it has collected ten times more light than the average DOM for an event. Only statistical uncertainties on photon count rates are included in the likelihoods of the reconstruction algorithms [36][37][38]. At the highest observed energies, bright DOM signals have very small statistic uncertainties and can In the signal (ν τ -induced double-cascade events) histogram (top), the region containing 95% of the expected signal is indicated with white dotted lines. In the background (all remaining events) histogram (bottom), the white vertical dotted lines mark the regions containing 68%, 90%, and 95% of the single-cascade induced background. In both histograms the two tau neutrino candidates are overlaid as white circles therefore lead to misreconstructions due to the lack of proper systematic uncertainty terms in the likelihood. For comparison of predicted photon counts for each hypothesis, the bright DOMs are displayed in Fig. 3. An event view of event #2, observed in 2014 and nicknamed "Double Double," is shown in Fig. 4. The two vertices of the cascades cannot be spatially resolved by eye, highlighting the need for the algorithmic topological classification employed in this work. Analogous to Fig. 3, collected photon counts as a function of time are displayed together with the predicted photon count distributions for single-and doublecascade hypotheses. The predicted photon count PDFs differ remarkably between the single-and double-cascade hypothesis, with the single-cascade hypothesis disfavored. Photoelectrons bright DOM bright DOM bright DOM Single, no brights Double, no brights Double, with brights Exp. Data Time Fig. 4 Double-cascade event #2 (2014). The reconstructed doublecascade vertex positions are indicated as grey circles, the direction indicated with a grey arrow. The size of the circles illustrates the relative deposited energy, the color encodes relative time (from red to blue). Bright DOMs are excluded from this analysis

Data from DOMs labeled as bright were excluded from the analysis , but are used for the comparison of predicted photon count PDFs in Fig. 4.

Figure 5 shows the distribution of the ratio of the doublecascade length L dc to reconstructed decay-cascade energy E 2 (top panel) and the energy asymmetry A E (bottom panel) of simulated events and data for the best-fit spectrum given in [19]. The distributions were not part of the topological classification chain. While the correlation between L dc and E tot is clear on average, there are large fluctuations in energy transfer from parent to daughter particle. Therefore, on the perevent basis, the more direct correlation between the doublecascade length L dc and the decay-cascade energy E 2 proves more informative. Event #1 has a length-to-energy ratio in a region where the ν τ contribution is larger than the background contribution, but outside of 90% of the simulated ν τinduced double cascades. Its high energy asymmetry is in a region with a background expectation which is on the order of the signal expectation. Event #2 has a length-to-energy ratio at the peak of the distribution for ν τ -induced double cascades and an energy asymmetry value in a highly signaldominated region. None of the classified double cascades are in a phase space greatly affected by the ice anisotropy.

To quantify the compatibility with a background hypothesis (i.e., not ν τ -induced) for the actual ν τ candidate events observed, a targeted MC simulation for each event was performed, consisting of simulation of ν e , ν μ , and ν τ interactions. In addition, for "Double Double," also atmospheric muons were simulated. However, none of the 1.2 • 10 10 generated muons passed the HESE veto undetected. See Table 3 Table 3 Parameter space for resimulated events. The upper value of the primary energy depends on the interaction type, reflecting the spread of visible energy losses typical of that interaction. The visible energy is the energy transformed into light, it equals the total energy deposited in the detector for electromagnetic showers and is lower for hadronic showers and events with final-state muons or neutrinos. rr evt is the two-dimensional distance in the x, y-plane.

The values in parentheses are for ν μ CC events Variable Event #1 Event #2 Primary energy > 1.5 PeV > 65 TeV Visible energy 1-3 PeV 60-300 TeV Vertex, rr evt 50 m 50 m Vertex, z -

for details on the restricted parameter space and Appendix A for a description of how this parameter space was chosen. Using targeted MC simulation for the analysis of exceptional events is a method often employed in IceCube [3,30,42,43]. These new MC events were filtered and reconstructed in the same way as the initial MC and data events. In total, ∼ 2•10 7 "Double-Double"-like events and ∼ 1 • 10 6 "Big-Bird"-like events from the targeted simulation pass the HESE selection criteria. A breakdown of simulated event types and their fractions passing the HESE double cascade selection criteria can be found in Appendix A. We define the tauness, P τ , as the posterior probability for each event to have originated from a ν τ interaction, which can be obtained with Bayes' theorem:

(

In the first line we have simply split the total probability of an event at the observed parameter space η evt into its ν τ and non-ν τ (written ν τ ) components in Bayes' theorem. In the second line we identify P(η evt |ν τ ) with the PDFs for ν τ , and express the prior probability P(ν τ ) = N ν τ /(N ν τ + N ν τ ) as the fraction of expected ν τ events evaluated at the observed parameter space of each event, η evt , obtaining the differential number of expected events N ν τ P ν τ (η evt ) (and analogous for the non-ν τ components indicated as ν τ ). For each tau neutrino candidate, the differential expected number of events at the point η evt , N ν τ P ν τ (η evt ) and N ν τ P ν τ (η evt ) is approximated from the targeted simulation sets using a multidimensional kernel density estimator (KDE) with a gaussian kernel and the Regularization Of Derivative Expectation Operator (rodeo) algorithm [44]. The rodeo algorithm provides an unbiased and computationally efficient way to find the optimal bandwidth in d dimensions for a d-dimensional set of n events. In the rodeo the bandwidth is reduced as long a the derivative of the kernel density estimate with respect to its bandwidth is large compared to its variance. The obtained optimal bandwidth for each considered dimension balances the relevance of the variable with the sparsity of the dataset at the evaluated point. The eight dimensions used in evaluating the tauness include the six dimensions (d = 6) of the restricted parameter space that the resimulation was carried out in: total deposited energy E tot , vertex position (x, y, z) and direction (θ, φ). Further, a region of interest is defined in the parameters not restricted during resimulation but used in the double-cascade classification: double-cascade length L dc and energy asymmetry A E [45]. The region of interest is obtained by slowly decreasing a two-dimensional box around the observed parameters as long as the statistical errors from the limited targeted MC stay below 10%. This procedure was established using the produced MC in a sideband region.

Having defined η evt = (E tot , x, y, z, θ, φ, L dc , A E ), and approximating

and

one obtains the tauness

Here, fν α (η evt , ĥν α ) is the density of ν α for the optimal bandwidth ĥν α determined by the rodeo algorithm in the region of interest. Originally developed for unweighted events, we extend the rodeo formalism to weighted events according to the procedure in [46]: Each of the n simulated events has a weight w i , with i = 1, . . . , n. We use the effective number of events n Eff = ( i w i ) 2 / i (w 2 i ), and their effective weight w Eff = i w 2 i / i w i . Note that the tauness is always evaluated under certain assumptions for the flux parameters. Computing the tauness for each of the events to originate from a ν τ interaction for the best-fit spectrum given in [19] with a 1/3:1/3:1/3 flavor composition yields P BB τ best fit ≈ 75% for "Big Bird," and P DD τ best fit 97% for "Double Double." For "Double Double," the statistics of the generated MC are not sufficient to evaluate the tauness to a higher precision. The tauness weakly depends on the astrophysical spectral index and decreases by ∼ 1% for a softening of γ astro by one unit.

We sample the posterior probability in the flavor composition, obtained by leaving the source flavor composition unconstrained and taking the uncertainties in the neutrino mixing parameters into account. When using the best-fit spectra given in [19] but varying the source flavor composition over the entire parameter space (i.e. ν e :ν μ :ν τ = a:b:1ab with 0 ≤ a, b ≤ 1 and a + b ≤ 1 at source), and the mixing parameters in the global fit NuFit4.1 [7] 3σ allowed range, the tauness is (97.5 +0. 3 -0.6 )% for "Double Double" and (76 +5 -7 )% for "Big Bird." "Double Double" is also identified as a candidate tau neutrino event in two complementary analyses using the double pulse method to search for tau neutrinos that have been performed while this analysis was ongoing [47,48].

A multi-component maximum likelihood fit is performed on the three topological subsamples using PDFs obtained from MC simulations. We account for the uncertainty due to limited MC statistics by using a variant of the effective likelihood L Eff , a generalized Poisson likelihood, presented in [46] and employed in [19]. This joint likelihood is composed of the contributions from the independent subsamples single cascades, double cascades, and tracks (SC, DC, and T, respectively):

where θ are the model parameters, j are the analysis bins, μ j is the expected number of events and the variance in the j-th bin with statistical uncertainty σ j , d j is the observed number of events in the j-th bin, and t = (SC, DC, T) are the event topologies. Each simulated event i has a weight w i which depends on the model parameters θ. The expected number of events is a product of the effective number of simulated events n Eff and the effective weight, w Eff introduced in Sect. 4.2:

For all topologies, the contributions from atmospheric and astrophysical neutrinos as well as atmospheric muons are taken into account in the likelihood analysis. The conventional atmospheric neutrino component is modeled according to the HKKMS calculation [49,50], the prompt atmospheric neutrino component is modeled following the BERSS [20] (for ν e , ν μ ) and MCEq [22] (for ν τ ) calculations. MCEq is using the SIBYLL-2.3c [51] model. The muon component is simulated using MUONGUN [52] which samples single muons from templates generated by CORSIKA [53] weighted to the Hillas-Gaisser-H4a cosmic-ray model [54] and employing the SIBYLL-2.1 hadronic interaction model [55] in the shower development. For the spectrum of the astrophysical neutrino flux Φ astro , a single power law with a common spectral index γ astro for all flavors is used,

where φ ν α is the astrophysical normalization of the ν +ν flux of flavor α at E 0 = 100 TeV. While for single cascades and tracks, atmospheric contributions pose the main background to the astrophysical signal, the main background to ν τ -induced double cascades arises from misclassified astrophysical ν e and ν μ . The background contributions from atmospheric neutrinos are small (0.2 events in 7.5 years expected), while those from penetrating atmospheric muons and prompt atmospheric ν τ are negligible (0 and 0.04 events in 7.5 years expected, respectively).

The systematic uncertainties are given in Table 5 found in Appendix C (reproduced from Table V of [19]), and are included in this analysis in the same way as in [19]. The main systematic uncertainty affecting the double-cascade reconstruction is the anisotropy of the light propagation in the ice [32,56].

While in [19,[57][58][59], the total likelihood is maximized assuming flavor equipartition, here we fit the three flavors' fractions f α of the overall astrophysical normalization Φ astro ,

To perform the flavor composition measurement using the multidimensional KDE, the likelihood is modified compared to the analyses in [19]. In the joint likelihood for the three topologies, L Eff = L SC Eff L T Eff L DC Eff [19], L DC Eff is replaced by the extended unbinned likelihood for the double-cascade events,

where c are the flux components used in the fit, c = ν astro,α , ν conv,α , ν prompt,α , μ atm for the flavors α = e, μ, τ . N c P c (η evt ) is computed using the rodeo algorithm introduced in Sect. 4.2. The aforementioned slight dependence on γ astro is parametrized in the extended double-cascade likelihood L DC Rodeo by evaluating N c P c (η evt ) as a function of γ astro .

The result of the flavor composition measurement is shown in Fig. 6. The fit yields

with a best-fit flavor composition of ν e :ν μ :ν τ = 0.20:0.39:0.42. (10) Comparing this result with previously published results of the flavor composition also shown in Fig. 6 clearly shows the advantages of the ternary topological classification. The best-fit point is non-zero in all flavor components for the first Fig. 6 Measured flavor composition of IceCube HESE events with ternary topology ID and extended multi-dimensional analysis of the double cascades (black). Contours show the 1σ and 2σ confidence intervals assuming Wilks' theorem [60] holds. The shaded regions show previously published results [27,61] without direct sensitivity to the tau neutrino component. Flavor compositions at source and after propagation expected from various astrophysical neutrino production mechanisms (see, e.g., [9]) are marked, and the entire accessible range of flavor compositions assuming standard 3-flavor mixing is shown time, and the degeneracy between the ν e and ν τ fraction is broken. The small sample size of 60 events in this analysis and the lower sensitivity of the HESE sample to ν μ than to ν e and ν τ flavors both lead to an increased uncertainty on the ν μ fraction as compared to [27,61].

The test statistic TS = -2 ln L(φ 0 ν τ )ln L(φ b.f. ν τ ) compares the likelihood of a fit with a ν τ flux normalization fixed at a value φ 0 ν τ to the free fit where φ ν τ assumes its best-fit value, φ b.f. ν τ . Evaluated at φ 0 ν τ = 0 and using Wilks' theorem, it gives the significance at which a vanishing astrophysical tau neutrino flux can be disfavored. The test statistic is expected to follow a half-χ 2 k distribution with k = 1 degree of freedom [62]. The validity of Wilks' theorem was tested with pseudo-MC trials as described in Appendix B. The observed test statistic is TS = 6.5, which translates to a significance of 2.8σ , or a p-value of 0.005. A one-dimensional scan of the astrophysical ν τ flux normalization is performed with all other components of the fit profiled over. The 1σ confidence intervals are defined by TS ≤ 1, and the astrophysical tau neutrino flux normalization is measured to

This constitutes the first indication for tau neutrinos in the astrophysical neutrino flux.

Seven and a half years of HESE events were analyzed with new analysis tools. The previously shown data set was reprocessed with improved detector calibration. A flavor composition measurement was performed using a ternary topological classification directly sensitive to tau neutrinos, which breaks the degeneracy between ν e and ν τ events that is present in a binary classification scheme (into tracks and cascades). This analysis found the first two double cascades, indicative of ν τ interactions, with an expectation of 1.5 ν τ -induced signal events and 0.8 ν e,μ -induced background events for the bestfit single-power-law spectrum with flavor equipartition [19].

The first event, "Big Bird," has an energy asymmetry at the boundary of the selected interval for double cascades. For the second event, "Double Double," the photon arrival pattern is well described with a double-cascade hypothesis, but not with a single-cascade hypothesis. A dedicated a posteriori analysis was performed to determine the compatibility of each of the events with a background hypothesis, based on targeted MC. The analysis confirms the compatibility of "Big Bird" with a single cascade, induced by a ν e interaction, at the 25% level. A "Big Bird"-like event is ∼ 3 (15) times more likely to be induced by a ν τ than a ν e (ν μ ), the result being only weakly dependent on the astrophysical spectral index. "Double Double" is ∼ 80 times more likely to be induced by a ν τ than either a ν e or a ν μ . All background interactions have a combined probability of ∼ 2%, almost independent of the spectral index of the astrophysical neutrino flux. Using a novel extended likelihood for double cascades, which allows for the incorporation of a multi-dimensional PDF as evaluated by a kernel density estimator, the flavor composition was measured. The best fit is ν e :ν μ :ν τ = 0.20:0.39:0.42, consistent with all previously published results by IceCube [27,61], as well as with the expectation for astrophysical neutrinos assuming standard 3-flavor mixing. The astrophysical tau neutrino flux is measured to:

A zero ν τ flux is disfavored with a significance of 2.8σ , or, p = 0.005. A limitation of the analysis presented here is the small sample size of 60 events. Merging the HESE selection with the contained cascades event selection [40] is expected to enhance the number of identifiable ν τ events by ∼ 40% [63]. Due to the small effective volume for ν μ -CC interactions of HESE, the ν μ fraction of the astrophysical neutrinos has large uncertainties. Work on updating the joint analysis of multiple event selections [27] is ongoing, where the strongest contribution to constraining the ν μ fraction is expected from through-going muons [41,64]. A few years from now, the Ice-Cube Upgrade [65] will greatly improve our knowledge and modeling of the optical properties of the South Pole ice sheet, which the ν τ -identification, via the double-cascade method, is sensitive to. The better modeling is expected to lead to a better distinction between single and double cascades around and below the length threshold of 10 m applied in this analysis. The planned IceCube-Gen2 facility [66] will provide an order-of-magnitude larger sample of astrophysical neutrinos and enable a precise measurement of their flavor composition, allowing to distinguish between neutrino production mechanisms [8][9][10][11] with high confidence.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permit-ted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.

The initial, untargeted simulation contains ∼ 7.4 • 10 5 events in the entire HESE analysis range and thus has insufficient statistics for events similar to the ones observed to calculate the probability for each double cascade to have been induced by a tau neutrino. Targeted MC sets were produced to obtain a large number of MC events with similar properties to the observed double cascade data. Such a simulation is computationally expensive, therefore the targeted MC was restricted to a parameter space around the reconstructed parameters of the observed events, as shown in Table 3.

The mapping between true and reconstructed quantities is not straightforward. The interaction vertex in the targeted simulation was restricted to a cylinder with radius 50 m and height 50 m, the direction of the incoming neutrino spans ±35 • in zenith and ±110 • in azimuth, centered on the reconstructed interaction vertex and direction of the events, respectively. For the zenith and azimuth angles, the resolution depends on the event topology. The azimuth region was chosen to cover a wide range to account for possible contributions from azimuthal regions affected by the ice anisotropy and due to the limited azimuthal resolution for single cascades. The zenith region was restricted more as the zenith resolution is better due to the much closer spacing of DOMs in the vertical direction. For event #1 simulated as ν μ CC interactions, the zenith and azimuth were restricted to ±17 • and ±40 • , respectively, to enhance the number of MC events with properties similar to the data, and reflecting the better angular resolution for tracks. In the case of the primary energy, the mapping depends on the neutrino spectrum and the interaction type, and is only well correlated to the reconstructed deposited energy for ν e CC interactions, as only in this case the neutrino deposits its entire energy in the form of visible energy in the detector. All other interactions have some non-visible energy losses -final state neutrinos, intrinsically darker hadronic cascades, muons leaving the detector -such that it is not a priori known what primary energy range will significantly contribute to the region around the reconstructed properties of the data events. The primary neutrino energy was restricted to cover the range of energies that can contribute to the observed reconstructed energies, which had to be determined by trial and error for each simulated interaction type. True quantities for the energy asymmetry and double cascade length are only defined for ν τ CC inter- actions. Those properties were therefore left unconstrained during the targeted simulation. Table 4 lists how many events were simulated for each of the interaction types and what fractions of the simulated events pass the visible energy requirements and HESE selection. These events (∼ 1 • 10 6 "Big Bird"-like and ∼ 21 • 10 6 "Double Double"-like events) are used in the tau neutrino probability assessment and the flavor composition analysis. For reference, the fraction of events classified as double cascades according to the procedure described in Table 1 are also given.

Figure 7 shows the flavor composition measurement using two-dimensional distributions for all three topologies (cos(θ z ) and E tot for single cascades and tracks, L dc and E tot for double cascades) that are employed in the analyses presented in [19,[57][58][59], and the fit using the extended likelihood and targeted simulation for double cascades as shown in Fig. 6 and used in this analysis.

The contours shown in Figs. 6 and 7 were obtained assuming Wilks' theorem holds. The validity of Wilks' theorem was tested for the untargeted MC, by generating pseudo-MC trails. For a one-dimensional fit as used for the astrophysical Fig. 7 Comparison of the flavor composition measurement using the HESE likelihood with ternary topology identification and 2D MC distributions for each topology (blue) with the measurement using the extended likelihood for double cascades as shown in Fig. 6 Fig. 8 Test statistic distribution of the one-dimensional pseudo Monte Carlo trials for an injected vanishing astrophysical tau neutrino flux. The expected half-χ 2 k=1 distribution (gray line) is followed. The χ 2 k=1 distribution is shown for reference (black line). The vertical lines that mark the events contributing to the 68% confidence interval (blue), and the expectation from the half-χ 2 k=1 distribution (gray) are almost perfectly aligned ν τ flux normalization measurement, Wilks' theorem holds over the entire available parameter space for the untargeted MC. The generated pseudo-MC trials are distributed according to a half-χ 2 k=1 distribution, as can be seen in Fig. 8. For the two-dimensional fit used for the flavor composition, Wilks' theorem holds within the flavor triangle and gives a conservative result at the boundary where one of the contributions vanishes.

The analysis model parameters are given in Table 5, modified from [19]. The atmospheric fluxes need to be carefully modeled. This is done via the parameters: φ conv scaling the overall conventional atmospheric neutrino flux normalization, φ prompt scaling the overall prompt atmospheric neutrino flux normalization, R K /π scaling the kaon-to-pion ratio, 2ν/(ν + ν) atmo providing the neutrino-to-antineutrino ratio, Δγ CR accounting for modifications in the cosmic ray spectral index, and Φ μ scaling the atmospheric muon flux normalization. The light yield of an optical module is affected by its overall efficiency and the propagation of photons through the ice to reach the module. Uncertainties in the former are parametrized by the DOM efficiency parameter, DOM , which describes changes in the total efficiency of the DOMs. Uncertainties on the ice properties are parameterized with the head-on parameter, which modifies the angular response of the DOM and depends on local ice properties of the refrozen ice surrounding the DOMs, and a s describing an azimuthal anisotropy of the photon propagation. Events have been simulated using variations of all of these parameters. With the exception of the anisotropy scale a s , the parameter uncertainties affect all topological classes of events in the same way, their effect and treatment is discussed in detail in [19].

The anisotropic photon propagation in the ice can affect the classification of events and needs careful attention. The ice model used in this analysis is called Spice3.2, and contains the South Pole ice sheet's optical properties (scattering and absorption coefficients) at each point in the detector and for each direction of photon propagation. Measurements with insitu IceCube calibration LEDs [56] have shown that the ice is not isotropic [32], i.e., the propagation of a photon depends on its direction. The anisotropy can be modeled as a sinusoidal modulation of the scattering coefficients in azimuth and zenith. Along the glacial flow direction, scattering is reduced by -10% while perpendicular to the glacial flow it is enhanced by +5%. Less scattering leads to photons traveling on straighter paths through the ice, which in turn leads to a cascade being elongated when aligned with the glacial flow. More scattering leads to photons traveling on more random paths, thus a cascade becomes compressed when its direction is perpendicular to the glacial flow. If uncorrected, this effect leads to a larger misclassification of true single cascades as double cascades due to their elongation along the glacial flow, and of true double cascades as single cascades due to their compression perpendicular to the glacial flow. The event reconstruction algorithm uses look-up tables Table 5 Analysis model parameters for the single power-law astrophysical model. Constraints for analysis parameters used in the analysis are shown. The mean and standard deviation are given for Gaussian con-straints, while constraint-free parameters are denoted with "-". Bounds are given for all parameters. Note that the f α parameters are unique to the analysis presented here and compares the received to the expected photon counts per receiving DOM. The look-up tables assume isotropy, and adding dimensions to incorporate the anisotropy would make the photon tables too large and their production as well as each event reconstruction computationally too expensive. As the expected light yield is looked up for each receiving DOM, a simple trick can be used to approximately correct for the anisotropy: the distance between source and receiving DOM is shifted and the expected light yield looked up for the effective distance which contains the effect of enhanced or inhibited photon propagation due to the anisotropy [35]. This first-order correction of the anisotropy of the photon propagation has been verified as sufficient, as the misclassification fraction of true single cascades as double cascades is now constant across the full azimuth range. To model uncertainties in the anisotropy scale, the scale parameter a s is used.

The length bias is defined as the mean difference of reconstructed lengths when the anisotropy is corrected for and when it is not, and is a function of the reconstructed zenith and azimuth. Under the assumption that the length bias scales linearly with the anisotropy scale parameter, the anisotropy scale can be fit from data. As the length bias is a function of reconstructed observables, the systematic uncertainties on the reconstructed length due to the anisotropy can be tested per event. Note that none of the classified double cascades is in a phase space greatly affected by the anisotropy. The uncertainty on the anisotropy scale in this analysis is 20%.

As the direction of the anisotropy is known with sub-degree precision, no uncertainty is assumed on it.