Measuring tau neutrino appearance probability via unitarity

We propose a {\em unitarity method} for determining $\tau$ neutrino appearance probability $P(\nu_{\mu} \rightarrow \nu_{\tau})$ in long-baseline (LBL) accelerator experiments and atmospheric neutrino observations. When simultaneous in situ measurements of $P(\nu_{\mu} \rightarrow \nu_{\mu})$ and $P(\nu_{\mu} \rightarrow \nu_{e})$ proceed, as is typical in the LBL experiments, one can use unitarity to"measure"$P(\nu_{\mu} \rightarrow \nu_{\tau})$. A theorists' toy analysis for the model-independent determination of $P(\nu_{\mu} \rightarrow \nu_{\mu})$ and $P(\nu_{\mu} \rightarrow \nu_{e})$ is presented by using the NOvA data. It is shown in our analysis that $\lsim$5\% (8\%) measurement of $\tau$ neutrino appearance probability in neutrino (antineutrino) mode is possible in the peak region $1.5 \lesssim E_\nu \lesssim 2.5$ GeV. The $\nu$SM-independent nature of determination of the probabilities is emphasized.


I. INTRODUCTION
There exists a prevailing feeling in our community that the third generation is special among the fundamental fermions. It is exemplified, in particular, by the exceedingly large mass of the top quark [1,2]. But, even before the top quark was discovered [3,4] signaling its exceptionally large mass, people examined, for example, the possibility that the Higgs field conceals its origin which comes from much higher energy scale represented by the tt condensation [5][6][7]. In more contemporary contexts, if the Higgs sector is the most likely place as portal of new physics beyond the Standard Model (SM) of particle physics, the third-generation fermions could be the best source for such information due to their strongest couplings to the Higgs boson [1,8].
Among the third-generation fermions tau neutrino may be the least studied one. From now on, our discussion anticipates understanding that the observed neutrino masses are embedded into the SM, a theory which will be dubbed as the "νSM". So far, only a handful ν τ events had been identified. τ neutrino has first been seen experimentally in an event-by-event basis by the DONuT Group [9]. With the use of the ν µ beam from CERN, the ν τ appearance events have first been identified by the OPERA experiment [10]. Then both experiments looked for the "kink" events characteristic to a τ decay in nuclear emulsion. By using the statistically enriched ν τ samples of the atmospheric neutrinos, the charged-current (CC) ν τ cross section has been measured by Super-Kamiokande (Super-K) with 21% uncertainty [11], while IceCube's Deep Core measured CC + NC (neutral current) cross section with about 50% uncertainty [12].
It should be emphasized that we are now in a very good, timely position: Soon we will have intense τ neutrino beams at the far detectors in the next-generation accelerator long-baseline (LBL) experiments, Tokai-to-Hyper-Kamiokande (T2HK) [13] and the Deep Underground Neutrino Experiment (DUNE) [14]. Thanks to the large mixing angle θ 23 ∼ π/4, about half of the µ neutrino beam from J-PARC and LBNF, respectively, will be transformed into the τ neutrino beam at the far detectors, Hyper-K and DUNE. Because of availability of CC production of τ leptons due to its higher beam energy, DUNE must offer the best place for exploring τ neutrino physics. Naturally, this fact is receiving keen interests in the community, see e.g., Refs. [15][16][17] and the references cited therein. We should also mention that large samples of the atmospheric neutrinos taken in these big detectors will also do the same physics, with the likely chance of much improving the existing Super-K results mentioned above.
To facilitate the τ neutrino physics in the far detectors DUNE and Hyper-K in their full strength we must resolve one particular problem. As indicated in Eq. (3), for example, the energy distribution of leptons produced by CC reactions depends on the two unknowns, the ν µ → ν α oscillation probability and the ν α nucleus CC cross sections.
This statement is made under the assumption that the initial ν µ flux is known. To measure the ν τ cross sections we have to know the probability P (ν µ → ν τ ), and vice versa.
One may argue that at the present stage one can use P (ν µ → ν τ ) calculated by the νSM three-flavor mixing scheme to measure the ν τ CC cross sections. It is a sensible attitude given the current large errors in ν τ cross sections. But, when we start searching for new physics beyond the νSM in the ν τ sector, much better accuracies would be required. In this era, we must keep in mind the possibility that it would show up both in the ν τ appearance probability and the ν τ induced CC reactions. Looking for new physics effects in the tau lepton production under the assumption of no new physics in the appearance probability P (ν µ → ν τ ) (or, vice versa) may miss the key features of the phenomenon.
In this paper we propose the "unitarity method" for model-independent determination of the appearance probability P (ν µ → ν τ ). The idea is very simple, use unitarity assuming νSM-independent measurements of the probabilities P (ν µ → ν e ) and P (ν µ → ν µ ), see section II. For an existence proof, we present an analysis for a model-independent extraction of the probabilities P (ν µ → ν e ) and P (ν µ → ν µ ) from the data. While it is certainly at the level of "theorists' toy" analysis, we hope that it triggers the experimentalists' interests in measuring P (ν µ → ν τ ), and eventually leads them to the real analysis. We emphasize that determination of everything by experimental measurements in a modelindependent manner must be the ultimate goal of experimental physics.
In this section we describe the unitarity method for determining the appearance probability P (ν µ → ν τ ). In most of the LBL accelerator ν µ beam experiments, including T2K [18], NOvA [19], DUNE [14], and T2HK [13], the experimental data are, and will be taken primarily in both the ν µ → ν µ and ν µ → ν e channels simultaneously. The same statement applies to the atmospheric neutrino observation even though the event characterization, e.g., identification of initial and final neutrino flavors, would be much more involved in some cases. Then, by using unitarity one can "measure" the ν τ appearance probability P (ν µ → ν τ ).
It is quite possible that many people thought about this or the similar ideas related to this. In that case, the present paper may add little new. But, to the best of our knowledge, the unitarity method has never been presented explicitly in an organized way, and it prompted us to write this paper. As we will learn in our discussions below, there exist many things to be understood in this method. We hope our considerations in this paper urge experimentalists to think about the unitarity method for measuring P (ν µ → ν τ ), because the real analysis can only be done by people inside the experimental collaborations.
A. Does assuming unitarity imply assuming the νSM?
One might ask: Is it not true that assuming unitarity is essentially equivalent to usage of the νSM expression of the probability P (ν µ → ν τ )? The answer is No, not at all. That is, unitarity is much more robust and it should generally hold. Toward having a clear cut discussion, we must first understand unitarity on generic ground: • If only the three active neutrinos span the complete state space of neutral leptons, unitarity must hold in a model-independent manner. There is no way to go outside the complete state space during propagation, assuming absence of inelastic scattering, absorption, etc., and hence neutrino evolution must be unitary. 1 Therefore, unitarity holds even in the case that neutrinos have additional interactions such as the non-standard interactions (NSI) [21][22][23][24] in propagation, which are under active search by various experimental methods which produced the numerous constraints [25][26][27].
B. How could non-unitarity come in?
Then, the question might be: In what circumstances can the three-neutrino unitarity be violated? The simplest answer would be existence of the fourth, sterile neutrino, which may be indicated by LSND and MiniBooNE [28,29].
For an overview see e.g., Refs. [30][31][32]. If this and the similar two or three sterile scenarios are confirmed experimentally by the various experimental methods, e.g., described in Refs [33,34], our unitarity method for P (ν µ → ν τ ) has to be revised.
Yet, all is not lost. Typically, there are two possibilities: (1) the unitarity method is still valid under the certain conditions, and (2) the unitarity method can be amended in such a way that it is valid as in the no sterile case. To discuss the first case (1), let us define the non-unitarity parameter ξ If the error of obtained P (ν µ → ν τ ) is larger than a few times ξ, we can ignore the issue of non-unitarity by the sterile neutrino for the moment, because the probability leaking to the sterile sector is smaller than the reachable accuracy for P (ν µ → ν τ ).
In the case (2) we assume that the sterile neutrino masses and the mixing parameters can be measured such that a modified unitarity relation P (ν µ → ν e ) + P (ν µ → ν µ ) + P (ν µ → ν τ ) + i P (ν µ → ν Si ) = 1 can be set up, where the sterile label i runs over a few sterile neutrinos. To the extent we know P (ν µ → ν Si ) well, the modified unitarity method should work better. The method works if the error in of P (ν µ → ν τ ) is comparable or larger than the estimated errors of i P (ν µ → ν Si ).

