Measurement of atmospheric tau neutrino appearance with IceCube DeepCore

We present a measurement of atmospheric tau neutrino appearance from oscillations with three years of data from the DeepCore subarray of the IceCube Neutrino Observatory. This analysis uses atmospheric neutrinos from the full sky with reconstructed energies between 5.6 and 56 GeV to search for a statistical excess of cascadelike neutrino events which are the signature of ν τ interactions. For CC þ NC (CC-only) interactions, we measure the tau neutrino normalization to be 0 . 73 þ 0 . 30 − 0 . 24 ( 0 . 57 þ 0 . 36 − 0 . 30 ) and exclude the absence of tau neutrino oscillations at a significance of 3 . 2 σ ( 2 . 0 σ ) These results are consistent with, and of similar precision to, a confirmatory IceCube analysis also presented, as well as measurements performed by other experiments.


I. INTRODUCTION
The phenomenon of neutrino flavor oscillations is now well-established experimentally, building on the discoveries of atmospheric neutrino oscillations by the Super-Kamiokande (SK) experiment [1] and solar neutrino oscillations by the Sudbury Neutrino Observatory (SNO) experiment [2,3].
These and most other neutrino oscillation experiments [4] are based on measuring the appearance or disappearance of electron neutrinos or muon neutrinos.In contrast, there are only two experiments with measurements of the appearance of tau neutrinos through neutrino oscillations, leaving the ν τ sector relatively underexplored.With the ν τ appearance measurement of the OPERA [5] experiment, using an accelerator-based beam of ν µ , the null hypothesis of no-ν τ appearance has been effectively ruled out.Additionally, a small excess of ν τ events has been measured by both OPERA (0.25σ) and SK [6] (1.47σ) relative to what is expected under the standard three-flavor oscillation paradigm.
The measured excess may be interpreted in a number of ways.The tau neutrino charged current cross section directly contributes to the total number of detected ν τ , with theoretical uncertainties [7] of O(10)% and much larger experimental uncertainties [8].These uncertainties can lead to an overall scaling of the number of observed ν τ interactions which can be measured by atmospheric oscillation experiments sensitive to tau neutrinos.This interpretation has been adopted in recent results from the SK collaboration, which recasts the excess in terms of a modification of the averaged tau neutrino charged current cross section.
Another potential interpretation for the observed excess in OPERA and SK would be the observation of non-unitarity in the neutrino sector.In the standard oscillation picture, the dominant appearance mode of ν µ → ν τ is given by ≈ cos 4 θ 13 sin 2 2θ 23 sin 2 ∆m 2 31 L 4E ν (2) with U denoting the elements of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) [9,10] mixing matrix (see also Eq. 3), ∆m 2 jk = m 2 j − m 2 k the mass-squared splittings, L the oscillation baseline, and E ν the neutrino energy.The angles θ 13 and θ 23 govern the amplitude of the mixing, while ∆m 2  31 drives the oscillations on the length and energy scales.The benefit of using trigonometric angles is that they conveniently preserve oscillation probabilities to be within 0 and 1 while also reducing the number of physics parameters to fit.But not all measurements of the same angle probe the same individual elements of the underlying PMNS mixing matrix.
Measurements of θ 23 from long baseline ν µ → ν µ disappearance probe |U µ3 | 2 , whereas measurements of θ 23 from ν µ → ν τ appearance probe |U µ3 | 2 and |U τ 3 | 2 , for further information see Supplementary Material of Ref. [11].Not only do different experimental measurements of the same angle probe different underlying elements, but the relation between the angles and the nine canonical matrix elements is only preserved if the PMNS matrix is 3 × 3 unitary.
A core aspect of any theoretically consistent neutrino mixing matrix is that the individual rows and columns preserve norms and rational probabilities, e.g.
While checks of unitarity across the entire matrix are important, the mixing elements for the third mass eigenstate are particularly interesting because it has been experimentally established that |U e3 | 2 + |U µ3 | 2 0.5, but it has only recently been confirmed by OPERA that |U τ 3 | 2 > 0 at 6.1σ significance [5].The only other evidence of |U τ 3 | 2 > 0 is from SK and reaches 4.6σ significance [6].Even with these two measurements, a global fit of leading oscillations results [11] illustrates that the current constraint of 0.2 < |U τ 3 | 2 < 0.61 at 3σ lacks the precision necessary to probe unitarity of the third mass eigenstate at even O(10)% precision.Unsurprisingly, the range of |U τ 3 | 2 from a global fit is not driven by the direct measurements of |U τ 3 | 2 , but rather that values outside that range would induce small deviations in the ν e and ν µ sectors that would exceed their 3 × 3 unitarity constraints.Using only current direct measurements, the allowed region is |U τ 3 | 2 > 0 but is otherwise weakly constrained.
A measurement of |U τ 3 | 2 differing from 0.5 would be further evidence for new physics beyond the Standard Model (SM), and would imply non-unitarity in the ν 3 mass eigenstate, i.e., |U e3 | 2 + |U µ3 | 2 + |U τ 3 | 2 = 1.The impact of such a deviation could indicate the existence of: • Non-standard interactions with the three active neutrinos in the SM.• At least one new neutrino (sterile) which has no SM gauge interaction with normal matter.
In principle, the three channels to measure |U τ 3 | 2 are ν e → ν τ , ν τ → ν τ , and ν µ → ν τ .But, the ν e → ν τ channel is unfavorable because 1) experimentally a ν e and ν τ interaction produce a similar signature in most detectors, and 2) the magnitude of the oscillation is low due to the flavor composition of the third mass eigenstate.The ν τ → ν τ channel probes |U τ 3 | 2 directly, but is also unfavorable because it requires a hitherto unrealized and experimentally challenging high-statistics focused ν τ beam.
Earlier IceCube neutrino oscillation measurements [17][18][19], and the measurement presented here, use atmospheric neutrinos arising mainly from the decay of pions and kaons produced in cosmic ray air showers in the Earth's atmosphere.The initial flux is dominated by ν e and ν µ , and contains negligible numbers of ν τ [20].The atmospheric neutrinos interacting in the DeepCore subarray of IceCube travel distances ranging from L ≈ 20 km (vertically downward-going) to L ∼ 1.3 × 10 4 km (vertically upward-going; the full diameter of the Earth).For vertically upward-going neutrinos, the first peak of maximal ν µ → ν τ oscillation probability occurs at roughly 25 GeV.This is comfortably above the E ν = 5 GeV threshold for the DeepCore neutrino reconstruction used in this analysis [21].The energy corresponding to the oscillation maximum is also above the kinematic energy threshold for charged current ν τ -nucleon interactions E ντ = 3.5 GeV, where for lower energies there is a complete suppression of the cross section due to the relatively high τ lepton mass as compared to the other charged leptons.Even so, there is still a suppression to the CC-ν τ cross section compared to CC-ν e,µ up until E ντ ≈ 10 TeV [7].
The identification of individual ν τ events at energies relevant for measuring atmospheric neutrino oscillation is precluded in DeepCore, as the outgoing tau lepton in CC interactions decays after ≈ 1 mm, far smaller than the meter-scale position resolution of DeepCore.The ν τ CC interactions mainly manifest as "cascades," similar to those from ν e CC and neutral current (NC) interactions of all neutrino flavors.Relative to the no-oscillation case, these ν τ -induced cascade events produce a distortion in the 2-D distribution of neutrino energy and direction (the zenith angle is directly related to the pathlength L in eqs.[1][2].This measurement is based on observing such oscillation-induced patterns between 5.6 GeV and 56 GeV in the atmospheric neutrino flux coming from all directions.
We present results based on two separate analyses that have different strategies for event selection and background estimation, but considerable overlap in their event selection variables and treatments of systematic uncertainties.The sample for our main analysis "A" targets a high acceptance of all-flavor neutrino events and its background estimation is simulation-driven.The sample for our confirmatory analysis "B" is optimized for higher rejection of non-neutrino events and its atmospheric muon background estimation is data-driven.

II. THE ICECUBE DEEPCORE DETECTOR
The in-ice array of the IceCube detector [22], buried in the South Pole glacial ice, comprises 5160 digital optical modules (DOMs).Each DOM houses a downward-facing 10" photomultiplier tube (PMT) [23] and its associated electronics [24] in a glass pressure sphere [23,24].The modules are arranged along 86 vertical strings with 60 DOMs per string (see Fig. 1).Of these strings, 78 are deployed in a nearly regular grid, with an inter-string distance of about 125 m and modules deployed between depths of 1.45 km and 2.45 km, instrumenting a total volume of roughly 1 km 3 .This part of the detector is optimized for neutrinos from 0.1 − 10 5 TeV, and for the analysis presented here primarily serves as an active veto against the downward-going atmospheric muon background.The remaining eight strings, situated at the bottom center of IceCube, form DeepCore [25].The PMTs on these strings have higher quantum efficiency and are primarily located below 2.1 km in the clearest instrumented ice.The DeepCore instrumented volume is roughly 10 7 m 3 with a module density about five times that of the surrounding IceCube array.
While the IceCube detector was fully commissioned in 2011, its noise rates were still stabilizing during the first year of operation.Therefore the data used here are limited to the period from April 2012 through May 2015.

A. Simulation
The simulation chain in IceCube involves three stages: generation, propagation, and detection in ice.Different software is involved at each stage depending on the particle type.

Neutrinos
Neutrino interactions in IceCube are generated following the flux calculation of Honda et al. [26] and using the interaction physics in GENIE 2.8.6 [27,28], which includes the nuclear model, cross sections, and hadronization process [29] based on KNO [30] and PYTHIA [31].The GRV98 [32]   in the DIS cross sections calculations.Muons created in ν µ CC interactions are propagated through the ice using PROPOSAL [33] for fast and precise modelling of the energy losses, while GEANT4 [34] is used to handle the direct propagation of tau leptons and their decay products, including muons, hadrons, and electromagnetic (EM) showers below 100 MeV.For events with EM showers above 100 MeV, shower-to-shower variations are small enough to use parametrizations [35] based on GEANT4 simulations.
The Cherenkov photons produced by the final state particles are then propagated through the ice using GPUbased software [36].This simulation takes into account the optical properties (scattering and absorption) of the ice.For the photons intersecting with a sensor module, the acceptance in terms of arrival angle and wavelength is then taken into account.For analysis B, a measure of the relative variation of optical efficiency among DOMs is included.Additional hits caused by thermal noise, decaying radioactive isotopes in the PMT and DOM glass, and scintillation are added.Finally, the PMT response and readout electronics are simulated and trigger algorithms are applied across the full detector in order to produce simulated neutrino events.

Atmospheric Muon Background
The generation of atmospheric muons is performed using a full CORSIKA [37] air-shower simulation with a hadronic interaction model from [38].The propagation of these background muons and the detection of the Cherenkov radiation are the same as those due to a secondary muon in a neutrino interaction.
At the final level of the event selections (see Sec. III B), the atmospheric muon background is reduced by roughly eight orders of magnitude.The standard simulation tools are too computationally inefficient to produce sufficient amounts of muon background surviving all the selection criteria.In order to estimate the muon background at final level, the two analyses use two distinct techniques: • Analysis A uses an atmospheric muon simulation employing a fast parametrized approach based on [39].This software targets the regions of the weakest background rejection: single low-energy muons aimed at the DeepCore fiducial volume, which make up approximately 75% of the final simulated muon sample.The simulation is approximately two orders of magnitude faster than one covering the entire IceCube detector.Unsimulated regions in zenith and energy are augmented by simulation produced with the CORSIKA simulation package.All simulated atmospheric muons are weighted using the H4a cosmic ray flux model [40].
• Analysis B follows an alternative, data-driven approach to estimate the shape of the remaining muon background.The method uses data side-bands consisting of events that would have been accepted in the final sample had they not included hits in DOMs in one of the corridor regions (see Fig. 1).

B. Selection
IceCube triggers on O(10 11 ) downward-going atmospheric muons, O(10 5 ) atmospheric neutrinos, and O (10) high-energy astrophysical neutrinos per year, placing stringent demands on background rejection efficiency for Ice-Cube analyses.At neutrino energies above about 50 GeV, standard techniques to accept neutrinos and reject atmospheric muon background in IceCube include selecting events which reconstruct as upward-going, have a starting vertex deep within the detector fiducial volume, fall within a very narrow temporal or directional window, or have a very high energy.For lower-energy neutrinos, however, only DeepCore's higher density of DOMs allows accurate reconstruction of these dimmer events, as described in the next section.The ν τ appearance analyses therefore focus on events that are contained within the DeepCore fiducial volume.Located at the bottom center of IceCube, DeepCore benefits from the exceptionally clear ice that has photon attenuation and absorption lengths of roughly 50 m and 150 m, respectively [41,42].An important benefit of DeepCore's location is the use of over 4500 IceCube DOMs as an active "veto region" to identify background muons for removal.
The selection of the final event sample is implemented in a series of "levels," the first three of which are very similar in Analyses A and B, while the subsequent ones differ.These differences primarily reflect the looser Analysis A selection criteria to prioritize the efficiency of selecting neutrino events, versus the tighter Analysis B criteria to prioritize the rejection of atmospheric muon background.Note that the Analysis B selection criteria were originally optimized to measure ν µ disappearance and follow closely the criteria used for that measurement in Ref. [19].Below we give a description of the selection criteria, highlighting important similarities and differences between the two analyses, and show distributions for some of the key variables central to the analyses.We provide a more detailed description of the selection criteria in Appendices B-D.

Common Selection Criteria
Analyses A and B share the first three levels of selection criteria, starting with the online triggering at the South Pole.Detected photons or "hits" are labeled "locally coincident" and included in the trigger if they occur within 1 µs of a hit on a nearby DOM on the same string.The trigger requires three or more locally coincident DeepCore DOMs to detect hits in a 2.5 µs time window.When this condition is met, the data acquisition system reads out all available data in the full detector, in a time window that extends 6 µs before and 6 µs after the dynamic trigger window (see Sec. 6.4.2 of [22] for more details).In Level 2, a filtering algorithm is used to reject any events consistent with a muon traveling at v c between the reconstructed interaction vertex within DeepCore and two or more more hit DOMs in the veto region [25].
After the application of the trigger and filter algorithms, a large number of background events are still present in the sample.Both analyses therefore perform a fast reconstruction at Level 2 that insures an adequate number of hit DOMs in IceCube consistent with either the track or cascade signature of a neutrino interaction.Both analyses then define a slightly enlarged fiducial volume, and require < 7 photoelectrons (p.e.) in the correspondingly smaller surrounding veto region.A set of criteria is also applied to remove low quality events with too many noise hits, too few DOMs with multiple hits, too much deposited charge, a reconstructed vertex in the upper region of the fiducial volume, too small a fraction of the event's total p.e. deposited at early times, or too large a fraction of DeepCore hits in the outer regions of the DeepCore fiducial volume.Detailed descriptions of these criteria, along with subtle differences between the two analyses, are discussed in Appendix B. In aggregate, these criteria remove events whose reconstruction is likely to be faulty, and those events that are likely to be downward-going atmospheric muon background.The event rates after each of these first three levels of the common event selection for Analyses A and B are shown in Table III of Appendix B.

Additional Selection Criteria: Analysis A
Event selection for Analysis A uses two boosted decision trees (BDTs) [43] to remove atmospheric muon background.The first BDT (Level 4) uses six different input variables adapted from [44]: three related to the charge measured by the PMTs, a simple vertex estimator, an event speed estimator, and a calculation of event shape.The resulting BDT output is shown in Fig. 2.
Accidental triggers due to random detector noise occur primarily in the DeepCore fiducial volume with few hit DOMs, appearing neutrino-like for this selection level.In order to limit the impact of these events, dedicated selection criteria, detailed in Appendix C, are introduced at later stages of the selection.
The second, subsequent BDT (Level 5) is used to further reduce the muon background based on six input variables: the time to accumulate charge, a vertex estimator, two variables using center-of-gravity calculations, a causal hit identifier, and a zenith angle estimation from a simple reconstruction.As an example, the distribution of this second BDT output for both simulation and data is shown in Fig. 3, and more distributions and information can be found in Appendix C.
The event rates after application of the Level 4 and 5 selection criteria are shown numerically in Table III of Appendix B and graphically in Fig. 4 below.After Level 5 the signal and background rates are roughly at parity.
Following the application of the two BDT-based selections, a series of individual event selection criteria are applied (Level 6).Requiring events to have a sufficient number of hits inconsistent with intrinsic DOM noise and to be spatially compact removes most remaining events caused purely by intrinsic noise hits.Removal of many of the remaining atmospheric muons is accomplished by requiring a likelihood-based vertex estimate to be well contained in the DeepCore fiducial region, and by reject-  and time-consuming reconstruction applied (Level 7).In addition, events are required to have a reconstructed deposited energy between 5.6 GeV and 56 GeV.

Additional Selection Criteria: Analysis B
The Analysis B sample applies several selection criteria at Level 4. These include requiring a sufficient number of p.e. deposited in the largest cluster of hits in the fiducial region, a minimum number of non-isolated hits in the fiducial region, an event vertex contained in the fiducial volume, a space-time interval between the first and fourth temporal quartiles of the hits in DeepCore consistent with v ≤ c, no more than 5 p.e. in the surrounding veto region, and no more than two hits in the veto region consistent with speed-of-light travel to the hit in DeepCore whose time is closest to the event trigger time.These criteria reject events caused by noise, reduce muon background, and favor the more cascade-like signature produced by most ν τ interactions.A BDT is then applied (Level 5), using 11 input variables, derived from the charge, time, and location of the hit DOMs, as well as reconstructed zenith angle and event speed using crude but fast track reconstructions.
Following the application of the BDT, events consistent with entering through corridor regions are rejected, and reconstructed events are further required to have starting and stopping positions in or near the DeepCore fiducial volume.These Level 6 criteria further reject atmospheric muon background.At this stage in the processing, the neutrino signal rate has been reduced by a factor of roughly 13 while the atmospheric muon background rate by a factor of 10 8 .A more detailed breakdown is provided in Table III of Appendix B.

C. Reconstruction
The reconstruction used in both Analyses A and B assumes that every event starts with an electromagnetic or hadronic shower followed by a finite, minimum ionizing muon at the same primary vertex.Due to the numerous charged particles in the shower, a cascade-like event is characterized by a localized Cherenkov light pattern centered at the interaction vertex.On the other hand, a track-like event involves a muon which deposits Cherenkov light uniformly along its trajectory, and travels much further than any non-muon particles produced in the primary shower.With the cascade plus track assumption, the reconstruction algorithm describes an event via eight parameters: the primary interaction vertex position (x, y, z) and time (t), the direction given by the zenith angle (θ ν ) and the azimuth angle (φ ν ) of the neutrino, the energy of the primary cascade (E cscd ), and the length of the track from the minimum ionizing muon (L µ ).
Based on the above hypothesis, a likelihood-based reconstruction method compares the observed pattern of photon counts from all active DOMs in an event to that predicted.The PMT measures a charge linearly related to the number of Cherenkov photons arriving at a DOM.Using the PMT charge as a proxy for photon counts, the number of photons arriving at the DOM is described by the time-binned PMT charge.The predicted pattern of charges from all DOMs in an event is then fitted to that of the observed event with the eight parameters in the event hypothesis allowed to vary freely.
To reduce computational complexity in running the reconstruction, energy deposition during an event is described using several independent light sources.In particular, the deposited energy from the primary cascade is treated independently of that from a muon track.Further, the energy deposition by the muon track is also discretized into segments with constant length.The total length of the track L µ is directly related to the energy of the track via E trck = L µ dEµ dx , where the differential energy loss of a minimum ionizing muon in ice, dE µ /dx, is fixed to 0.22 GeV/m.
The energy deposition of a muon along its track is not constant in reality nor in our simulation.This simplification is only used for reconstruction of low energy events and yields a good approximation at the O(10 GeV) scale.The approximation begins to break down above about 50 GeV when stochastic losses along the muon track become non-negligible [45].
The expected charge q i (t) at the i th DOM at time t is estimated by the charge due to energy depositions by the cascade E cscd and by the track E trck plus a timeindependent noise term n i .The expected charge can be expressed as, where Λ cscd represents the charge expectation for a 1 GeV cascade and Λ track for a minimal ionizing muon of one segment length.A linear relation between Λ cscd in Eq. 4 and the deposited cascade energy is assumed.To obtain the values of Λ cscd and Λ trck , large sets of look-up tables are generated from simulations of photon propagation in the ice [46].These tables, used with the assumption that the number of Cherenkov photons emitted is directly proportional to the deposited energy of the particle, allow for the calculation of the expected charge from an arbitrary cascade or track.
The process of finding the maximum likelihood hypothesis for an event is an eight-dimensional optimization problem, and the likelihood space is typically non-convex,i.e.populated with local maxima.To cope with these challenges the MultiNest algorithm [47] is used to find the best-fit hypothesis.
Both presented analyses follow the above reconstruction algorithm but with two main differences.First, each track length segment in Analysis A is 5 m long, whereas analysis B uses coarser 15 m long segments.Second, the reconstruction used in analysis A ignores the observed charges, instead implementing a binary response of 0 p.e or 1 p.e. per 45 ns in each DOM individually, while Analysis B uses the observed charge in each DOM.The treatment of charge in Analysis A reduces the impact of observed discrepancies observed between the distributions of the average charge per DOM in data and simulation, which affect mainly the stochastic nature of charge depositions in events with a small number of hit DOMs.

E
reco (GeV) Despite the differences between the two analyses, the energy and cosine of zenith angle resolutions of the two analyses are similar, as shown in Fig. 5 and Fig. 6.

D. Classification
A ν µ CC neutrino interaction often produces an event with an identifiable track, whereas events from ν e CC and all-flavor ν µ,τ,e NC have a cascade-like topology.Most ν τ CC interactions also produce cascade-like events, with the short-lived τ lepton decaying roughly 83% of the time to non-muon modes [4].The ≈ 17% muonic decay mode is τ − → µ − νµ ν τ (and charge conjugate), where the daughter muon may have sufficient energy to create a visible track indistinguishable from a ν µ CC event causing it to be identified as a track-like.To improve the sensitivity of ν τ measurement, both analyses divide their samples into cascade-and track-like subsets, enhancing the purity of ν τ events in the cascade channel.
To determine if an event is cascade-like or track-like Analysis A relies on the reconstructed track length L µ .Events with a track length between 0 m and 50 m are considered cascade-like, and events with track lengths longer than 50 m are considered track-like.For Analysis B, an additional reconstruction is performed with the track length forced to 0 m.Events are then classified based on the log-likelihood difference between the cascade-and-track hypothesis and that of cascade-only; ∆LLH reco = ln L cascade+track − ln L cascade .Events with ∆LLH reco > 2 are considered as track-like, while events with −3 < ∆LLH reco < 2 are cascade-like.The cascade only reconstruction should in principle never yield a likelihood that is better than the track+cascade one, but due to finite precision of the minimization process, negative ∆LLH reco do occur.We allow events with a negative ∆LLH reco as low as −3; the remaining events are removed from the analysis due to their bad reconstruction quality.As shown in Fig. 7, the cascade and track separation powers from the two analyses are similar above 20 GeV.

IV. ANALYSIS
Our tau neutrino appearance analyses yield two distinct quantities: the level at which the null hypothesis of no ν τ appearance is rejected and the measurement of the ν τ normalization, which is defined as the ratio of the measured ν τ flux to that expected assuming bestfit oscillation and other nuisance parameters for that ν τ normalization.These best-fit nuisance parameters are obtained simultaneously with the best-fit tau normalization during the optimization process, meaning that the expected distribution of tau neutrino events can be understood as: (ν τ normalization) × ((baseline ν τ expectation) + (nuisance parameter)×(ν τ systematic change)).
Since DeepCore cannot distinguish between ν τ CC and NC interactions, our analyses benefit from treating them on an equal footing by applying the ν τ normalization to both CC and NC tau neutrino interactions.However, to facilitate comparisons with results from other experiments, a second set of measurements are also performed applying the ν τ normalization only to the ν τ CC component.In this second case, the ν τ NC component is unaffected by the value of the ν τ normalization.In both the CC+NC and CC-only cases, there is a separate uncertainty assigned to all neutral current events.
In Analysis A, data is binned into a 3-d histogram with eight reconstructed energy bins spaced logarithmically between 5.6 GeV and 56 GeV, 10 reconstructed cosine zenith bins spaced linearly between −1 and 1, and two reconstructed track length bins for particle identification (PID).Track-like events in Analysis A have a reconstructed energy of at least 10 GeV associated with the minimum track length of 50 m.Therefore, the first two energy bins for track events are empty by construction and not included in the analysis.Figure 8 shows the S/ √ B as a figure of merit, where S and B are the number of signal and background events, respectively.The figure indicates that upward-going cascade events with reconstructed energies around 20 GeV dominate the measurement.With the same energy binning as Analysis A, Analysis B covers the same cosine zenith range with eight bins instead of 10 and uses ∆LLH reco for PID instead of reconstructed track length.
In each of Analyses A and B, a χ 2 minimization is performed on the binned data as a function of the ν τ normalization and nuisance parameters associated with the relevant systematic uncertainties, see Sec.V.The χ 2 function is defined as where N exp i is the number of total events expected from the signal and all background events in the i th bin, and N obs i is the number of events observed in the i th bin.For both Analyses A and B, the denominator consists of the standard Poisson variance N exp i and the uncertainty in the prediction of the number of expected events σ exp i of the i th bin.Analysis A uses Monte Carlo simulation for the prediction of all event types, and the term σ exp is the sum of uncertainties due to finite statistics of MC simulation from each event type.In Analysis B, the term σ exp encompasses both the uncertainty due to finite MC statistics as well as the uncertainty in the data-driven muon background estimate described in the Sec.V F. The second term of Eq. ( 5) is the sum of penalty terms for nuisance parameters that have prior constraints imposed, where s j is the central value of j th systematic parameter, ŝj is its maximum likelihood estimator, and σ 2 sj is the prior's Gaussian standard deviation.
For both analyses the uncertainty due to limited MC statistics is small for signal neutrinos, as the effective livetime for simulation is an order of magnitude higher than that of the acquired data.The situation is different for the muon background predictions: for Analysis A the uncertainty arises from simulation with less effective livetime than the actual data and for Analysis B from a data side-band, in both cases resulting in larger uncertainties than for signal neutrinos.However, any ensuing variations are predominantly constrained to the track-like and downward-going region of the event sample which is away from the cascade-like and upward-region region associated with our targeted signal events.
While both analyses use data from the same operating period of April 2012 through May 2015, minor differences in the event selection criteria lead to a total livetime of 1006 days for Analysis A and 1022 days for Analysis B. Table I shows the expected number of events at the best fit point for each neutrino flavor and interaction type, and for atmospheric muons and noise-triggered backgrounds.The effect of systematic uncertainties is included in the analyses with nuisance parameters that impact the shape and normalization of the expected event distributions.The uncertainties considered can be broadly grouped in categories according to their origin: the initial unoscillated flux of atmospheric neutrinos, neutrino-nucleon cross sections, neutrino flavor oscillation parameters, detector response, and atmospheric muon background estimates.The associated parameters, together with their best-fit values, are summarized in Table II.Each category of uncertainties will be discussed in turn.
To quantify the impact of each systematic uncertainty, the 1σ confidence interval of the expected tau neutrino normalization measurement was calculated while fixing one parameter at a time.The resulting change in the confidence interval is shown in Fig. 9.Of the fitted systematic uncertainties, the neutrino mass splitting provides the strongest impact on the final confidence interval.The impact of each category of systematic uncertainties was also tested in a similar way.When entire categories of systematic uncertainties are fixed at the same time, the largest impact comes from the detector uncertainties, which account for 41% (36%) of the NC+CC (CC) measurement FIG. 9.The relative impact from each systematic uncertainty and each group on the final 1σ confidence interval width in Analysis A. Each systematic uncertainty is fixed to the bestfit value in turn and the change in the interval is measured.The most important systematic uncertainty is ∆m 2  31 , with a 14% (16%) impact on the NC+CC (CC) measurement.The detector uncertainties show degeneracies that limit the impact of individual parameters, but together account for 41% (36%) of the uncertainty in the NC+CC (CC) measurement in Analysis A.
in Analysis A. This is due to individual systematic variations being correlated, especially the ones in the detector uncertainty group.
In addition to the systematic uncertainties mentioned above and included in the analysis, we have studied different optical models for the glacial ice as well as a newly available charge calibration for the detector.In both cases, the impact on the final result was found to be negligible, and they were thus omitted from the fit and the error calculation.

A. Atmospheric neutrino flux
The measurement presented in this work is extracted from an observed distortion of the flux of atmospheric neutrinos.Our nominal model is the calculation of Honda et al. [26].The calculation covers the energy range 100 MeV to 10 TeV, and was produced specifically for a detector situated at the geographic South Pole, so local geomag-netic effects are included.The cosmic rays that contribute the most to the neutrino production at the energies of interest, between 5.6-56 GeV, are protons and helium.Honda et al. model the energy spectrum of each of these incident particles using a single power law, fitting the flux to data from satellite and balloon experiments.In this calculation, interactions of cosmic rays with the Earth's atmosphere are simulated using a combination of the JAM interaction model [48] and a modified version of DPMJET-III [49].The modifications, discussed in [50], are changes to the yields of π and K mesons to reach a better agreement with muon measurements from the BESS experiment [51].The atmospheric conditions such as temperature and column density are taken from the NRLMSISE-00 model [52], whose authors estimate the resulting calculation has an uncertainty on the neutrino flux of ≤ 15%.Seasonal variations are included in the flux calculations, one-year averaged values are used in our analyses.
In both Analyses A and B, a detailed modification of the neutrino flux prediction as a function of energy, zenith angle, and particle species has been used.The basis of this modification is the work of Barr et al. [53], who have performed a detailed study of the uncertainties on neutrino flux predictions by systematically modifying the inputs required to perform the calculation.Their work suggests that, for the energies that are of interest here, the flux calculation is mostly affected by the uncertainties on the spectral index assumed when modeling the cosmic ray fluxes, and the lack of measurements on the production of π and K mesons with energies above 500 GeV and 30 GeV, respectively, and where the secondary particle contains > 10% of the incident particle energy.
A modification to the spectral index on the cosmic rays translates into a very similar modification of the neutrino flux.We therefore account for this uncertainty by modifying the neutrino flux using the function E ∆γν , which only depends on neutrino energy.Modifying the yields of pions and kaons in hadronic interactions produces changes in the neutrino flux, not only as a function of energy, but also incoming zenith angle for each of the particle species in it.In [53] a summary of these modifications is shown for the ν/ν flux ratio as function of energy and as function of zenith angle for three energy regions, and the upward-going to horizontal ν ratio as a function of neutrino energy.We use that information to build a model able to reproduce the effects described as function of both energy and zenith angle.
In summary, four effective parameters account for the uncertainties considered on the atmospheric neutrino flux.These are a modification of the spectral index (∆γ ν ), the ratio of ν e to ν µ fluxes ("ν e /ν µ ratio"), the ratio of the ν to ν fluxes as function of zenith angle and energy ("ν/ν ratio"), and an additional parameter for the remaining uncertainty in the upward-going vs. horizontal flux of electron neutrinos ("up/hor ratio").All parameters are introduced assuming that they are uncorrelated.A 5% uncertainty is assumed for the ν e /ν µ flux ratio.The two parameters that modify ν/ν and ν e up/hor receive an uncertainty such that, when both are evaluated at 1σ, the results from [53] are reproduced.They roughly correspond to a 10% energy-dependent change to the neutrino flux with a 3% zenith-dependent modulation.The top two panels of Fig. 10 demonstrate the effect of these parameters in the reconstructed final sample.The error assigned to ∆γ is discussed in the next section.
Sources of uncertainty that result in a global scaling of the neutrino flux, independent of energy or zenith angle, are not considered in this work as the normalization is left free in the fit (scaled by the effective livetime parameter).

B. Atmospheric muon flux
While the sources of uncertainties discussed above are included in both Analyses A and B, an additional uncertainty related to neutrino-muon coincidence is taken into account in Analysis A. An extra simulation set was produced, in which every neutrino event is contaminated by an atmospheric muon resulting from an independent air shower.Together with the baseline neutrino sets with no muon contamination, the event count is parametrized per bin as a function of coincident fraction.Because previous high-energy analyses using the IceCube volume found less than 10% contamination due to coincident muons, a one-sided Gaussian prior centered at 0 with a width of 10% is applied to the coincident fraction for Analysis A. The effect from neutrino-muon coincidences is normalized to leave the total event rate unchanged.
Analysis A also considers an uncertainty related to the cosmic ray spectral index in the atmospheric muon flux.Atmospheric background muons in Analysis A are produced in air showers of energies 1 TeV to 1 PeV.These shower energies are higher than the expected energies from the atmospheric neutrinos making it into the final analysis.To be conservative, the effect of a change in the cosmic ray spectral index is treated independently between neutrinos and muons to account for the separate energy regimes probed.
Measurement uncertainties from a fit to cosmic ray experimental data [54] are used to obtain an estimate for the uncertainty on the spectral indices associated with proton and helium cosmic ray primaries.Based on the error bars from the experiments, the deviation from the central fit value is determined as a function of primary energy using CORSIKA simulations.This change in the flux weighting for atmospheric muons is parametrized as a function of true energy and zenith angle and applied to the final simulated atmospheric muon sample.A Gaussian prior is applied to the spectral index uncertainty, with a 1σ deviation in the parameter corresponding to a 1σ change in the cosmic ray spectral index.

C. Neutrino-nucleon interactions
Deep-inelastic-scattering (DIS) interactions make up the bulk of neutrino interactions visible in DeepCore.The uncertainties associated with these interactions were investigated in the final samples.
The first studies were on the parameters used in the Bodek-Yang model to allow the parton distribution functions used in the calculation of cross sections to be extended to the lower Q 2 region [55].These DIS events were re-weighted on an event-by-event basis in response to changes in the higher-twist parameters and valence quark corrections using the reweighting scheme included in the GENIE generator [27].Though this did have a small impact on the final analysis, they were fully degenerate with either the overall neutrino scaling provided by the neutrino event rate (via the "effective livetime" parameter) or the energy dependent scaling provided by the spectral index parameter ∆γ ν .Since these two systematics fully absorb the effect of the uncertainty in the Bodek-Yang model, no additional parameter was included in the final analysis.
We further investigated the impact of both high-and low-W averaged charged hadronization multiplicity, a systematic uncertainty also related to DIS interactions [56].These studies were done by modifying PYTHIA to change the multiplicity of outgoing charged particles to be within the range observed by bubble chamber experiments [57][58][59].These changes were then propagated through GENIE to evaluate the effect on the final sample.It was found this has less than 0.1% impact on events at the final level, with the change being energy dependent.Due to the small size of this effect and its shape being degenerate with that of spectral index changes (∆γ ν ), we did not include this as an additional parameter in the final fit.
The final DIS uncertainty studied was its differential cross section.The approach here was to modify the structure function as a function of the Bjorken-x within the uncertainties measured by NuTeV [60].This resulted in a change at final level of less than 1% up to 3% at 200 GeV.As with the studies on hadron multiplicity, these changes are degenerate with a change in the spectral index uncertainty and so are not included in the final fit.
Many cross section systematic uncertainties were tested, but the only two which were not already degenerate with other systematic uncertainties were the axial mass form factors for charged current quasi-elastic (M CCQE A ) and resonant (M res A ) events.Both of these are included in the final analysis and change the expected number of CCQE or resonant events seen in the event sample.The systematic is implemented so that a change in M CCQE(res) A will result in a change to each CCQE (resonant) event weight on an event-by-event basis using GENIE's re-weighting capabilities.
The nominal value used for M CCQE A is 0.99 GeV, with an uncertainty of (−0.1485, +0.2475) GeV used as a prior; for M res A we used 1.12 GeV and ±0.22 GeV.These are the same values used as GENIE's default model and reweighting scheme, respectively.The last row of Fig. 10 shows that the impact of the M res A uncertainty on the event distribution is energy dependent, with the largest impact at lower energies where the majority of resonant events are expected.The effect of M CCQE A follows a similar shape with even smaller changes, as the quasielastic events are peaked at lower energies.The axial mass uncertainties have little impact as a function of cos θ ν .
Measurements of the ν τ cross section exist from only a few experiments, with DONUT providing a ratio of σ(ν τ )/σ(ν e,µ ) of 1.37±0.35±0.77[8].Uncertainties on the ν τ CC cross section in the energy range of interest in this analysis [7] differ primarily by a factor degenerate with the ν τ CC normalization tested.Indeed, this degeneracy is used by SK to reinterpret their best-fit ν τ normalization as a modification to the ν τ neutrino CC cross section [6].Due to this degeneracy, we do not include any nuisance parameters specifically modifying the ν τ CC cross section.

D. Oscillation Parameters
The model in this analysis assumes three-flavor oscillations and hence relies on three mixing angles, two mass-squared splittings, and a CP violating phase.We use the Prob3++ [61] software which incorporates matter effects for full three-flavor oscillations calculations.The earth is approximated with 12 radial layers of constant density [62].For earth crossing neutrinos, matter effects start to significantly alter the ν e ↔ ν µ transition probabilities only at energies of around 6 GeV and below, hence the effect is very small for these analyses.
With atmospheric neutrinos we are not sensitive to the solar parameters, so we fix the mass splitting ∆m 2 21 to 7.5 × 10 −5 eV 2 and the mixing angle θ 12 to 33.48 • .The reactor angle θ 13 is treated as a systematic uncertainty in Analysis B and is assigned a Gaussian prior with a central value of 8.5 • and an uncertainty of ±0.21 • .All of the above values are taken from [63].
No prior constraints are used for the two atmospheric parameters ∆m 2  31 and θ 23 which vary freely in the fit.Since this analysis is insensitive to δ CP it is fixed to 0 • .Also, since the neutrino mass ordering is not yet known, we check both normal and inverted orderings in the fit and accept the one yielding the better likelihood.To avoid any bias in the fitted value of θ 23 , we fit its value in both octants (sin 2 θ 23 < 0.5 and > 0.5) and accept the value yielding the maximum likelihood.

E. Detector Uncertainties
Systematic uncertainties related to the response of the detector itself play an important role in the analyses.The impact of these uncertainties is complex, depending upon the properties of the detector, on the impacts in event TABLE II.Nuisance parameters along with their associated priors where applicable and the best fit values from Analysis A when fitting the charged and neutral current ντ normalization combined (NC+CC) and the charged current alone (CC), and the same for Analysis B. Priors are given as central value together with the ±1σ ranges when a Gaussian prior is imposed, while "-" denotes that no external prior constraint (i.e.flat prior) is used.
Analysis selection, and on the effect in the reconstruction used to estimate particle properties.In order to account for these complexities, separate simulations for different settings of the detector response were produced and propagated through each step of the event selection and reconstruction as described in Sec.III.Each simulation set includes a change to at least one detector uncertainty parameter.The change in the number of expected events for each of the analysis bins relative to the baseline simulation set is used to estimate the effective impact of each systematic uncertainty for each simulated discrete point of parameter settings.
To arrive at a continuous description, the effects are approximated using a function with linear dependencies on the nuisance parameters.For N linear parameters, we use N -dimensional "hyperplanes" as given in the following equation for each bin k in the analysis histogram: with the nuisance parameters p i , the fitted hyperplane slopes a i , and the common offset b.Thus for N parameters N + 1 values are fitted.Such parameterizations are obtained independently for every analysis bin, separately for each of the three neutrino flavors in CC interactions, and combined for all NC interactions.These relative changes of event rates are then applied as scale factors to the event weights during the analysis.In Analysis A, detector response uncertainties of simulated atmospheric muons are also parametrized in a similar way to the neutrino uncertainties.Variations in the overall efficiency of the optical modules and the absorption result in particularly strong changes in the observed light yields in the veto region from muon tracks, leading to large changes in the atmospheric muon event rates after selection.These simulated muon rates are not well-modeled with linear parametrizations.In these two cases, an exponential form is instead used, giving the form for each bin k as a jk e −b jk pj +c k (7) where N parameters describe the lateral and head-on optical efficiency as well as the scattering of the glacial ice and M parameters cover the overall efficiency and the absorption.
The values f k give the fractional change for each histogram bin given the values of the detector nuisance parameters p.This is applied as a multiplicative reweighting factor for each bin of the analysis histogram.
Both analyses incorporate six nuisance parameters to account for detector uncertainties.Each nuisance parameter is modeled by 2-5 additional simulation sets for each neutrino flavor and, in the case of Analysis A, atmospheric muons.Using the obtained parametrizations, we obtain an average χ 2 /expected degrees of freedom, per flavor and bin, of 13.1/13 across the included neutrino simulation sets and 6.0/6 for background muon sets in Analysis A. Similarly, a χ 2 distribution with 24.0/25 degrees of freedom is obtained from the neutrino simulation sets in Analysis B.
The transparency of the ice in our fiducial volume was calibrated using remotely-controlled light-emitting diodes (LEDs) inside every deployed DOM.The optical properties affect the light yield and temporal arrival distributions of photons that are produced from events seen by the DOMs.The parameters in the model-scattering and absorption coefficients as a function of depth-were determined as a function of location within the detector as described in [41,42].Both coefficients have associated uncertainties of ±10% and are included as systematic uncertainties in this measurement.Additional MC sets were produced with enhanced scattering (+10%), enhanced absorption (+10%), and diminished scattering and absorption (−7%, −7%) to estimate the effects.
The overall photon detection efficiency of the IceCube DOMs depends on both individual PMTs as well as properties of the glass housing and nearby cables.Dedicated measurements of the efficiency of the DOMs yield a relative uncertainty of 10% [22].This effect is modeled by changing the light collection efficiency of the DOMs in simulation, with the efficiency of all modules scaled simultaneously by a common factor.Simulated data sets ranging from 88% to 112% of the nominal optical efficiency were used to parametrize the effect of the DOM efficiency uncertainty and a Gaussian prior with a width of 10% was applied to the overall photon collection efficiency for these analyses.
In addition to modifying the absolute efficiency, any bubbles in the refrozen ice in the borehole ("hole ice") near the DOMs can cause increased scattering of Cherenkov photons.The effect of the refrozen ice column is modeled by two effective parameters controlling the shape of the DOM angular acceptance curve (see Fig. 11).The lateral parameter controls the relative sensitivity between photons traveling roughly 20 • above and below the horizontal.The uncertainty on this parameter is constrained by LED calibration data [41].
Simulated data sets were generated covering the ±1σ uncertainty range and a Gaussian prior based on the calibration data is used for this parameter.The head-on parameter modulates the sensitivity for photons traveling upwards and arriving near the DOM's lower face.This is a region that is poorly constrained by the string-to-string LED calibration because no bright, upward-pointing LEDs were deployed.To account for this uncertainty, the acceptance curve is altered using a dimensionless parameter ranging from −5 (corresponding to a bubble column completely obscuring the DOM's lower face for vertically incident photons) to 2 (no obscuration).Simulated data sets covering the range from −5 to 2 were used to parameterize this effect.No prior is imposed on this parameter due to lack of information from calibration data.Modelling the hole ice via the angular acceptance curve is an approximation, as it only truly holds in the far field.In addition, it can only model hole ice radii significantly larger than the DOM radius as no azimuthal dependence is incorporated.
An additional model of the hole ice has also been tested in Analysis B, incorporating an explicit simulation of the bubble column consisting of ice with enhanced scattering located in the refrozen holes [64].In addition, photons arriving at a DOM are not accepted based on their incident angle, but by requiring that they impact the DOM's lower hemisphere.Although in principle more realistic than the angular acceptance model, the tuning of all parameters involved in such a simulation is a challenge.Various MC sets for a range of different settings (optical properties of bubble column ice, column radius) and using the best knowledge of the position of the column with respect to each DOM were produced.For comparison, the fraction of photons arriving at several DOMs as a function of the arrival direction is also shown in Fig. 11.A fourth parameter (the local ice model) is introduced in Analysis B to account for differences not covered by the angular acceptance model.A value of zero corresponds to the purely angular acceptance base simulations, while a value of one is assigned to the explicit bubble column simulations.Since this model is disfavoured by the data of Analysis B, Analysis A only incorporates the angular acceptance model.

