Physics with Beam Tau-Neutrino Appearance at DUNE

We explore the capabilities of the upcoming Deep Underground Neutrino Experiment (DUNE) to measure $\nu_\tau$ charged-current interactions and the associated oscillation probability $P(\nu_\mu \to \nu_\tau)$ at its far detector, concentrating on how such results can be used to probe neutrino properties and interactions. DUNE has the potential to identify significantly more $\nu_\tau$ events than all existing experiments and can use this data sample to nontrivially test the three-massive-neutrinos paradigm by providing complementary measurements to those from the $\nu_e$ appearance and $\nu_\mu$ disappearance channels. We further discuss the sensitivity of the $\nu_\tau$ appearance channel to several hypotheses for the physics that may lurk beyond the three-massive-neutrinos paradigm: a non-unitary lepton mixing matrix, the $3+1$ light neutrinos hypothesis, and the existence of non-standard neutral-current neutrino interactions. Throughout, we also consider the relative benefits of the proposed high-energy tune of the Long-Baseline Neutrino Facility (LBNF) beam-line.


I. INTRODUCTION
Over the last twenty years, our understanding of neutrino-flavor change as a function of the neutrino proper time has improved exponentially. The old solar and atmospheric neutrino anomalies evolved into the very robust three-massive-neutrinos paradigm, which postulates that neutrinos, while still interacting with ordinary matter as prescribed by the standard model of particle physics (SM), have distinct, non-zero masses and that the neutrino flavor-eigenstates ν α , α = e, µ, τ , are non-trivial linear superpositions of the neutrino mass-eigenstates ν i , with masses m i , i = 1, 2, 3: ν α = U αi ν i , where U αi are the coefficients of the 3 × 3 unitary lepton mixing matrix U .
Assuming the three-massive-neutrinos paradigm is correct, the oscillation parameters -those that determine U , along with the neutrino mass-squared differences ∆m 2 ij ≡ m 2 i − m 2 j , i, j = 1, 2, 3 -are well constrained by existing neutrino data, with a few exceptions, including the neutrino mass ordering (or the sign of ∆m 2 31 ) and the strength of leptonic CP-invariance violation. When it comes to both confirming the three-massive-neutrinos paradigm and determining the oscillation parameters, all current statistical power comes from studies of ν µ and ν e (plus antineutrinos) disappearance and ν µ → ν e (plus antineutrinos) appearance. While there is definitive evidence for ν µ → ν τ appearance from both atmospheric [1][2][3] and beam [4] neutrino experiments, the quantitative impact of the current ν τ appearance data is, at best, sub-dominant. * Among the goals of the next-generation long-baseline neutrino experiments -the Long-Baseline Neutrino facility to the Deep Underground Neutrino Experiment (LBNF-DUNE) in the United States and the Tokai to Hyper-Kamiokande (T2HK) experiment in Japan -are the precision measurements of neutrino oscillation parameters, exploring CP-invariance violation in the neutrino sector, and, perhaps most important, testing the validity of the three-massive-neutrinos paradigm and looking for more new physics in the neutrino sector. There are several studies in the literature, including those pursued by the collaborations, of the physics reach of LBNF-DUNE and T2HK, including their sensitivities to a variety of hypothetical new neutrino physics scenarios. All of these explore the ν µ disappearance and the ν µ → ν e appearance (plus antineutrinos) channels.
Here, instead, we explore the physics reach of the beam ν τ appearance data that is accessible to LBNF-DUNE. Direct measurements of ν τ appearance, both in atmospheric and beam neutrino experiments, are very challenging, for a variety of reasons. All neutrino experiments are of the fixed-target type and the large τ mass translates into relatively large thresholds for charged-current ν τ scattering off of ordinary matter (E ντ 3.35 GeV for ν τ + N → τ + N , where N is a nucleon, and E ντ 3.1 TeV for ν τ + e → τ + ν e ). Given the beam energies of LBNF-DUNE and T2HK, driven by the requirement that large ∆m 2 31 -driven oscillation effects are observed in the far detector, only LBNF-DUNE is expected to observe beam events with neutrino energies above the τ -threshold. Even at LBNF-DUNE, phase-space-suppression effects are large and the ν τ appearance event sample is expected to contain, for the lifetime of the experiment, between 100 and 1,000 events. Identifying and reconstructing τ -leptons in neutrino detectors is also very challenging. The tracking resolution of liquid argon detectors, around several millimeters, is such that τ -leptons decay promptly and hence must be identified via their decay products. Furthermore, all τ leptons decay channels involve missing energy -τ → ν τ + something else. These imply that the ν τ appearance channel is "dirty" and that the reconstruction of the ν τ energy is a larger challenge than that of the ν µ and ν e energies. We discuss our simulations of the ν τ sample at LBNF-DUNE in Sec. II, including background estimates and the challenges of neutrino-energy reconstruction.
In spite of all the challenges, LBNF-DUNE is expected to collect an unprecedented number of reconstructed beam ν τ events. Here, we proceed to understand what nontrivial particle physics information one should be able to extract. We explore, assuming the three-massive-neutrinos paradigm, how well one can measure the ν τ charged-current scattering cross section, and, assuming the standard model expectation for the ν τ charged-current scattering cross section, how well the ν τ appearance data can constrain oscillation parameters. As far as the latter exercise is concerned, we compare, in different ways, the sensitivity of the ν τ appearance sample with that of the ν e appearance and ν µ disappearance samples. These results are presented and discussed in detail in Sec. III A.
While, not surprisingly, the reach of the ν τ appearance sample is comparatively weak, the value of using different oscillation channels to measure the oscillation parameters cannot be overstated. Indeed, comparing the results of different oscillation channels consists of one of the most robust tests of the threemassive-neutrinos paradigm. In order to further pursue how the ν τ appearance sample complements the search of physics beyond the three-massive-neutrinos paradigm, we explore different concrete scenarios, including the existence of new heavy and light neutrino degrees of freedom -Sec. III B and Sec. III C, respectively -and the existence of new neutral-current neutrino-matter interactions (Sec. III D).
In Sec. IV, we present some concluding remarks and highlight a few future research directions one can pursue with ν τ -appearance at LBNF-DUNE, including searches in the near detector facility and the atmospheric ν τ sample.