C. More generic scenarios for non-unitarity
More generically, if there exists an extra sector which is somehow isolated from the νSM one, but has a contact by having a weak coupling or small mixing with the three neutrino species, our three neutrino system is approximately unitary, but not perfectly. The most well known example for the extra sector is the heavy right-handed (RH) neutrinos in the seesaw model of neutrino masses [35][36][37][38][39]. In the three active plus three RH neutrino model, the 6 × 6 flavor mixing matrix is unitary, but 3 × 3 sub-matrix for the active neutrinos is not. But in the original scenario, since the RH neutrinos are so heavy, m RH ∼ 10 15 GeV, the non-unitarity of the flavor mixing matrix for the three active neutrinos is practically undetectable.
In fact, much less model-dependent descriptions for more generic unitarity violation (UV) scenarios exist for beyond νSM new physics both at high scale m W [40] and low scale m W [41]. For the terminology of high-and low-scale UV see Ref. [41]. In high-scale UV, the UV effect has to be small because the prevailing weak SU (2) × U (1) structure allows us to use the charged lepton probe to constrain UV effect in the neutrino sector [40]. It indeed entails the severe bounds on UV [42][43][44][45]. If the UV effects are parametrized by the α parameters [43], they are constrained to be < ∼ 10 −3 , or smaller [45], and we should obtain the ξ parameter bound of the similar order. Thus, our unitary method works in the presence of high-scale UV.
The low-scale UV scenarios may be described by using the system of three active plus generic N sterile neutrinos, the model known since the early days, see e.g. Refs. [46,47]. Within this framework, a sterile sector model-independent description of low-scale UV is attempted [41,48]. In such scenarios the bounds on the α parameters are milder by at least one order of magnitude, and even more milder for α τ τ . See e.g., Refs. [45,49,50]. A rough estimation in Appendix A reveals the current upper bound on ξ of about 0.1 or less. Clearly, we need the better bound on ξ to ensure the validity of our unitarity method.