F. Atmospheric Muon Uncertainties
The last nuisance parameter pertains to the amount of atmospheric muon contamination in the final data sample, where Analysis A is based on Monte Carlo simulation and parameterizations while Analysis B is data-driven.For Analysis A, uncertainties due to atmospheric muons flux include the uncertainties associated with the cosmic ray spectral index in Section V B based on [54].Additional uncertainties due to detector response are treated the same way as the case of neutrinos, where additional sets are produced and a hyperplane fit is performed per bin.
For Analysis B, a data-driven method is used to estimate the shape of this background as described in Section III B (see Fig. 12).With the absolute efficiency for tagging background events not amenable to direct measurement, the normalization of the muon contribution is left unconstrained in the fit.Its nominal value is set to match the expected rate from simulated atmospheric neutrinos, and error terms are calculated with respect to this nominal value.In addition, we account for uncertainties in these background templates arising from shape changes when modifying the selection cuts.Two samples are obtained by requiring more than one hit and more than two hits in the muon veto regions, with the latter being a more muon-rich sample.The difference in shape between the two (ignoring normalization differences) is added in quadrature, together with the limited statistics term, to the uncorrelated uncertainties σ exp in Eq. 5.The output shape and uncertainty are in agreement with muon simulations.

VI. CONFIRMATORY MEASUREMENT OF ATMOSPHERIC NEUTRINO OSCILLATION PARAMETERS
Under the assumption of a unitary PMNS mixing matrix, the atmospheric neutrino oscillation parameters ∆m  Previous measurements of the atmospheric neutrino oscillation parameters have been performed using the Ice-Cube detector, including a measurement of atmospheric muon neutrino disappearance performed using the event sample from Analysis B, here referred to as "IC2017" [19].The IC2017 analysis included a subset of systematic uncertainties from Analysis B found to be significant for the disappearance measurement.Detector systematics related to the optical efficiency were included, but used a different parametrization of the detector uncertainties than that described in Sec.V E.
Figure 13 shows the 90% allowed region of atmospheric neutrino oscillation parameters for the analyses based on Analysis A and IC2017, and the allowed regions reported by other experiments.Overall, the 90% allowed regions from Analysis A and IC2017 are statistically consistent, and both results compare favorably with the latest published 90% contours from other neutrino experiments [65][66][67][68].
The shift between the two IceCube contours in both ∆m 2  23 and sin 2 θ 23 can be explained by statistical fluctuations alone, as 65% of events in Analysis A are unique with respect to Analysis B (while 48% of events in Analysis B are unique).Furthermore, detailed investigation showed that differences in the analyses-namely differences in the parametrization of detector effects (hyperplane), inclusion of bulk ice uncertainties, and the differences in the event selection and reconstruction as described in Sec.III-can lead to small (< 0.5σ) systematic shifts in the result as well.
Separate analyses based on the same event samples as presented here, but testing the neutrino mass orderings, were performed using only up-going events.Results including best-fit oscillation parameters are reported in [69], and are also compatible with the values reported in this section.

