A Measurement of the Tau Neutrino Cross Section in Atmospheric Neutrino Oscillations with Super-Kamiokande

Using 5,326 days of atmospheric neutrino data, a search for atmospheric tau neutrino appearance has been performed in the Super-Kamiokande experiment. Super-Kamiokande measures the tau normalization to be 1.47$\pm$0.32 under the assumption of normal neutrino hierarchy, relative to the expectation of unity with neutrino oscillation. The result excludes the hypothesis of no-tau-appearance with a significance level of 4.6$\sigma$. The inclusive charged-current tau neutrino cross section averaged by the tau neutrino flux at Super-Kamiokande is measured to be $(0.94\pm0.20)\times 10^{-38}$ cm$^{2}$. The measurement is consistent with the Standard Model prediction, agreeing to within 1.5$\sigma$.


I. INTRODUCTION
In the three-flavor neutrino framework, the three neutrino flavor states (ν e , ν μ , ν τ ) are superpositions of three neutrino mass states (ν 1 , ν 2 , ν 3 ). The oscillation parameters in the framework have been measured in atmospheric neutrino experiments [1][2][3], solar neutrino experiments [4][5][6], reactor neutrino experiments [7][8][9][10], and long-baseline neutrino experiments [11][12][13]. Atmospheric neutrino observations are characterized by a large deficit of muon events. The deficit is largely explained by the quantum mechanical mixing of the propagating mass states such that weak interaction at the detector is comprised of a mixture of muon and tau flavors, whereas most of the tau neutrino flux has energy below the tau lepton production threshold. The objective of this paper is to observe those tau neutrino interactions that are above that threshold. A direct detection of tau neutrinos from neutrino oscillation is important for an unambiguous confirmation of three-flavor neutrino oscillations.
However, the detection of tau neutrino appearance is challenging. Charged-current neutrino interactions are required to determine the flavor in neutrino detection. Charged-current tau lepton appearance has an energy threshold of 3.5 GeV, and the charged-current tau neutrino cross section is greatly suppressed at low energies due to the large mass of the tau lepton relative to the electron or muon. The DONUT experiment first directly observed the tau neutrino by measuring charged-current interactions using a high-energy neutrino beam that contained tau neutrinos [14]. Long-baseline experiments tuned for maximum oscillation have the bulk of their neutrinos below this energy. In addition, the tau lepton has an extremely short lifetime, making a direct detection very difficult. Nevertheless, the long-baseline neutrino experiment OPERA measured tau neutrino appearance in a highenergy muon neutrino beam by observing five ν τ events with a background expectation of 0.25 events [15].
Atmospheric neutrinos are mostly electron or muon neutrinos at production [16]. Tau neutrino appearance is expected in the atmospheric neutrinos from neutrino oscillations. In three-flavor neutrino oscillation in the vacuum, the probability of ν τ appearance can be approximately expressed as P ν μ →ν τ ≃ cos 2 θ 13 sin 2 ð2θ 23  where Δm 2 32 ≡ m 2 3 − m 2 2 is the mass splitting in eV 2 , θ ij is a mixing angle in the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, L is neutrino path length in km, and E is neutrino energy in GeV. Atmospheric neutrinos have energies spanning many orders of magnitude from 10 MeV to more than 1 TeV; the high energy component of the atmospheric neutrinos have enough energy for charged-current tau neutrino interactions. Super-Kamiokande is expected to detect roughly one charged-current tau neutrino interaction per kiloton of water per year. Super-Kamiokande previously published a measurement of atmospheric tau neutrino appearance consistent with three-flavor neutrino oscillation with data collected in SK-I through SK-III [17]. This analysis has been updated with data collected in SK-IV between 2008 and 2016, and the simulation and reconstruction have been improved. Using the measured charged-current tau neutrino events, Super-Kamiokande also measures the charged-current tau neutrino cross section.
This paper proceeds as follows: Section II describes some basic features of the Super-Kamiokande experiment (Super-K, SK). Section III describes Monte Carlo simulations of both the charge-current tau neutrino signal and the atmospheric neutrino background. Section IV describes standard data selection and reconstruction algorithms used in the analysis. Section V describes a neural network algorithm developed to select the tau signal. Section VI describes a search for atmospheric tau neutrino appearance, and a measurement of charge-current tau neutrino cross section. Section VII presents our results and conclusion.