III. DETERMINATION OF OSCILLATION PROBABILITY WITHOUT νSM ANSATZ
In the rest of this paper, we proceed with assumption of no unitarity violation in the three active neutrino space until Appendix A. To put the unitarity method for measuring P (ν µ → ν τ ) into practice we need to determine the neutrino oscillation probabilities P (ν µ → ν µ ) and P (ν µ → ν e ). As we have learnt in section II our unitarity method does not necessitate the three-flavor νSM ansatz, we want to carry this task out in the theoretical-model independent way, as much as possible.
To give this general idea a concrete shape, we present a toy analysis in this section assuming the experimental setting of the LBL accelerator neutrino experiment with muon neutrino beam. Since analyses of the atmospheric neutrino data are quite involved we focus on accelerator LBL experiment from now on. Among the two ongoing LBL experiments, T2K [18] and NOvA [19], we focus on the latter because of its higher energy beam. While exploration of τ neutrino physics using the CC τ production may require a higher energy neutrino beam, we should wait for LBNF nominal, or preferably its τ -optimized configurations [16] for this purpose. We will be merited by the fact that NOvA has the functionally identical near and far detectors: A large fraction of the systematic errors would cancel between the two detectors.
To show the basic idea of our toy analysis, we assume the quasi-elastic CC reactions ν e N → e − N and ν µ N → µ − N for detection of ν e and ν µ at both the near and far detectors. The choice, where N and N denote, respectively, the target and produced nuclei, enables us to reconstruct the initial neutrino energy E ν via the two-body kinematics. Nonetheless, by using the data in which the events with four hadronic energy-fraction quartiles [51,52] are added, purity of the quasi-elastic nature of the CC events sample may be slightly harmed. Yet, we hope that the major part of this and the related problems are taken care of by the resultant relatively large error bars possessed by the results of P (ν µ → ν µ ) and P (ν µ → ν e ) obtained by our method.
To describe the principle of our analysis, we hereafter discuss explicitly only the neutrino channels, but the way how the antineutrino channels can be handled should be obvious from the neutrino channel discussion. After a brief description of event number distribution via the quasi-elastic CC reactions in section III A, we carry out our analyses for P (ν µ → ν µ ) and P (ν µ → ν e ) in sections III B and III C, respectively. Then, we obtain P (ν µ → ν τ ) by our unitarity method in section IV. The similar analyses for the antineutrino channel probabilities will be repeated in section V.

