Leptonic Unitarity Triangles

We present a comprehensive analysis of leptonic unitarity triangles, using both current neutrino oscillation data and projections of next-generation oscillation measurements. Future experiments, sensitive to the degree of CP violation in the lepton sector, will enable the construction of precise triangles. We show how unitarity violation could manifest in the triangles and discuss how they serve as unitarity tests. We also propose the use of Jarlskog factors as a complementary means of probing unitarity. This analysis highlights the importance of testing the unitarity of the leptonic mixing matrix, an understanding of which is crucial for deciphering the nature of the neutrino sector.

We present a comprehensive analysis of leptonic unitarity triangles, using both current neutrino oscillation data and projections of next-generation oscillation measurements. Future experiments, sensitive to the degree of CP violation in the lepton sector, will enable the construction of precise triangles. We show how unitarity violation could manifest in the triangles and discuss how they serve as unitarity tests. We also propose the use of Jarlskog factors as a complementary means of probing unitarity. This analysis highlights the importance of testing the unitarity of the leptonic mixing matrix, an understanding of which is crucial for deciphering the nature of the neutrino sector.
Introduction. -The discovery of neutrino oscillations confirmed that lepton flavor and mass eigenstates are distinct entities. The mixing of the known lepton flavor and mass eigenstates is canonically parameterized by the 3 × 3 unitary Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [1,2], in direct analogy to the Cabibbo-Kobayashi-Maskawa (CKM) matrix [3,4] in the quark sector. A crucial difference between the leptonic and quark sectors, however, is the origin of such mixing, since the Standard Model as formulated does not allow for a renormalizable neutrino-Higgs Yukawa interaction. Misalignment of the flavor and mass eigenstates in the lepton sector, i.e., the origin of the PMNS matrix, remains an open question. Predictions of neutrino masses, for example those involving right-handed neutrinos, often lead to a non-unitary 3 × 3 leptonic mixing matrix (LMM), which is a sub-matrix of a larger unitary matrix. Searches for deviation from unitarity of the LMM are directly probing our fundamental understanding of neutrino masses. Throughout this work, we will refer to a general 3 × 3 LMM as U LMM and one that is assumed to be unitary as U PMNS .
The matrix U PMNS is parameterized by three mixing angles, θ 12 , θ 13 , θ 23 , and a phase δ CP . 1 The degree to which the combination of charge C and parity P symmetry, CP, is violated is proportional to the Jarlskog invariant [5], J PMNS ≡ c 12 c 2 13 c 23 s 12 s 13 s 23 sin δ CP , where s ij = sin θ ij and c ij = cos θ ij . This quantity is relevant to understanding the origin of the baryon asymmetry of the universe, as CP-violation is one of the required conditions for such an asymmetry to exist [6]. Furthering our understanding of the PMNS matrix is therefore important for probing fundamental questions.
In the quark sector, many experimental tests of the unitarity of the mixing matrix have been carried out, e.g., meson mixing and decays [7]. Their results are often visualized using unitarity triangles [8][9][10][11][12][13][14][15]. In such triangles, many measurements would meet at a point in the complex plane if the matrix is unitary. The areas of the triangles are proportional to the Jarlskog invariant of the CKM matrix.
In this work, we present a comprehensive analysis of leptonic unitarity triangles using current neutrino oscillation data and projections of future experiments. Our main results are shown in Fig. 1, where we present six unitarity triangles, and in Fig. 3, where we show nine Jarlskog factors and the PMNS Jarlskog invariant. This set of triangles and Jarlskog factors allows for a complete understanding of the LMM in a compact form. In our companion paper [16], we discuss how well current and next-generation experiments are testing all unitarity conditions of U LMM . Our results show the importance of separately analyzing appearance and disappearance data, and demonstrate the power of future oscillation measurements to constrain fundamental physics related to the neutrino sector.
The Leptonic Mixing Matrix and Unitarity Triangles. -The LMM describing the mixing of leptonic flavor and mass eigenstates, is defined by its appearance in the charged current interaction. In full generality, this interaction can be written as where U αi is an M × N matrix coupling M charged leptons to N neutral leptons, of which the LMM would then be a non-unitary 3 × 3 sub-matrix. If only 3 generations of leptons exist, then U αi ≡ U LMM ≡ U PMNS .
This mixing of flavor and mass eigenstates leads to neutrino oscillation, the precise measurement of which is the focus of current and future detection efforts. The primary measurements in oscillation experiments are the so-called appearance and disappearance probabilities, P αβ and P αα , where α = β and α, β ∈ e, µ, τ . The general transition probability for any α, β is where ∆ ji = ∆m 2 ji L/2E ν (with ∆m 2 ji ≡ m 2 j − m 2 i , and L and E ν are the neutrino's propagation distance and energy, respectively). Appearance and disappearance probabilities thus directly measure a combination of PMNS matrix elements and mass-squared differences.
For a 3 × 3 matrix, a set of six unitarity triangles may be defined as arising from the conditions U † U = U U † = I. These triangles correspond to the closure of the products of cross-terms of a unitary matrix: The elements U αi are in general complex, so the above conditions can be shown as closed triangles in complex planes. The unitarity triangle is constructed by dividing out one of the three terms in the sums of Eq. (5), and defining a vertex of the triangle as ρ xy + iη xy , (x = y), such that a closed triangle is defined as having vertices at the origin, (1, 0), and (ρ xy , η xy ). There is ambiguity in the choice of the denominator, and hence (ρ xy , η xy ) for a given row/column. We explain our choices, made in an attempt to cover measurements of all parameters of U PMNS , and show the definitions of (ρ xy , η xy ) in Appendix A. Data Analysis and Methodology. -We take global neutrino data and recast their joint measurements onto leptonic triangles. For current measurements, our goal is not to do the most precise, comprehensive global fit on mixing parameters, so we interpret the majority of experimental results as rate-only measurements. More concretely, this means we assume a given experiment measures an oscillation probability P αβ at a fixed L and E ν . We find that this reasonably reproduces the reported experimental results, and apply it to T2K [23] and NOvA [21] (which are sensitive to δ CP ), as well as Daya Bay [18], SNO [19,20], and OPERA [22]. For KamLAND, we include a more detailed analysis that utilizes the measured energy spectrum of neutrinos [17]. We do not include other solar neutrino observations [26][27][28][29][30] nor any atmospheric measurements [31][32][33] as they do not contribute significantly to our results. Appendix C provides further details.
Using a given combination of data sets, we construct a likelihood function depending on a set of oscillation parameters. Our fits in Figs. 1 and 2 depend on six parameters: sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 , δ CP , ∆m 2 21 , ∆m 2 31 . We include Gaussian priors on the mass-squared splittings from respective experimental results when analyzing current data. Appendix C explains in more detail. We use the Bayesian inference tool pyMultinest [58,59] to construct credible regions in this parameter space. The posterior distributions are then projected onto any combination of parameters of interest.
Results. - Figure 1 shows the 95% and 99% credible regions of the six unitarity triangles with all current data (green contours), and with the addition of JUNO and DUNE projections (red). Note that these results inherently assume that the LMM is unitary and follow the PMNS parameterization. We discuss the implication of assuming unitarity a priori, and how the results would differ without this assumption in Appendix B. Figure 1 also shows the projections of future sensitivity to neutrino disappearance (P αα , blue contours) and appearance (P αβ , orange) probabilities separately. This distinction between the two measurements is crucial, both for understanding how the best fit regions arise and for determining how non-unitarity manifests itself.
Let us discuss the e-µ plane ( Fig. 1 top-left), as it is amenable to dissection. This triangle can be expressed in terms of PMNS parameters ρ eµ = c 2 12 +cos δ CP s 12 c 12 c 23 s 13 s 23 , η eµ = sin δ CP s 12 c 12 c 23 s 13 s 23 .
(7) Precise measurements of the disappearance channels P ee and P µµ allow for the determination of all three mixing angles: θ 12 , θ 13 , and θ 23 , which determine the center and the radius of a circle. Disappearance measurements of P µµ in long-baseline contexts, e.g., DUNE, are only sensitive to sin 2 (2θ 23 ), leading to an octant degeneracy. This produces ambiguity in the radii of the circles in the e-µ and e-τ planes, and results in the genus-two structure in the 1-2 plane. Appearance channel measurements, such as P µe for DUNE, can determine the value of δ CP , and therefore determine a preferred direction in the (ρ, η) plane. The definitions of (ρ, η) can significantly impact the appearance of these triangles, but the general features of disappearance measurements giving ring-type structures and appearance selecting out a direction are universal (For comparisons to other choices, see Refs. [45][46][47] and Appendix A). Tau appearance measurements, present or projected, are insufficiently precise to provide further discriminatory power. In addition, no planned measurement of P µτ is actually sensitive to δ CP , even with improved precision [36,37]. In contrast with the CKM fits, where multiple observables can independently measure the CP-violating phase, in the neutrino sector there is only one present or near-future observable sensitive to δ CP , namely ν e appearance.
One might wonder how unitarity violation would appear in these leptonic triangles. To test this, we simulate data with injected non-unitarity, and analyze it assuming the LMM is unitary. Note that in order to include effects of non-unitarity, we must adopt a parameterization beyond PMNS, which requires 13 parameters, as described in detail in Appendix B. We inject non-unitarity by making i U ei U * µi = 0.01+0.02i and α U α2 U * α3 = −0.004+0.006i, allowed by current data. Figure 2 shows the constructed triangles in the e-µ and 2-3 planes, where the fit (incorrectly) assumes the LMM is unitary. With a joint disappearance and appearance analysis, there is no indication of unitarity violation as the joint contours appear similar in shape and size to those in Fig. 1. However, individual channel measurements reveal tension. The appearance-only sensitivity slightly improves with visible shifts in the central values. Most notably, disappearance and appearance measurements disagree at nearly 2σ in the 2-3 plane. Importantly, no tension is present (even at 1σ) when viewed in terms of measurements of the PMNS parameters. This demonstrates that such analyses of disappearance and appearance measurements can be a powerful probe of unitarity violation in the lepton sector, and complementary to sterile neutrino searches [63][64][65][66][67] described in Appendix C.
Testing whether certain subsets of data agree on these planes is one general test of unitarity. However, constructing the triangles assumes the LMM is unitary. To numerically constrain non-unitarity, one needs to discard such an assumption and include sterile neutrino searches, which contribute significantly. This leads to non-intuitive connections between unitarity triangles and numerical results, e.g., the unitarity violation in e-µ plane for Fig. 2 will be excluded by future searches at ∼2σ, contrary to what the figure suggests, when sterile search results are included. In what follows, we develop another intuitive means of using oscillation measurements to probe unitarity in the lepton sector, which accounts for information loss when triangles are constructed under the assumption of unitarity. Specifically, we will consider the consistency between measurements of the degree of CP violation in the LMM.
Jarlskog Factors and the Jarlskog Invariant. -For a unitary LMM, the Jarlskog invariant (Eq. (1)) is a measure of CP violation and is related to the area of the unitarity triangles in Fig. 1. The area of triangles as defined in Eq. (5) would be J PMNS /2. However, by constructing the unitarity triangles as in Eq. (6), their area becomes (J PMNS /2)/(|U αi | 2 |U βj | 2 ). As long as δ CP is not 0 or π, J PMNS , and therefore the triangle areas, are non-zero. In general, J αi can be calculated by taking the cofactor of the (α, i) element of the LMM, modulo taking the complex conjugate of two elements. If the LMM is not unitary, there is no guarantee that these nine different cofactors are the same (Ref. [68] demonstrated this potential in a 3 + 1 sterile neutrino scenario explicitly). We define these Jarlskog Factors J αi , related to each cofactor, using such that the nine possible |J αi | are the same and equal to |J PMNS | when the matrix is unitary. This condition is necessary, but not sufficient, for U LMM to be unitary. Without assuming unitarity when constructing Fig. 1, choosing six triangles can provide information related to at most six different J αi (see Appendix A for further details). Therefore, to obtain a full characterization of the LMM and its potential non-unitarity, six unitarity triangles do not suffice. Indeed, one option is to construct nine unitarity triangles corresponding to the nine Jarlskog factors. However, as there are only six closure relations (Eq. 5), this would lead to some redundancy. Therefore, a compact manner of representing all six closure relations, as well as characterizing the full LMM is to to show six unitarity triangles assuming unitarity, and nine Jarlskogs without assuming unitarity.
To compare the Jarlskog invariant and Jarlskog factors, we perform a fit to all current and current plus future data while (not) assuming unitarity of the LMM to construct J PMNS (J αi ). When not assuming unitarity, we adopt the parameterization explained in Appendix B. Figure 3 shows the result. Our current measurement of J PMNS (bottom row, purple) is consistent with the results of other more detailed fits [60]. Each independent J αi measurement agrees, consistent with the expectation that the LMM is unitary. Furthermore, the data prefer J αi = 0 at nearly 1σ, indicating that CP is violated in the lepton sector. If unitarity is assumed, J PMNS = 0 is disfavored at roughly 2σ. The inclusion of future JUNO and DUNE data strengthens these claims, disfavoring J PMNS = 0, and therefore CPpreservation, at over 5σ. Certain Jarlskog Factors, particularly those involving knowledge of |U τ i |, will remain difficult to measure when unitarity is not assumed.
Conclusions. -We have presented a comprehensive analysis of leptonic unitarity triangles using all current and future neutrino oscillation experiment data. Figure 1 displays how the closure of six unitarity triangles is currently constrained, and how these constraints will improve with future experiments. By virtue of the nature of these measure-ments, it is noteworthy that, unlike for the CKM matrix, it is difficult to construct intersections of many measurements of PMNS matrix parameters. Non-unitarity can nevertheless explicitly manifest itself as shown in Fig. 2, though observation of non-unitarity requires distinguishing between appearance and disappearance data sets when constructing unitarity triangles. Figure 3 presents an alternative and complementary visualization of constraints on the unitarity of the LMM in terms of Jarlskog factors. The allowed ranges of the J αi are consistent with each other and with non-zero CP violation in the lepton sector. If the leptonic mixing matrix was not unitary, these measurements would yield different J αi .
The Standard Model demands new physics to explain the origin of neutrino masses and therefore oscillations. The impending era of precision experiments will enable us to further understand the structure of the leptonic mixing matrix, and constraints of the matrix's unitarity serve as a probe of the mechanism of neutrino masses. Meanwhile, the origin of the baryon asymmetry of the universe will be better understood through studies of the degree of CP violation in the lepton sector. Performing detailed studies of the PMNS unitarity triangles therefore bears directly on both of these problems in the Standard Model.

Supplemental Material
In these appendices we provide additional details on the analyses described in the main text. In Appendix A we explain our choices for the set of ρ xy and η xy that we use and express them in terms of PMNS elements. Appendix B explains the parameterization we adopt when we do not assume the leptonic mixing matrix (LMM) is unitary. We show how the injected unitarity violation in the analysis surrounding Fig. 2 is obtained, and we also show how measurements of the unitarity triangle parameters are modified when unitarity is not assumed. Finally, Appendices C and D list, in detail, the current and future experiments that we consider in our analyses, respectively.
From this, one can write down six real equations, requiring the normalization of rows and columns of matrix-elementssquared |U αi | 2 to be 1: One can also write six complex equations corresponding to "closures" between two different rows (α and β) or two different columns (i and j): From these closure relations, one can construct the familiar unitarity triangles in the (ρ, η) plane by dividing each term in the closure by one of the sides. For a given row/column, there are 3 different triangles one could define that are not related to one another by a simple inversion. For the rows, the normalized triangle definitions can be chosen as: • e − µ triangle: • e − τ triangle: • µ − τ triangle: Recalling that Jarlskog factors are related to the areas of the un-normalized triangles, we can naturally relate the normalized triangle definitions to the Jarlskog factors as well. In the cases above, we obtain the following relations: For the columns, the triangle definitions can be chosen as • 1 − 3 triangle: • 2 − 3 triangle: Repeating the same analysis for the column triangles as we performed for the rows, we can derive the following relations between the triangles and Jarlskog factors: The information contained in the triangles defined in Eqs. (A5-A7), (A11-A13) is duplicated by taking the reciprocal triangles (we abuse notation and refer to the reciprocal triangle T αβ ) −1 for k = 1, 2, 3), so one need only consider one set. If we do not assume unitarity when performing the triangle analysis, we can see that 9 unitarity triangles must be constructed in order to fully measure the elements and phases of the LMM. Looking only at the imaginary parts of the triangles, this corresponds to a set of 9 triangles measuring all 9 Jarlskog factors once in combination with all matrix element norms twice.
If we assume unitarity when constructing triangles, we would find that all J-factors are identical by definition, and equal to the Jarlskog invariant in the PMNS parameterization. Thus constructing triangles allows us to measure J PMNS and products of norms. In this case, to display all the information contained in the 18 possible triangles, we must pick 6 triangles to cover all 9 norms with minimal redundancy. Any set of 6 triangles which includes 9 separate norms actually contains 12 norms, such that there are always three norms which are measured twice. An example of a set of triangles which would encapsulate all possible information would be T ⊃ T eτ . With this set, |U ei |, i = 1, 2, 3 would be repeated twice. A more "flavor-democratic" approach to choosing a set of six triangles, and the one we use in the main text, is to choose one triangle from each row and column: In this way, we measure all information in the LMM matrix under the assumption of unitarity, while only repeating measurements of |U e3 |, |U µ2 | and |U τ 1 |. This specific set of choices is further motivated by the discussion around Fig. 2 of how to use unitarity triangles to observe unitarity violation. Given the degree of non-unitarity allowed by current measurements, which is then used to construct Fig. 2, it was determined that the above set of 6 triangles was most instructive for observing tension between the appearance and disappearance data in the various (ρ, η) planes.
Other choices of triangles have been made previously in Refs. [45][46][47]. For the 1-3 triangle, their choice of (ρ, η) is This choice corresponds to T (−1) 13 as defined above, and is therefore measuring J τ 2 . Note that this definition is the leptonic equivalent of the d − b CKM triangle that is commonly shown. When we adopt this definition, our joint-fit region from current experiments is consistent with the result of Refs. [46,47], which can be broken down to a disappearance circle that is centered at (0, 0) and an appearance region that is more visibly radially oriented.
In order to fully characterize a potentially non-unitary LMM in an economical and intuitive fashion, it is clear that we must account for the fact that assuming unitarity sets all J-factors equal in the triangles. Thus we must separately measure all nine possible J-factors. By showing an appropriate set of 6 unitarity triangles when assuming unitarity, 9 J-factors computed without assuming unitarity, and J PMNS , we graphically represent all possible information in the LMM. As discussed in the main text, there is an added benefit to computing J-factors separately, as they include information obtained from sterile searches that is not otherwise visible in the triangle planes.
The chosen set of (ρ, η), under the assumption that the LMM is unitary, can be expressed as: 3 ρ eµ = c 2 12 + cos δ CP where the 13 free parameters (ignoring the potentially physical Majorana phases) are necessary to describe a 3 × 3 complex matrix after accounting for rephasing of the charged lepton fields. Alternative parameterizations are commonly adopted in the literature [69][70][71], all with the same number of free parameters. We motivate the use of our parameterization, and discuss maps between this and others in the literature in Ref. [16]. While the complex phases in U LMM appear on different matrix elements than in U PMNS , the two parameterizations are related (if U LMM satisfies the unitarity conditions) by rephasing of the neutrino fields. For any set of PMNS parameters, an equivalent set of LMM parameters may be determined uniquely. We obtain best-fit values for the 13 LMM parameters by making use of the observation that the LMM fit must match the global fit for the PMNS parameters when the LMM is unitary. Analyzing current data when assuming unitarity yields the maximum-likelihood values of sin 2 θ 12 = 0.316, sin 2 θ 13 = 0.0220, sin 2 θ 23 = 0.565, and δ CP = 223 • . These four best-fit values may then be used in conjunction with the 9 constraints applicable to a 3 × 3 unitary matrix to solve for the 13 LMM elements, yielding When we study the nine separate Jarlskog factors J αi cf. Fig. 3, we use this set of 13 parameters, projecting down to relevant combinations for the allowed regions of different J αi . We also use this parameterization to simulate future data for JUNO and DUNE assuming U LMM is not unitary, i.e., U † U = I. In generating the data that are analyzed for Fig. 2, we modify the input values of φ e2 and φ e3 from those in Eq. (B2) by ∆φ e2 = −3.44 • and ∆φ e3 = 1.72 • . This preserves the unitarity constraint that the individual rows and columns of U LMM are properly normalized, i |U αi | 2 = 1, α |U αi | 2 = 1, but causes non-closure of the triangles in several planes.
One non-closure is in the e-µ plane, with i U ei U * µi = 0.01 + 0.02i. Sterile neutrino searches that look for zero-distance neutrino oscillation (as described in Appendix C) are sensitive to the absolute value squared of the non-closure, and this level is at the upper end of what is currently allowed by data. In addition, there is also non-closure in the 2-3 plane, with α U α2 U * α3 = −0.004 + 0.006i. This is shown in Fig. 2. In the results shown in Fig. 1 and Fig. 2, fits to existing/future data were performed with the PMNS mixing angles as free parameters (sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 , δ CP , and mass-squared splittings), such that unitarity was explicitly assumed. Confidence level contours were then constructed by mapping these parameters onto those that enter each unitarity triangle, using the PMNS parameterization. Here we perform a fit using the parameterization described above, where the unitarity of U LMM is not guaranteed. The difference between these two fits is most apparent in triangles that depend on mixing matrix elements that are not powerfully measured on their own in experiments, but can be inferred by other measurements if unitarity is assumed. Specifically, that is the case for the mixing matrix elements U τ i . For instance, in the PMNS parameterization, U τ 3 = cos θ 23 cos θ 13 , which can be constrained fairly well by atmospheric and reactor experiments. Without unitarity, the strongest measurement of U τ 3 comes from OPERA's ν µ → ν τ appearance, a significantly less precise measurement.
We show this difference in Fig. A1, focusing on two different triangle planes, (ρ eµ , η eµ ), where the differences are small, and (ρ 23 , η 23 ), where the differences are the largest. All contours shown are 95% (dark contours) and 99% confidence level (light). Here, the purple and green regions correspond to current data fits, where the purple (green) region is the result of the fit without (with) assuming unitarity. Likewise, light blue (unitarity not assumed) and red (unitarity assumed) regions are from fits including current and future data. The green and red regions here correspond with those of the same color shown in the appropriate panels of Fig. 1. In each of the two panels, the filled-in star denotes the best-fit point in this parameter space obtained by analyzing all current data when unitarity is assumed (green dataset), where the open star indicates the best fit point when unitarity is not assumed (purple dataset).
Most notable here is the difference in the size of allowed regions between when unitarity is or is not assumed for (ρ 23 , η 23 ). As mentioned above, this is largely due to the uncertainty regarding the magnitude of the elements |U τ 2 | and |U τ 3 |. We also see that the current data prefer a much larger triangle in this plane if unitarity is not assumed -this is due to the preference for large |U τ 3 | from the OPERA experiment [21]. Less dramatic is the difference in the left panel of Fig. A1, where an additional region of parameter space is viable in the future analysis when unitarity is not assumed (pale blue contours) compared to when it is (red ones). This comes from not being able to precisely determine whether |U µ3 | 2 is less than or greater than 0.5 when unitarity is not assumed, analogous to the octant degeneracy for θ 23 .