II. THE SUPER-KAMIOKANDE DETECTOR
Super-Kamiokande is a 50 kiloton cylindrical water-Cherenkov detector located in the Kamioka mine under about 1 km rock overburden (2.7 km water equivalent) at the Ikenoyama mountain in Japan [18,19]. The detector is arranged into two optically separated regions: the inner detector (ID) and the outer detector (OD). The ID is instrumented with 11 129 20-inch inward-facing photomultiplier tubes (PMTs) and the OD is instrumented with 1,885 8-inch outward-facing PMTs. The PMTs collect Cherenkov light produced in the ultrapure water in the detector. A fiducial volume of the ID is defined as the cylindrical volume 2 meters inward from the ID wall, and has a mass of 22.5 kilotons [18].
Super-K has been in operation since 1996, and has had four data-taking periods.  [20]. The period after the upgrade is referred to as SK-IV. In this paper, SK-IV data are used up to March 2016, totaling 2519.9 live-days. The complete SK-I through SK-IV data set comprises a total exposure of 5,326 live-days.

III. SIMULATION
In order to predict the rate of both tau signal and atmospheric neutrino background, a full Monte Carlo (MC) simulation is used to model both the neutrino interactions and the detector response of Super-K. Since the four Super-K periods have different detector configurations, separate sets of MC for both tau signal and atmospheric neutrino background are generated for each period.
Atmospheric neutrinos are produced from the decays of charged mesons and muons in the cosmic-ray induced atmospheric showers, and are mostly ν μ and ν e at production. The intrinsic tau neutrinos in the atmospheric neutrino flux are negligible for this analysis [21]. Three-dimensional neutrino fluxes of ν μ and ν e are modeled from the calculation of Honda et al. [16]. The calculation predicts the fluxes of electron and muon neutrinos as a function of neutrino direction and neutrino energy at the Super-K site.

A. Neutrino fluxes
Although atmospheric neutrinos consist of ν μ and ν e at production, ν τ are expected to appear due to neutrino oscillations. The probabilities of ν τ appearance from neutrino oscillations of ν μ or ν e in vacuum are shown in Eq. (1).
However since neutrinos coming from below travel through the Earth, the oscillation probabilities are altered by the matter effect. Therefore, a custom code [22] is used to calculate the oscillation probabilities, which takes into account the effect of neutrino types, path lengths, neutrino energies and the matter effect. The oscillation parameters used are Δm 2 32 ¼ 2.1 × 10 −3 eV 2 , Δm 2 21 ¼ 7.6 × 10 −5 eV 2 , sin 2 2θ 23 ¼ 1.0, sin 2 2θ 13 ¼ 0.099, δ CP ¼ 0 [23]. A method from [24] is used to account for the matter effect in the calculation of oscillation probabilities based on the matter density structure of the Earth in Ref. [25]. Figure 1 illustrates the tau neutrino appearance probabilities from muon neutrino or electron neutrino oscillations in three-flavor neutrino oscillation under the assumption of the normal hierarchy. For a neutrino with a given energy and path length, a muon neutrino has a larger probability than a electron neutrino to be detected as a tau neutrino. Following the oscillation calculation, we can predict the atmospheric tau neutrino flux at Super-K. Figure 2 shows the expected atmospheric tau neutrino fluxes from neutrino oscillations at Super-K.