VII. RESULTS AND CONCLUSION
The distributions of the reconstructed neutrino energy, the reconstructed zenith angle, and the event class for the best fit tau neutrino hypothesis for Analysis A are shown in Fig. 14, overlaid with background-subtracted data, i.e., cosmic-ray muons and all non-ν τ neutrinos subtracted.Figure 15 shows all events projected onto the L/E axis for the best fit expectations overlaid with the observed data for both analyses separately.The excellent agreement of the model with the data can be seen qualitatively in the figure .Using the actual measurement bins and setting all parameters to their best fit values, the model agrees well quantitatively in Analysis A (B) with the observed data with a total χ 2 of 127.6 (113.3),corresponding to a p-value of 55% (20.3%), estimated via pseudo-data trials.The corresponding values for the nuisance parameters can be found in Table II.
Figure 16 shows the expected and observed ∆χ 2 values for a ν τ normalization ranging from 0 to 2.0.The band of expected values assumes standard oscillations with a ν τ normalization of 1.0.Our main result for the CC+NC measurement has a best fit value of 0.73 with the 68% confidence interval (C.I.) covering the range (0.49, 1.07) and the 90% C.I. covering (0.34, 1.30).For the CC-only normalization, we observe the best fit at 0.57 with the 68% C.I. (0.30, 0.98) and the 90% C.I. (0.11, 1.25).
These measured values are compatible with corresponding values obtained from Analysis B within less than 1σ   All values are also compatible within the 90% confidence interval with expectations assuming the three-flavor neutrino oscillation paradigm (i.e., ν τ normalization = 1.0) and the assumed ν τ CC cross sections.The significance at which we can reject the null hypothesis of no ν τ appearance is 3.2 σ and 2.1 σ for the CC+NC and the CConly case for Analysis A, respectively.The confirmatory Analysis B yields slightly weaker limits of 2.5σ (1.4σ).
The confidence intervals for the measurements presented here, shown in Fig. 17, are calculated using the approach of Feldman and Cousins [70] to ensure proper coverage.
The presented results are of a comparable precision to those of SK and OPERA (see Fig. 17), and complementary to those measurements in terms of energy scale, L/E range, systematic uncertainties, and statistics.Specifically, the SK measurement is based on lower-energy events where roughly 50% interact via CC quasi-elastic or resonant scattering, while the IceCube data are dominated by higher-energy events that interact primarily via the deep inelastic scattering interaction and are thus subject to different sources of neutrino interaction uncertainties [71].Additionally, the event samples used here are considerably larger than both OPERA and SK, with an estimated 1804 CC and 556 NC ν τ events for the final sample in Analysis A and 934 CC and 445 NC ν τ events in the final sample in Analysis B.
Determining the impact on tests of PMNS matrix unitarity requires global fits incorporating results from other experiments, as our result is only sensitive to the two elements U µ3 and U τ 3 of the matrix, while unitarity tests involve elements from a full row or column of the matrix.Also, as noted earlier, one could also use the measured ν τ normalization reported here along with the previously reported results from OPERA and SK to better constrain the CC ν τ cross section.
The measurement is limited by systematic uncertainties, in particular uncertainties in the initial flux of atmo-   spheric neutrinos and uncertainties in our detector model.Nevertheless, our result will improve with more statistics, as the aforementioned uncertainties are constrained by the data in the measurement itself-the increased sample size from more data allows us to control various detector effects and other sources of systematic uncertainties at a higher precision.
This defines a clear path forward towards a higher precision tau neutrino appearance measurement: more data, extended event selection and better control of detector uncertainties.With ten years of DeepCore data we expect an analysis similar to the one presented here to attain a precision of 15%.Better reconstruction algorithms-currently under development-promise to improve the precision, as do approved detector upgrades [72].The upgrades will include advanced calibration devices to improve our understanding of detector-related uncertainties, and the additional optical modules will be better and more efficient at identifying and reconstructing low energy neutrinos.These improvements will yield an anticipated precision of the tau neutrino normalization of better than 10% with a single year of operation.is expected to leave early hits in the veto region.Therefore, the variable VertexGuessZ provides a quick guess at an early stage in the event selection process for the z position of the interaction point.