A. Theoretical expression of the event number distribution
Muon neutrinos ν µ of energy E ν in the neutrino beam reach a detector at distance L from the production point as ν µ or ν e with the probabilities P (ν µ → ν µ : E ν , L) and P (ν µ → ν e : E ν , L), respectively. The event number distribution at the detector by the CC reaction ν µ N → − α N , where α (α = e, µ, τ ) are SU (2) L doublet, can be written as a function of neutrino energy E ν as where N T denotes the number of target particles, Φ νµ (E ν , L) is the ν µ flux at the distance L from the source, and dσ dE α is the cross section of the CC reaction ν µ N → − α N . which produces α lepton of energy E α . (E α ) denotes energy-dependent efficiency of identifying α lepton. In the last line of Eq. (3) S α (E ν , L) and B α (E ν , L) denote, respectively, the contributions from signal events without oscillation, and from background events.
As it stands, the expression in Eq. (3) does not fully respect the experimental reality. The energy of neutrinos that undergo the CC reactions must be reconstructed using the reaction kinematics, and E ν in Eq. (3) must be understood as the reconstructed energy. In this process the various issues, e.g., the detector energy resolution and the effect of Fermi motion (as the target nucleus is in nuclei) have to be taken into account. Equation (3) assumes that the error associated with this reconstruction process is small compared to the genuine neutrino energy. The assumption seems to be supported by the result of simulation which reports less than 10% error in the reconstructed energy [52]. See section VI A for a brief description of how Eq. (3) may be justified. A final comment on Eq. (3) is that the sum over the various CC reactions, a a (E µ ) dσa dEµ where a denotes indices for the varying reaction channels, must be introduced with varying efficiencies. As it can be done without affecting the validity of our following discussion, we keep our simple expression Eq. (3) as it is, with understanding that the summing over the CC reactions is always meant.
Despite these and possibly other drawbacks, we use the expression in Eq. (3) as the toy model for the event number distribution as a function of reconstructed neutrino energy. Despite that we do not write down the explicit expressions of the similar formulas in the antineutrino channels, they are easily obtained in an analogous way as in the neutrino channels.
A few comments on the NOvA data used in this paper: From start to almost the end of our analysis we have consulted and used the information given in Ref. [19], which is then updated in Ref. [53]. Very recently a new paper appeared from the NOvA collaboration [54] which reports all the available data to date in the neutrino and antineutrino channels. It appears that the Monte Carlo analysis code is completely renewed. In each period, the data are conveniently made available at the NOvA data release [55], and we utilize the most recent version of it in our analysis.
B. Determination of disappearance probability P (νµ → νµ) Now, we describe a method for extracting the survival (or disappearance) probability P (ν µ → ν µ ). Quite conveniently for our purpose, the experimental groups not only provide the experimental data of dNe dEν (ν µ N → µ − N : L far ), the left-hand-side (LHS) of Eq. (3) (α = µ), but also Monte Carlo expectation of the same quantity without oscillation. If we take the ratio of these two quantities at the far detector distance, we obtain where we have defined the background to signal ratio The right-hand side (RHS) of Eq. (4) is almost the probability, apart from the r µ terms, because all the factors other than these cancel out between the numerator and the denominator. This cancellation takes place even in the case that sum over the varying reaction channels are introduced in the CC reactions to produce muons, as mentioned earlier.
It may be relevant for higher hadronic energy-fraction quartiles [51,52].
Thanks to the experimental group the information of the background is also provided [55], and hence we can obtain the disappearance probability P (ν µ → ν µ ). The background for ν µ disappearance CC events is about 4% level for neutrino and 3% level for antineutrino channels, respectively [19]. Notice that in plotting the event number distribution as a function of reconstructed neutrino energy, the effects of energy smearing through the event reconstruction process as well as by the Fermi motion are taken care of by the experimental group. The same comment applies to the plot for extracting P (ν µ → ν e ) in section III C.
Therefore, the determination of P (ν µ → ν µ ) through Eq. (4) would be the cleanest way among the methods for determining the oscillation probability we discuss in this paper. Notice that our method is a data-driven way, and we do not rely on the expression of the oscillation probability calculated by the νSM standard three-flavor oscillation. The obtained result for P (ν µ → ν µ ) is presented in Fig. 1 with the black histogram and its 1σ error band as the shaded gray region.  Table IV in Ref. [19]. In the legend "NOvA BF" stands for the "NOvA best fit".
The blue line in Fig. 1 is the νSM three-neutrino expression of P (ν µ → ν µ ) with the mixing parameters used by the NOvA group, in Table IV in Ref. [19]. For simplicity and brevity we call this parameter set as the "NOvA best fit". Figure 1 indicates that the obtained result of P (ν µ → ν µ ) is consistent with the standard three neutrino oscillation. In fact, the νSM line passes through the 1σ uncertainty band of the obtained histogram in almost all the bins. As mentioned above, and will be further discussed in section III D, our method for determining P (ν µ → ν µ ) does not rely on the νSM. Therefore, we do not judge whether our method is successful or not by how close our result is to the νSM.
In certain limited energy regions in Fig. 1, the 1σ error band of P (ν µ → ν µ ) penetrates into the unphysical regions of P (ν µ → ν µ ) < 0. Similarly, later in section V we will see that P (ν µ →ν µ ) expands into the region > 1. See Fig. 4. We expect that these features will disappear as the better statistics of events is accumulated.
C. Determination of appearance probability P (νµ → νe) Now, we discuss determination of P (ν µ → ν e : E ν , L far ). In our simple-minded experimental setting of the LBL neutrino experiment we assume that the near detector is placed at a location so close to the neutrino production point such that one can safely assume that P (ν µ → ν µ : L near ) = 1. This is a good approximation for the NOvA experiment because L near /L far ≈ 10 −3 . Then, the relevant ratio of the event number distributions is given, by using Eq. (3), as 2 where r e and r µ are defined in Eq. (5), and f eµ is defined by The ratio f eµ is the far-to-near flux ratio weighted by (i) the detector volumes and (ii) the efficiencies averaged over the event number distributions. Since the LHS of Eq. (6), the both numerator and denominator, is given by the experimental group, we can determine the appearance probability P (ν µ → ν e : E ν , L far ) if we know f eµ , r e (L far ), and r µ (L near ). Plotted by the thick black histogram is the appearance probability P (ν µ → ν e : E ν , L far ) calculated by using (9), and the shaded gray region is its 1σ error. The blue line shows the νSM oscillation probability calculated with the mixing parameters given in Table IV in Ref. [19]. In the legend "NOvA BF" stands for the "NOvA best fit".
Despite that the experimental group keeps the information on f eµ not public, the result of Monte Carlo calculation is given for the event number distribution of electrons at the far detector in Slide 23 of Ref. [53] (see Fig.4 in Ref. [54]): where dNµ dEν (ν µ N → µ − N : L near ) is given by the NOvA experimental measurement and P (ν µ → ν e : E ν , L far )| M C is calculated by using the "NOvA best fit". Quite conveniently, the both quantities are included in the NOvA data release [55]. Then, one can solve Eq. (8) for f eµ . Now, we must note that f eµ defined in (7) contains only the information on the signal events, not background. Therefore, to evaluate f eµ by using Eq. (8) we must restrict both LHS and RHS of Eq. (8) the information of the signal events only. One can easily satisfy this condition for the quantities obtained by MC, but not the muon number distribution, the last factor in Eq. (8) because it is the data. But, this problem is easily avoided if one insert everything into Eq. (6), which entails The obtained result of P (ν µ → ν e : E ν , L far ) is presented in Fig. 2. Again the νSM blue line is consistent with P (ν µ → ν e ) in Fig. 2.
As in the case of P (ν µ → ν µ ) in Fig. 1, P (ν µ → ν e : E ν , L far ) goes into unphysical regions. Even the think black line (3.0 ≤ E ν ≤ 3.5 GeV) as well as the lower end of gray-shaded region (2.5 ≤ E ν ≤ 4.0 GeV) go down into minus. We expect that these features will disappear as the measurements further proceed.