II. THE ν τ SAMPLE AT DUNE
The production of τ leptons by charged-current ν τ -nucleus scattering requires neutrino energies E ν 3.4 GeV. Furthermore, the prompt decay of the τ combined with the fact that all τ decays contain ν τ in the final state prove a challenge for both identifying a scattering event as a ν τ charged-current interaction and reconstructing the initial neutrino energy. Nonetheless, the LBNF beam will have a significant portion of its flux above the τ production energy threshold and liquid argon detectors are excellent at reconstructing final-state particles, allowing for at least a modest but unique ν τ data sample. Here, we discuss our expectations regarding the capability of DUNE to identify ν τ interactions and measure their energies.
We assume that the DUNE far detector, 1300 km from the neutrino source at Fermilab, consists of 40 kton of liquid argon (fiducial mass) and that the LBNF beam will deliver 1.1 × 10 21 protons on target (POT) per year. We consider three different modes of beam operation -forward horn current (which we will refer to as "neutrino mode"), reverse horn current ("antineutrino mode"), and the currently-underconsideration tau-optimized flux [7] ("high-energy mode"). The high-energy mode contains significantly higher-energy neutrinos than the other two modes, with a large fraction of the neutrinos above the τ threshold. For neutrino and antineutrino modes, we simulate event yields using the "CP-optimized fluxes" from Refs. [7,8]. In our analyses we consider two different hypotheses for DUNE's data-collection strategy. We (a) assume 3.5 years of operation each in neutrino and antineutrino modes, for total of seven years. We do not perform any analyses considering only the high-energy mode, but (b) we perform analyses combining 3 + 3 years of neutrino and antineutrino mode with an additional year of data taking in the high-energy mode, also for a total of seven years. We refer to these two data collection strategies as "3.5 + 3.5" and "3 + 3 + 1", respectively.
When a ν τ interacts via a charged-current interaction, producing a τ lepton, the decay length of the τ is significantly smaller than the resolution of the DUNE detector, meaning that one must reconstruct the τ decay products in order to classify the incoming neutrino as a ν τ . The authors of Ref. [9] investigated the capability of liquid argon detectors to identify ν τ events using hadronic τ decays, which also make up the largest branching fraction for the τ (around 65% [10]). This approach takes advantage of the detector's capability to identify pions and kaons and to measure their kinematic properties. Understanding the usefulness of the leptonic τ decays, which are, naively, heavily contaminated by charged-current ν µ or ν e scattering, is beyond the scope of this manuscript. Building on Ref. [9], efforts within the DUNE collaboration [11] have optimistically estimated that one can isolate a ν τ -rich event sample where 30% of hadronically-decaying τ events are successfully identified while only 0.5% of neutral-current background events remain. More sophisticated analyses are currently under intense investigation by the collaboration. We restrict ourselves to this rather optimistic scenario throughout this work, and hope that the results discussed here will help motivate more detailed studies of the ν τ -reconstruction capabilities of DUNE and other liquid argon detectors.
When performing an oscillation analysis, it is paramount to measure the energy of the incoming neutrino. This is relatively straightforward † for charged-current events induced by ν e and ν µ , where the outgoing charged lepton is well-reconstructed. However, τ leptons are impossible to reconstruct perfectly, as their decays necessarily involve at least one outgoing neutrino. We have constructed a migration matrix mapping the true neutrino energy E true ν to the reconstructed neutrino energy E reco ν , parameterizing it with a bias in the preferred reconstructed energy as a function of the true energy, as well as an uncertainty on the reconstructed energy. This is done by simulating final-state hadronic τ -decays using MadGraph [14] for different τ energies and computing the assumed-to-be observed energy in the hadronic system. Our simulations are consistent with the following simple picture. For a given E true ν , the reconstructed energy follows a Gaussian distribution with a mean value bE true ν and width σ E = rE true ν , where b is the bias and r is the resolution of the measurement. Our simulations point to b ≈ 45% and r ≈ 25%. While the bias does not have a significant impact on the quantitative results presented in the following section, the resolution does; more detailed studies indicate that the DUNE collaboration will be able to achieve 8% r 25%, and we choose the conservative upper limit [15]. Fig. 1 depicts the migration matrix we have obtained and will use for all forthcoming analyses. We have also performed a more realistic simulation using the Genie3.0 software [16], and the results we obtained agree with Fig. 1.
Virtually no ν τ are produced at the neutrino source. In the DUNE detector, however, the three-massiveneutrinos paradigm predicts a healthy ν τ flux, mostly from ν µ → ν τ oscillations over the 1300 km baseline. Assuming the three-massive-neutrinos paradigm, the SM ν τ scattering cross section, and incorporating the expected identification capability for hadronic τ events, we calculate the expected number of reconstructed events as a function of neutrino energy for each mode of beam operation. Fig. 2 depicts the expected number of signal events for the neutrino mode (left), the antineutrino mode (center), and the high-energy mode (right). Here, we assumed the following values for the oscillation parameters [10], in agreement with the most recent global-fit results obtained by the NuFit collaboration [17]: sin 2 θ 12 = 0.310, sin 2 θ 13 = 0.02240, sin 2 θ 23 = 0.582, δ CP = 217 • = −2.50 rad, † Albeit not at all trivial and still under intense investigation, see, for example, Refs. [12,13]. Unless otherwise stated, we will assume these oscillation parameters to be the true values in our analyses going forward. We discuss the details on the ν µ → ν τ oscillation probability in the next section. Fig. 2 depicts both "smeared" (solid histograms) and "unsmeared" (dashed histograms) event yields using the energy migration matrix discussed above. We divide the simulated data in energy bins of constant width ∆E ν = 0.5 GeV, between 0 and 20 GeV, for our analyses. In each panel, we show the contribution due to neutrinos in green and antineutrinos in orange. Each distribution has been normalized to the expected runtime in each mode, 3.5 years for neutrino and antineutrino modes and 1 year for high-energy mode. Fig. 3 depicts the expected event yields -stacked -taking neutral-current backgrounds into account, along with the energy smearing. In seven years of running, we anticipate a healthy, relatively clean ν τappearance sample. In the 3.5 + 3.5 case, the sample includes over 200 events, while in the 3 + 3 + 1 running scheme we expect over 300 events. In both cases, neutral-current backgrounds make up less than 50% of all events [11]. We display stacked histograms of background (grey), events from ν τ producing τ + leptons (orange), and events from ν τ producing τ − (green). Each distribution has been normalized to the expected runtime in each mode, 3.5 years for neutrino and antineutrino modes and 1 year for high-energy mode.
From time to time, we will compare the results obtained using our simulated ν τ data sample with those obtained using samples of ν µ → ν µ disappearance and ν µ → ν e appearance (and their CP-conjugated channels). The details of those simulations are provided in Refs. [18][19][20][21][22] with fluxes updated to match those in Refs. [7,8]. The expected yields are presented in Appendix A.
In our analyses, we include a 25% normalization uncertainty on the number of signal events, a conservative assumption related to systematic uncertainties regarding the neutrino flux, cross section, etc. We find, in practice, that this uncertainty does not have a strong impact on the results discussed in the next section.
Assuming the three-massive-neutrinos paradigm is correct, one can use the ν τ sample to measure the ν τ charged-current cross section on argon. Excluding the 25% normalization uncertainty mentioned in the last paragraph ‡ and holding oscillation parameters fixed, we find that DUNE can measure the fluxaveraged ν τ charged-current cross section on argon to be σ = σ 0 (1 ± 0.07) (or σ = σ 0 (1 +0. 15 −0.14 ) at the two sigma level), where σ 0 is the SM expectation, after seven years of data taking (3.5 years in neutrino and antineutrino mode, respectively). This can be compared to current results from OPERA, and Super-Kamiokande, σ = σ 0 (1.47 ± 0.32) [2]. Note that these are not apples-to-apples comparisons. The expected flux-averaged neutrino energy for DUNE is E ν 4 GeV, compared to the one for OPERA E ν 10 GeV and Super-Kamiokande E ν 7 GeV. ‡ When measuring the charged-current cross section, one should not include systematic effects related to the cross section. Here, for simplicity, we are assuming that other normalization-related uncertainties will be constrained in a variety of ways, including measurements of νµ disappearance and νe appearance.