Charge Variables
The event selection of Analyses A and B use the following charge information from a cleaned hit series to identify downward-going atmospheric muon events.
First, NAbove200 is the integrated charge from hit DOMs above z = −200 m and 2 µs before the Deep-Core trigger.Compared to an upward-going neutrino, a downward-going atmospheric muon is more likely to hit DOMs in the upper part of the detector and deposit a larger total charge in the veto region before trigger time.
The second charge variable is the charge ratio (QR) and has several distinct implementations.The variable QR6 is the fraction of accumulated charges from cleaned hits during the first 600 ns after the DeepCore trigger with respect to the total accumulated charge from all cleaned hits; that is, for 0 ns < t(Q i ) < 600 ns.A contained neutrino event deposits more charge in a shorter time scale than a throughgoing muon, so a QR6 value closer to zero indicates a potential atmospheric muon event.Similarly, the variable QR3 is calculated with a tighter time window of 300 ns instead of 600 ns.Further, to reduce the impact from noise hits, the charge ratio variables C2QR6 and C2QR3 are calculated as described above for QR6 and QR3, respectively, but with the two initial hits in the cleaned hit series removed.

Veto Regions
Veto variables are used in both analyses to identify and remove atmospheric muon backgrounds.These variables define veto regions, event-by-event, based on the positions and/or photon arrival times of the hit DOMs.
First used in [73] and optimized for DeepCore atmospheric oscillation searches in [74], the Veto Identified Causal Hits (VICH ) algorithm uses the cleaned hit closest to the trigger time as a reference hit and determines a veto region in which the hits may be causally connected with the reference hit.Fig. 18 shows the causal veto region (in red), which is defined by where ∆r and ∆t are the distance and photon arrival time between a given hit and the reference hit, respectively, and c is the speed of light in vacuum.The width of the veto region accounts for a reasonable amount of scattering in the ice combined with the typical time scale of a GeVscale neutrino event in IceCube (1 µs).Restrictions of the red region close to the time and position of the reference hit are made to allow variations between the interaction time, the first trigger, and subsequent hits that could occur in a neutrino event.With the defined veto region, the VICH variable is the integrated charge from cleaned hits that lie in the veto region.The more charge from cleaned hits deposited inside the veto region, the more likely that the event is caused by an atmospheric muon.
Another effective way to identify atmospheric muon background is to count the number of corridor DOMs.As shown in Fig. 1, the IceCube and DeepCore strings are arranged in a roughly triangular lattice in the horizontal plane.An atmospheric muon interacting in DeepCore can potentially come from a corridor (shown as the purple arrow in the figure) leaving no detected light in the veto region.To identify these muons, known corridor regions are studied.Given an interaction vertex, the corridor algorithm first finds direct hits from a cleaned hit series.Direct hits are hits due to minimally-scattered photons in the ice between their emission point and their detection in the DOM (the procedure to identify direct hits can be found in [18]).The algorithm then looks through the closest IceCube strings along the known corridors and counts the number of direct hits on those strings, which is defined as the number of corridor DOMs.

