Current and Future Neutrino Oscillation Constraints on Leptonic Unitarity

The unitarity of the lepton mixing matrix is a critical assumption underlying the standard neutrino-mixing paradigm. However, many models seeking to explain the as-yet-unknown origin of neutrino masses predict deviations from unitarity in the mixing of the active neutrino states. Motivated by the prospect that future experiments may provide a precise measurement of the lepton mixing matrix, we revisit current constraints on unitarity violation from oscillation measurements and project how next-generation experiments will improve our current knowledge. With the next-generation data, the normalizations of all rows and columns of the lepton mixing matrix will be constrained to $\lesssim$10\% precision, with the $e$-row best measured at $\lesssim$1\% and the $\tau$-row worst measured at ${\sim}10\%$ precision. The measurements of the mixing matrix elements themselves will be improved on average by a factor of $3$. We highlight the complementarity of DUNE, T2HK, JUNO, and IceCube Upgrade for these improvements, as well as the importance of $\nu_\tau$ appearance measurements and sterile neutrino searches for tests of leptonic unitarity.


Introduction
With the discovery that neutrinos oscillate came a new understanding of the standard model (SM) of particle physics -neutrinos have mass and leptons mix. Many experiments have since been performed, with more planned, to deepen our understanding of the nature and origin of neutrino masses and their mixing. A coherent picture is forming regarding leptonic mixing and the three-massive-neutrinos paradigm through the experimental data gathered to date. However, open questions regarding the dynamics of the neutrino sector remain, with substantial room for new physics to provide answers.
Unitarity, the requirement that the matrix governing the transformation between two eigenbases satisfies U † U = U U † = I, forms the basis of our understanding of SM fermion mixing [1][2][3][4]. This theoretical paradigm has been thoroughly tested to great acclaim in the quark sector [5][6][7][8][9][10][11][12]. However, our understanding of the corresponding leptonic mixing matrix (LMM) remains limited [13][14][15]. The phenomenon of neutrino mixing predicates nonzero neutrino masses, and yet the SM does not provide a mechanism for such masses to exist. As a result, a plethora of models has been postulated to explain the origin of neutrino masses, and hence oscillations, involving new physics beyond the standard model (BSM). A key feature of many such models is that they predict the existence of new neutrino eigenstates, leading to non-unitarity of the active neutrino LMM used to characterize neutrino oscillations [16][17][18][19][20][21][22][23][24][25].
Many studies have been undertaken to study the effect of LMM non-unitarity and to determine existing and projected constraints on non-unitarity [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40]. Such constraints can be derived from analyzing a multitude of processes, such as decays involving leptons, and, crucially, neutrino oscillations. The latter are among the most theoretically clean probes of LMM unitarity. With this in mind, and given that future neutrino oscillation experiments will be capable of precise measurements, we revisit current constraints and project future constraints on the unitarity of the LMM from oscillation experiments. A previous exploration of oscillation constraints on LMM unitarity was performed in 2015 in Ref. [14], utilizing contemporary data. Experimental precision has since improved, with better precision expected in near-future experiments, motivating our in-depth study.
In this work, expanding on the set up in Ref. [14], we offer a more comprehensive perspective on leptonic unitarity. We explore a set of reasonable assumptions regarding the possible origin of unitarity violation and discuss how they can affect tests of unitarity. We also break down how different subsets of experiments contribute to the constraints on specific rows and columns of the LMM, highlighting the importance of sterile neutrino searches and the uniqueness of τ -appearance searches. While we are interested specifically in oscillation-based constraints on unitarity, we discuss other probes, and their model-dependence as well. We include all existing oscillation measurements that make major contributions to unitarity constraints, as well as projections for oscillation-based constraints through the next decade. These include the planned IceCube Upgrade, Jiangmen Underground Neutrino Observatory (JUNO), Deep Underground Neutrino Experiment (DUNE), and Tokai to Hyper-Kamiokande (T2HK) experiments. In our companion paper [15], we explored this combination of current and future data to address the unitarity constraints and CP violation present in the LMM through unitarity triangles, an approach familiarized by studies of the quark mixing matrix.
This manuscript is organized as follows. Section 2 introduces the formalisms we adopt when computing neutrino oscillations, including the theoretical assumptions one can adopt when performing an analysis of non-unitarity, and how these assumptions impact results. In Sections 3 and 4, we explain the current and future datasets included in our analyses, respectively. In Section 5, we present the primary results of our analyses in a number of ways, resulting in our constraints on the unitarity conditions U U † = U † U = I in Section 5.3. We consider some alternate assumptions that impact the results, and present the results in light of these alternate assumptions, in Section 6. Finally, in Section 7 we provide discussion on our results and conclude.
We also wish to highlight the results that are included in our appendices. In Appendix A, we discuss how non-oscillation probes, such as rare charged-lepton decays, can be used in certain scenarios to constrain the unitarity of the LMM. Appendix B derives neutrino oscillation probabilities (both for appearance and disappearance/survival) in vacuum when unitarity is not assumed. Appendix C demonstrates how, with subsets of data, we can constrain the elements of the LMM in the electron and muon rows separately. In Appendix D, we discuss the Bayesian approach used in many of our analyses, and the priors that enter this type of analysis. Appendix E includes the measurement of the phases present in the LMM, a parameterizationdependent measurement. Lastly, Appendix F offers some discussion regarding the LSND and MiniBooNE anomalies, whether they can be resolved in this framework, and how they may be tested in next-generation experiments.

Neutrino Oscillations and the Leptonic Mixing Matrix
In this section, we summarize the phenomenon of neutrino oscillations, and how the structure of the LMM enters the calculations for oscillation probabilities. We introduce the formalism we use throughout our analyses, which allows for the possibility that the LMM is not unitary. Given that we allow this possibility, we discuss the possible origins of the unitarity violation and different theoretical assumptions that map on to these different origins. These different theoretical assumptions will affect our analyses, and so we will spend considerable time discussing their effects.

Unitarity of the Leptonic Mixing Matrix
Neutrino oscillation studies are generally carried out assuming a unitary 3 × 3 mixing matrix for rotating between eigenstates of flavor and mass. However, this assumption only strictly holds in a rather limited number of models for neutrino masses, some of which suffer from fine-tuning issues.
In many models for neutrino masses, while there is non-unitarity of the lepton mixing matrix, it is expected to be small. For example, in a generic type-I see-saw scenario, the non-unitarity of the light-neutrino mixing comes from the mixing (angle squared) between the light-heavy states, which is proportional to the mass ratio between the light and heavy states (see Appendix A for details). We expect the deviation from unitarity to be at most O(10 −13 ) even for an O(TeV) seesaw scale.
There are nevertheless abundant examples of neutrino mass models that can lead to large non-unitarity. It has been shown by various groups that for non-trivial neutrino Yukawa textures, the see-saw mechanism can lead to substantial deviations from unitarity (see e.g. [41][42][43][44][45]). In addition, mass models invoking symmetry arguments may also produce large unitarity violation [19][20][21][22][23]. It is therefore important to test the unitarity of the 3 × 3 lepton mixing matrix with experimental data.
To test unitarity with oscillation data, one could adopt mathematical assumptions on the 3 × 3 matrix corresponding to different theoretical assumptions on the origin of the unitarity violation. To clearly describe our choices, let us first look at how a neutrino state is defined. A flavor eigenstate neutrino field ν α (x) can be written as a linear combination of mass eigenstate fields ν k (x): (2.1) The flavor index α ∈ e, µ, τ, . . . m includes the usual SM flavor fields, along with m − 3 possible additional right-handed (sterile) fields, while the mass index k ∈ 1, 2, 3, . . . n, allowing for n − 3 additional mass eigenstate fields. U αk is thus an m × n matrix. In an oscillation experiment, neutrinos are produced and detected via weak interactions. The produced state, neglecting the effect of neutrino masses, can be defined as [30,46] 1 Only active flavors participate in weak interactions, so the flavor index α ∈ e, µ, τ . However, it is important to note that the mass eigenstate index k does not run to n, the total number of mass eigenstates. Instead, the sum over k is only performed over the total number of kinematically accessible mass eigenstates [32,35,37]. For pion-decay sources which are used for most experiments, a conservative cutoff is to include all mass states below 140 MeV. 2 For models where additional sterile neutrinos are heavier than the electroweak scale, this sum truncates at k = 3. The normalization factor in Eq. (2.2) guarantees that the flavor states are always properly normalized: ν α |ν α = 1. Note that if U is not unitary, the flavor states are not necessarily orthogonal: ν α |ν β = δ αβ . In this work, we focus on understanding the structure of the 3 × 3 (sub-)matrix of the full m × n mixing matrix U. We will refer to the 3 × 3 mixing matrix as the LMM or U LMM , or simply U αk , where now α = e, µ, τ and k = 1, 2, 3. We will parameterize and derive oscillation formulae using only U . An implicit assumption we adopt at this step is that there exist no additional sterile states with masses below 140 MeV. In this case, oscillation measurements provide a direct test of unitarity (see Appendix A for discussions on the model-dependence of charged lepton decay searches). For sterile neutrinos with masses between 0.1-10 eV, dedicated searches for spectral distortions in oscillation experiments are very sensitive [49][50][51]. For heavier steriles that are still kinematically accessible, see Refs. [38,39] for detailed discussions on their oscillation signatures. Sterile neutrinos in other mass ranges can be probed e.g. via beta decay, meson decay and neutrino-less double beta decay. See Ref. [52] for comprehensive discussions of these experimental searches.
We organize our discussion of oscillations around the following three cases: 1. The "standard" case, where m = n = 3, U = U.
2. The "sub-matrix" case, where m = n > 3. U is unitary, and U is not.
3. The "agnostic" case, m 3 and/or n 3. U αk is not assumed to be unitary, and U is not.
The standard case is the most commonly adopted in oscillation studies. It is worth pointing out that even when this is phenomenological applicable, it nevertheless involves fields beyond the SM, highlighting the need for BSM physics to fully understand the neutrino sector. The sub-matrix scenario is the most commonly adopted in unitarity-violation studies. It applies to all cases when the unitarity violation is induced by the existence of new particles. The agnostic case is a peculiar one, where the full m × n mixing matrix is not assumed to be unitary. This is of course a difficult case to realize, as unitarity is one of the fundamental principles upon which theories are typically built. However, it is useful to consider this possibility so as to verify whether experimental data support the theoretical bias that U αk should be unitary. Additionally, Refs. [37,53] explored scenarios in which non-standard neutrino interactions during neutrino propagation through matter may be mapped on to the effects of a non-unitary mixing matrix. While this is a specific scenario, it is one in which the agnostic case applies, and provides motivation for adopting this case to allow for generality in the form of U αk . For the rest of this paper, we adopt the agnostic assumption as our default scenario, aimed to be the most conservative with our bounds on non-unitarity. We distinguish the sub-matrix and agnostic cases because the sub-matrix assumption imposes additional criteria on the structure of U and hence leads to more stringent bounds. We discuss what these criteria are and how they improve certain bounds throughout our analysis.