III. PHYSICS WITH THE BEAM ν τ SAMPLE
In this section, we simulate and analyze the ν τ sample that is expected to be collected at DUNE after either 3.5 years of running in the neutrino mode plus 3.5 years of running in the antineutrino mode (3.5 + 3.5) or 3 years of running in both the neutrino and antineutrino modes plus 1 year of running in the high energy mode (3 + 1 + 1), as discussed in the previous section. The simulated data are consistent with the three-massive-neutrinos paradigm, where the input oscillation parameters are set to the values listed in Eq. (II.1); the associated event yields for ν τ -appearance are depicted in Figs. 2 and 3 while the ones for disappearance and ν e -appearance are depicted in Appendix A. When discussing oscillation parameters, following tradition, we will refer to ∆m 2 21 and sin 2 θ 12 as the 'solar parameters,' while ∆m 2 31 and sin 2 θ 23 are referred to as the 'atmospheric parameters. ' We will discuss four distinct hypotheses: the three-massive-neutrinos paradigm (Sec. III A), a nonunitary lepton mixing matrix, in the context of heavy § new neutrino states (Sec. III B), the existence of a fourth, light neutrino (Sec. III C), and non-standard neutral-current neutrino-matter interactions (Sec. III D). We do not advocate for or favor any particular hypothesis for the physics that may exist beyond the three-massive-neutrinos paradigm. Instead, these are to be viewed as different, quantifiable modifications to the neutrino oscillation probabilities that serve as proxies for what may ultimately turn out to be the new physics. All hypotheses have been heavily scrutinized in the past so it is easy to make comparisons between different neutrino oscillation experiments and different oscillation channels.