D. ν Standard Model independence of our method and its significance in wider contexts
Now some of the readers may argue that by using the NOvA Monte Carlo simulation results in Eq. (9) our analysis depend on the standard three-flavor model of oscillation. If so, we can no longer claim that it is independent of the νSM paradigm. Fortunately, this is not the case. Notice that dNe dEν (ν e N → e − N : L far )| M C scales as P (ν µ → ν e : E ν , L far )| M C apart from the small background contributions. Then, the dependence on P (ν µ → ν e : E ν , L far )| M C cancels between the numerator and the denominator in Eq. (9), allowing us to remain essentially in the νSM independent analysis.
Extraction of the oscillation probabilities P (ν µ → ν e ) and P (ν µ → ν µ ) in a model-independent way may be important in much wider contexts beyond the unitarity method for P (ν µ → ν τ ). Currently the experimental results are reported by showing the best-fit values of the mixing angles and the CP phase by assuming the νSM parametrization. While it is a valid way, the result is of course νSM dependent. Instead, a model-independent extraction of the oscillation probability itself could directly signal effects outside the νSM. It can be done immediately with the currently available data, but it would become an indispensable alternative in high-statistics experiments like T2HK and DUNE.

IV. DETERMINATION OF P (νµ → ντ )
Given our estimates of P (ν µ → ν µ ) and P (ν µ → ν e ) in Figs. 1 and 2 in sections III B and III C, respectively, it is now straightforward to obtain P (ν µ → ν τ ) by using unitarity (1). The result is given in Fig. 3. The errors of P (ν µ → ν µ ) and P (ν µ → ν e ) are added in quadrature. As before, the blue line shows the νSM oscillation probability calculated with the mixing parameters given in Table IV in Ref. [19], the "NOvA best fit".
We need to make some comments on Fig. 3, because we have presented the two panels. They differ in the binning, mainly at low energies E < ∼ 3 GeV. The issue is that while P (ν µ → ν µ ) is determined with finer bins as seen in Fig. 1, P (ν µ → ν e ) has coarse bins as in Fig. 2. If we use the coarse bins for the both P (ν µ → ν µ ) and P (ν µ → ν e ), the result in the upper panel is obtained. But, since P (ν µ → ν µ ) is much larger than P (ν µ → ν e ) in most bins, we could combine P (ν µ → ν µ ) and P (ν µ → ν e ) in such a way that the respective binning of P (ν µ → ν µ ) and P (ν µ → ν e ) are kept as they are. If we take this attitude the obtained result of P (ν µ → ν τ ) is presented in the lower panel.
A problem in our treatment for the lower panel P (ν µ → ν τ ) is, of course, we have to assume that P (ν µ → ν e ) is constant over the energy regions of e.g., E = 1.0 − 1.5 GeV, or E = 1.5 − 2 GeV, whereas P (ν µ → ν µ ) changes in the region. Nonetheless, P (ν µ → ν µ ) significantly varies in the region E = 1.0 − 2.0 GeV, so that keeping the information with finer bins would make sense. These are the reasonings for which we wind up to present the two panels in Fig. 3. The blue line for the νSM three-neutrino expression of P (ν µ → ν τ ) reasonably fit to our results both in the upper and lower panels. The τ neutrino appearance probability P (ν µ → ν τ : E ν , L far ) is plotted by the thick black histogram, which is calculated by using unitarity with the ν µ → ν µ and ν µ → ν e probabilities in Figs. 1 and 2, respectively. The gray shaded area is its 1σ error. In the upper panel, to combine P (ν µ → ν µ ) and P (ν µ → ν e ), we take the coarse bins for the both channels. In the lower panel, we have kept the original bin sizes of the both P (ν µ → ν µ ) and P (ν µ → ν e ) as in Figs. 1 and 2, respectively. See the text for more details. The blue line shows the νSM oscillation probability calculated with the "NOvA best fit".
The < ∼ 5% measurement of P (ν µ → ν τ ) around the peak region 1.5 < E < 2.5 GeV reported in Fig. 3 is certainly intriguing. But, we postpone our comment to section V A where we make comparison between the results of P (ν µ → ν τ ) in the neutrino and antineutrino (ν →ν) channels.