Appendix C. Current Experimental Results Included
Here we list the current experimental results that we include in our data analysis, and specify to which parameters each experiment is most sensitive. We also show the results of our data analysis (when unitarity is assumed) of all of the current data included in measuring the combination of the parameters δ CP and sin 2 θ 23 , validating our approach.
Mass-Squared Splittings When mentioned in the following list, we allow the mass-squared splittings ∆m 2 21 and ∆m 2 31 to vary independently, allowing the possibility of both mass orderings. Based on our methods of incorporating existing measurements, KamLAND, Daya Bay, T2K, NOvA, and OPERA are sensitive to mass-squared splittings. For each experiment, we include a Gaussian prior on the relevant mass-squared splitting from the experimentally reported 1σ range. The resulting fit region of the two mass-squared-splittings that we obtain is consistent with those from more sophisticated global fits [60][61][62]. While tensions exist between different measurements of ∆m 2 21 , we find that the analyses we perform do not change whether we also include a measurement of ∆m 2 21 from Super-Kamiokande or not [30]. Normalization Many of the current (and future) experiments we consider infer a neutrino oscillation probability by measuring a far-detector-to-near-detector ratio, i.e., they measure the neutrino flux times cross sections of one flavor α at the near detector and one flavor β at the far detector. The ratio of these two, up to cross section and flux effects, gives the oscillation probability P αβ . If the LMM is not unitary, however, this is not completely true -zero distance effects lead to P αα not being 1 at the near detector, but (C1) When performing an analysis that does not assume unitarity, like those surrounding Figs. 3 and A1, the inferred oscillation probability of an experiment with a near detector must be normalized by the factor in Eq. (C1). This normalization factor has the same impact on an analysis as including an overall systematic normalization uncertainty in an analysis, so as long as the normalization uncertainty of a given experiment is larger than the uncertainty on the quantity in Eq. (C1), these effects are unimportant. This is the case for all experiments considered here with a muon-neutrino source (T2K, NOvA, OPERA, DUNE), so we neglect this effect for those experiments. As described below, we include this effect for Daya Bay but the results are largely insensitive to it. KamLAND The reactor antineutrino experiment KamLAND measures the oscillation probability P (ν e → ν e ) over a wide range of baseline length L divided by neutrino energy E ν . Appendix B of Ref. [17] provides measurements of this oscillation probability for different values of L/E ν . We use these measurements, which take into account matter effects, as well as either the standard, three-neutrino oscillation probability or a modified one to account for non-unitary mixing, to place constraints on mostly sin 2 θ 12 . If the matrix is not unitary, KamLAND is mostly sensitive to the product |U e1 | 2 |U e2 | 2 . A more recent analysis is Ref. [25], but it does not contain enough information for us to reasonably reproduce its results using this approach.
SNO is also sensitive to neutral current scattering, which observes the effective oscillation probability P NC = i |U ei | 2 prod. |U ei | 2 det. . This becomes This measurement is limited by theoretical uncertainties associated with the Standard Solar Model [20], so we conservatively assume that SNO measures it at the 25% level.
T2K For the electron-neutrino appearance channels P (ν µ → ν e ) and P (ν µ → ν e ) measured at T2K, we assume that the experiment measures this probability for a fixed energy E T2K = 600 MeV (the mean energy of the J-PARC beam) at a distance of L = 295 km. We also assume a constant matter density of ρ = 2.6 g/cm 3 along the path of propagation. While this approach is an oversimplification and does not include systematic uncertainties from T2K, we find that it reproduces the results of Ref. [23] well. We use the predicted signal and background event rates from Ref. [23] for the ν e chargedcurrent quasielastic (CCQE), single-pion, and ν e CCQE samples. We allow ∆m 2 31 to vary (including the prior mentioned above) in this portion of the analysis. T2K also measures ν µ and ν µ disappearance. We interpret this measurement as information on the quantity 4|U µ3 | 2 (|U µ1 | 2 + |U µ2 | 2 ) to agree with the results of Ref. [23]. Matter effects are much smaller in this channel, so we ignore them. If U LMM is unitary, this translates effectively into a measurement of sin 2 (2θ 23 ). We assume T2K measures 4|U µ3 | 2 (|U µ1 | 2 + |U µ2 | 2 ) = 1.00 ± 0.03. We find that including disappearance information in this way reproduces the measurement capability of the experiment from Refs. [23,24] better than assuming a fixed length, fixed energy measurement in this channel. More recent results exist in Ref. [24], however, we use Ref. [23] as it allows us to extract signal and background rates for the appearance channels for any combination of oscillation parameters.
NOvA Our methodology for NOvA is very similar to that of T2K: we assume a fixed energy for the electron-neutrino appearance measurements of E NOvA = 1.9 GeV and L = 810 km (as well as a constant matter density of 2.84 g/cm 3 ). Expected signal event rates from Ref. [21] allow us to construct a log-likelihood as we vary oscillation parameters. Like with T2K, we allow ∆m 2 31 to vary within its prior for the NOvA measurements. Since NOvA prefers a value of sin 2 θ 23 slightly away from maximal, we include its disappearance channel measurements in our fit by assuming it measures 4|U µ3 | 2 (|U µ1 | 2 + |U µ2 | 2 ) = 0.99 ± 0.02. We find good agreement between our simplified analysis and the results of Ref. [21] for all different oscillation parameters of interest.
OPERA We include the OPERA collaboration's measurement of tau neutrino appearance via P (ν µ → ν τ ), where 10 ν τ signal events were observed with an expected background of 2.0 ± 0.4 events. Assuming sin 2 (2θ 23 ) = 1 and ∆m 2 32 = 2.5 × 10 −3 eV 2 , OPERA expected 6.8 ± 1.4 signal events. We include this information, assuming a mono-energetic measurement at E OPERA = 17 GeV and L = 730 km, giving results consistent with those from OPERA [22]. Matter effects are included, even though they are small for ν τ appearance oscillation probabilities.
Sterile Neutrino Searches When unitarity is not assumed, sterile neutrino searches provide additional constraints on the unitarity of the LMM. Specifically, we include results from sterile neutrino searches at experiments in regions where such new oscillations would have "averaged out" in the experiment's detector. This corresponds to the high ∆m 2 41 mass-squared splitting region of experimental sensitivities/exclusions and utilizes the zero-distance effects. Searches for anomalous appearance of a different neutrino flavor constrain triangle closure information, where searches for anomalous disappearance constrain row normalizations. References [16,42] provide further explanation of these effects. We do not include sterile search constraints where mass-squared-splitting information is utilized because they do not apply to a generic unitarity violation scenario where we are agnostic about the mass scale of the violation.
Neutrino Appearance Searches: Experiments that search for anomalous appearance (such as NOMAD searching for anomalous ν µ → ν e or ν µ → ν τ ) are sensitive to these zero-distance effects, which, if the LMM is not unitarity, correspond to the non-closure of a unitarity trianglei U αi U * βi = 0. We include results from KARMEN [63], NOMAD [64,65], and CHORUS [66] in these analyses. The LSND [72,73] and MiniBooNE [74] experiments have famously observed an excess of electron-like events in a ν µ beam, which can be interpreted as observing a non-closure between the electron and muon rows of i U ei U * µi 2 ≈ 2.6 × 10 −3 [74]. We do not include information from MiniBooNE and LSND in our analyses due to the tension between these appearance searches and disappearance searches.
Neutrino Disappearance Searches: For muon neutrino disappearance, we include results from the MINOS/MINOS+ experiments [67] that constrain the muon row normalization. In addition, hints for the existence of a sterile neutrino with a new mass-squared splitting around ∆m 2 41 ≈ 1 eV 2 have been observed in a variety of reactor antineutrino experiments (see Refs. [75][76][77][78] for reviews of these), which could point to the electron row being not properly normalized, i.e., i |U ei | 2 = 1. However, in order to interpret these results in terms of a constraint on unitarity, we must go to the averaged-out regime of these experimental sensitivities, which is limited by the predicted reactor antineutrino fluxes [55,56]. We therefore do not include these results due to the uncertainty regarding reactor antineutrino flux predictions.
Fit Results: To demonstrate the validity of our methods, Fig. A2 displays the result of our fit (when unitarity is assumed) to all of the current data discussed above. We show the fit as a measurement of the parameters δ CP and sin 2 θ 23 , where the other, unseen parameters have been marginalized. Such results can be compared against those of NOvA [21] and T2K [23,24], and we find that our results are consistent with other global fits [60][61][62].