B. Neutrino interactions
The NEUT code [26] is used to model the neutrino nucleon interactions including quasielastic scattering, single meson production, coherent pion production, and deep-inelastic scattering (DIS). In the simulation of atmospheric neutrino background, all ν μ and ν e interactions are included. All flavors of neutrinos interact with neutral-current (NC)  interactions, and hence are unaffected by oscillations. Atmospheric neutrino neutral current events are simulated based on the total neutrino flux. The simulation of tau signal contains only charged-current (CC) ν τ interactions whose cross sections are calculated following the same models as those used for ν μ and ν e . The relatively large mass of the tau lepton produced in the interactions greatly suppresses the cross section of charged-current tau neutrino interactions at low energies and results in an energy threshold of 3.5 GeV. Figure 3 shows the total cross section of charged-current interactions for ν τ andν τ in the simulations. Tau leptons produced in the CC tau neutrino interactions are polarized, and the polarization affects the distributions of its decay particles. Therefore, a polarization model from Ref. [27] is also included in the simulation. Figure 4 shows the polarization of τ − =τ þ in the simulation for interactions of neutrinos with energy of 10 GeV. This analysis selects events at relatively high neutrino energies, at which the CC interactions contain a high percentage of DIS (45%) in the background, with CC ν τ events containing 60% DIS. The GRV98 [28] parton distribution functions are used in the calculation of the DIS cross sections. In order to smoothly match the DIS cross sections with the resonance region, an additional correction developed by Bodek and Yang [29] is also applied.
The tau lepton has a mean lifetime of 2.9 × 10 −13 s and it decays very quickly after production in the detector. The decays of tau lepton are divided into leptonic decays and hadronic decays based on the particles produced. The kinematics of tau lepton decay is simulated with TAUOLA version 2.6 [30].
The particles produced in both atmospheric neutrino background interactions and CC ν τ interactions are input to a custom detector simulation based on GEANT3 [31]. The code simulates the propagation and Cherenkov light emission of the particles and the Super-K detector [19].

IV. REDUCTION AND RECONSTRUCTION
This analysis only uses fully contained (FC) multi-GeV events in the fiducial volume. Fully contained events are defined as events which only have activity in the ID, and FC events in the fiducial volume are selected by requiring the reconstructed event vertex be at least 2 meters away from the ID wall. In addition, the events are required to have more than 1.3 GeV of visible energy (E vis ), which is defined as the energy to produce the observed light in the event if it were produced by a single electron. As shown in Fig. 5, the E vis cut selects the majority of the tau signal but rejects the bulk of low-energy atmospheric neutrino background events. The selection efficiencies for this set of cuts are 86% for the ν τ CC signal and 23% for the background in four Super-K periods. The same selection is applied to the simulations and the observed data.  The selected events are passed through a reconstruction program to determine the event vertex, the number of Cherenkov rings, the particle type and momentum of each ring, and the number of Michel electrons. The same reconstruction algorithms are applied to the MCs and the observed data. Events are assumed to originate from a single vertex, and the reconstruction uses the distribution of observed charge and the PMT timing to find the vertex and the brightest Cherenkov ring. A Hough transformation method [32] is used to find additional rings. Each ring candidate is tested using a likelihood method to remove fake rings and determine the final number of rings. A likelihood method based on the ring pattern and ring opening angle is used to identify each ring as e-like (showering type from e AE or γ) or μ-like (nonshowering type). Michel electrons from stopping muons are tagged by searching for clusters of hits after the primary event. The time window for such clusters extends to 20 μs after the primary event. In the SK-I to SK-III periods, there was an impedance mismatch in the electronics which caused signal reflection around 1000 ns after the main event. Therefore, the time period 800-1200 ns after the main event was excluded. The improved SK-IV electronics avoids such signal mismatch, thus no exclusion is required. As a result, the tagging efficiency was improved from 80% to 96% for μ þ decays and 63% to 83% for μ − decays between SK-I-II-III and SK-IV. More details of the reconstruction can be found in Refs. [33,34].