Center of Gravity
The Center of Gravity (CoG) is a parameter that measures space-time correlations between assumed signal hits in DeepCore and likely veto hits in the surrounding Ice-Cube DOMs.For a total of N hits, the average position x CoG is given by where x i is the position of the i th hit DOM relative to the center of IceCube.Then, an average CoG time t CoG is calculated by assuming a simple cascade hypothesis with light propagating from the CoG without scattering; Here, t i is the photon arrival time at the i th hit DOM, and c ice is the speed of light in ice.The same calculations can be applied to obtain the average position and time from a specific group of cleaned hits.
CoG is often used to check for causality.The CoG position and time are calculated from the cleaned hits inside the DeepCore fiducial volume.Then, for each cleaned hit in the veto region, a causally connected veto hit is identified if its x and t satisfy This requirement is used to reject muon tracks entering the fiducial volume after leaving hits in the veto region.
CoG is also used to evaluate event topology by dividing the time-sorted and cleaned hit series equally into four quartiles.Quartile 1 (Q1) consists of the earliest hits; quartile 4 (Q4) the latest hits.The CoG position and time are calculated for the hits in each quartile.The following three variables, based on these CoG quartiles, are used to identify background events.
First, to identify atmospheric muons, a z-travel variable is defined as the vertical distance between the CoG z position of Q1 and that from all cleaned hits in the fiducial volume; that is, z-travel ≡ z all − z Q1 .For a downwardgoing atmospheric muon event, its earlier hits tend to have a z Q1 position above the average z all , resulting in a negative z-travel.Similarly, an upward-going neutrino event usually has a positive z-travel.
Second, the spatial separation between Q1 and Q4 is also used for rejecting muon background events.The Q1-Q4 separation is defined as | x Q4 − x Q1 |.Because atmospheric muons usually travel long distances across the detector, they often have a larger spatial separation between their earlier and later hits compared to neutrino events.
Third, to discriminate noise triggers from physics triggers, the space-time interval ∆s 2 between Q1 and Q4 is also used.By definition, ∆s 2  2 .For an event caused by random detector noise, its ∆s 2 will appear either too time-like or too space-like compared to an event due to a neutrino interaction.
In addition to hit causality and event topology, a chargeweighted CoG can also estimate the size of an event by determining the standard deviations of z position (σ z ) and photon arrival time (σ t ) from all cleaned hits.An atmospheric muon tends to produce hits across the detector for a longer period of time, so its σ z and σ t are typically larger than for a neutrino event.