V. DETERMINATION OF PROBABILITIES IN THE ANTINEUTRINO CHANNELS
In this section, we repeat the same exercise for the antineutrino channels, the ones we have carried out in sections III and IV for the neutrino channels. The antineutrino channels are important to obtain the information on CP violation in combination with the neutrino channel. The blue line shows the νSM oscillation probability calculated with the "NOvA best fit".
In Fig. 4, plotted is the disappearance probability P (ν µ →ν µ : E ν , L far ) (upper panel), and the appearance probability P (ν µ →ν e : E ν , L far ) (lower panel), which are calculated by using the antineutrino versions of Eqs. (4) and (9), respectively, and their 1σ errors. The νSM oscillation probability calculated with the "NOvA best-fit" is also shown. Roughly speaking, the uncertainties in determination of P (ν µ →ν µ : E ν , L far ) and P (ν µ →ν e : E ν , L far ) are comparable to each other. However, in the disappearance channels bin-to-bin fluctuations look somewhat larger The antineutrino appearance probability P (ν µ →ν τ : E ν , L far ) and its 1σ error region are presented with the same style as in Fig. 3 for the neutrino version. The blue line shows the νSM oscillation probability calculated with the "NOvA best fit".
in the antineutrino channel with a few vanishing number of event bins at around the oscillation maximum, i.e., the oscillation minimum in the disappearance channels P (ν µ →ν µ ).
Probably due to lack of statistics the probability exceeds unity in a few low and high energy bins of P (ν µ →ν µ : E ν , L far ). Similarly, P (ν µ →ν e : E ν , L far ) goes into minus at the similar low and high energy bins.
A. Result of P (νµ →ντ ) and its comparison with P (νµ → ντ ) In Fig. 5, theν τ appearance probability P (ν µ →ν τ : E ν , L far ) and its 1σ error are plotted. The style of presentation and line symbols are the same as before. They are calculated by using the antineutrino version of unitarity (1). The upper panel is for the case of common coarse bin as used in P (ν µ →ν e : E ν , L far ), while the lower panel is for use of different binning, the finer bin for P (ν µ →ν µ : E ν , L far ), and the coarse bin for P (ν µ →ν e : E ν , L far ), as done in the lower panel of Fig. 3 in the neutrino channel.
By comparing between the obtained ν τ andν τ appearance probabilities in Fig. 3 and Fig. 5, one can say that (1) The uncertainties of the appearance probabilities are comparable but slightly larger in the antineutrino channel. (2) More visibly, the bin by bin fluctuations are larger in the antineutrino channel. The accuracy of P (ν µ → ν τ ) itself is quite good with less than 5% (8%) error in the peak region 1.5 < E < 2.5 GeV in the neutrino (antineutrino) channel. One may say that the accuracy of 5%, or 8%, is a superb performance, but it is basically achieved by the experimental measurement of P (ν µ → ν µ ) and P (ν µ → ν e ), and what is done by our analysis is to translate the accuracies to P (ν µ → ν τ ). Whereas the central value as well as the error of the antineutrino probability considerably fluctuate bin by bin, but at less than ±10% level for the central value in the same peak region as above. In the above we are referring the finer bin versions of P (ν µ → ν µ ) and P (ν µ →ν µ ).
Here, we note a possible mechanism of error reduction for P (ν µ → ν τ ), which is characteristic to our unitarity method. First of all the effect of P (ν µ → ν e ) is relatively minor, and hence we disregard it in this discussion. In the disappearance channels the statistics is high, and we could assume that the errors are well characterized as a relative, percent error. The peak region of the appearance channel ν µ → ν τ corresponds to the region where P (ν µ → ν µ ) is small, so that the error of P (ν µ → ν µ ) is also small. The small error, through unitarity, leads to the small error of P (ν µ → ν τ ) in its peak region. Since P (ν µ → ν τ ) is large in the peak region, its percent error is even smaller. If this is the qualitatively correct explanation, it is a new merit of the unitarity method. When much higher statistics are gained, a smaller percent error of P (ν µ → ν τ ) than P (ν µ → ν µ )'s would manifest in regions where the both probabilities are large.
The smallness of the error might also be because the experimental errors are not taken into account to a sufficient level or the error correlations is important. On the other hand, large fluctuations in P (ν µ →ν τ ) seem to tell us that accumulating a better event statistics is necessary, which is not easy to achieve in the LBL neutrino experiments. Even though we have included the T2K data to our analysis, it would not improve so much the accuracy of our unitarity-reconstructed P (ν µ → ν τ ) because the T2K events mostly span lower energy region than NOvA's, E ν < ∼ 1 GeV. In this sense these two LBL experiments are complementary with each other by covering the different energy regions. 3