A. Three-Massive-Neutrinos Paradigm
Assuming the three-massive-neutrinos paradigm, in vacuum, for the LBNF-DUNE baseline L = 1300 km and neutrino energies above τ threshold (E ν 3.4 GeV), such that all reconstructed ν τ appearance events occur above the first oscillation maximum. The oscillation amplitude is 4|U τ 3 | 2 |U µ3 | 2 = cos 4 θ 13 sin 2 2θ 23 ∼ 0.95, using the NuFit results listed in Eqs. (II.1). The subleading terms include the "solar" oscillations and the interference term. Here, the subleading terms are indeed subleading: ∆m 2 21 /∆m 2 31 ∼ 0.03 and, unlike the case of ν µ → ν e oscillations, all relevant elements of the lepton mixing matrix are large. ¶ More quantitatively, for the energies of interest, we expect the interference term to be around a few percent of the leading term while the solar term is at the per-mille level. Matter effects, which we include in all numerical computations, are small and do not modify this picture. We provide more details in Appendix B.
In summary, as far as relevant DUNE energies and baseline are concerned, P (ν µ → ν τ ) depends, at leading order, only on two effective oscillation parameters: ∆m 2 31 and Furthermore, there is no access to even the first oscillation maximum. Under these circumstances, it is virtually impossible, using ν µ → ν τ , to determine the octant of θ 23 (whether θ 23 < π/4 or θ 23 > π/4), since P (ν µ → ν τ ) is, at leading order, invariant upon |U µ3 | 2 ↔ |U τ 3 | 2 . The same is true for the neutrino mass ordering since resolving the neutrino mass ordering with neutrino oscillations requires visible matter effects or the ability to resolve more than one mass-squared difference. Fig. 4 depicts the ν µ → ν τ oscillation probability. All oscillation parameters are fixed to the values listed in Eqs. (II.1), except for sin 2 θ 23 . The shaded region corresponds to sin 2 2θ µτ ∈ [0.5, 1]. We analyze the simulated 3.5 + 3.5 ν τ appearance data. Fig. 5(center) depicts the preferred regions of the sin 2 2θ µτ × ∆m 2 31 -plane at the one (solid) and three (dashed) sigma level. The star indicates the best-fit point. In the analysis, the solar oscillation parameters are held fixed, along with the neutrino mass ordering (normal), while the other mixing parameters orthogonal to sin 2 (2θ µτ ) -in practice, sin 2 θ 13 and δ CP -are marginalized. We do not include any other external information on the mixing parameters, including the current very precise measurement of sin 2 θ 13 by reactor neutrino experiments. We emphasize that only the ν τ appearance data from DUNE is used in the fit -the ν µ disappearance and the ν e appearance data are not included in this fit -so this is, except for the solar parameters, a ν τ -appearance-only measurement of neutrino oscillation parameters.
This result is easy to understand. Since the accessible neutrino energies are all above the first oscillation maximum, the ν τ data, to zeroth order, measures the combination sin 2 2θ µτ × (∆m 2 31 ) 2 . This behavior cuts-off at ∆m 2 31 ∼ 5 × 10 −3 eV 2 , when the oscillation maximum shifts significantly to energy values above τ -threshold.
We perform similar analyses of the simulated ν e -appearance and ν µ -disappearance samples. The results are depicted in Fig. 5(left) and Fig. 5(right), respectively. For the sake of comparison, we include these preferred regions in Fig. 5(center), making it clear that they correspond to different observables in the horizontal axis. The effective mixing angles are generalizations of Eq. (III.3),  The ν e and ν µ 3σ results are shown as inserts in the center panel. In these analyses, we fix sin 2 θ 12 and ∆m 2 21 to their best-fit values, fix the neutrino mass ordering to normal, and marginalize over sin 2 θ 13 and δ CP . In each panel, we show results assuming 3.5 years each of data collection in neutrino and antineutrino modes in dark colors, and results assuming 3 years in neutrino/antineutrino modes and 1 year in high-energy mode in light colors.
In both analyses the solar parameters and the neutrino mass ordering are held fixed, and we marginalize over all other orthogonal oscillation parameters. As with the ν τ -appearance analysis, we do not include other external information on the mixing parameters. Fig. 5 reveals that, in some sense, and as expected, the ν τ -appearance sample is much less precise than the ν e -appearance or ν µ -disappearance samples. On the other hand, the three data sets are virtually independent and, strictly speaking, measure different phenomena. The information they provide is complementary. Assuming the three-massive-neutrinos paradigm, however, the different sin 2 2θ µα , α = e, µ, τ are not independent but satisfy sin 2 2θ µe + sin 2 2θ µτ = sin 2 2θ µµ . This assumption would allow one to determine sin 2 2θ µτ at (better than) the percent level with DUNE ν e -appearance and ν µ -disappearance data only. We return to tests of the unitarity of the mixing matrix in the next subsection. Fig. 5 also depicts the results obtained in the 3 + 3 + 1 case (faint colors). These are virtually indistinguishable from those obtained in the 3.5 + 3.5 case (bold colors). In the ν τ -appearance channel, the larger sample size leads to marginally tighter allowed regions than in the 3.5 + 3.5 case.
Within the three-massive-neutrinos paradigm, ν τ -appearance can also be used to constrain some of the other oscillation parameters. Fig. 6 depicts the allowed regions of parameters space obtained after 3.5 + 3.5 years of DUNE data taking, including information from ν τ -appearance (green), ν e -appearance (blue), and ν µ -disappearance (red). Contours represent the one (dashed) and three (solid) sigma preferred regions in the sin 2 θ 13 × sin 2 θ 23 , sin 2 θ 13 × ∆m 2 31 , and sin 2 θ 23 × ∆m 2 31 planes. We do not depict twodimensional projections of the parameter space that include the δ CP -direction since the ν τ -appearance channel is virtually insensitive to δ CP . In all analyses, the solar parameters and the mass ordering are held fixed at their current best-fit values, while we marginalize over all other oscillation parameters not depicted. The stars indicate the best-fit points. The well-known complementarity of the ν e -appearance and ν µ -disappearance channels is apparent, while the impact of the ν τ -appearance channel is marginal when it comes to measuring oscillation parameters within the three-massive-neutrinos paradigm. We repeated this exercise in the 3 + 3 + 1 case and obtained results that are very similar to those of the 3.5 + 3.5 case. The 3 + 3 + 1 results are not depicted here.