Quick Track Reconstructions
At lower selection levels with high event rates, computationally-inexpensive reconstruction algorithms are often used to provide a rough estimate of the event parameters related to the interacting particle.These parameters include the particle's speed, direction, and point of interaction.The following two quick algorithms assume a track event hypothesis and are used in both Analyses A and B.
The improved LineFit, or iLineFit, is based on the LineFit reconstruction.Assuming an infinitely-long muon track, the simple LineFit algorithm analytically minimizes a least-squares fit of cleaned hits with respect to the event parameters.The iLineFit then takes into account effects such as random detector noise and the scattering and absorption properties of the ice.The fitting procedure is described in [75].
The second reconstruction is a likelihood-based single photoelectron fit with eleven seeds (SPEFit11 ).Given a track-like event with a set of event parameters, the likelihood between the expected and observed photon arrival times is determined for each DOM.The total likelihood from all DOMs is minimized with respect to the event parameters.This fit runs iteratively from eleven different starting orientations to avoid falling into local minima.More information is found in [76].event selection uses a boosted decision tree (BDT) [43] to further reduce atmospheric muon background at Level 4. Six variables are used in training the BDT, including QR6, C2QR6, and NAbove200 discussed in Appendix A 2, as well as VertexGuessZ in Appendix A 1.
The remaining two variables are the estimated speed from a quick iLinefit reconstruction (see Appendix A 5) and the rough event topology.The concept of moment of inertia in classical mechanics is adapted to describe the overall shape of an event with N hits, each of which has a charge of q i .The diagonal elements of the moment of inertia are then given by where x i , y i , z i are the distance of the i th hit from String 36 at z = 0 m.For a spherically-shaped cascade-like event, the numerical values of I x,y,z are similar, while an elongated track-like muon background is more likely to have a smaller lateral distribution of hits than longitudinal.Therefore, a variable based on the tensor of inertia eigenvalue ratio, (ToIEVal ), is defined as: ToIEVal is included in training the Level 4 BDT, and its distribution is shown in Fig. 19.
The BDT score distribution is shown in Fig. 2, and a score cut is applied to keep events with a score above 0.04.Compared to Level 3, the muon background from MC estimates after the BDT score cut is reduced from 970 mHz to 51 mHz, while ≈ 55% of all neutrino events are kept.