Appendix D. Future Experiment Simulations
In our analyses, we include simulations of both the DUNE and JUNO experiments. In this appendix, we briefly describe how these simulations are included and some of their details.
DUNE The upcoming Deep Underground Neutrino Experiment (DUNE) is planned to collect data from beam neutrinos in 2026. Its beam-based goals are driven by a wide-energy ν µ beam, with energies between roughly 0.5 and 10 GeV, with a baseline distance of 1300 km. This will allow DUNE to study both ν µ disappearance and ν e appearance, and allow for a powerful measurement of δ CP . Reference [36] demonstrated DUNE's ability to use its beam to study ν τ appearance as well. We include all three of these channels in our simulations, assuming 7 years of data collection (equal operation in neutrino and antineutrino modes) with a 1.2 MW beam and 40 kt of far detector fiducial mass.
Our simulation of both ν µ disappearance and ν e appearance channels follows those developed for Refs. [50][51][52], and the analyses are designed to match the official collaboration sensitivities and expected signal and background event yields [34,35]. We take the neutrino fluxes from Ref. [48] and neutrino cross sections from Ref. [49]. We include energy uncertainty by using migration matrices, assuming that the energy resolution is σ E = 7% × E ν (GeV) + 3.5% E ν (GeV), consistent with the analyses performed in Refs. [34,35]. We include all of the different background channels discussed in Refs. [34,35] for the ν µ disappearance and ν e appearance channels -the largest of which are due to neutrino neutral-current scattering and 0 π/2 π 3π/2 2π δ CP beam contamination of opposite sign or different flavor neutrinos. Efficiencies for both signal identification and background rejection are both taken to be constant as a function of neutrino energy, where we normalize our expected signal and background event rates to those from Ref. [34].
Our simulation for ν τ appearance channel follows Ref. [36]. For a given true neutrino energy E true ν , we assume that the reconstructed energy follows a Gaussian distribution with a mean energy bE true ν and an uncertainty σ E = rE true ν , where b = 45% and r = 25% [36]. We include a 25% systematic normalization uncertainty to account for uncertainties associated with the ν τ charged-current cross section. We also assume a 30% signal identification efficiency for all hadronically-decaying τ ± events, and that 0.5% of neutral current events will contribute to backgrounds in this search.
For all channels, we include a correlated systematic normalization uncertainty on the muon neutrino flux for each beam mode (separate nuisance parameters for neutrino and antineutrino modes) of 5%. As discussed in Appendix C, if we do not assume the unitarity of the LMM and an experiment has a near detector, the inferred measurement of an oscillation probability must be normalized by the appropriate channel, i.e., the normalization of the muon row of the LMM in this case. DUNE is one such experiment, however, since we include a 5% normalization uncertainty on the muon neutrino flux, this effect is negligible in our analyses. This is because the MINOS/MINOS+ sterile neutrino search constrains this normalization effect to the 2.5% level, so the systematic uncertainty of 5% covers any impact of this type.
Since the measurements obtained using P (ν µ → ν τ ) are much less powerful than those for ν µ disappearance and ν e appearance, we do not include any of the tau appearance simulations in the analyses that assume that the LMM U is unitary -those surrounding Figs. 1 and 2. The ν τ appearance data provide useful information when unitarity is not assumed, however, specifically in measuring the quantity 4|U µ3 | 2 |U τ 3 | 2 , which no other channel has access to. This channel provides important statistical power in measuring the J αi of Fig. 3, specifically those for α = e, µ. Likewise, we do not include any simulation of DUNE's solar neutrino capabilities [38] as we expect JUNO to provide more powerful measurements of the relevant parameters. We discuss this in further detail in Ref. [16].
JUNO The Jiangmen Underground Neutrino Observatory (JUNO) is planned to collect data from reactor neutrinos in 2021. It measures the oscillation of reactor ν e of 2-8 MeV at a propagation distance of 53 km. Matter effects can cause O(1%) level modifications to this oscillation probability for these energies and propagation distance. However, they do not impact measurement sensitivity of the parameters of interest, so we do not include them. This enables a measurement of the neutrino mass hierarchy, the primary science goal for JUNO, as well as of θ 12 and ∆m 2 21 . Our analyses are designed to match the official collaboration sensitivities on sin 2 θ 12 [40]. For reactor flux calculation, we follow the strategy in Ref. [53], taking the fission isotope fractions from Ref. [54], 235 U, 239 Pu, and 241 Pu spectra from Ref. [55], and 238 U spectrum from Ref. [56], and this leads to the following total reactor neutrino flux: For each event, the detected energy, which is smeared with an energy resolution of 3% E(MeV) [40], is the total energy of the positron plus its rest mass. We do not consider detector efficiencies and backgrounds, and only match the total sample size to the CDR nominal choice, 120k events (six years). For systematics, we include a correlated flux uncertainty of 2%, an uncorrelated flux uncertainty of 0.8%, the spectrum shape uncertainty of 1%, and the energy scale uncertainty of 1% [40].
Similarly to the case for DUNE, the oscillation probability that JUNO measures depends on what one assumes of a near detector. There is no official JUNO statement on whether there will be a near detector yet, and hence, we perform our analysis assuming there will not be a near detector.