B. Non-Unitary Mixing Matrix
One of the fundamental basic questions one can ask of neutrino oscillations is whether the lepton mixing matrix is unitary. Unitarity tests of the quark mixing matrix, for example, were a key component of the FIG. 6: Measurement potential of seven years of data collection at DUNE, assuming 3.5 years each in neutrino and antineutrino modes. We show measurements assuming separate analyses of the three oscillation channels: ν e Appearance (blue), ν τ Appearance (green), and ν µ Disappearance (red). We show 1σ and 3σ CL regions of these measurements. We fix sin 2 θ 12 and ∆m 2 21 to their best-fit values, assume the mass-ordering is known to be normal, and marginalize over δ CP . Results assuming 3 years each in neutrino/antineutrino modes and 1 year in high-energy mode are qualitatively similar. experimental and phenomenological quark flavor-physics program and stand as one of the landmarks of contemporary particle physics. Unitarity tests of the lepton mixing matrix are still rather unimpressive [24,25] and the fact that there is very little direct information concerning ν τ appearance -plus there is literally no information on ν τ -disappearance -contributes significantly to the current state of affairs.
The results discussed in the last section allow one to perform a simple, mostly model-independent, unitarity test. As discussed earlier, unitarity relates the effective mixing angles that govern, for the most part, the three oscillation channels accessible to DUNE and can be translated into the following sum rule: sin 2 2θ µτ + sin 2 2θ µe − sin 2 2θ µµ = 0. Fig. 5 reveals this sum rule can be tested at, roughly, the 20% level. We performed a fit of our simulated data, assuming it is consitent with the three-massive-neutrinos paradigm, for |U e3 | 2 + |U µ3 | 2 + |U τ 3 | 2 , for fixed values of the solar parameters, and the mass ordering, and marginalizing over the CP-odd phase δ CP and the mass-squared difference ∆m 2 31 . We find that DUNE should be able to measure, in the 3.5 + 3.5 case, a significant improvement over the current constraints [25]. The precision is dominated entirely by the ν τ -appearance data. If one allows for different values of ∆m 2 31 for the different appearance channels, we estimate DUNE would measure |U e3 | 2 + |U µ3 | 2 + |U τ 3 | 2 = 1 +0.07 −0.20 (1 σ) or 1 +0.15 −0.32 (3 σ). In order to assess more quantitively how well the unitarity of the mixing matrix can be probed, we adopt the following parameterization for the not-necessarily-unitary lepton mixing matrix [26]: where U is a unitary matrix. The diagonal α ii (i = 1, 2, 3) are real and the off-diagonal α ij ≡ |α ij |e −φ ij (i = j = 1, 2, 3) are complex, totaling nine additional free parameters. A unitary lepton mixing matrix corresponds to α ij = δ ij . The most conservative constraints on the off-diagonal parameters are of order |α ij | 10 −2 , while the diagonal parameters are constrained to be ( [27,28]. In certain model-dependent situations, more stringent constraints apply [29][30][31][32][33].
Simply stating that the lepton mixing matrix is not unitary, is not, strictly speaking, sufficient when it comes to computing the effects of nontrivial α ij on neutrino oscillations. Here, we will restrict the discussion to the case when the 3 × 3 mixing matrix is not unitary because there are more than three neutrino species and the new neutrino mass eigenstates are heavy enough not to be kinematically accessible in the LBNF beamline. This happens, roughly, for heavy neutrino masses heavier than the kaon mass. We also assume the new neutrino flavor eigenstates are sterile, i.e., they do participate in the electroweak interactions.
Under these circumstances, it is convenient to express the unitary n×n lepton mixing matrix (n = 3+k, where k is the number of extra neutrino species) in block form where V ≡ N U is a 3 × 3 matrix. In the mass basis, integrating out the heavy neutrino degrees of freedom (and the heavy gauge bosons), where i, j = 1, 2, 3, α, β = e, µ, τ , G F is the Fermi constant, and P L is the left-chiral projection operator. J µ CC and J µ NC are the relevant hadronic currents for the charged-current and neutral-current weak interactions, respectively. The fact that V is not unitary will impact neutrino production and lead to, for example, zero-baseline "flavor-change," and modify neutrino propagation through matter. Here, we are most interested in the latter and will comment briefly on the former. As far as matter effects are concerned, the modified effective Hamiltonian used to describe neutrino flavor-evolution as a function of the baseline is, in the mass-eigenstate basis, where A CC and A N C are the standard charged-current and neutral-current matter potentials, respectively. Keeping all this in mind, we analyze our simulated data samples allowing for the hypothesis that the N matrix is non-trivial. As in the previous subsection, we fix the neutrino mass ordering and the solar parameters, and marginalize over sin 2 θ 13 , sin 2 θ 23 , ∆m 2 31 , and δ CP . The oscillation probability P (ν µ → ν τ ) is most sensitive to the lowest row of N , the parameters α 31 , α 32 , and α 33 . For oscillations in vacuum and in the limit sin θ 13 → 0, there is, however, no sensitivity to the parameter α 31 in the ν τ appearance channel. We do not, therefore, expect sensitivity to |α 31 | competitive with the existing upper bound of a few ×10 −2 , and fix α 31 to zero in our analysis. We also allow the phase of α 32 to vary. Fig. 7 depicts (green contour) the expected sensitivity of the ν τ -appearance channel at DUNE in the |α 32 | × (1 − α 33 )-plane, in the 3.5 + 3.5 case. It also depicts the sensitivity of the ν µ disappearance channel (red contour), and the reach of a combined analysis of all channels (solid black contour). We also show expected sensitivity of the combined analysis assuming 3 + 3 + 1 years of data collection (dashed black contour). In the parameterization of non-unitarity using α ij , the parameters are not completely independent: the off-diagonal elements |α ij | must be less than the product 2 (1 − α ii )(1 − α jj ). The existing constraint (1 − α 22 ) < 2.2 × 10 −2 therefore excludes combinations of |α 32 | and (1 − α 33 ) in the grey shaded region of Fig. 7. The expected limit on |α 32 | is considerably weaker than existing limits, while that on (1 − α 33 ) is comparable to the O(10 −1 ) existing limit, even after one includes information from α 22 . For higher neutrino energies, the sensitivity to the unitarity-violating parameters is expected to increase. This is due to the increased importance of matter effects for higher neutrino energies and the fact that the matter potential depends non-trivially on the unitarity-violating parameters. Fig. 7 depicts (dashed black contour) the expected combined sensitivity in the α 32 × (1 − α 33 )-plane, in the 3 + 3 + 1 case. It is easy to see that in this case the sensitivity is significantly stronger than the one in the 3.5 + 3.5 case. Fig. 7 reveals that, in the 3 + 3 + 1 running scheme, we are sensitive to (1 − α 33 ) > 6 × 10 −2 at the 3σ confidence level.
As mentioned earlier, a non-unitary lepton mixing matrix, assuming it arises as hypothesized here, leads to non-trivial flavor-change in the limit where the baseline vanishes. This is a consequence of the fact that, for example, the neutrino produced in π + → µ + + ν is not orthogonal to the one that participates in ν + n → p + e − . As far as ν τ appearance is concerned, the search for ν τ candidates in the DUNE near detector complex would be sensitive to |α 33 | 2 |α 32 | 2 . The large neutrino flux, combined with the fact that there are no beam-related backgrounds for ν τ appearance in the DUNE near detector site, should translate into exquisite sensitivity to |α 33 | 2 |α 32 | 2 . The pursuit of this and other near-detector-related question, is beyond the aspirations of this manuscript. See Refs. [34][35][36] for discussions of non-unitarity effects at near detectors not involving ν τ appearance.

C. The Three Plus One Neutrino Hypothesis
In the previous subsection, we discussed the sensitivity of ν τ -appearance to the hypothesis that there are heavy neutrino states. If, instead, the new neutrino states are light -much lighter than the pionthey cannot be integrated out but must be included among the kinematically accessible neutrino mass eigenstates. The so-called 3+1 scenario has been the subject of immense phenomenological scrutiny, and remains a possible solution, albeit under significant stress, to the so-called short-baseline anomalies. In this subsection, we discuss how well the ν τ -appearance data sample to be collected at DUNE can constrain the 3+1 neutrino hypothesis, and how its information complements that which can be obtained via the ν e -appearance or the ν µ -disappearance channels.
In more detail, we add a fourth mass-eigenstate, with mass m 4 , to the three existing neutrino mass eigenstates. We further assume that the new interaction eigenstate is sterile, i.e., does not participate in charged-current or neutral-current weak interactions. In this case, the mixing matrix is 4 × 4, and, in vacuum and in the limit ∆m 2 12 → 0, where δ 24 is the new CP -violating phase associated with U µ4 (we use the parameterization from Ref. [18]). The first term in Eq. (III.10) is the standard "atmospheric" term, Eq. (III.1), while the extra contributions are proportional to |U µ4 U τ 4 | = sin θ 24 sin θ 34 cos 2 θ 14 cos θ 24 sin θ 24 sin θ 34 in the limit where all "new" mixing angles are small.
We simulate data as in the previous subsection and analyze it assuming the 3+1 neutrino hypothesis. Here, however, we include priors on ∆m 2 31 and |U µ3 | 2 from the T2K Experiment [37], |U e3 | 2 from the Daya Bay experiment [38], and the solar parameters (|U e2 | 2 and ∆m 2 21 ) from the solar experiments and KamLAND [17]. * Note that the priors are imposed on the magnitudes of some elements of the mixing matrix, not on the values of the mixing angles. For a more detailed discussion on how these bounds are imposed, see [18]. The reason for the priors is as follows. We want to gauge the impact of ν τ -appearance at DUNE without the benefit of information from the other oscillation channels. On the other hand, the ν τ -appearance channel, as we argued earlier, provides only very limited information regarding the dominant oscillation frequency and, for the most part, is unable to distinguish one from several oscillation frequencies.
The priors impose the constraint that there are two known oscillation frequencies -proportional to ∆m 2 31 * The reason for using a subset of existing bounds is one of convenience. These are the same priors we use in the next subsection, when we are interested in the subset of experiments that are not expected to be signifciantly impacted by non-standard neutrino-matter neutral-current interactions. and ∆m 2 21 and that both ν µ and ν e "see" these two frequencies. Both are well established facts. This way, we infer what we believe is a more meaningful estimate of the sensitivity. Fig. 8 (left) depicts the sensitivity, in the 3.5 + 3.5 case, of DUNE, in the 4|U µ4 | 2 |U τ 4 | 2 × ∆m 2 41 -plane. We fixed the "123" neutrino mass ordering -normal -and the solar parameters, and marginalize over all other oscillation parameters. We restrict the parameter space to m 2 4 > m 2 1 and restrict our analysis to values of ∆m 2 41 less than 0.1 eV 2 . For larger values of ∆m 2 41 , we expect ν τ -appearance searches in the near detector -beyond the scope of this manuscript -to play a decisive role. For ∆m 2 41 10 −4 eV 2 (or lower), the sensitivity of DUNE does not depend on ∆m 2 41 .
] ν τ 95% CL FIG. 8: Expected 95% CL sensitivity to the combination of mixing matrix elements 4|U µ4 | 2 |U τ 4 | 2 vs. ∆m 2 41 using only the ν τ appearance channel, assuming 3.5 years of data collection each in neutrino and antineutrino modes. This analysis includes priors on the solar parameters (|U e2 | 2 and ∆m 2 21 ) from the solar experiments and KamLAND [17], priors on ∆m 2 31 and |U µ3 | 2 from the T2K Experiment [37], and |U e3 | 2 from the Daya Bay experiment [38], as discussed in the text. It also fixes the mass ordering to be normal. All other parameters, including the additional mixing angles and CP -violating phases are marginalized over in our analysis.
Current data constrain |U τ 4 | only rather poorly and mostly indirectly. There are, however, significant constraints on |U µ4 |, mostly from ν µ disappearance searches. The MINOS and MINOS+ collaborations [39] constrain |U µ4 | 2 10 −2 for ∆m 2 41 10 −3 eV 2 and |U µ4 | 2 0.5 for ∆m 2 41 10 −3 eV 2 . Significantly better sensitivity is expected from the disappearance channel at DUNE [18,[40][41][42]. Matter effects also lead to some sensitivity to |U τ 4 | in the ν µ -disappearance and ν e -appearance channels [30,43,44]. Inspecting Fig. 8, one would be tempted to believe, naively, that the current (and future, assuming no discovery) bounds on |U µ4 | 2 would render the bounds on |U τ 4 | 2 from ν τ -appearance at DUNE trivial, especially for ∆m 2 41 10 −3 eV 2 . This is not correct. Even in the limit |U µ4 | 2 → 0, the ν τ -appearance data, assuming it is consistent with the three-massive-neutrinos paradigm, can rule out large |U τ 4 | 2 via the (unitarity) sum rule |U τ 4 | 2 ≤ 1 − |U τ 3 | 2 . The lower bound on |U τ 3 | 2 translates into a robust bound on |U τ 4 | 2 . Figure 9 depicts, in the 3.5 + 3.5 case, the sensitivity of DUNE to the fourth light-neutrino hypothesis, this time in the sin 2 θ 24 ×sin 2 θ 34 plane, for fixed values of the new mass-squared difference. In the left panel, ∆m 2 41 is low enough that oscillations driven by the new mass-squared difference have not yet developed at the DUNE far detector. In the middle panel, ∆m 2 41 is close to ∆m 2 31 , meaning oscillations are relevant at DUNE. In the right panel, the new oscillations have averaged out at the far detector. Information from the ν τ -appearance channel (green) complements that from the combined ν µ -disappearance and ν e -appearance channels (purple). The sensitivity of combined analyses is depicted in black. In the case of small sin 2 θ 24 , constraints from ν τ -appearance on sin 2 θ 34 always surpass those from the other oscillation channels, for small, intermediate, or large values of ∆m 2 41 . We note that the results depicted in Figure 9 are less sensitive to the priors on ∆m 2 31 , |U e3 | 2 , and |U µ3 | 2 than those depicted in Fig. 8. . We compare sensitivity using ν τ appearance data only (green) with a combined ν e appearance ν µ disappearance analysis (purple), as well as a joint analysis (black). Solar parameters are fixed to their best-fit values, the mass-ordering is fixed to normal, and all other unseen parameters have been marginalized over.
We repeated this exercise in the 3 + 3 + 1 case and obtained results -not presented here for the sake of conciseness -that are very similar to those of the 3.5 + 3.5 case.
In order to further illustrate the impact of the ν τ -appearance data on the DUNE's capability to test the 3+1 hypothesis, we addressed a slightly different question: how does the sensitivity of DUNE change once ν τ -appearance data are added to the more traditional ν e -appearance and ν µ -disapparance samples? We concentrate our discussion on |U τ 4 | 2 ∝ sin 2 θ 34 . Fig. 10 depicts, in the 3.5 + 3.5 case, the sensitivy of DUNE in the sin 2 θ 34 × ∆m 2 41 -plane when analyzing only ν e -appearance and ν µ -disappearance (purple) or further including information from the ν τ -appearance data (black). We draw attention to our use of a linear scale for sin 2 θ 34 . Here, unlike the analysis that led to Figs. 8, and 9, we only make use of external priors on |U e2 | 2 and ∆m 2 21 [17]. The impact of the ν τ -appearance data is nontrivial.

D. Non-standard Neutrino Interactions
New interactions between neutrinos and ordinary matter can also modify neutrino propagation in a way that is testable at long-baseline experiments. In this subsection, we consider additional neutral-currentlike interactions between neutrinos and matter (electrons, up quarks, and down quarks) -non-standard neutrino interactions (NSI) -that modify the matter potential for neutrino oscillations. The modified matter potential is traditionally written, in the flavor basis, as ] ν e + ν µ ν e + ν µ + ν τ FIG. 10: Expected 95% CL sensitivity in the sin 2 θ 34 × ∆m 2 41 -plane, analyzing only ν e -appearance and ν µdisappearance data (purple) or further including the ν τ -appearance channel (black). Note the linear scaling on sin 2 θ 34 . We include external priors on the solar parameters (|U e2 | 2 and ∆m 2 21 ) from the solar experiments and Kam-LAND [17] and fix the mass ordering to be normal. All other parameters, including the additional mixing angles and CP -violating phases are marginalized over.
propagation close to the surface of the Earth, A 6 × 10 −4 eV 2 /GeV. This is to be compared to the different oscillation frequencies, proportional to ∆m 2 21 /E and ∆m 2 31 /E. The NSI parameters are complex, while the flavor-propagation Hamiltonian is Hermitian, leading to a total of nine free parameters. Neutrino oscillations are insensitive to an overall shift to all the eigenvalues of the propagation Hamiltonian, so we may subtract off any NSI parameter times the identity -in practice this is typically µµ , as external constraints on this parameter are quite strong [52]. Hence, limits from oscillations are usually derived for ee − µµ and τ τ − µµ . In practice, for the sake of DUNE, the existing strong limits on µµ imply that there is effectively no difference between the limits on αα − µµ and αα (α = e, τ ), so we will make no distinction between the two in what follows. For oscillation-based probes † , the existing constraints on NSI parameters are at the level of O(10 −1 ) − O(10), the weakest being on the ατ parameters.
The sensitivity of an experiment to NSI parameters is dependent on which oscillation probability is being measured and, in general, the channel P (ν α → ν β ) is most strongly dependent on the αβ parameters (α = e, µ, τ ). Since the ν e appearance measurements at DUNE are quite sensitive to matter effectsthe 1 in the ee-component of Eq. (III.11) -we expect DUNE to be sensitive to O(0.1) values of the NSI parameters. Refs. [19,58] explored the ability of the ν e appearance and ν µ disappearance channels to probe these NSI parameters in great detail.
Writing the oscillation probabilities in matter with nonzero αβ is a nontrivial exercise, and many perturbative approaches have been developed to express the probabilities analytically [53,59,60]. In general, for the baselines and energies of interest at DUNE, the probability P (ν µ → ν τ ) depends predominantly on the parameters µτ and τ τ . Moreover, the oscillation probability is more sensitive to changes to the real part of µτ than the imaginary part, so we expect sensitivity to be weakest when φ µτ , the phase associated with µτ is π/2 or 3π/2. In the ν τ appearance channel, we expect sensitivity to eµ and eτ only if the parameters are large ( 1). Additionally, in the ee 1 limit, ν e becomes an eigenstate of the flavor-evolution Hamiltonian, and all sensitivity to eµ and eτ will vanish in the ν τ appearance channel. For this reason, we set ee = 0 going forward. Fig. 11 depicts the expected 3σ sensitivity of the DUNE ν τ appearance channel (in green) to the NSI parameters | eµ |, | eτ |, | µτ |, and τ τ in the 3.5 + 3.5 case. Degeneracies between NSI parameters and the three-massive-neutrinos-paradigm parameters have been extensively studied in the literature [19,21,33,46,47,49,51,58,[61][62][63][64][65][66][67][68]. To reduce the effect of such degeneracies on this analysis, we have included as priors information from existing neutrino oscillation searches that are generally insensitive to NSI. As in previous sections, we simply fix the solar neutrino parameters to their best-fit points from NuFit and assume the mass ordering is normal. Here, however, we also include a prior on sin 2 θ 13 = 0.02240 ± 0.00066, since the results of the Daya Bay experiment are effectively insensitive to matter effects [38]. We also include priors from the T2K experiment on the measurements of sin 2 θ 23 = 0.526±0.036 and ∆m 2 31 = (2.537±0.071)×10 −3 eV 2 [37]. In our analysis, we keep the phases φ αβ of the off-diagonal NSI parameters free, as well as sin 2 θ 13 , sin 2 θ 23 , δ CP , ∆m 2 31 , along with the parameters depicted in Fig. 11. For comparison, Fig. 11 also depicts the sensitivity of the ν e appearance channel (blue) or the ν µ disappearance channel (red) to the NSI parameters, subject to the same priors as in the ν τ channel. In black, we display the results of a combined analysis using all three channels. We see that, generically, these channels are more powerful than the ν τ channel, however sensitivities of individual channels can be comparable (for instance, | µτ | vs. τ τ , comparing ν e and ν τ appearance channels). As with other probes, we stress the importance of performing independent analyses of the diferent channels as a way of cross-checking results.
We also repeated this analysis assuming 3 years of neutrino beam, 3 years of antineutrino beam, and 1 year of the high-energy, ν τ optimized beam, and found that the sensitivity was comparable to what is displayed in Fig. 11. The only noticeable difference is the 3 + 3 + 1 analysis is slightly more sensitive to | µτ | than the 3.5 + 3.5 one. This is due to the large number of events in the ν µ disappearance channel at large neutrino energy, where effects of | µτ | are more pronounced. We choose not to display the results of this analysis due to their similarity to those depicted in Fig. 11.

IV. DISCUSSION & CONCLUSIONS
As of today, only a handful of ν τ events have been directly observed via the standard model chargedcurrent weak interactions. The LBNF-DUNE experimental setup is expected to collect and isolate a ν τ -enriched sample in the DUNE far detector that dwarfs all existing data samples in both its size and purity. Here we explored the future physics sensitivity of seven years of ν τ -appearance searches in the DUNE far detector.
The collection of a ν τ -enriched sample is a nontrivial task. The high τ -production threshold, a consequence of the heavy τ -lepton mass, implies relatively very small statistics, and the fact that the τ leptons decay promptly and only semi-visibly implies that reconstructing the ν τ energy is difficult and that the rejection of neutral-current backgrounds is challenging. These challenges are reflected in the results presented here. The sensitivity of the ν τ -appearance channel, whenever comparisons are meaningful, is markedly inferior to that of the ν e -appearance or the ν µ -disappearance ones as one can readily confirm in Figs. 6 or 11. On the other hand, we emphasize the complementarity of the different oscillation channels. Different hypotheses about the physics responsible for neutrino oscillations, including the very successful three-massive-neutrinos paradigm, imply different correlations among oscillation channels and, in many cases, the ν τ -appearance channel provides information that cannot be accessed in other ways.
Another interesting side effect of the large τ -production threshold is that, in the case of the LBNF-DUNE setup, all ν τ -appearance events are above the first oscillation maximum, see Fig. 4 2 31 = (+2.537 ± 0.071) × 10 −3 eV 2 (free, ordering fixed) δ CP = −2.496 rad = 217 • (free) ee = µµ = 0 (fixed) eµ = eτ = µτ = τ τ = 0 (free) Phases φ αβ free FIG. 11: Expected sensitivity using ν τ appearance data (green), ν e appearance data (blue), ν µ disappearance data (red), and a combined analysis of all three (black) at DUNE assuming 3.5 years each of neutrino and antineutrino modes. As shown in the figure and explained in the text, we have included priors on sin 2 θ 13 , sin 2 θ 23 , and ∆m 2 31 , have fixed sin 2 θ 12 and ∆m 2 21 to their best-fit values, and assumed the mass ordering is known to be normal. The NSI phases have also been marginalized over. that one is not able to explore the oscillatory behavior of the oscillation probability, a fact that translates into the strong correlation between ∆m 2 31 and sin 2 2θ µτ in Fig. 5. This fact also partially explains why, in the case of the three-massive-neutrinos paradigm, data taking in the high energy configuration does not lead to qualitatively better results. To better illustrate this point, we repeated the exercise that led to Fig. 5 assuming that the LBNF-DUNE baseline is L = 2000 km. The associated measurement of sin 2 2θ µτ and ∆m 2 31 is depicted in Fig. 12 for both 3.5 years of running in the neutrino mode and 3.5 years of running in the antineutrino mode, and for 3 years in both neutrino and antineutrino modes and one year in the high energy beam configuration. Even though the neutrino fluxes are smaller at 2000 km (relative to 1300 km), the "measurement" of the oscillation parameters is more precise for the longer baseline. Furthermore, the 3+3+1 case allows for a slightly more precise measurement than the 3.5+3.5 case. We explored the impact of the ν τ -appearance channel on testing different hypotheses: the three-massive neutrinos paradigm, the existence of new heavy neutrino degrees of freedom (nonunitarity), the existence of new light neutrino degrees of freedom (3+1 neutrino oscillations), and the existence of new neutrino-matter neutral-current interactions (NSI). Our results are summarized in Figs. 6, 7, 8, 9, 10, and 11. In many cases we also depicted the sensitivity of the ν e -appearance or the ν µ -disappearance channels. It is always instructive to appreciate the limitations of the different oscillation channels and how they complement one another. We also compared two different running schemes for LBNF-DUNE -the 3.5+3.5 case and the 3+3+1 case -in order to gauge the benefits of running LBNF in the so-called high energy configuration. In some cases, see, for example, Fig. 7, the larger ν τ -appearance sample leads to better sensitivity, while in other cases the two running strategies are practically identical. Overall, our study does not unambiguously indicate that running part of the time in high-energy mode is qualitatively better than running in the so-called CP-optimized mode [7,8].
While there is a guarantee flux of ν τ in the DUNE far detector, one can also search for ν τ in the DUNE near detector. In this case, the neutrino flux -combining all flavors -is huge, while the beam-related backgrounds are zero: virtually no ν τ are produced as a direct consequence of the original beam-target collision. This allows one to perform "background free" searches for ν τ -appearance and, we expect, probe some of the scenarios discussed here with great sensitivity. Similar searches are currently being pursued by the NOνA collaboration [69]. This subject is beyond the aspirations of this manuscript. When it comes to ν τ -appearance, there are also neutrino sources for the DUNE far detector other than the LBNF beam. The atmospheric ν τ flux is large and spans a few energy decades. Some of the challenges of observing atmospheric ν τ in large liquid argon detectors, including the unknown neutrino direction and the absence of timing information, were explored in [9], and the subject is currently under investigation by the DUNE collaboration.
Here we present the details of our simulations for the ν µ → ν e appearance and ν µ → ν µ disappearance channel analyses, updated from Refs. [18][19][20][21][22]. As in the case of the expected ν τ -yields, discussed in Sec. II, the LBNF-DUNE neutrino and antineutrino fluxes match those available in Refs. [7,8] and when computing the effects of neutrino oscillations we have used the oscillation parameters listed in Eq. (II.1). Fig. 13 depicts signal (blue) and background (black) event yields for the ν µ → ν e appearance channel as a function of the reconstructed neutrino energy, assuming 3.5 years of data collection in neutrino mode (left), 3.5 years antineutrino mode (center), and 1 year of high-energy mode (right). The data are organized into energy bins of width 0.125 GeV. The signal events are stacked on top of backgrounds, which consist of the following: opposite-sign background (ν µ →ν e ); intrinsic ν e beam-contamination; misidentified ν µ chargedcurrent events where a muon is misidentified as an electron; misidentified ν τ charged-current events; and neutral-current events, misidentified as ν e charged-current events. Likewise, in Fig. 14 depicts the expected signal (red) and background (black) event yields for the ν µ → ν µ disappearance channel for the same beam configurations as in Fig. 13. The data are also organized into energy bins of width 0.125 GeV. Here, the backgrounds considered are the opposite-sign background (ν µ →ν µ ); misidentified ν τ charged-current events where the τ is misidentified as µ; and neutral-current events misindentified as muon charged-current events. The distributions in Figs. 13 and 14 are as a function of reconstructed neutrino energy: we consider an energy resolution of σ E = 7% Eν 1 GeV + 3.5% Eν 1 GeV , consistent with the goals of the DUNE experiment [8]. When performing analyses including the appearance and disappearance channels, we include a normalization uncertainty of 5% as a nuisance parameter, marginalized over when presenting results. and hence A > ∆ 31 for neutrino energies above τ -threshold. It is easy to show that P (ν µ → ν τ ) = sin 2 2θ 23 sin ∆ 31 L 2 (B.8) In the limit A ∆ 31 (when ∆ 31 L 1, see Eq. (III.2)), P (ν µ → ν τ ) in matter is very well approximated by the vacuum expression, Eq. (III.1). We have verified numerically that matter effects are also negligible close to τ -threshold, thanks in no small part to the fact that sin 2 θ 13 is very small.
The inclusion of solar effects -i.e., including nonzero ∆m 2 21 -is much more cumbersome -see for example, [70] -and does not modify the conclusions reached here, as we verified both semi-analytically, using the results of [70], and numerically.