Sample A, Level 5
A second BDT at Level 5 is trained to further reduce background muon contamination.After Level 4 cut is applied, over fifteen times more atmospheric muons and three times more pure-noise events remain compared to the integrated number of events from all neutrino flavors.Six variables are used for training the BDT, four of which are the radial position ρ of the FirstHLC (Appendix A 1), the Q1-Q4 separation and z-travel from the CoG algorithm (Appendix A 4), and the VICH veto variable (Appendix A 3).The distribution of VICH is shown in Fig. 20 as an example.The remaining two variables are the accumulated time and a simple zenith angle estimation.Because atmospheric muons are likely to travel across the detector, more time is needed, compared to neutrino events, for the hit DOMs to detect photons from the light source.Thus, the time to accumulate 75% of the charge from a cleaned hit series is included in the Level 5 BDT training.Further, because a majority of muon background events are down going, a fast reconstruction SPEFit11 (see Appendix A 5) is performed on a cleaned hit series to provide a rough estimate on the zenith angle of the interacting particle.
The Level 5 BDT score distribution is shown in Fig. 3, and a BDT score cut is applied to keep events with a score above 0.04.Compared to Level 4, the background muon rate after the Level 5 BDT score cut is reduced to ≈ 4 mHz, which is slightly more than a factor of 10.Moreover, ≈ 85% of noise-triggered events are rejected.

Sample A, Level 6
Previous selection criteria have reduced the background rate to ≈ 6 mHz, which is comparable to the total neutrino rate of > 2 mHz.At this stage, most obvious background muons are rejected, and the remaining atmospheric muons are more difficult to be identified.Therefore, Level 6 includes two straight cuts to reject events caused by random detector noise and two more cuts to identify sneaky atmospheric muons.
To identify events caused by random detector noise, an algorithm called Fill-Ratio is performed to look for the topology of hits in an event around an estimated vertex position given by FirstHLC (see Appendix A 1).A sphere is defined around the vertex with a radius 60% larger than the average distance from all cleaned hits.The pattern of hits in an event can be estimated with a Fill-Ratio variable defined by the ratio of number of hit DOMs inside the sphere to the total number of DOMs inside the sphere.Since hits in a pure-noise event are randomly scattered across the detector, such event tends to have a longer mean distance away from its vertex and less hit DOMs inside the sphere.Therefore, pure-noise events often have smaller values of Fill-Ratio compared to physics events, as shown in Fig. 21.A straight cut is placed at 0.05, events above which are kept.
A second cut placed on NChannel also help further remove noise-triggered events.NChannel is the number of hit DOMs in a cleaned hit series, and its distribution is shown in Fig. 22.Most noise-triggered events have less than eight hit DOMs.Further, the computationallyintensive, likelihood-based reconstruction algorithm described in Sec.III C fits eight parameters.In order to have enough degrees of freedom for a reasonable fit, a direct cut is applied on NChannel to remove events with less than eight hit DOMs, reducing the contribution from pure-noise events by a factor of 20.
For atmospheric muons, two straight cuts are applied.The first cut is applied on the number of corridor DOMs as explained in Appendix A 3. Events with two or more direct hits found along those corridors are rejected.The second cut involves a likelihood-based algorithm known as FiniteReco.Given an infinite track event hypothesis, FiniteReco reconstructs an event based on the probabilities of the individual DOMs to see a hit or not [74].It provides a relatively quick estimate on the interaction vertex of an event.Figure 23 shows the fractional 2D distribution of radial ρ and depth z positions from FiniteReco Color axis represents the fraction of atmospheric muon with respect to the total expected MC events per bin.Red lines represent the cut, within which events are kept.
96% of atmospheric muon and noise-triggered events are rejected compared to Level 5. Neutrinos contribute more than 66% of the sample, which opens the opportunity to use a computational-expensive but more comprehensive reconstruction algorithm discussed in Sec.III C.

Sample A, Level 7
The comprehensive reconstruction discussed in Sec.III C is performed between Level 6 and 7. To reduce the remaining background events and to improve agreement between data and MC, final selection cuts are applied at Level 7.
The first selection criteria is similar to the containment requirement from the FiniteReco vertex positions at Level 6.With a more sophisticated reconstruction method described in Sec.III C, the fitted interaction vertex of an event is better than the estimates from FiniteReco. Figure 24 shows the fractional 2D distributions of radial ρ and depth z positions obtained from the final reconstruction.An extra containment condition is added to exclude events below z position of −500 m, and events outside the red lines are rejected.
Two final cuts are applied during the development of the A selection.The first cut is based on a 2D distribution of reconstructed energy per number of hit DOMs and the root mean square (RMS) of photon arrival times from a cleaned hit series.The second cut is related to 'flaring' DOMs, which emit light sporadically.The extra light is not simulated in MC, and a disagreement between data and MC is shown from the distribution of normalized RMS of total charges (see Fig. 25).Therefore, a cut is placed to remove events in which the normalized RMS of total charges are above 85%.The above two cuts only remove 5% of the total events and do not alter any physics results.Nonetheless, they are applied.
Finally, only events within the analysis histogram ranges stated in Sec.IV are used for the oscillation analyses.This restricts to events from all sky with a reconstructed energy between 5.6 GeV to 56 GeV and a reconstructed track length less than 1,000 m.In summary, the event rates as a function of event type and selection level is shown in Table III.The neutrino rates are the combination of the NC+CC channels and use the atmospheric neutrino flux predictions from [26] with values of θ 23 and ∆m 2 from [78].At the final level, 92% of the A sample is neutrino events, while the contamination from atmospheric muon and noise-triggered events are 8% and 0.1% respectively.

Appendix D: Higher Level Selection for Sample B
This section focuses on the event selection method for analysis B. It features a set of straight cuts and a boosted decision tree to improve the purity of neutrino events at the final selection level.Based on simulations, about 40,000 neutrinos are expected given three years of detector exposure.The resultant sample is also used for the most recent published measurements of atmospheric neutrino oscillation parameters from the IceCube Collaboration [19].After the common Level 3 filtering discussed in Appendix B 2, a set of straight cuts is applied at Level 4 to further remove events due to random detector noise and atmospheric muons.These selection requirements rely on hit information from a cleaned hit series, charge information in veto regions, and an estimated interaction vertex position of an event.
To ensure that enough information are detected in an event, the first two selection variables are NChannel and RTFiducialQ.NChannel is the number of hit DOMs in a cleaned pulseseries, and only events with at least eight cleaned hits are kept.Further, RTFiducialQ is a chargerelated variable from an algorithm, which searches for clusters of cleaned hits in the DeepCore fiducial region that satisfy space-time correlations.Events with a minimum of one cluster with at least 7 p.e. are kept.
Three additional selection criteria are placed based on variables from the Center of Gravity (CoG) algorithm discussed in Appendix A 4. First, the space-time interval ∆s 2 has some power to distinguish events caused by random detector noise from physics events.Therefore, a cut is applied to keep events where ∆s 2 is between −(400 m) 2 and 0 m 2 .The remaining two cuts depend on the size of an event, which is estimated by the charge-weighted spread of vertex z position (σ z ) and photon arrival time (σ t ) from a cleaned hit series.To be specific, only events with σ t ≤ 1,000 ns and 7 m ≤ σ z ≤ 100 m are kept.
Because an atmospheric muon event can be identified using the charge information in the veto regions, two veto requirements are placed.The first variable counts veto charges using the CoG algorithm (see Appendix A 4).This cut is very similar to the total veto charge requirement at Level 3 where the veto region is outside the extended DeepCore volume.At Level 4, the same algorithm is applied to the veto region outside the standard DeepCore fiducial volume defined in Fig. 1.With a slightly larger veto volume, a tighter cut is applied at Level 4 to remove events with a total veto charge greater than 5 p.e.The second veto charge requirement is applied on the VICH variable discussed in Appendix A 3. With a veto volume defined by the estimated point of interaction, VICH looks for potentially causally-related hits in the event-by-event veto region.A cut is placed at 7 p.e. to reject potential atmospheric muon events.
To increase the purity of ν µ CC events, a cut is applied on the number of direct hits.Direct hits are hits due to photons that experience minimal scattering in the ice between its emission point and detection in the DOM.When a muon track passes next to a string, the intersection of its Cherenkov cone with a string results forms a hyperbolic pattern as a function of the photon direct arrival depth and times.The orientation of this pattern is determined by the angle between the string and the passing muon track.The procedure to identify direct hits is explained in [18], and only events with at least three direct hits are kept.The last four BDT variables are based on the CoG algorithm discussed in Appendix A 4. One of them is the separation between the first and the last quartiles of a cleaned hit series, or Q1-Q4 separation.As shown in Fig. 26, atmospheric muon background tends to have a longer spatial distance between the two quartiles, compared to neutrino events.The remaining three variables are reused from the previous level.They are the chargeweighted spread of vertex z position (σ z ) from all cleaned hits and the estimated radial ρ Q1 and depth z Q1 positions of the interaction vertex from the first quartile (Q1) of CoG.
Figure 27 shows the BDT score distribution, and a cut is applied to accept events with a score above 0.2.Compared to Level 4, 99.9% of the atmospheric muon background is removed, whereas 58% of all neutrinos is kept after the BDT score cut is placed.