V. A NEURAL NETWORK ALGORITHM FOR TAU NEUTRINO SEARCH
As described in Sec. III, tau leptons produced in CC ν τ interactions decay quickly to secondary particles. Because of the short lifetime of tau lepton, it is not possible to directly detect them in Super-K. The decay modes of the tau lepton are classified into leptonic and hadronic decay based on the secondary particles in the decay. The leptonic decays produce neutrinos and an electron or a muon. These events look quite similar to the atmospheric CC ν e or ν μ background. The hadronic decays of the tau are dominant and produce one or more pions plus a neutrino. The existence of extra pions in the hadronic decays of tau allows the separation of the CC ν τ signal from CC ν μ , CC ν e and NC background. As shown in Fig. 6, CC ν τ events typically produce multiple rings in the detector. Multiplering events are relatively easy to separate from single-ring atmospheric neutrino events. However, the multiring background events, resulting from multipion/DIS atmospheric neutrino interactions, are difficult to distinguish from the tau signal. Simple selection criteria based on kinematic variables do not identify CC ν τ events efficiently. In order to statistically identify events with the expected characteristics that differentiate signal and background, a multivariate method is applied in this analysis. Specifically, a multilayer perceptron (MLP) method is used. It is implemented in the ROOT-based TMVA [35] library, and was also used in our previously published ν τ search [17].
A multilayer perceptron is a feed-forward artificial neural network (NN), which maps between a set of inputs and a set of outputs. It is typically organized in layers of interconnected neurons with one or more hidden layers between the input and output layer. Neurons in the input layer receive inputs, then normalize the inputs and forward them to the neurons in the first hidden layer. Each layer is fully linked to neurons in the next one with weighted connections. The output of a neuron is scaled by the connecting weights and fed forward to the neurons in the next layer. A MLP has the ability to learn through training, during which the weights in the network are adjusted. Once trained with representative training data, the MLP can be applied to new, unseen data.
A MLP is used in this analysis, which has seven inputs, ten neurons in one hidden layer, and one output. It takes seven input variables for both the CC ν τ signal and atmospheric neutrino background to produce a single discriminating output variable that separates signal and background. To prepare the MLP algorithm for event identification, three stages are required. They are referred to as training, testing, and analysis. Separate sets of signal and background MC are used in each stage. We describe the MLP that we implemented for this analysis below.
Seven variables are used as inputs to the MLP based on the expected separation between signal and background in these variables. visible energy than background. The E vis spectrum of the CC ν τ signal peaks around 4 GeV, as shown in Fig. 7. By contrast, the E vis spectrum of the background falls with increasing E vis . (2) The particle identification likelihood parameter of the ring with maximum energy. Tau leptons decay quickly to daughter particles after production through leptonic and hadronic decays. Except for the leptonic decay to a muon, most decay channels have at least one showering particle. A showering particle has a negative value in the definition of particle identification likelihood, compared with a positive value for a nonshowering particle. The particle identification of the most energetic ring for the signal has a distribution mostly in the negative region, while the background has a broad distribution in both negative and positive regions. (3) The number of decay electron candidates in the event. Naively, we expect more decay electrons for signal from pion decays which are produced in hadronic tau decays. This variable does not depend on ring reconstruction, so it is relatively independent of most other variables. (4) The maximum distance between the primary interaction point and any decay electron from a pion or muon decay. Energetic muons can travel a long distance in water. Therefore, CC ν μ background involving a high energy neutrino is expected to have a large distance between the primary interaction point and the decay electron from the muon. In comparison, the pions from hadronic tau decay are expected to have smaller momentum, resulting in a smaller value of the variable. (5) The clustered sphericity of the event in the center of mass system. Sphericity is a measure of how spherical an event is. A perfectly isotropic event has sphericity 1, while a perfectly one-directional event has sphericity 0. We follow the definition from [36], defining the spherical tensor as where α,β ¼ 1, 2, 3 are three Cartesian momentum vectors pointing to binned photoelectric charge in the event. Sphericity is then constructed by finding the eigenvalues, λ 1 > λ 2 > λ 3 , of the tensor: The hadronic decay of the heavy tau lepton is more isotropic than a typical ν μ or ν e background. The spectrum of sphericity is centered near 0.8 for signal, The background event with multirings has a similar pattern to the signal event, and requires more effort to statistically distinguish.
while the spectrum for background has an almost flat distribution between 0.1 and 0.8. (6) The number of possible Cherenkov ring fragments.
In the ring reconstruction, these ring candidates are formed using a method based on a Hough transformation to find rings. We expect more ring candidates for signal because of the multiple charged particles and pions in hadronic tau decay. This variable is sensitive to even partial ring fragments. (7) The fraction of the total number of photoelectrons carried by the most-energetic ring in an event. This variable is calculated from the number of photoelectric charge in each PMT (q i ) and the reconstructed vertex and direction of an event as rfrac ¼ where θ i is the angle between the reconstructed direction of the first ring and the direction of the reconstructed vertex to a PMT. The variable calculates the ratio of charge within 48°of the direction of the first ring in the event. The variable is expected to be smaller for the signal because energy is carried by multiple particles in the hadronic decay of the tau. Because the oscillation-induced tau events come from below, the downward sample is expected to have no tau events. Therefore, the data in the downward sample can be used to study the extent to which the atmospheric neutrino simulation for background events agrees with data. Figure 7 shows the seven variables for data and MC in the downward sample, along with the expected tau signal. The data and background MC have good agreement.
The neural network is trained with a set of signal and background MC with the target of output ¼ 1 for signal and 0 for background. The weights in the MLP are adjusted iteratively during the training such that the difference between the actual output of the MLP and the target is minimized.
The MC is generated with a Honda flux that calculates the fluxes of neutrinos at production in the atmosphere [16]. Therefore, the events need to be weighted with oscillation probabilities before being fed to the MLP. Tau neutrinos from oscillations are expected to mostly come from below because the oscillation length of neutrinos in excess of the tau threshold is at least 4,100 km. Used naively, the oscillation weight described in Sec. III will encode this up-down asymmetry into the signal MC. To avoid training the neural network to select signal events based on event direction, the training sample is instead weighted based on the average probability at its visible energy, as shown in Fig. 8. In other words, instead of weighting each MC event with its oscillation probability, the weight is calculated as an average of the oscillation probabilities for MC events in each bin of log 10 ðE vis Þ. In this way, the upward and downward samples are treated the same in the training process, while the overall weight is still correct. Moreover, since the weights of the downward signal simulation are not set to zero, the whole of signal MC statistics are preserved for training.
During training, a testing data set is used as validation to avoid overtraining. Figure 9 shows the neural network output for background and signal with the training and testing samples. The clear separation of signal and background in both samples demonstrates that the MLP learned to separate signal from background. Also, the good agreement between training and testing samples shows that it is properly trained.
The testing sample is also used to plot the efficiencies of signal selection and background rejection by cutting on NN output, as shown in Fig. 10. By varying the cut on NN output, the efficiencies of signal selection and background rejection can be changed. When selecting tau-like events from the events after reduction by requiring the NN output be greater than 0.5, 76% of the signal events and only 28% of the background remain. Table I summarizes the breakdown of the interaction modes in different samples, including the fraction for tau and non-tau-like samples. Table II summarizes the decay modes of the largest branching fractions and the fraction of tau-like events in each mode. These efficiencies are only shown to assess the  performance of the neural network in selecting tau events.
No cut is used in the following analysis. The analysis sample is finally processed with the trained MLP. Unlike the training and testing processes, no information is given to the MLP regarding the composition of the analysis samples as either signal or background. The analysis sample processed with the trained MLP is used in this analysis.