VI. A FEW FINAL REMARKS
We have described our unitarity method for determining ν τ appearance probability P (ν µ → ν τ ), and examined performance of the method by taking the concrete case of the NOvA experiment. We believe that our analysis method is reasonably set up to allow model-independent determination of the probabilities P (ν µ → ν α ) (α = e, µ, τ ), and the results are indeed sensible. But, there are limitations inherent to our method.

A. Assumptions and limitations of our analysis
The most important approximation we have made in deriving our basic equation (3) is that the error in reconstructing the neutrino energy is much smaller than the genuine neutrino energy. Without this assumption we cannot factorize the oscillation probabilities as in Eq. (3). The point may be illuminated by a toy-model expression of the event number distribution as a function of the reconstructed neutrino energy E rec in the reaction ν µ is oscillated to ν α , and ν α undergoes CC reaction ν α + N → α + N , where we have assumed the Gaussian shape of E ν reconstruction error function. Under the limit σ E ν Eq. (10) reproduces Eq. (3). Fortunately, the detailed study in Ref. [52] assures the smallness of the error in neutrino energy reconstruction to be less than 10%, which is indeed small but not vanishingly small.
Further "limitation discussions" on our analysis would entail an endless list. For example, mistreatment of error correlations, or double counting of the errors, etc. Or, one could raise the possibility of analysis without binning. We are reluctant to enter into the detailed discussions of these or the other points here. It is because, we believe, improving our toy analysis is not the right way to proceed. What is really needed is the real analysis by the experimental group.
B. Improving the bound on non-unitarity Improving the constraints on non-unitarity, in our case on the ξ parameter, is important to strengthen the basis of our unitarity method for P (ν µ → ν τ ). In more generic context including the ξ bound, we expect that the better constraints on UV which improve the current ones [45,49,50] will be obtained before DUNE starts to do τ neutrino physics. It will be done, for example, by the ongoing and upcoming experiments such as SBN program at Fermilab [33], JSNS2 [34], T2K [18], NOvA [19], Super-K [56], IceCube [57,58], KM3NeT [59], JUNO [41,60], and possibly Hyper-K [13]. These are the case of low-scale UV (or low mass sterile leptons) and the bound is already much severer in high-scale UV case, < ∼ 10 −3 [45].