Sample B, Level 6
In between Level 5 and 6, the comprehensive reconstruction discussed in Sec.III C is performed.At Level 6, two final selection requirements are placed to further improve the quality of the final sample.
First, final containment criteria is required based on

AbsorptionFIG. 1 .
FIG.1.Top and side views of IceCube indicating the positions of DeepCore DOMs with red circles and surrounding IceCube DOMs with green circles.The DeepCore fiducial region is shown as a green box at the bottom center.The Deep-Core DOMs were deployed mostly > 2100 m below the surface (shown highlighted in green) with some DeepCore DOMs also deployed around 1800 m below the surface (shown highlighted in red) to aid in rejection of atmospheric muons.The bottom left of the plot shows the absorption length for Cherenkov light vs. depth.The purple arrow in the top view shows one example of a "corridor" path along which atmospheric muons can circumvent the simple veto cuts, as they may not leave a clearly detectable track signature (see Sec. III B for details).The gray band indicates the dust layer, a region of higher scattering and absorption.

FIG. 5 .
FIG.5.Reconstructed energy vs. true energy for each neutrino flavor separately (CC interactions) and all flavors combined (NC interactions).The red and blue solid lines are the resolutions from Analyses A and B respectively, and the dashed lines represent the 68% ranges.The solid black lines are the references indicating perfect reconstruction.For ντ CC and ν NC events the final state ensemble of out-going particles include at least one "invisible" neutrino which manifests as missing energy when comparing Etrue to Ereco.

FIG. 6 .
FIG.6.Difference between reconstructed and true cos θ vs. true energy for each neutrino flavor separately (CC interactions) and all flavors combined (NC interactions).The red and blue solid lines are the resolutions from Analyses A and B respectively, and the dashed lines represent the 68% ranges.The solid black lines are the references indicating perfect reconstruction.

FIG. 7 .
FIG.7.Fraction of track-like events as a function of true neutrino energy for each neutrino event type in Analyses A (left) and B (right).Differences in particle classification lead to different fractions of track-like events at lower energies.

FIG. 8 .
FIG. 8. Expected signal ντ (CC+NC) divided by the squareroot of the expected background (νe, νµ, atmospheric µ, and noise-triggered) events as a function of reconstructed cosine of the zenith angle and reconstructed energy.Cascade-like events are shown on the left and track-like events on the right.The plots include both neutrinos and anti-neutrinos.

FIG. 10 .
FIG. 10.Effect of selected systematic uncertainties on the nominal event distribution shown as a percentage change of the expectation per bin.With cascade-like events on the left and track-like events on the right, shown from top to bottom are: νe/νµ flux ratio at +1σ, ν/ν flux ratio at +1σ, head-on optical efficiency at +1, ∆m 2 32 at 2.778 × 10 −3 eV 2 instead of 2.526×10 −3 eV 2 , and M res A at +1σ. (See text for definitions of these parameters.)A complete collection of plots is provided in Appendix E.

FIG. 11 .
FIG. 11.Relative acceptance of photons versus photon arrival angle for different optical models of the ice.Zenith angles θ with cos θ = −1.0indicate vertically downward-going photons (hitting the top of a DOM), cos θ = 0.0 horizontal photons, and cos θ = −1.0vertically upward-going photons.The black lineshows the angular photon sensitivity of a module as measured in the laboratory.The green line and surrounding green band show the angular acceptance used and its uncertainty, respectively, and is based on two parameters, lateral and head-on sensitivity.The head-on area has a large associated uncertainty.Data points obtained from the direct simulation of a bubble column (not based on angular acceptance) are overlaid in blue.

FIG. 12 .
FIG.12.Event distributions of the atmospheric muon background for Analysis A (top row) obtained from the best-fit simulation, and for Analysis B (bottom row) obtained from the data sideband.

2 23 and sin 2 θ
23 are measured as a cross-check of the validity of Analysis A presented earlier.With the ν τ normalization fixed to 1, all sources of systematic uncertainties listed in Table II are taken into account.With 140 non-zero bins and 133 effective degrees of freedom, a χ 2 defined in Eq. 5 of 129.4 is obtained when letting all 16 nuisance and two oscillation parameters float.The best fit values of ∆m 2 23 and sin 2 θ 23 are 2.55 +0.12 −0.11 × 10 −3 eV 2 and 0.58 +0.04 −0.13 , respectively.

2 FIG. 13 .
FIG. 13.The 90% allowed region using the data sample from Analysis A in blue compared to other experiments[19, 65-68].The best fit point from Analysis A is shown as the blue cross mark.The IceCube 2017 result [19], represented in black, uses the data sample from Analysis B. The top and right plots are the 1-d ∆χ 2 profiles of the measured oscillation parameters.

FIG. 14 .
FIG.14.Distributions of the data with best-fit neutrino and muon backgrounds subtracted and signal simulation.Statistical errors are shown for the data.The best-fit hypothesis shows good agreement in the reconstructed energy axis (left), the cosine of the reconstructed zenith angle (middle) and PID categories (right) for Analysis A.

FIG. 15 .
FIG.15.Distribution of the data as a function of reconstructed L/E, overlaid with the best fit neutrino and cosmic-ray muon histograms for Analysis A (top) and B (bottom).The bottom portion of each shows the ratio of the data to the predicted distribution at the best fit point, with black points representing data and the height of the shaded band the uncertainty of the best fit (statistical errors only).

FIG. 16 .
FIG. 16.Observed ∆χ 2 from the best fit CC+NC (CC) ντ normalization of 0.75 (0.62) as a function of the ντ normalization (black lines).Shaded bands show the 68% ranges of the expected distribution of ∆χ 2 values obtained from pseudo-experiments assuming nominal values for oscillation parameters and a ντ normalization of 1.0.
ment 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 clus-ter of the RWTH Aachen; Sweden -Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; 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, Danish National Research Foundation (DNRF), Carlsberg Foundation; 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).

FIG. 18 .
FIG.18.Definition of veto region based on the VICH algorithm.The four lines surrounding the red region define the causal veto volume, event-by-event, based on the hit closest to the trigger time.The more total charge from cleaned hits deposited inside the veto region, the more likely that the event is caused by an atmospheric muon.

FIG. 22 .FIG. 23 .FIG. 24 .
FIG. 22. (top) NChannel distribution at Level 6 (before cut applied).Each shaded color represents the stacked histogram from each event type.Black dots represents the data distribution.The vertical line is the cut value of 8, events below which are rejected.MC events are weighted by world averaged best fit oscillation parameters.(bottom) Ratio of distribution from data to that from MC. Black error bars are the statistical fluctuation from data, whereas shaded red areas are the uncertainties from limited MC statistics.

5 FIG. 25 .
FIG. 25. (top) distribution of normalized RMS of total charge at Level 7 (before cut applied).Each shaded color represents the stacked histogram from each event type.Black dots represent the data distribution.The vertical line is the cut value of 0.85, events above which are rejected.MC events are weighted by world averaged best fit oscillation parameters.(bottom) Ratio of distribution from data to that from MC. Black error bars are the statistical fluctuation from data, whereas shaded red areas are the uncertainties from limited MC statistics.

FIG. 27
FIG. 27. (top)BDT score distribution at Level 5 cut applied).Each shaded color represents the stacked histogram from each event type.Black dots represents the data distribution.The blue vertical line is the cut value of 0.2, events below which are rejected.MC events are weighted by world averaged best fit oscillation parameters.(bottom) Ratio of distribution from data to that from MC. Black error bars are the statistical fluctuation from data, whereas shaded red areas are the uncertainties from limited MC statistics.
parton distribution functions are used

TABLE I .
Expected number of events at the NC+CC best fit point, grouped by flavor and interaction type, and including atmospheric muons.The observed counts from the data are shown in the last row.Associated ±1σ uncertainties due to limited simulation statistics are also shown (the uncertainty showed on the observed count is just the Poisson error).
FIG. 17.The measured values for CC+NC and CC-only results in both analyses.Also shown are previous best-fit values of the CC-only ντ normalization from OPERA and SK, which were performed with different energy ranges and fluxes and a different definition of the ντ normalization from those used in IceCube.All measurements of tau neutrinos are consistent with standard oscillations (ντ normalization of 1.0), with the two analyses presented here showing excellent internal agreement.

TABLE III .
The event rate in mHz for the common filtering and the subsequent event selection levels for Analyses A and B, respectively.After the final selection level, the analyses only include events with energies in the region from 5.6 GeV to 56 GeV which is denoted as "LE".
FillRatio distribution at Level 6 (before cut applied).Each shaded color represents the stacked histogram from each event type.Black dots represent the data distribution.The vertical line is the cut value of 0.05, events below which are rejected.MC events are weighted by world averaged best fit oscillation parameters.(bottom) Ratio of distribution from data to that from MC. Black error bars are the statistical fluctuation from data, whereas shaded red areas are the uncertainties from limited MC statistics.