A. Search for atmospheric tau neutrino appearance
To search for atmospheric tau neutrino appearance, the data is fit to a combination of the expected tau signal resulting from neutrino oscillations and atmospheric neutrino background with neutrino oscillations. In order to extract maximum information from the sample, the analysis uses a two-dimensional unbinned maximum likelihood fit implemented in ROOFIT [38]. Using two-dimensional histograms of the neural network output and the reconstructed zenith angle of the events, two-dimensional probability distribution functions (PDFs) are built for background and tau signal. The probability density follows the normalized bin contents in the histograms. Figure 11 is an example of the 2D distributions for oscillated signal on the top and background on the bottom. The horizontal axis of the plots is the cosine of the reconstructed zenith angle as determined by the energy-weighted sum of the ring directions in the event. The vertical axis is the NN output, in which tau-like events have a value close to 1 and non-tau-like events have a value close to 0. The signal events (top panel) are primarily tau-like and come from below (cosine of the zenith angle, Θ, less than zero), while the background (bottom panel) is more non-taulike and come from all directions. The amount of signal and background events can be adjusted by varying the normalization of the distributions. Figure 12 shows a combination of signal PDFs and background PDFs for SK-I to SK-IV with both tau normalization and background normalization equal to 1, with the data overlaid on the combined PDF.
The systematic errors used in this analysis are selected from the systematic errors in the Super-K atmospheric  neutrino analysis [17]. Only systematic errors for which the maximum bin content changes in the two-dimensional signal/background histograms after shifting that systematic error by 1σ is larger than 2.5% are considered. After the reevaluation, 28 systematic errors from atmospheric neutrino analysis are considered in the fit. The systematic errors are summarized in Tables IV and V. In order to simultaneously fit the systematic errors, a set of PDFs are built for each systematic error with the same structure as the PDFs for the signal and the background in Fig. 11. The construction of PDFs for systematic errors is based on the three-flavor oscillation framework in Ref. [39]. The framework is capable of estimating the resulting change in a given event distribution after changing in a single systematic error. A two-dimensional histogram of NN output and reconstructed zenith angle is built which results from a 1σ change in each systematic error. An example of a twodimensional distribution for a 1σ change is shown in Fig. 13. The size of a systematic error determines the normalization of the distribution, which adjusts the size of the systematic error in the fit.
Since the uncertainties in oscillation parameters can also change the measured results and significance, a set of PDFs is built for each oscillation parameter. The values of sin 2 2θ 23 and Δm 2 32 are based on the Super-K atmospheric neutrino oscillation analysis result (Δm 2 32 ¼ 2.10 þ0.12 −0.18 × 10 −3 eV 2 , sin 2 θ 23 ¼ 0.5 AE 0.13) [17]. The value of sin 2 ð2θ 13 Þ based on the combined Daya Bay [8] and RENO [10] measurements of sin 2 2θ 13 ¼ 0.099 AE 0.014. For this analysis, the value of δ CP is set to be zero. Varying the value of δ CP from 0 to 2π in the three-flavor oscillation formula results in less than 1% change in the number of fitted tau events. The analysis is performed for both normal and inverted hierarchy. The data is fitted to the sum of background PDF, signal PDF and systematic PDFs varying the normalizations using a ROOFIT-based [38] unbinned likelihood fit algorithm as TABLE V. Systematic errors used in the tau neutrino appearance search that are dependent on Super-K run periods. The systematic errors are ordered by the maximum fractional change in the bins of the two-dimensional event distribution after shifting the systematic error by 1σ, and only systematic errors with the maximum fractional change larger than 2.5% are shown. The estimated 1σ error size is shown in percentage.
Multiring e-like background 12.1 11.1 11.4 11.6 Multiring PID Multi-GeV e-like −2.9 −3.9 2.7 −1.6 multi-GeV μ-like 6.5 9.7 −4.  IV. Systematic errors used in the tau neutrino appearance search that are common to all Super-K run periods. The systematic errors are ordered by the maximum fractional change in the bins of the two-dimensional event distribution after shifting the systematic error by 1σ, and only systematic errors with the maximum fractional change larger than 2.5% are shown. The estimated 1σ error size is shown in percentage.
where α is the normalization of the tau signal, with 0 meaning no-tau appearance and 1 meaning the expected tau appearance based on the neutrino oscillation parameters assumed in the simulation. The size of the ith systematic error in the fit, ϵ i , has a Gaussian univariate constraint in the fit. The PDFs of systematic errors are built separately for signal and background, but are combined together in the fit because the normalizations of both PDFs are adjusted by the same normalization factor ϵ i simultaneously. The fit is performed jointly with all data periods being fit at the same time. First, we present the results of the fit assuming the normal hierarchy of neutrino mass splitting. Relative to the expectation of unity, the tau normalization is found to be 1.47 AE 0.32 (stat þ syst) in the joint fit. This corresponds to a significance level of 4.6σ of rejection the hypothesis of no-tau-appearance. To estimate the statistical uncertainty of the fitted tau normalization, the systematic errors are excluded from the fit, and the tau normalization is found to be 1.41 AE 0.28. Therefore, the total uncertainty is dominated by the statistical uncertainty. The measured significance is larger than the expected significance of 3.3σ because more events are measured than expected. The number of tau events observed is evaluated after the fit by adding the tau events in the signal PDF rescaled by the fitted tau normalization and tau events in the systematic PDFs rescaled by the fitted values of systematic errors. The number of tau events is found to be 290.8 in the sample selected for this analysis. After correcting for efficiency, the observed number of fitted CC ν τ events over the entire running periods is calculated to be 338.1 AE 72.7 (stat þ syst), compared with an expectation of 224.5 AE 57.2 (syst) interactions.
The fit is repeated with the inverted hierarchy when calculating the oscillation probabilities, resulting in a reduction in the expected number of θ 13 induced upward-going electron neutrino. Under the assumption of inverted hierarchy, the fit results in a higher fitted value of tau normalization, 1.57 AE 0.31 and a correspondingly higher significance of 5.0σ. The higher fitted tau normalization is due to the reduction in θ 13 -induced upward-going electron neutrinos when calculating the oscillation probabilities under the assumption of the inverted hierarchy.
In order to test the stability of our analysis to changes in measured oscillation parameters, we also repeated our analysis procedure with the values of oscillation parameters in [37] and 2017 update [Δm 2 32 ¼ ð2.45 AE 0.05Þ × 10 −3 eV 2 , sin 2 ðθ 23 Þ ¼ 0.51 AE 0.04, and sin 2 θ 13 ¼ 0.0210 AE 0.0011]. Under the assumption of the normal hierarchy, the tau normalization is fitted to be 1.42 AE 0.31. The result of fitted tau normalization is basically unchanged with the updated oscillation parameter values.
The results of the final combined fit are examined graphically by plotting the binned projections of the fitted results. Figure 14 demonstrates the projection in zenith angle for both tau-like (NN output > 0.5) and non-tau-like (NN output < 0.5) events, along with the projections in NN output for both upward-going (cos Θ < −0.2) and downward-going (cos Θ > 0.2) events. In these plots, the signal PDFs have been rescaled to the fitted normalization values, and PDFs of systematic errors for signal and background have been rescaled by the fitted magnitudes of systematic errors and added to the signal and the background respectively. The fitted tau signal is shown in gray. All distributions have good agreement between data and MC simulations. In these plots, the PDFs and data from all of the run periods are combined.