Mixing Matrix Parameterizations
In the standard scenario, U LMM is a 3 × 3 unitary matrix. It is well-known that in order to parameterize such a matrix, three angles and three complex phases are required. Two of the (Majorana) phases are irrelevant for neutrino oscillations, and are unphysical if neutrinos are Dirac particles. The standard parameterization employs three mixing angles, θ 12 , θ 13 , and θ 23 , and one complex phase δ CP . Often referred to as the PMNS [1,3] or PDG [9] parameterization, this form of the LMM is where s ij ≡ sin θ ij and c ij ≡ cos θ ij . The mixing angles are often referred to by the regime of neutrino oscillations in which they have been studied in the most detail: solar (θ 12 ), reactor (θ 13 ), and atmospheric (θ 23 ). A number of global fit efforts in the three-flavor hypothesis have been performed, leading to relatively precise understanding of the mixing angles under this hypothesis [54][55][56][57]. More generally, a complex 3 × 3 matrix U can be described by eighteen real parameters. There are 9 conditions for relating a generic complex matrix for leptonic mixing to a unitary one. These conditions can be obtained from the requirement that a unitary matrix satisfies U † U = I. This is equivalent to requiring that all columns of the matrix are normalized to one: as well as requiring that the column unitarity triangles close: Note that these are nine real constraints as t kl can be complex. Because U † U is Hermitian, the unitarity condition can equivalently be written as U U † = I, which can be translated to row normalization conditions: (2.6) and the closure of row unitarity triangles: For the general case where U is a non-unitary 3 × 3 matrix, the number of real parameters needed to describe the matrix for neutrino oscillation is 18 − 3 − 2 = 13, where 3 phases can be absorbed by charged lepton fields and 2 Majorana phases do not participate in oscillations. Equivalently, one can see that 13 parameters are required, as a unitary LMM would have 4 parameters, and the extension to include potential non-unitarity involves relaxing 9 unitarity conditions. The magnitudes of the elements of the mixing matrix are parameterization-independent, therefore we choose to adopt the following parameterization: (2.8) Here, we have nine magnitudes and four CP-violating phases. 3 Going forward, we refer to the parameterization given in Eq. (2.8) as the Magnitudes & Phases (MP) parameterization. Note that in this case, the row and column normalizations can be larger than 1, and the 13 parameters are completely independent of each other. This parameterization applies straightforwardly to the agnostic case described above.
When U LMM is unitary, we can relate the MP parameterization in Eq. (2.8) to the PMNS parameterization, with a straightforwardly obtainable correspondence between parameters. The phases can be related to the PMNS parameterization by using Jarlskog factors J αi , which are defined as where the ε are Levi-Civita tensors. It is straightforward to see that If U LMM is unitary, all constructible Jarlskog factors must be equal to each other and equal to the PMNS matrix Jarlskog invariant: 4 J PMNS = c 12 s 12 c 23 s 23 c 2 13 s 13 sin δ CP . (2.12) (2.11) allows for the simple derivation of a relation between the phases φ e2 , φ e3 , φ τ 2 , and φ τ 3 and PMNS parameters when U LMM = U PMNS . Finally, we briefly discuss the sub-matrix case, where U is a 3 × 3 sub-matrix of a larger unitary matrix U. This introduces two additional constraints on the structure of U : the row and column normalizations of LMM must not exceed unity: Further, by applying the Cauchy-Schwarz inequality on the vectors U αk , where α or k runs from 4, 5... n, we obtain the following inequalities [14] 5 : In this work, the bulk of our results will be presented under the minimal set of theoretical assumptions, corresponding to the agnostic case. Where we discuss sub-matrix case results, the conditions of Eqs. (2.13) and (2.15) are imposed on U LMM , and the comparison with the agnostic case will be analyzed. Such comparisons will appear throughout our analysis, as well as in Section 6. The most commonly adopted parameterization in the sub-matrix case [37,[58][59][60] is the following: When there are three active neutrinos and any number of sterile neutrinos, one can express α kl in terms of mixing angles and phases between active and sterile neutrino mixing. For the full expressions, see Ref [59].
Here, unitarity is achieved in the limit that α kl → δ kl . The off-diagonal α kl may be complex, so there are nine free parameters corresponding to the nine constraints discussed above. We note here that this parameterization is useful in that unitarity is obtained in a relatively simple limit, i.e., α kk → 1 and α k =l → 0, compared to the MP parameterization. However, it is not straightforward to map between the "α" 4 This can be used to test the unitarity of the LMM. For details, see Ref. [15]. 5 One can directly apply Cauchy-Schwarz inequality on the matrix U , which leads to the weaker conditions (2.14) These inequalities hold for both the sterile and agnostic cases. However, we know that LMM is at least very close to unitary, so these conditions are met for all viable parameter space.
parameterization and the constraints of unitarity, specifically the normalizations and closures of columns. As constructed, the map between the α kl and the normalization of rows and closures between two different rows is relatively simple: On the other hand, for the normalizations of the columns N k and the closures of the triangles between different columns t kl , such a mapping depends on both the α kl as well as the mixing angles and δ CP from the PMNS parameterization.

Oscillation Probabilities
We now review the oscillation formulae for different parameterizations. See Refs. [30,35] for further discussions. In vacuum, the mass states |ν k are eigenstates of the Hamiltonian, such that they form an orthogonal basis, ν k |ν l = δ kl , and evolve in time as A flavor state is created as Eq. (2.2). After traveling for time t, the flavor state is evolved to The oscillation probability for a neutrino of energy E produced as ν α and detected as ν β after propagating L = ct is therefore (2.25) Here, E k = m 2 k /2E ν . In order to better understand the behavior of the oscillation probability, we separate the discussion here into two cases, α = β (disappearance/survival probability) and α = β (appearance probability). Defining ∆ ij ≡ ∆m 2 ij L/2E ν , the disappearance/survival probability may be written as For the appearance probability, we define ϕ αβ k ≡ arg (U * αk U βk ) and ϕ αβ ≡ arg (t αβ ). The appearance probability may be written as When U is unitary, Eqs. (2.26) and (2.27) reduce to the standard oscillation formulae. In most neutrino experiments, one or two mass terms dominate the oscillation behavior due to the large hierarchy between ∆m 2 21 and ∆m 2 32 , such that Eqs. (2.26) and (2.27) can be simplified. These oscillation formulae are derived in Appendix B. We also provide the simplified expressions when either ∆m 2 21 or ∆m 2 31 is the dominant term of interest in Appendix B and comment on which experiments belong in each of these regimes.
Equation (2.27) reveals an interesting consequence of a non-unitary mixing matrix, namely the zerodistance effect. When L → 0 (or in an experiment where ∆ 21 and ∆ 31 1), the disappearance probability P (ν α → ν α ) → 1, but the appearance probability becomes This implies that searches for short-baseline anomalous appearance from one flavor eigenstate α to another β provides a direct constraint on the closure between the α and β rows. We discuss how these searches, typically interpreted in the context of searches for light, coherently-oscillating sterile neutrinos, may be applied to our scenario in Section 3.7.
Matter effects in the context of non-unitarity: Even though it is a good approximation to use the vacuum oscillation probabilities in Eqs. (2.26) and (2.27) for many oscillation experiments of interest (as well as providing useful analytic interpretation of results), matter effects are important for several existing experiments, and crucial for the future DUNE experiment. Interactions of neutrinos with matter as they traverse can be included by adding a potential to the Hamiltonian that governs the time-evolution of the neutrino states, which is diagonal in the flavor basis: where n e and n n are the electron and neutron density in the medium, respectively. This potential is rotated by U LMM into the mass basis, and combined with the (mass-basis-diagonal) energy values ∆ k1 /(2E ν ). Typically, the neutron density is removed as its contribution to the total Hamiltonian is proportional to the identity matrix, and represents a phase common to the propagation of all three neutrino states. However, since U † U = I is not assumed in our analysis, n n must be included in our calculation, as the phase is no longer common to all three propagating states. We note here that the inclusion of matter effects in light of nonunitary mixing depends strongly on the assumptions regarding any new neutrino states (for instance, whether they interact with matter via the standard weak interactions or any other new interactions, and what their masses are). The form of V αβ above is that obtained in the minimal-unitarity-violation context, in which the new physics scale is assumed to be much higher than the electroweak scale [30], and we adopt it for the remainder of our work.

Normalization Effects
In this subsection we discuss further non-trivial consequences of a non-unitary LMM. Much of this discussion is adapted from Ref. [30]. Without the assumption of the LMM being unitary, the expected flux of a given neutrino flavor (e.g. produced by pion decay-in-flight) will be modified from the unitary expectation. The same will be true of charged-current (CC) scattering cross sections, where non-unitarity of the LMM leads to deviations in the rate of charged leptons produced from interactions involving neutrinos of the corresponding flavor. The flux of neutrinos of flavor α, and the corresponding CC cross section, may be expressed (relative to their unitary, SM expectations) as The neutral-current cross section is also modified, (2.31) Neutrino oscillation experiments infer oscillation probabilities by measuring event rates and spectra, which are a convolution of fluxes, cross sections, and efficiencies of detection. The number of detected neutrinos of flavor β, n β , is given by Thus, if an experimental analysis is performed assuming SM predictions of fluxes and cross-sections as truth, the measurement of n β corresponds to an inferred measurement ofP (ν α → ν β ). Recalling the vacuum oscillation formula of Eq. (2.27), this leads to the conclusion that experiments making such assumptions are inferring the oscillation probability given bŷ such that these measurements are not sensitive to the normalization factors N α and N β . Experiments often use a near detector to measure the neutrino flux. This is the case in, e.g., DUNE. Normalization effects will therefore manifest themselves differently for appearance and disappearance probabilities, as well as for near-detector-only measurements compared with near-to-far ratio measurements. These effects are both taken into account in our analysis.
Let us first consider disappearance measurements at near detectors, i.e., sterile neutrino searches. These experiments measure an energy spectrum of events n ND α (E), where α is the flavor label: At L = 0, the oscillation probability is P αα (E; L = 0) = 1. Experiments that measure disappearance spectra at near detectors to constrain an oscillation probability rely on understanding of the (SM) predictions of Φ α and σ α , so the measured spectrum can be expressed as Therefore, with a sufficiently precise measurement of the spectrum, as well as understanding of the SMexpected flux and cross section, a constraint on N 2 α can be placed. These results are usually reported in the context of a limit on disappearance probabilities, i.e. P αα (E; L = 0) is close to 1 with some degree of confidence. In Section 3.7 we discuss how these reported limits map onto constraints on N α . When considering searches for anomalous short-baseline appearance of a flavor α to a flavor β, Eq. (2.35) must be modified accordingly, and we find that these searches are sensitive to thereby providing constraints on the closure of row unitarity triangles.
Let us now consider near-to-far-ratio measurements. These experiments infer a far-detector oscillation probability P αβ ≡ P (ν α → ν β ; E, L) by measuring a near detector spectrum n ND α as above, as well as a far detector spectrum n FD β , We may express the measured ratio of oscillation probabilities as For disappearance measurements where α = β, the measured probabilities are the true probabilities, e.g., Eq. (2.27) for the vacuum case. For appearance measurements where α = β, the measured probabilitiesP αβ are the true probabilities P αβ multiplied by a factor N β /N α .

Current Experiments and Statistical Treatment
Before turning to upcoming experiments, we first consider the most powerful set of existing experimental results that are sensitive to quantities of interest in the LMM. Our analysis is performed only over the experimental data which, when combined, provides the dominant sensitivity to the corresponding parameters in the mixing matrix. As such, it does not constitute a complete set of experimental results, but nevertheless captures our current knowledge of neutrino oscillation parameters. This section serves to describe the inputs to our fit, as well as the translation from experimental measurements in the PMNS parameterization to the MP parameterization we use for our analysis.
For most experiments, we take the reported event spectra and fit for LMM elements, i.e., mixing angles in the PMNS parameterization. Because our focus is on the LMM matrix elements and not the mass-squared splittings, we include the best measurement of the latter for each experiment as a Gaussian prior. We will specify throughout the values we take from each experiment. Included in this, we marginalize over the mass ordering (i.e., the sign of ∆m 2 31 ). While there is a long-standing tension between solar and reactor experiment measurements of ∆m 2 21 , we find that allowing the mass-squared splittings to vary does not impact the results of measuring the elements of the mixing matrix. Table 1 summarizes the experiments that enter our analysis. We show the dominant quantity to which each experiment is sensitive both when unitarity is assumed (the middle column, labeled "PMNS Quantity"), and when unitarity is not assumed (the right column, labeled "LMM Quantity").

Solar Neutrino Measurements
The Sudbury Neutrino Experiment (SNO) [61] measures solar neutrinos via neutral-current interactions. The oscillation probabilities for solar neutrinos can, to very good approximation, be calculated by considering an incoherent sum over neutrino mass eigenstates of their production in the sun and their scattering cross sections in a detector. Critically, this includes their journey from production to exiting the sun. The effective probability that a neutrino begins as an electron flavor eigenstate and exits as a mass eigenstate |ν k is related to the effective matrix-element-squared in propagation in matter, which we call | U ek | 2 . Since all of the results we consider are in the regime where matter effects dominate, we focus on this region, where we can express these mixing angles as | U ek | 2 = 0, (|U e1 | 2 + |U e2 | 2 )/N e , |U e3 | 2 /N e . Table 1. Quantities to which each experiment is sensitive: using the PMNS parameterization when unitarity is assumed (center column), using the MP parametrization when unitarity is not assumed (right column).

Experiment
PMNS Quantity LMM Quantity Solar Neutral Current 1 Solar Charged Current sin 2 θ 12 cos 4 θ 13 + sin 4 θ 13 Following discussion similar to that of Sec. 2.4, we can write the detected number of neutral current events as where we defineP The terms of the last line in Eq. (3.6) are small relative to those in the line above it, so we can express the measured oscillation probability at leading order asP NC (|U e1 | 2 + |U e2 | 2 )N 2 2 + |U e3 | 2 N 2 3 . This measurement is limited by uncertainties associated with the 8 B neutrino flux prediction from the standard solar model [62]. We conservatively assume that this probability is measured at the 25% level,P NC = 1 ± 0.25.
In addition, SNO, and Super-Kamiokande (Super-K) [63] measure solar neutrinos via CC interactions. The oscillation probability is entirely dominated by matter effects. When assuming unitarity, the measured survival probability can be expressed as P ee = sin 2 θ 12 cos 4 θ 13 + sin 4 θ 13 . We take the results from a preliminary joint SNO and Super-K analysis [64], which reports sin 2 θ 12 = 0.306 ± 0.014 with a constraint of sin 2 θ 13 = 0.0219±0.0014, and interpret it as a measurement of the survival probability P ee = 0.2932±0.0134. 6 If unitarity is not assumed, the survival probability is We also include the measured value of ∆m 2 21 = (6.11 ± 1.21) × 10 −5 eV 2 from the joint analysis [64].

KamLAND
The Kamioka Liquid Scintillator Antineutrino Detector (KamLAND) was an experiment that observed the oscillation of reactor electron antineutrinos ν e (with energies between 2-10 MeV) at distances of roughly 180 km. This allows KamLAND to be sensitive to the solar mass-squared splitting, measuring ∆m 2 21 = (7.50 ± 0.20) × 10 −5 eV 2 . We use a slightly older measurement from Ref. [65] to be consistent with the corresponding measurement of the oscillation probability. A more recent analysis measures this mass-squared splitting slightly more precisely [66], but this does not affect our results.
Appendix B of Ref. [65] gives weighted measurements and uncertainties on the oscillation probability P (ν e → ν e ) for different values of x ≡ sin 2 (∆ 21 /2) , where the averaging is performed over effective mixing angles in matter to incorporate matter effects in their calculations. If unitarity is assumed, this oscillation probability (since oscillations associated with ∆m 2 31 have averaged out in KamLAND) is written as [65] P ee = (cos 4 θ 13 + sin 4 θ 13 ) − cos 4 θ 13 sin 2 (2θ 12 ) x.
When not assuming unitarity, this becomes KamLAND does not use a near detector in its analysis, so we need to re-scale this oscillation probability by N 2 e . The probability measured by KamLAND is thereforê

Daya Bay
The Daya Bay experiment observes the disappearance of electron antineutrinos P (ν e → ν e ). Daya Bay operates in the regime of ∆ 21 1 and the coefficient of the dominant oscillation term is given by sin 2 (2θ 13 ) in the PMNS parameterization. Given the derivations in Appendix B, we apply Eq. (B.8), which shows that the disappearance probability (without assuming unitarity, and including near detector normalization effects) is The most recent measurement from Daya Bay is sin 2 (2θ 13 ) = 0.0856 ± 0.0029 [67], which we map on to a measurement of the coefficient in Eq. (3.11) in our fit. Daya Bay's measurement of the mass-squared-splitting ∆m 2 31 = (2.471 ± 0.070) × 10 −3 eV 2 is also included in our analysis [67].

OPERA
The Oscillation Project with Emulsion-tRacking Apparatus (OPERA) collaboration has completed data collection and has provided results of searches for ν µ → ν τ oscillations, where the neutrinos have a mean energy of 17 GeV and a baseline of 730 km [68]. Across four channels, the experiment observed 10 ν τ signal events with an expectation of 2.0 ± 0.4 background events, and 6.8 ± 1.4 signal events (assuming sin 2 (2θ 23 ) = 1 and ∆m 2 32 = 2.5 × 10 −3 eV 2 ). This measurement predominantly constrains the quantity |t µτ | 2 and the product |U µ3 | 2 |U τ 3 | 2 . For our analysis, we assumed that OPERA measures the oscillation probability P (ν µ → ν τ ) for a fixed energy of 17 GeV and baseline of 730 km. We compute the oscillation probability numerically (including matter effects) and multiply it by N µ N τ to account for the normalization effects discussed in Section 2.4. We find that this approximation reproduces the results reasonably well. Finally, we include OPERA's measurement of the mass-squared-splitting ∆m 2 32 = (2.7 ± 0.7) × 10 −3 eV 2 [68].

T2K
The Tokai to Kamioka (T2K) experiment has performed searches for both electron (anti-)neutrino appearance and muon (anti-)neutrino disappearance in a mostly muon (anti-)neutrino beam. We include all searches possible, making some simplifications. We use the preliminary results reported in Ref. [69], using the figures therein to extract the expected signal and background rates for ν e and ν e events as a function of oscillation parameters. 7 In Ref. [69], expected signal plus background rates are shown for the ν µ → ν e and ν µ → ν e channels for values of δ CP of −π/2, 0, π/2, and π. These are given for both the normal and inverted mass orderings, assuming the other oscillation parameters are fixed. We make the simple assumption that T2K measures these event rates for a fixed energy E T2K = 600 MeV at a fixed baseline length L T2K = 295 km, with a constant matter density of ρ T2K = 2.6 g/cm 3 [70]. Using the figures provided in Ref. [69], we arrive at the following expected event rates: Here, the values 21.90 and 10.66 are our extracted background rates for each of the two channels. In order to extract these predictions, we make the assumption that the background rates are independent of the other oscillation parameters. The pre-factors 1282.84, 179.59, can be interpreted as a weighted flux times cross-section for this particular energy, translating an oscillation probability into an expected number of signal events. Figure 1 presents the expected number of signal plus background events for these two different channels, given by the formulae in Eqs. (3.12)-(3.13), along with stars that indicate the values given by the figure in Ref. [69]. Note that, despite assuming a mono-energetic measurement, our curves intersect the stars nearly perfectly. Also shown in each panel of Fig. 1 is the observed number of events in each channel (again, from Ref. [69]), and its ±1σ statistical range. Since the statistical uncertainties are large, we only include them (and no other systematic uncertainties) in this T2K electron-neutrino appearance measurement. We compute the oscillation probabilities including matter effects numerically and multiply it by N e /N µ to account for the use of a near detector for T2K. For T2K's measurement of muon-neutrino and muon-antineutrino disappearance, we find that instead of assuming a mono-energetic measurement, that we obtain results more compatible with those of the collaboration if we assume a fixed measurement of the relevant coefficient of the disappearance probability. The disappearance probability is such that the relevant coefficient is We include T2K's combined ν µ and ν µ disappearance searches by assuming a measurement of C µµ = 1.0 ± 0.03, and find that it gives us results consistent with Refs. [69][70][71]. Finally, when analyzing results from T2K (when studying its appearance channels, disappearance channels, or both), we include a Gaussian prior on the mass squared splitting ∆m 2 32 = (2.49 ± 0.082) × 10 −3 eV 2 [69].

NOvA
The last long-baseline experiment we include is the NuMI Off-Axis ν e Appearance (NOvA) experiment. It operates at similar values of L/E ν as T2K, but at longer distance/higher energy. Like our T2K analysis, we assume that NOvA measures event rates at a fixed energy and baseline -we assume this to be E NOvA = 1.9 GeV, L NOvA = 810 km, with a constant matter density along the path of propagation of ρ NOvA = 2.84 g/cm 3 [72]. Using the preliminary results of Ref. [73], we extract expected signal and background rates for different sets of oscillation parameters, using our mono-energetic assumption and assuming that the backgrounds are independent of neutrino oscillations. 8 Reference [73] reports observed event rates for both neutrino and antineutrino appearance, as well as expected background and signal rates for a set of oscillation parameters. We take these expected rates and our mono-energetic assumption to infer our expected event rates, Instead of presenting the expected event rates for NOvA as a function of δ CP in separate panels, we choose to show a "bi-event" plot, which shows the two expected rates simultaneously, in Fig. 2. We have fixed ∆m 2 21 , sin 2 θ 12 , and sin 2 θ 13 as given by the values in the top-right of the panel. Each ellipse is generated by fixing (sin 2 θ 23 , ∆m 2 31 ) (given by the legend), and varying δ CP , while using Eqs. The four different ellipses correspond to different combinations of (sin 2 θ 23 , ∆m 2 31 ) given in the legend where we vary δ CP between −π and π. The black cross indicates the number of observed signal neutrino and antineutrino events with statistical uncertainties only. Colored markers on the solid green ellipse correspond to the expected values of event rates presented in Ref. [73] for δ CP = ±π (star), −π/2 (square), 0 (circle), and π/2 (triangle). mass ordering ∆m 2 31 > 0 corresponds to the best-fit according to Ref. [73]. Since corresponding values for NOvA's preferred values in the inverted mass ordering and/or lower octant of θ 23 are not provided, we simply choose values such that sin 2 (2θ 23 ) and |∆m 32 | 2 are constant for these choices. For the normal ordering, upper octant (green solid line) choice, we display as markers the four expected event rates corresponding to the figure shown in Ref. [73], displaying how well our results agree with the official collaboration ones. Finally, the black cross indicates the observed event rate of N NOvA νe = 82, N NOvA νe = 33 with statistical uncertainties. As with T2K, we do not include any systematic uncertainties in this portion of the analysis.
For NOvA's measurement of ν µ and ν µ disappearance, we find, similar to our analysis of T2K, that the reported results are more realistically matched if we simply assume a fixed measurement of the disappearance coefficient given in Eq. (3.15). Since NOvA measures this to be slightly away from maximal, we include this as a measurement C Dis. µµ = 0.99 ± 0.02, which replicates the results from Ref. [73,74] fairly well. When we analyze NOvA results, either appearance data alone, disappearance data alone, or combined, we include its measurement ∆m 2 32

Sterile Neutrino Searches
When unitarity is not assumed, there could be additional zero-distance effects. Given our current knowledge of the mass-squared-splittings ∆m 2 21 and ∆m 2 31 and assuming that there are no additional neutrinos, any experiment, including near detectors, that operates at L/E ν such that should see no neutrino oscillations, and therefore is sensitive to these zero-distance effects. In this regime of L/E ν , oscillation probabilities without assuming unitarity are Due to the nature of such experiments, none of the sterile neutrino searches employ a supplementary "near" detector. Based on our discussions in Section 2.4, the measured oscillation probabilities must be rescaled due to Monte Carlo predictions aŝ Sterile neutrino searches are typically carried out in the following way. If a fourth neutrino exists with a mass-squared splitting ∆m 2 41 ∆m 2 31 , and an experiment operates in the L/E ν regime given by Eq. (3.18), then it can search for oscillations given by the probabilities Limits or potential observations from these searches are presented in terms of sin 2 (2θ αβ ) and ∆m 2 41 (see, for example, Refs. [49][50][51]) For any particular experiment situated at a baseline L and measuring oscillations for a specific energy range, the oscillations associated with ∆m 2 41 become very rapid (as a function of E ν ) as ∆m 2 41 becomes larger and larger. Eventually, these oscillations average out, and the term sin 2 ∆m 2 41 L/4E ν → 1/2. In this averaged-out regime, the oscillation probabilities become These have the same lack of energy-dependence as searches for unitarity violation. Equating Eq. (3.25) with Eq. (3.21) therefore allows us to map constraints on sin 2 (2θ αα ) from sterile neutrino searches in the averaged-out regime onto constraints of N 2 α , while comparing Eq. (3.26) with Eq. (3.22) allows us to map sin 2 (2θ αβ ) onto |t αβ | 2 for α = β. Table 2 summarizes the null sterile neutrino searches included in our analysis: KARMEN, NOMAD, CHORUS, and MINOS/MINOS+.
We do not consider any sterile neutrino searches for P (ν e → ν e ) from reactor antineutrino experiments. The averaged-out regime of these searches, which is required to perform this mapping, depends on flux and cross section uncertainties to be well understood in a disappearance search. The overall flux of these searches is notoriously difficult to constrain [80][81][82][83], so these experiments do not place robust, strong limits in the high-∆m 2 41 regime. Both the Liquid Scintillator Neutrino Detector (LSND) [84,85] and MiniBooNE [86] experiments have observed an excess of electron-like events in the presence of a beam that is mostly ν µ (or ν µ ), which can be interpreted as a short-baseline oscillation with P (ν µ → ν e ) ≈ 2.6 × 10 −3 . A combined study of these two experiments favors P (ν µ → ν e ) = 0 at roughly 6σ. When analyzed in the context of a light sterile neutrino, the preferred parameter space is compatible with the averaged-out regime, however that is not where their best-fit point lies. We inspect the effect of including the favored |t µe | 2 = 0 preference from LSND/MiniBooNE in Appendix F. Table 2. Sterile neutrino searches included in our analysis, and the associated 90% CL limit on the effective mixing angle from the given experimental search.

Search
90% CL Limit Angle Constrained Unitarity Constraint

IceCube Upgrade
The IceCube experiment is capable of detecting atmospheric neutrinos over a broad range of energies, 1 GeV E ν 100 GeV. By measuring track-like (ν µ ) and cascade-like (ν e and ν τ ) events, IceCube is sensitive to effects of the mass-squared splitting ∆m 2 31 and oscillations associated with that splitting. Of importance for constraining leptonic unitarity is IceCube's ability to constrain the normalization of ν µ → ν τ appearance in its data sample. Currently, this measurement has a precision of ∼40% [93]. However, in the coming years, O(10%) precision will be attainable by considering either eight years of IceCube DeepCore data or one year of IceCube Upgrade data [94]. Given the regime of oscillations that IceCube measures, for ν µ → ν τ appearance it is dominantly sensitive to 4|U µ3 | 2 |U τ 3 | 2 /N 2 µ . We do not include any information from IceCube on this quantity in our current fits, but include a 10% measurement of that quantity in our future projections.

JUNO
JUNO is an upcoming reactor-based neutrino experiment (scheduled to start operation in 2022 [95]), where neutrinos from 10 different nuclear reactors travel for baselines of 53 km before reaching a 20 kiloton detector (with 2 additional reactors at a baseline of 200 km), comprised of liquid scintillator. JUNO is primarily sensitive to the atmospheric mass-squared splitting ∆m 2 31 , as well as being sensitive to the solar mass-squared splitting ∆m 2 21 . Its design goal is to measure ∆m 2 31 precisely enough to determine the mass ordering of the neutrinos, i.e. whether m 1 < m 2 < m 3 or m 3 < m 1 < m 2 . Because it operates in neither the regime ∆ 21 1 nor ∆ 31 1, when considering JUNO we must use the full oscillation probability in Eq. (B.6). The effects of the matter potential on the oscillation probability P (ν e → ν e ) are negligible for 9 With the recent update of oscillation data from T2K [69] and NOvA [73], the preference for the normal mass ordering (∆m 2 31 > 0) over the inverted mass ordering (∆m 2 31 < 0) has diminished [57,92]. We choose the best-fit point according to the normal ordering, given data not included in our fit, specifically Super-Kamiokande's atmospheric neutrino sample. sensitivity studies on the solar sector mixing [88], which is the main contribution from JUNO dataset to the global unitarity constraints, so we employ the form for said probability given in Eq. (B.6).
We simulate the expected event rate assuming six years of data collection, corresponding to a total of 1.2×10 5 signal events from ν e inverse beta-decay scattering in the JUNO detector [88,96]. Following Ref. [88], we include the following sources of systematic uncertainties in our simulation: correlated (among different reactors) flux uncertainty of 2%, uncorrelated flux uncertainty of 0.8% for each reactor, spectrum shape uncertainty of 1%. We do not include backgrounds in our simulation. Our projected sensitivity to sin 2 θ 12 in the standard three-flavor oscillation scenario is 0.42%, compared to the official collaboration sensitivity of 0.54%. We also simulate the JUNO experiment assuming it does not have a near detector, as current plans for the experiment do not incorporate one.
JUNO will in principle be sensitive to each element of the electron row of the LMM, |U e1 | 2 , |U e2 | 2 , and |U e3 | 2 . As we will show in Section 5, JUNO will allow for more precise measurements on |U e1 | 2 and |U e2 | 2 than current experiments, but will not improve on the precision of |U e3 | 2 achieved by Daya Bay.

DUNE
DUNE is a future beam-based neutrino experiment, scheduled to begin data collection in the late 2020s. It consists of an O(GeV) neutrino beam, consisting primarily of ν µ when operating in neutrino mode and ν µ when operating in anti-neutrino mode [97], and a 40-kton liquid argon far detector. The baseline length is 1300 km, meaning that DUNE operates near the regime ∆ 21 1. Nevertheless, it will be sensitive to the effects of ∆m 2 21 , allowing for a precise measurement of the CP-violating phase δ CP in the PMNS matrix, among other goals.
We consider three different beam-related channels when simulating DUNE, each effectively measuring a different neutrino oscillation probability: the electron-neutrino appearance channel, sensitive to P (ν µ → ν e ) (and its CP-conjugate); muon-neutrino disappearance/survival P (ν µ → ν µ ) (and its CP-conjugate); and tau-neutrino appearance, sensitive to P (ν µ → ν τ ) (and its CP-conjugate). To simulate the expected event rates for these different channels, we employ simulation code developed with Refs. [98][99][100][101][102] for ν e appearance and ν µ disappearance and Ref. [103] for ν τ appearance. For all of our simulations, we assume seven years of data collection with DUNE, divided evenly between operation in neutrino mode and antineutrino mode. We include signal and background normalization uncertainties (5% for the ν e appearance and ν µ disappearance channels, and 25% for ν τ appearance). As we will show in Section 5, DUNE will allow improve on the precision of the measurements made by NOvA and T2K for both ν µ disappearance and ν e appearance, as well as OPERA for ν τ appearance.
Finally, DUNE is capable of improving on existing measurements of the 8 B solar neutrino flux using ν e CC and elastic electron scattering [104]. We simplify the analysis by assuming that DUNE will be able to measure P ee ≡ |U e2 | 2 |U e1 | 2 + |U e2 | 2 + |U e3 | 4 at the 3% level, consistent with the more complete analysis of Ref. [104].

T2HK
In the next decade, T2K will be upgraded with a larger waterČerenkov detector and begin operating as T2HK. It will operate in a similar region of L/E ν as DUNE, albeit at a lower length and energy. This, along with the different detection mechanism, allows for tests between the results of the two experiments, and further power in validating (or discovering new physics beyond) the three-massive-neutrinos paradigm. T2HK will also collect a very large sample of atmospheric neutrinos, which we do not include in our analysis. Its beam-based program plans to collect data in a 1:3 ratio between neutrino and antineutrino modes. While T2HK intends to operate for ten years or longer, we rescale all of our expected signal and background rates to a data collection period of seven years to be consistent with our DUNE projections.
We perform simulations of the T2HK expected yields for ν µ → ν e appearance and ν µ → ν µ disappearance, consistent with Refs. [91,105], developed from Refs. [102,106]. As with our DUNE simulation, this has been modified to allow for a non-unitary LMM. We include expected signal and background yields in our simulations, along with energy resolutions discussed in Refs. [102,106], and signal and background normalization uncertainties of 5%.

Primary Results
Throughout this section, we present the results of analyzing the datasets described in Sections 3 and 4. We begin in Section 5.1 with consistency checks where unitarity is assumed, and determine whether different datasets agree on their measurements of different parameters. These consistency checks can serve as a simple test of unitarity. Subsequently, we abandon the unitary assumption and consider the agnostic case, adopting the MP parameterization, and present results for our full analysis in Sections 5.2 and 5.3. Section 5.2 presents how well we can currently constrain the LMM matrix element magnitudes |U αk | 2 and how much we can improve with future data. Section 5.3 demonstrates how well we can determine the normalization of the rows/columns, and closure of the row/column triangles of the LMM. To further understand how subsets of data contribute to the sensitivity in the electron and muon rows, see Appendix C.
Before presenting our results, we clarify the statistical approaches taken in our analyses. In Section 5.1 and Appendix C, the analyses only rely on a handful of parameters. We perform frequentist analyses, scanning a given likelihood function over the parameters of interest, and determining the confidence levels (CL) in these parameter spaces. In Sections 5.2 and 5.3, the analyses depend on 15 parameters, and we use the Bayesian inference tool pyMultiNest [107][108][109][110] to construct credible regions (CR) based on the posterior likelihood density. For further details, see Appendix D.

Simple Unitarity Constraints & Consistency Checks
One straightforward way to test whether the LMM is unitary is by analyzing different experimental measurements separately and checking for consistency. This is demonstrated conceptually in Fig. 3 for two pairs of mixing angles, sin 2 θ 13 vs. sin 2 θ 12 and sin 2 θ 13 vs. sin 2 θ 23 . 10 We assume the LMM is unitary, and interpret experimental measurements as combinations of PMNS mixing angles (see Table 1). If the LMM is indeed unitary, all of these measurements should meet at a single point in the sin 2 θ jk -sin 2 θ nl planes, while if it not unitary, an intersection is not guaranteed.
Before continuing, we note that the analytical expressions in Table 1 are good approximations of measurements near the best-fit regions of mixing parameters. If one deviates too far from best-fit regions, these approximations break down. For example, current data indicates small |U e3 | 2 , or sin 2 θ 13 if unitarity is assumed. In what follows, when analyzing different oscillation coefficients, we will generically allow sin 2 θ 13 to be large. This very likely produces results at large sin 2 θ 13 that are inconsistent with real experimental observations. Nevertheless, it is instructive to see how the results of these analytic estimates change over the full, potential range of the mixing angles. We indicate the region where these expressions likely break down, sin 2 θ 13 0.2, with dashed lines.
In the left panel of Fig. 3, we compare how three types of measurements can pin down sin 2 θ 13 vs. sin 2 θ 12 given infinite experimental precision, under the assumption that the full oscillation probabilities are dominated by the coefficients listed in Table 1. The three measurements we incorporate are: • P SBL ee : short-baseline measurement of P (ν e → ν e ) from a reactor neutrino experiment, e.g., Daya Bay. • P Solar ee : solar neutrinos measured via CC interactions.
• P LBL ee : long-baseline measurement of P (ν e → ν e ) from a reactor neutrino experiment, e.g., KamLAND. Whether all three measurements intersect is a test of the e-row normalization, i.e., whether N e = 1.
In the right panel of Fig. 3, we show similar hypothetical infinite-precision measurements of oscillation probabilities that are sensitive to a combination of sin 2 θ 13 and sin 2 θ 23 . They are the following: • P SBL ee as above.
• P LBL µµ : long-baseline measurement of P (ν µ → ν µ ) as performed by the current T2K and NOvA, or future DUNE and T2HK experiments.
• P LBL µτ : long-baseline measurement of P (ν µ → ν τ ) as performed by the current OPERA experiment, upcoming IceCube results, or future DUNE experiment.
From the right panel of Fig. 3, we see that long-baseline experiments alone do not suffice to resolve the octant degeneracy, i.e., whether sin 2 θ 23 is larger or smaller than 1/2, as the green, red, and orange curves meet in two locations. Reactor data providing a precise measurement of sin 2 θ 13 is needed to break the degeneracy. Also, we see that an infinitely precise measurement of the ν τ appearance probability is useful, but not necessary, to precisely determine the oscillation parameters. Figure 3 (right) allows for a test of the normalization of the third column of the LMM, i.e., whether N 3 = 1 when one relaxes the unitarity assumption. Referring to the right column of Table 1, given measurements from KamLAND, solar CC experiments, and Daya Bay, we see that |U e3 | 2 can be determined fairly robustly. Next, long-baseline ν µ → ν e appearance allows us to measure 4|U e3 If additional information (for instance from MINOS/MINOS+) on N 2 µ is obtained, this combination allows us to determine |U µ3 | 2 precisely. Finally, long-baseline ν µ → ν τ appearance is sensitive to 4|U µ3 | 2 |U τ 3 | 2 /N 2 µ , from which we can extract |U τ 3 | 2 . In tandem then, we can measure each of the elements |U α3 | 2 , which allows the placing of a constraint on N 3 . Note that the measurement capability of N 3 will be mostly limited by one's measurement of long-baseline ν µ → ν τ appearance, as it is the least well-constrained of these measurements currently. See the discussion around Fig. 6 for an illustration of this.
We now analyze how well the combination sin 2 θ 13 vs. sin 2 θ 12 is measured by current experiments, as well as prospects for near-future experiments. This result is presented in Fig. 4, where the left (right) panel demonstrates our current (expected future) knowledge of the two parameters. For easy comparison, the right panel also includes a joint fit of current data in black.     The shapes of the measurement regions in Fig. 4 match those in the left panel of Fig. 3, as expected. We see that JUNO will significantly improve the precision on measuring sin 2 θ 12 . However, it will not measure sin 2 θ 13 as precisely as Daya Bay. If we allow sin 2 θ 12 > 1/2, JUNO has an allowed region (analogous to KamLAND) that cannot be resolved by JUNO alone. In addition, DUNE will modestly improve on existing solar neutrino measurements.
Similarly, Fig. 5 demonstrates our knowledge of the combination sin 2 θ 13 vs. sin 2 θ 23 . In the left panel, because the T2K and NOvA experiments both measure ν µ → ν µ disappearance and ν µ → ν e appearance, and their measurements are qualitatively similar, we combine the two in each analysis, but separate by the two channels (all unseen parameters, including δ CP , are marginalized over in these analyses). Like with the current set of measurements shown in the left panel of Fig. 4, these measurements all agree at the 1σ level. DU NE/ T2H K P µµ DUNE/ T2HK P µe DUNE P µτ , U † U = I Figure 6. Projected measurements of sin 2 θ 13 vs. sin 2 θ 23 when unitarity is violated (N 3 ≈ 2). For DUNE's longbaseline measurement of P µτ (green), we simulate data assuming the underlying mixing matrix is non-unitary, and extract the measurement of these parameters assuming the matrix is unitary.
NOvA and T2K both see modestly higher event rates than expected in ν e appearance, driving the orange contours up and to the right relative to where the blue and red contours overlap, however, any tension is modest at best. DUNE and T2HK (like with NOvA and T2K, we combine these two experiments because of their similar sensitivity) will improve on the NOvA, T2K, and OPERA measurements in the right panel of Fig. 5, leading to the tighter contours of the right panel. If non-unitarity is present, then these regions may not overlap at high CL, as they do today.
Finally, we revisit the point raised when discussing the right panel of Fig. 3, regarding how one might use such a set of measurements to determine whether the LMM is unitary through their sensitivity to N 3 = 1. Here, we inject unitarity violation by making |U τ 3 | 2 significantly larger than it should be, making N 3 ≈ 2, beyond the edge of the currently allowed 3σ range. We then fit the simulated data using DUNE's ν τ appearance measurement while still assuming the mixing matrix is unitary. The resulting fit is shown in green in Fig. 6 -note the lack of overlap between the green (DUNE P µτ ), orange (DUNE/T2HK P µe ), and red (DUNE/T2HK P µµ ) measurements -indicating the need for a more thorough test of unitarity violation. If such a discrepancy arises, IceCube's measurement of 4|U µ3 | 2 |U τ 3 | 2 /N 2 µ could allow for verification of this arising from unitarity violation.
We caution the reader that the approach taken in Fig. 6 was to illustrate how channel combinations test unitarity. However, this is not the most robust way of testing unitarity, as the sensitivities of different measurements to unitarity violation are not easily disentangled. Furthermore, this framework does not accommodate sterile neutrino searches. An alternate example of how to test unitarity when analyzing data in the PMNS paradigm can be found in Ref. [15], where it was demonstrated how unitarity triangles ρ xy +iη xy can be used. In Ref. [15], we showed that, like here, separating analysis results by different oscillation channels can lead to inconsistent fits. In the following subsections, we carry out a global fit to the LMM that can directly test unitarity.  Results are generated using pyMultiNest [107][108][109][110]. All results here are obtained under the agnostic assumption.

Joint Measurement of All Matrix-Elements-Squared
In this subsection, we present the current constraints on the parameterization-independent |U αk | 2 , and project how well these will be constrained once future data from DUNE/T2HK, JUNO, and IceCube Upgrade are included. In Appendix E, we show how well the parameterization-dependent phases {φ e2 , φ e3 , φ τ 2 , φ τ 3 } are and will be constrained. Figure 7 displays the individual measurements of |U αk | 2 including current data (blue) and current data with future data (red). 11 Each panel displays the one-dimensional ∆χ 2 measurement of the element after marginalizing the 15-parameter fit down to the individual element. Here, we define ∆χ 2 = −2∆L, where L is the posterior likelihood obtained in our analysis. The improvement in |U e1 | 2 and |U e2 | 2 is driven predominantly by JUNO, while the improvement in the muon row is driven by DUNE and T2HK ν µ → ν µ measurements.The tau row improvements are driven by ν µ → ν τ measurements in IceCube and DUNE. The IceCube measurement will reach a higher precision than DUNE because of the large systematic uncertainties on neutrino-nucleus cross sections at DUNE's beam energies.
Using these results, we can determine the allowed 3σ CR for each of the nine mixing-matrix-elements squared according to the current data, as well as projections to including future data. The allowed ranges for these are As is evident in Fig. 7, the improvement is noticeable especially for the elements |U e1 | 2 , |U e2 | 2 , |U µ1 | 2 , and the τ row elements. This analysis is performed under the agnostic case -we compare these results with those obtained under the sub-matrix case in Section 6.2.

Constraining the Normalization and Closure Conditions with Current and Future Data
In this subsection, we check the consistency of data with the requisite conditions to determine whether the LMM is unitary. Specifically, we measure the row/column normalizations N α and N k and triangle closures t αβ (between two rows) and t kl (between two columns), using the same analyses as in the previous subsection. The left panel of Fig. 8 displays the results of this analysis, projecting down to two-dimensional CR measuring the row normalizations N e , N µ , and N τ at 95% and 99% credibility. We see that the analysis of all current data is consistent with unitarity for these values. Future data will lead to a modest improvement in the constraint on N µ , some improvement in N e , and significant improvement in N τ .
Similarly, the right panel of Fig. 8 presents the current constraints, as well as projected future ones, on the column normalizations N 1 , N 2 , and N 2 , at 95% and 99% credibility. The correlation between measurements of each pair of column normalizations is due to the fact that these constraints are limited by the measurement of the tau-row elements, |U τ k | 2 . Future data will improve the constraint on each column normalization by a factor of roughly 3. Table 3 summarizes the current and expected future measurements of the row and column normalizations of the LMM. Here, we give the current best-fit (maximum likelihood point) value of each normalization, as well as the extents of its current 3σ CR. We also show the projected future 3σ CR, assuming a true value of N X = 1, demonstrating the improvement attributable to future data. We note here that we obtain a 0.8% 1σ constraint on N e when we project to including future data including JUNO, DUNE, and T2HK. This is in contrast with Refs. [13,88] which project a 1.2% measurement of N e once JUNO is considered   Figure 8. Left: constraints (and projected constraints) on the row normalizations N e , N µ , and N τ at 95% (dark) and 99% (faint colors) credibility. Right: constraints on the column normalizations N 1 , N 2 , and N 3 at 95% (dark) and 99% (faint colors) credibility. All results here are obtained under the agnostic assumption. with current data. We find that if we repeat the future analysis here without DUNE and T2HK, we obtain a 1% measurement of N e , which is roughly 20% stronger than the result shown in Refs. [13,88], and consistent with our statement in Section 4.2 that our JUNO simulation is roughly 20% stronger than the official collaboration projections of Ref. [88]. Figure 9 presents the results on the closures of different triangles t αβ and t kl . Each panel in this figure presents constraints on the real and imaginary part of t αβ (top row) or t kl (bottom row) at 95% credibility (dark) and 99% credibility (faint). We draw circles corresponding to constant values of the magnitude of |t αβ | 2 and |t kl | 2 as labeled, where each successive inward circle is an order of magnitude smaller. When constraining t eτ and t µτ , the expected future constraints are nearly degenerate with the current ones -constraints here are dominated by the sterile neutrino searches discussed in Section 3, specifically the NOMAD and CHORUS results discussed in Table 2. Constraints on t µe will improve modestly once information from DUNE and JUNO are incorporated. In contrast, measurements of the closures of the different pairs of columns will improve significantly with future data. Currently, each of these can be constrained |t kl | 2 10 −1 at 95% credibility. With future data, this will improve to roughly |t kl | 10 −2 for each of the three triangles. We summarize the current and future 3σ credibility upper limits on the triangle closures in Table 4. The analysis yielding Figs. 8 and 9 was conducted assuming the agnostic case of Section 2.2, whereby the matrix of which the LMM is a sub-matrix need not be unitary. The sub-matrix approach was taken in Ref. [14], where it was pointed out that assuming unitarity of the larger matrix leads to strong constraints from Cauchy-Schwarz inequalities. By remaining agnostic about the larger matrix, the improved measurement capability of future data is more apparent. An analysis assuming the larger matrix is unitary is contained in Section 6.2.

Secondary Results with Alternate Assumptions
As discussed throughout this work, different assumptions regarding the origin of unitarity violation, as well as which datasets are included in the analysis, can have significant impact on the resulting constraints on the unitarity of the LMM. The primary results of our work, where we analyzed all possible data under the agnostic case, were shown in Section 5. In this section, we explore two alternate assumptions. In Section 6.1, we repeat our analysis without including any short-baseline sterile neutrino searches (discussed in Section 3.7 and Table 2). In Section 6.2, we conduct an analysis in the sub-matrix case of Section 2.2, comparing the results with those obtained in the agnostic case presented above.

Impact of Short-Baseline Sterile Neutrino Searches
In the bulk of the analyses performed in our work, we have included results of short-baseline sterile neutrino searches, with results adapted from these sterile neutrino searches reinterpreted as limits on unitarity violation (see Table 2 for a summary of these results). To better understand how unitarity constraints rely on sterile neutrino searches, we repeat the analyses of the main text surrounding Fig. 9 without short-baseline results. Figure 10 shows the results. Here we note that the ranges on each of the panels in Fig. 10 measuring t αβ and t kl are much larger than the corresponding ranges in Fig. 9. However, it is apparent that in the absence of sterile searches, future data from IceCube, DUNE, JUNO, and T2HK would nevertheless allow us to understand the closure of all triangles of the LMM considerably better than current data allow. As in Fig. 9, we draw lines of constant |t αβ | 2 and |t kl | 2 = 10 −1 and 10 −2 in each panel, where the outer (inner) dashed line corresponds to a constant 10 −1 (10 −2 ) in these planes.   Tables 4 and 5, the improvement in the absence of sterile searches is much more dramatic, highlighting the importance of such experiments. We have not included any additional short-baseline searches in our future projections, as we do not expect any upcoming experiments to provide stronger sensitivity in the "averaged-out" regime [98,106,111,112] (as discussed in Section 3.7) than those summarized in Table 2.
Comparing the measurements of the individual matrix-elements-squared |U αk | 2 , as well as the row and column normalization conditions N α and N k , is difficult in this scenario. Short-baseline sterile neutrino searches, particularly the information from ν τ appearance that |t eτ | 2 and |t µτ | 2 are small, provide significant information on the elements |U τ k | 2 . Additionally, the constraint from MINOS/MINOS+ that N µ ≈ 1 is very important for determining the muon elements |U µk | 2 . If this information is discarded, every other probe of |U µk | 2 we consider is subject to a rescaling degeneracy. This is a direct result of the discussion of normalization effects throughout Section 2.1. Again, this highlights the importance of short-baseline sterile neutrino searches, such as MINOS/MINOS+ ν µ disappearance, for precise tests of leptonic unitarity. This analysis without short-baseline measurements results in a lower limit on N µ comparable to the one given in Table 3, however, N µ can be as large as ≈ 2.

Impact of the Sub-matrix Case Assumption
As discussed in Section 2.2, the sub-matrix case places certain Cauchy-Schwarz restrictions on the elements of U LMM , specifically requiring N α 1, N k 1, (6.1) For our main analyses in Sections 5.2 and 5.3, we assumed the agnostic case. In this subsection, we compare the results analyzed under the agnostic and sub-matrix hypotheses. First, we repeat the process that generates Fig. 7 for the current data analyzed, under the two different case assumptions. The results of this procedure are shown in Fig. 11.
We note several effects of the sub-matrix case in Fig. 11. First, the measurement of the electron row elements is largely unchanged -this is because the combination of KamLAND, solar neutrino measurements, and Daya Bay measure these elements very precisely regardless of the sub-matrix or agnostic assumptions. We discuss this further in Appendix C. We note that the measurement |U e1 | 2 results in a slightly lower preferred value under the sub-matrix hypothesis than the agnostic one -this arises from the Bayesian approach where we place a prior forbidding N e > 1 in the sub-matrix analysis. In looking at the muon row elements, we see marginal improvement in the measurement capability in the sub-matrix case than the agnostic one. This is again unsurprising, as restricting N µ 1 allows for improved measurements of these parameters. Finally, the largest difference between the two cases is in the tau row elements. Under the agnostic assumption, due |U τ 3 | 2 Figure 11. Measurements (and projections) of |U αk | 2 with current experimental data under the agnostic assumption (blue) and the sub-matrix case (green). Similar to Fig. 7, these are one-dimensional ∆χ 2 measurements of each parameter after marginalizing a 15 parameter fit down to each individual one. to the mild excess over expectation of OPERA ν τ appearance events, the current data prefer N τ > 1 (at very low significance). When we analyze the data under the sub-matrix assumption, this solution is forbidden, and the |U τ k | 2 elements are both preferred to be lower in magnitude, and we end up with smaller measurement regions.
In Section 5.3, we analyzed how current and future data can constrain the normalizations of the LMM rows and columns to be close to 1, asking how well we measure N α and N k . If we were to repeat this analysis under the sub-matrix assumption, the resulting figure analogous to Fig. 8 would look very similar, modulo the regions N α > 1 and N k > 1 being forbidden.
For completeness, Fig. 12 presents the projected future measurements of |U αk | 2 under the agnostic (red) and sub-matrix (purple) assumptions. The red lines are identical to those in Fig. 7. We see that the sub-matrix hypothesis improves constraints on |U µk | 2 somewhat, and |U τ k | 2 significantly.
Finally, we repeat the procedure that generated Fig. 9, which determined how well we can constrain the closure of the six different pairs of rows/columns focusing on current data only, we compare the results of this process when data are analyzed under the agnostic or sub-matrix case, in Fig. 13.
Again, we note several features of this result. First, when looking at the closure between two rows |t αβ | 2 (the top three panels of Fig. 13), we see that the resulting constraint on |t αβ | 2 is largely unchanged. This |U τ 3 | 2 Figure 12. Projections of future measurements of |U αk | 2 under the agnostic case assumption (red) and the submatrix case (purple). Similar to Figs. 7 and 11, these are one-dimensional ∆χ 2 measurements of each parameter after marginalizing a 15 parameter fit down to each individual one. Table 6. Summary of current and expected future constraints on the row closures |t αβ | and column closures |t kl |, under the sub-matrix assumption.

Current 3σ Upper Limit Future 3σ Upper Limit
is because both analyses include the searches for short-baseline neutrino appearance discussed in Section 3.7 which directly constrain the closures |t eµ | 2 (from ν µ → ν e appearance), |t eτ | 2 (from ν e → ν τ appearance), and |t µτ | 2 (from ν µ → ν τ appearance). The mild improvement seen in each of these panels comes from the Cauchy-Schwarz constraint |t αβ | 2 (1 − N α )(1 − N β ), where the normalization constraints assist in these planes. For the bottom panels, the closure of triangles formed between two columns |t kl | 2 , the difference between the agnostic and sub-matrix analyses is more drastic. Here, the Cauchy-Schwarz constraints are of the form |t kl | 2 (1 − N k )(1 − N l ), and because there are no direct experimental constraints on the closure, these Cauchy-Schwarz inequalities play a much more significant role. This feature was observed in Ref. [14], where they analyzed data under the sub-matrix case and noted that these inequalities place the strongest constraints on the closures between two columns. Likewise, Fig. 14 presents projected future constraints on |t αβ | 2 and |t kl | 2 under these two hypotheses, similar to Fig. 13. Moderate improvement on each parameter going from the agnostic to the sub-matrix cases is present in each panel here. Table 6 summarizes the results of these analyses -current and projected constraints on |t αβ | and |t kl | when operating under the sub-matrix hypothesis. Comparing Tables 4 and 6, we see the improved constraints that are obtained under the sub-matrix hypothesis compared to the agnostic case, coming from the Cauchy-Schwarz constraints and forbidding N α , N k > 1 in the fits. We do not reproduce an analogous version of Table 3 here -the same 3σ lower limits on N α and N k apply here, where the upper limits are now simply 1 from the theoretical constraints of the sub-matrix case.

Discussion & Conclusions
We have comprehensively analyzed how current and future neutrino oscillation data can be used to constrain non-unitarity of the LMM. Neutrino oscillation probabilities are functions of the LMM elements, the masssquared differences of the neutrino mass eigenstates, and the interaction potential of the neutrinos in matter. Measurements of the oscillation probabilities are therefore particularly useful for probing the structure of the LMM, as they are generally less sensitive to other new physics than non-oscillation probes (e.g. lepton universality, lepton-flavor violation, etc.). In our analysis we have characterized our ignorance of the LMM structure by separating our analysis according to three hypotheses: 1. The "standard" case, with only 3 neutrino flavors, and the LMM is the PMNS matrix, and therefore unitary.
2. The "sub-matrix" case, where there are n > 3 neutrino flavors, and the 3 × 3 LMM is non-unitary, but the greater n × n matrix is unitary.
3. The "agnostic" case where we assume nothing regarding the unitarity of a possibly m × n neutrino mixing matrix.
In the bulk of our analysis, we choose not to impose any a priori theoretical biases, and therefore compute constraints on non-unitarity in the agnostic case. However, we provide comparisons and sanity checks with the other two cases as necessary as validation of this approach.
When performing fits to current data, we include a representative sample of experiments that, when combined, provide the strongest set of constraints on the unitarity of the LMM possible. This set of data is as follows: solar neutrino experiments (SNO and Super-K), reactor antineutrino experiments (KamLAND and Daya Bay), long-baseline muon-neutrino disappearance (T2K and NOvA), long-baseline electron-neutrino appearance (T2K and NOvA), and long-baseline tau-neutrino appearance (OPERA). Where possible, up to being able to adapt experimental results to account for non-unitary mixing, we include the most recent results from each experiment. We also include results from short-baseline searches for anomalous neutrino disappearance/appearance from KARMEN, NOMAD, CHORUS, and MINOS/MINOS+.
Projecting to the next decade or so, we include the upcoming experiments that will qualitatively improve our understanding of leptonic unitarity. First is JUNO, where we simulate its reactor antineutrino capabilities. Second, we include T2HK's long-baseline beam-based searches for muon-neutrino disappearance and electronneutrino appearance. Next, we include near-future projections for IceCube's DeepCore and/or Upgrade measurements of tau-neutrino appearance. Finally, we consider all DUNE searches currently under studyits long-baseline beam-based searches for muon-neutrino disappearance, electron-neutrino appearance, and tau-neutrino appearance; as well as its solar neutrino capabilities.
Our main results are visually represented in Figs. 7, 8 and 9 for matrix elements-squared, column/row normalizations (Eqs. (2.4) and (2.6)) and column/row closures (Eqs. (2.5) and (2.7)) respectively, all in the agnostic case. The 3σ constraints on |U µ1 | 2 will improve by almost an order of magnitude, and further significant improvements are expected for |U e1 | 2 , |U e2 | 2 and |U µ2 | 2 . These are driven primarily by the improved precision of JUNO and DUNE over current experiments. The τ -row elements will be measured better by a factor of 2, owing to the expected sensitivity of DUNE and IceCube's τ -appearance search. Of note is that while the constraints on most substructures of the LMM will improve with future experiments, some do not. In particular, |U e3 | 2 has been extremely well-measured by the Daya Bay experiment, and the NOMAD constraints on |t eτ | 2 and |t µτ | 2 will not be improved on in the near future.
Our results highlight the improvements achievable by long-baseline oscillation experiments. However, they also emphasize that many constraints on the LMM are dominated by sterile neutrino searches capable of measuring triangle closure and row/column normalizations, in some cases extremely precisely. Indeed, the results of Fig. 8 demonstrate that much of the power in constraining N µ arises from the MINOS/MINOS+ search for ν µ disappearance. This is suggestive that dedicated sterile searches measuring ν e and ν τ disappearance precisely could be important for future improvements in our understanding of the LMM. A dedicated τ beam in particular would be extremely useful in understanding many of the poorly-constrained components of the LMM. An electron disappearance sterile search (specifically one free of large systematics that allows for stronger constraints in the "averaged-out" regime), meanwhile, would provide an alternative handle on the improvements expected from long-baseline experiments, and would therefore provide greater understanding of possible systematic effects in the two types of measurements.
Obtaining a detailed understanding of the LMM is critical to developing a theory of neutrino masses, as LMM unitarity or non-unitarity can be the result of non-intersecting subsets of neutrino mass mechanisms. Our results highlight the power of neutrino oscillation measurements to provide theoretically clean constraints on non-unitarity of the LMM, and therefore act as probes of these mechanisms. The neutrino puzzle is rapidly being untangled, as evidenced by the significant improvements expected between the current and future generation of neutrino oscillation experiments. Continuing improvement of both our theoretical and experimental understanding of neutrino oscillations, and therefore of the LMM, is critical to solving fundamental questions left unanswered by the Standard Model. the partial width for the decay must be expressed in terms of mass eigenstates: µ → eν i ν j , and the partial width is summed over these unmeasurable indices: When N µ and N e are not 1, the value of G F inferred from muon decay would be different from other processes in which it is measured. Measurements such as leptonic universality, the weak mixing angle, Z decays, and rare charged lepton decays place constraints on this type of scenario. We direct the reader to Refs. [30,34,36,52] for a thorough discussion of these types of constraints, as well as the model-dependent cases in which they apply.
We will discuss the model-dependence of such probes through the example of rare charged lepton decays. These processes would be forbidden in the SM if neutrinos were massless. The most constrained of these rare decays are loop-mediated processes such as µ → eγ, µ → e conversion and µ → 3e. 12 Unlike in cases where the neutrinos are in the final state, in these loop-mediated processes, neutrinos only appear within the loop. These processes are therefore in principle able to probe contributions from all contributing neutrino states, including weakly-coupled sterile states, as well as the structure of the full lepton mixing matrix U αk . However, as we explain below, the interpretation of such processes as constraints on unitarity are somewhat subtle and model-dependent.
The measurement of loop-mediated lepton flavor-violating decays is often expressed as a branching ratio between α → β γ and α → α νν. The branching ratio for this decay is [113][114][115][116] where F (x k ) is a loop function of the ratio x k ≡ (m ν k /m W ) 2 , whose relevant limits are Note that the matrix in Eq. (A.2) is not the 3 × 3 sub-matrix U αk , but the full matrix including all sterile states U αk , where k runs through 1, 2, ... n. From here we may proceed to interpret constraints from rare leptonic decays according to the three cases we discussed in Section 2.1.
In the standard case with three neutrinos and U ≡ U PMNS , by construction, any observation of charged lepton flavor violation (CLFV) should be interpreted as arising due to new physics that does not affect the unitarity of the LMM.
In the sub-matrix case, let us first examine the expected level of unitarity violation in a canonical type-I seesaw scenario, where there are three sterile neutrinos with masses much larger than the weak scale. Equation (A.2) can consquently be written as [117]: where X k = (m ν H,k /m W ) 2 and U αk is the unitary PMNS matrix. The newly introduced mixing angle is defined as tan 2θ k ∼ m D k /m M in terms of the Dirac masses m D k and Majorana mass scale m M of the canonical type-I seesaw mechanism. This manner of writing the branching ratio is dependent on the assumption that the Majorana mass matrix is diagonal with identical entries in the diagonal. A more complex ultraviolet structure will undoubtedly lead to variations on the estimates we present here, but should not affect the conclusion as long as the assumptions we make hold to a good approximation (e.g. small off-diagonal terms in the Majorana mass matrix). Since by construction, U αk is unitary, the leading non-zero terms arising due to the inclusion of heavy neutrinos are O(θ 4 k ) and O(x k θ 2 k ). Because θ 2 k ∼ m ν L /m ν H , and we assumed m ν H m W , both these terms lead to unobservably small CLFV.
One could interpret CLFV in the sub-matrix case without model-dependence, and therefore allow CLFV to impose a constraint on the LMM unitarity. This is often referred to as the Minimal Unitarity Violating approach, see e.g., Ref [34]. Remaining agnostic about the scale of neutrino masses, but imposing that all the sterile states are in the same limit compared to the weak scale, the leading contribution to Eq. (A.2) can be written as where c 1 = 4/3, 17/6, 10/3 depending on the mass scale of the sterile neutrinos, and c 2 = c 1 − 10/3. 13 In going from the first line to the second above, we make use of the fact that by virtue of studying the sterile case, the overall mixing matrix U must be unitary. This also results in the first and third terms in Eq. (A.5) being zero due to row closure of the full U matrix. Then, In this case, we can see that the non-observation of CLFV can indeed be viewed as a constraint on the closure of the α-β triangle, and therefore a test of unitarity. Note that our analysis of this scenario assumes that all steriles have masses in the same limit with respect to the weak scale. This would therefore not apply to e.g., models with one light and two heavy steriles. Finally, in the agnostic case, we can see that owing to the fact that the leading terms of F (x k ) are all constants, CLFV is a clear test of unitarity of the full U αk matrix. Below we quote the constraints on closure of the α − β rows in terms of the full matrix in the agnostic case. It is then straightforward to map this onto a constraint on t αβ when assuming Minimal Unitarity Violation for the sub-matrix case as above.
For the µ − e row, the strongest such constraint comes from the MEG collaboration, and leads to the requirement that | k U µk U * ek | 10 −5 [118]. Future measurements of µ → eγ will improve this limit by roughly one order of magnitude [119]. Planned searches for the related process µ → 3e, which are expected to have a similar sensitivity [120], while searches for µ → e conversion in nuclei can improve on this future constraint by a further order of magnitude [121].
For the τ − e and τ − µ rows, the strongest current constraints come from the BaBar collaboration and require that | k U τ k U * (e,µ)k | 10 −3 [122]. These measurements will be improved at B-factories, leading to a factor of ∼ 3 improvement in the τ − e row, and an order of magnitude improvement in the τ − µ row [123].

B Derivation of Vacuum Oscillation Probabilities Without Assuming Unitarity
In this appendix, we provide derivations for three-neutrino oscillations in vacuum where unitarity is not assumed. We analyze these oscillation probabilities in different distance and neutrino energy regimes that are appropriate for a variety of experiments. Experiments operating in particular length/energy limits are most sensitive to only certain matrix elements, as we will explain below.
In Section 2.3 we introduced the formalism for time-evolving a neutrino flavor eigenstate |ν α when its mixing is not unitary, beginning with Eq. (2.24). We can then project this onto a flavor eigenstate |ν β to determine the transition amplitude for ν α → ν β , which we express as l , E ν is the neutrino energy, and N α is defined in Eq. (2.6). The overall phase, exp(−im 2 1 L/2E ν ), does not enter oscillation probabilities so we will drop it henceforth. We separate our discussion into oscillation probabilities for disappearance/survival channels (α = β) and appearance channels (α = β).

B.1 Disappearance/Survival Probabilities
This can be rewritten as To obtain the oscillation probability P αα ≡ P (ν α → ν α ) ≡ |A αα | 2 , we square the transition amplitude. After rearranging terms, we obtain After minor rearrangements, this becomes oscillations associated with ∆m 2 21 ≈ 7.5 × 10 −5 eV 2 have yet to develop. Taking the limit ∆ 21 1 in Eq. (B.6), we obtain the oscillation probability In the limit of unitarity, this has a familiar form oscillations associated with ∆m 2 31 ≈ 2.5 × 10 −3 eV 2 will have averaged out over the energy resolution of an experiment. In this regime, the terms proportional to sin 2 (∆ 31 /2) will average out to 1/2. The term with sin(∆ 31 /2) cos(∆ 32 /2) will not, as one might naively expect, average out to 0, but Putting this into the complete oscillation probability in Eq.(B.6), we arrive at The ∆ 31 1 approximation is valid for the KamLAND experiment's measurement of P ee . The proposed JUNO experiment, which aims to measure the mass-ordering of neutrinos, will not operate in this regime. Instead, JUNO will operate in a regime between the two limiting cases discussed here, and as such we must adopt the full oscillation probability.

B.2 Appearance Probabilities
Now we discuss the appearance oscillation probabilities, where α = β. We define the phases associated with each quantity as ϕ αβ , ϕ αβ 2 , and ϕ αβ 3 with t αβ ≡ |t αβ |e iϕ αβ and U * αk U βk ≡ |U αk ||U βk |e iϕ αβ k . In analogy to Eqs. (B.3) and (B.4), we write the transition amplitude as Squaring this, we arrive at the oscillation probability, In the limit of unitarity, |t αβ | 2 → 0 and only the first two lines of Eq. (B.16) (sans the first term) remain. ∆ 21 1 Regime: Here the second term in Eq. (B.16) goes to zero, as do the second and third lines: This regime serves as a decent approximation for long-baseline electron-neutrino and tau-neutrino appearance oscillation probabilities (T2K, NOvA, DUNE, and OPERA). Note that unlike for the case of disappearance/survival probabilities, matter effects make a significant impact on these expressions, and we account for that in our simulations. ∆ 31 1 Regime: If the ∆m 2 31 -driven oscillations are averaged out, the oscillation probability becomes The appearance channel oscillation probability in the ∆ 31 1 regime has various terms capable of probing the second column of the LMM, and interference between terms from the third and second columns, as well as row triangle closure. Unfortunately, no experiments probing neutrino appearance in this distance/neutrino energy regime exist, or are planned.

C Measurements of the Electron and Muon Rows
In Sections 3 and 4, we discussed the sensitivity of each experiment to the various LMM parameters when unitarity is not assumed. Here, we show the results for the square of matrix elements in the electron and muon rows from a fit to certain subsets of the experiments we consider. Showing these subsets of experiments illustrates how combinations of different measurements affect our understanding of the various |U αk | 2 . This is performed for the electron and muon rows, where individual experimental measurements are only sensitive to those elements. We do not do this for the tau row, because there is not nearly as much experimental information for it, and the oscillation probability P τ τ has never been measured.
Electron Row Only: Figure 15 (left) displays the current knowledge of the electron row |U ek | 2 for a subset of the existing experiments we discussed in Section 3. Each panel displays two-dimensional projections of the test statistic ∆χ 2 for these CL, after marginalizing over the third, unseen parameter.  Figure 15. Left: Current measurements coming from experiments that are only sensitive to electron-row parameters |U e1 | 2 , |U e2 | 2 , and |U e3 | 2 . Here we compare measurements from KamLAND, Solar CC, and Daya Bay, and a joint fit region (black). Note that the measurements are combined prior to marginalizing over any of the three parameters, hence the joint measurement appears stronger than expected. Right: Future projections of these parameter measurements by JUNO and DUNE's solar neutrino capabilities, compared against the joint fit of current data (black).
We show results individually from KamLAND, SNO and Super-K measurements of solar neutrinos from CC interactions, and Daya Bay, in addition to a joint fit to these three sets of results. Note that due to not marginalizing over any parameters before performing the joint fit, the resulting two-dimensional ∆χ 2 contours do not follow the naïve expectation as a result of degeneracies in the parameter space. For example, Daya Bay measures the combination 4 |U e3 | 2 (|U e1 | 2 + |U e2 | 2 )/N 2 e , which is degenerate under |U e1 | 2 ↔ |U e2 | 2 , such that it appears that Daya Bay places no constraint in the upper-left panel of Fig. 15 (left), showing constraints in the |U e1 | 2 − |U e2 | 2 plane. However, Daya Bay constraints on |U e3 | 2 as a function of |U e1 | 2 , |U e2 | 2 combine with solar measurements in the |U e2 | 2 − |U e3 | 2 plane to bound |U e2 | 2 and |U e1 | 2 , such that the best-fit regions are the small black elliptical contours shown in the figure. Note that in this procedure the only constraint we impose is that we require each of the parameters satisfy 0 |U ek | 2 1. If one were to impose the constraint that the sum not exceed one, as would apply in the sub-matrix case we discussed in Section 2.1, the upperright triangular half of the |U e1 | 2 vs. |U e2 | 2 panel would be forbidden, somewhat limiting the KamLAND, solar CC, and joint fit contours. The stars in Fig. 15 represent the best-fit point in this parameter space given by this combination of datasets: |U e1 | 2 = 0.670, |U e2 | 2 = 0.306, and |U e3 | 2 = 0.022.
The upcoming JUNO experiment will measure a combination of these parameters as well, with the oscillation probability given by Eq. (B.6). JUNO will operate in the regime of L/E ν in which the effects of both mass-squared splittings are relevant, and therefore, with enough statistical power, can be independently  Here we compare measurements from T2K ν µ and ν µ disappearance, NOvA ν µ and ν µ disappearance, and the MINOS/MINOS+ sterile neutrino search, as well as a combined fit of these (black). Right: Comparison between the current, joint fit (black) and the future DUNE and T2HK ν µ and ν µ measurements of these parameters. sensitive to each mixing element |U ek | 2 . JUNO's capacity to measure each of these three parameters is depicted in Fig. 15, alongside the combined current measurement from Fig. 15 (left) in black. On its own, JUNO suffers the same degeneracy discussed above for Daya Bay under the interchange |U e1 | 2 ↔ |U e2 | 2 , and requires solar neutrino experiments to break the degeneracy. However, this interchange |U e1 | 2 ↔ |U e2 | 2 also requires changing the neutrino mass ordering (or the sign of ∆m 2 31 ). If the mass ordering can be determined independently of JUNO at high enough significance (for instance, by DUNE, T2HK) then the solution where |U e2 | 2 > |U e1 | 2 may be eliminated. The availability of redundant but independent data to select the right (|U e1 | 2 , |U e2 | 2 ) is a powerful tool to test new physics scenarios such as possible non-standard interactions of neutrinos.
In Fig. 15 (right), we do not perform a joint analysis of current and future data, and simply note that once future data from JUNO are included, the measurements of |U e1 | 2 and |U e2 | 2 will be dominated by JUNO, whereas the measurements of |U e3 | 2 will be dominated by current experiments, specifically Daya Bay. JUNO will measure |U e3 | 2 at slightly worse precision than Daya Bay. This measurement, since it is performed at a significantly different baseline from Daya Bay, serves as an important cross-check.
Muon Row Only: We perform a similar procedure focusing on the elements |U µk | 2 in Fig. 16. In the left panel of Fig. 16, we include the MINOS/MINOS+ sterile neutrino search, NOvA ν µ and ν µ disappearance, and T2K ν µ and ν µ disappearance. The combined fit of these data sets is shown in black. Instead of presenting these measurements in terms of the elements |U µ1 | 2 , |U µ2 | 2 , and |U µ3 | 2 , we present the measurements in terms of the combinations |U µ1 | 2 ±|U µ2 | 2 and |U µ3 | 2 . This is because the combination |U µ1 | 2 +|U µ2 | 2 is more precisely measured by MINOS/MINOS+, NOvA, and T2K, where their difference is not well-constrained by current data.
Again, parameter degeneracies cause the combined fit to appear stronger than naïve expectations. We see that |U µ3 | 2 is measured the most precisely of these three parameters, and the combination |U µ1 | 2 + |U µ2 | 2 is well measured by the combination of datasets in the bottom-left panel. The difference |U µ1 | 2 − |U µ2 | 2 is not well-measured, and is predominantly constrained by the requirement that both |U µ1 | 2 and |U µ2 | 2 are both between 0 and 1. This forces the difference |U µ1 | 2 − |U µ2 | 2 to be less in magnitude than the sum |U µ1 | 2 + |U µ2 | 2 . The stars in Fig. 16 (left) correspond to the best-fit points of these elements from the full analysis discussed in the main text: |U µ1 | 2 = 0.096, |U µ2 | 2 = 0.352, |U µ3 | 2 = 0.552. Figure 16 (right) displays our future projections on the muon row element measurements, where we compare the current joint fit (black) to projections of DUNE (green) and T2HK (red) ν µ and ν µ disappearance. 15 DUNE and T2HK measure oscillations over a wider range of L/E ν than their predecessors, thus they are sensitive to more than just the "dominant" term in the disappearance channel oscillation probability, namely the prefactor of sin 2 (∆ 31 /2) in Eq. (B.6). Indeed, these experiments have some sensitivity to the interference term, namely the final term of Eq. (B.6). This allows DUNE and T2HK to be sensitive to |U µ1 | 2 − |U µ2 | 2 unlike the current data considered in Fig. 16 (left). We see that these experiments can both demonstrate |U µ2 | 2 > |U µ1 | 2 at high significance.

D Bayesian Priors for Fifteen Parameter Fits
In many of our results, we use the Bayesian inference tool pyMultiNest [107][108][109][110] in order to analyze current and future measurements of the parameters of interest. We include all fifteen parameters in the following -nine matrix-elements-squared, four phases, and two mass-squared splittings. In this appendix, we explain the assumptions made in our Bayesian analysis and how they impact the applicable results.
We allow the mass-squared splittings to vary within the ranges with flat priors in those ranges on the two splittings. The current measurements pull the posterior likelihood of the fit to the current best-fit values of ∆m 2 21 ≈ 7.5 × 10 −5 eV 2 and ∆m 2 32 ≈ 2.45 × 10 −3 eV 2 . For the matrix-elements-squared |U αk | 2 , we require that they lie in the range [0, 1]. We include flat priors on each of these elements-squared between 0 and 1. For the electron-and muon-row elements |U ek | 2 and |U µk | 2 , we find that the current data are powerful enough that the posterior likelihood is largely independent of the choice of prior (unless anything atypical is adopted), whereas the tau-row elements |U τ k | 2 still do not have strong information from the OPERA experiment and are mildly sensitive to the choice of prior. Once future data are included, this point becomes unimportant.
Lastly, the phases {φ e2 , φ e3 , φ τ 2 , φ τ 3 } are defined to span [0, 2π] and the prior included is flat in terms of the phase. The relative sizes of some of the matrix element magnitude combinations require the phases to take on certain values. For instance, the phase φ e2 is quite constrained by the closure of the e-µ triangle, so it is insensitive to the choice of the prior. On the other hand, since any current constraint on φ e3 is relatively weak, its posterior likelihood can be sensitive to the choice of prior. As with the case of the matrix-elements-squared above, this issue is unimportant once future data are included in the analysis.

E Measurement of LMM Phases
In the main text (Section 5.2), we presented the results of our analysis for the current and future measurements of the elements-squared |U αk | 2 , which are parameterization-independent quantities. Here, we discuss the current and future measurements of the phases {φ e2 , φ e3 , φ τ 2 , φ τ 3 }, the four phases we consider in our analyses. These measurements are parameterization-dependent and only apply in the MP parameterization we employ. . Current and future measurements of the four phases φ e2 , φ e3 , φ τ 2 , and φ τ 3 in the MP parameterization. We show the measurements given current data under the agnostic case assumption (blue) and the sub-matrix case assumption (green), compared with future data analyzed under the agnostic case (red) and the sub-matrix case (purple). All contours shown are 95% (dark colors) and 99% (fainter colors) credibility. Figure 17 presents each set of two-dimensional measurements at 95% and 99% credibility of these four phases. We show this for four sets of data/assumptions -all current data are analyzed in blue and green, where the blue (green) contours correspond to data analyzed under the agnostic (sub-matrix) case. Future data are analyzed under the agnostic (sub-matrix) case in red (purple).
We see in Fig. 17 that, even with current data, φ e2 and φ τ 2 are constrained to be close to π, and φ τ 3 must be close to π. These requirements come from the current constraints on the closures between different rows |t αβ | 2 and the relative sizes of different products of magnitudes of elements. For instance, the relative size of the legs of the e-µ triangle are |U e1 ||U µ1 | ≈ 0.25, |U e2 ||U µ2 | ≈ 0.33, and |U e3 ||U µ3 | ≈ 0.11. This requires φ e2 to be near π if the triangle is to close. Because the product |U e3 ||U µ3 | is relatively small, the constraints on φ e3 from the closure of this triangle are relatively weak. As precision experiments that are sensitive to CP-violation begin collecting data, all of the phases will be measured more precisely -with DUNE and JUNO data included, φ e3 will be measured to much higher precision. We also see that the assumptions of sub-matrix or agnostic do not impact the measurement precisions of the phases, with current or future data.
In Section 3.7 we discussed how various short-baseline searches for anomalous neutrino disappearance/appearance are used to constrain the unitarity of the 3 × 3 LMM. We briefly commented on the LSND and MiniBooNE anomalies before ultimately not including them in our more complete analysis, due to tensions between their anomalous ν µ → ν e appearance results with searches for ν µ and ν e disappearance.
If we assume that the anomalous appearance of ν µ → ν e in LSND and MiniBooNE is instead due to non-unitarity, we must use the "averaged-out" region of their sterile neutrino parameter space to draw connections. Their combined result points to a preference for nonzero |t eµ | 2 ≈ 2.6×10 −3 , excluding |t eµ | 2 = 0 at roughly 6σ confidence. Under the sub-matrix assumption, this closure is constrained using