C. Absolute neutrino flux
If our purpose is restricted to determine the oscillation probability only, the necessity of knowing the precise muon neutrino flux may be relaxed because the near-far detector comparison basically does the job. The fact that NOvA has the functionally identical near and far detectors certainly helps. However, to measure ν τ CC cross sections with comparable accuracy with ν µ 's, and to study possible new physics effects in the ν τ induced reactions, we would need to know the absolute neutrino flux, hopefully to the accuracy better than what are achieved for the ongoing projects [61,62]. A method of using ν µ − e scattering is suggested [63] based on the measurement in MINERvA [64].

VII. CONCLUDING REMARKS
In this paper we have described a way of determining τ neutrino appearance probability P (ν µ → ν τ ) using unitarity in the νSM independent way. Despite our analysis is at the level of theorists' toy exercise, we hope, we were able to demonstrate the "in principle feasibility" of the unitarity method for measuring P (ν µ → ν τ ). Of course, the experimentalist-level real analysis must be performed to give the idea a realistic shape. If this paper acts as a trigger for this, it would be the most successful outcome of this paper.
Once DUNE and Hyper-K turn on in the near future, we will enjoy the rich prospects for τ neutrino physics. They will carry out simultaneous measurement of P (ν µ → ν µ ) and P (ν µ → ν e ), and ν τ appearance probability P (ν µ → ν τ ) can be determined by the unitarity method. Then, the promising expectation is that the oscillation-produced intense ν τ neutrino beam in DUNE can be used to investigate the properties of τ neutrino CC reactions. It is worth to note that all these processes take place in the DUNE experiment in an in situ manner. This feature would allow the reduction of systematic uncertainties by comparing between near and far detectors, and in mutual simultaneous analyses of the three observables. From the viewpoint of unitarity measurement of P (ν µ → ν τ ) from low to high energies (say, 400 MeV to ∼10 GeV), T2HK and DUNE will play complementary role as analogous to T2K-NOvA complementarity, but at much higher level of the accuracies.
The method for measuring P (ν µ → ν τ ) via unitarity may be applicable to the atmospheric neutrino observation, because extraction of P (ν µ → ν µ ) and P (ν µ → ν e ) from the data should be possible in the analyses. 4 If it works in Super-K it will allow DUNE to enjoy the knowledge of P (ν µ → ν τ ) in the energy region of 1 − 10 GeV from the first day of its operation. Notice that Super-K will be able to accumulate the atmospheric neutrino data for 30 years at the DUNE turn-on, which provide a rare, valuable chance of the international collaboration for tau neutrino physics.