B. Charged-current tau neutrino cross section measurement
This sample of CC ν τ interactions observed in Super-K offers the opportunity to measure the CC ν τ cross section. By scaling the theoretical cross section in the MC simulations to match the data, we can measure the inclusive charged-current tau neutrino cross section in water: where S τ is the factor that is used to scale the theoretical cross section to match simulations and data. For this analysis, S τ is the tau normalization measured in the search for tau neutrino appearance in Sec. VI A. Therefore, the measured CC ν τ cross section is expressed as: hσ theory i is the flux-averaged theoretical charged-current tau neutrino cross sections used in the NEUT code. To calculate the flux-averaged theoretical cross section, the differential CC ν τ cross section as a function of neutrino energy is weighted with the energy spectrum of atmospheric tau neutrinos from neutrino oscillations. Because CC ν τ interactions are not distinguishable from CCν τ interactions in Super-K, the theoretical cross section is a flux average of ν τ andν τ cross sections. The flux-averaged theoretical cross section, hσ theory i, is calculated as where dΦðE ν Þ dE ν is the differential flux of tau neutrinos as a function of neutrino energy as shown in Fig. 2, and σðE ν Þ is the differential charged-current tau neutrino cross sections used in NEUT code as seen in Fig. 3. The range of the integral is determined to be between 3.5 and 70 GeV from the tau neutrino energies in the simulation. As shown in Fig. 15, the neutrinos have energies more than 3.5 GeV in the CC ν τ interactions because of the energy threshold, and the expectation of CC ν τ interactions with more than 70 GeV is less than one in the entire run period.
The flux-averaged theoretical charged-current tau neutrino cross section is calculated to be 0.64 × 10 −38 cm 2 between 3.5 and 70 GeV, and thus the measured fluxaveraged charged current tau neutrino cross section: The measured cross section is shown together with the theoretical cross sections and the MC simulations in Fig. 15. The measured and theoretical cross section values are consistent at the 1.5σ level. C. Comparisons of charged-current tau neutrino cross section measurement with previously reported results Because of the difficulties in tau neutrino production and detection, charged-current tau neutrino cross sections have not been well measured. DONUT [14] and OPERA [15] are the only two experiments that have directly observed charged-current tau neutrino interactions, and DONUT is the only experiment that reported a measurement of the cross section. The DONUT measurement was based on nine observed charged-current tau neutrino events with an estimated background of 1.5 events. In DONUT, 800 GeV protons from the Fermilab Tevatron were used to produce neutrino beam by colliding with a beam dump, and tau neutrinos were produced via decays of charm mesons. The mean energy of the detected tau neutrino interactions was estimated to be 111 GeV, an energy at which deep inelastic interactions are dominant. Assuming that the DIS charged-current tau neutrino cross section had a linear dependence on neutrino energy, DONUT measured the energy-independent slope of the cross section, σ const , after correcting for the kinematic effect of tau lepton mass from the standard model calculation: where σðEÞ is the charged-current cross section per nucleon as a function of neutrino energy, σ const is the asymptotic slope which is constant in σ=E for deep inelastic scattering, and KðEÞ is the kinematic effect of tau lepton mass. DONUT measured σ const to be ð0.39 AE 0.13 AE 0.13Þ × 10 −38 cm 2 GeV −1 in their final results paper [40]. DONUT was incapable of distinguishing the charge of the τ lepton, therefore, the measurement is an average of the ν τ andν τ cross sections assuming equal number of ν τ andν τ in the neutrino flux. We wish to compare the ν τ cross section measured with atmospheric neutrinos by Super-K at relatively low energies to that measured by DONUT with a neutrino beam at higher energies. We recalculate the DONUT value of σðEÞ from Eq. (10) with the kinematic correction KðEÞ integrated over neutrino energies between 3.5 and 70 GeV and weighted to the world average ratio of cross sections between ν μ andν μ [37]. The calculated DONUT value of σðEÞ is then further weighted by the predicted ν τ andν τ flux ratio of 1.11 for atmospheric neutrino tau appearance at Super-K. The resulting σðEÞ is shown in Fig. 16. The charged-current tau neutrino DIS cross section inferred from the DONUT published number and reweighted to lower energy is ð0.37 AE 0.18Þ × 10 −38 cm 2 . This is smaller than our measurement of ð0.94 AE 0.20Þ × 10 −38 cm 2 , but the measurements are not yet directly comparable. DONUT measured the cross section with a neutrino beam that had a much higher average energy than that of the tau neutrinos in the atmospheric neutrino flux at Super-K. Quasielastic scattering and resonant pion production is a small component of the DONUT measurement and was neglected in their calculations. However, the tau neutrino flux at Super-K has a large component of neutrinos below 10 GeV, where charged-current quasi-elastic (CCQE) and resonant pion production makes a significant contribution to the detected event rate. We complete the comparison using the predicted CC DIS fraction in the Super-K sample. According to our Monte Carlo simulation, the fraction of DIS events in Super-K CC tau neutrino sample is estimated to be 41%. Therefore, the ν τ DIS-only cross section determined by Super-K atmospheric neutrinos is found to be ð0.40 AE 0.08Þ × 10 −38 cm 2 by scaling the measured cross section in Eq. (9) by 41%. This resulting DIS-only cross section is comparable and consistent with the DONUT measurement of the DIS ν τ cross section extrapolated to lower neutrino energy

VII. CONCLUSION
Using 5326 days of atmospheric neutrino data in SK-I through SK-IV, Super-K measured the tau normalization to be 1.47 AE 0.32, excluding the hypothesis of no-tau appearance with a significance of 4.6σ. A flux-averaged charged current tau neutrino cross section is measured to be ð0.94 AE 0.20Þ × 10 −38 cm 2 for neutrino energy between 3.5 and 70 GeV in Super-K, to be compared with the fluxaveraged theoretical cross section of 0.64 × 10 −38 cm 2 . Our result is consistent with the previous DONUT result, and is consistent with the Standard Model prediction to within 1.5σ.

ACKNOWLEDGMENTS
We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting Company. The Super-Kamiokande experiment has been built and operated from with σðEÞ inferred from DONUT (dashed lines). The DONUT cross section considers only DIS, and is digitized from [40].