11institutetext: SLAC National Accelerator Laboratory, 2575 Sand Hill Road,
Menlo Park, CA, USA22institutetext: Theoretical Physics Department, Fermilab,
P.O. Box 500, Batavia, IL 60510, USA

# Current and Future Neutrino Oscillation Constraints on Leptonic Unitarity

Sebastian A. R. Ellis 2 Kevin J. Kelly 1 Shirley Weishi Li
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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.

††preprint: FERMILAB-PUB-20-364-T, SLAC-PUB-17552

##  1 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^{\dagger}U=UU^{\dagger}=\mathbb{I}$, forms
the basis of our understanding of SM fermion mixing Maki:1962mu ;
Cabibbo:1963yz ; Pontecorvo:1967fh ; Kobayashi:1973fv . This theoretical
paradigm has been thoroughly tested to great acclaim in the quark sector
Hocker:2001xe ; Bona:2006ah ; CKMfitter ; UTfit ; Tanabashi:2018oca ;
Wolfenstein:1983yz ; Buras:1994ec ; Charles:2004jd . However, our
understanding of the corresponding leptonic mixing matrix (LMM) remains
limited Qian:2013ora ; Parke:2015goa ; Ellis:2020ehi . 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)
Minkowski:1977sc ; Mohapatra:1979ia ; Schechter:1980gr ; Mohapatra:1980yp ;
Cheng:1980qt ; Zee:1980ai ; Gelmini:1980re ; Ma:2006km . 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 Minkowski:1977sc ; Weinberg:1979sa ; Mohapatra:1979ia ;
Wyler:1982dd ; Langacker:1988up ; Hewett:1988xc ; Nardi:1993ag ;
Tommasini:1995ii ; Gluza:2002vs ; Abazajian:2012ys .

Many studies have been undertaken to study the effect of LMM non-unitarity and
to determine existing and projected constraints on non-unitarity
Bilenky:1992wv ; Bergmann:1998rg ; Czakon:2001em ; Bekman:2002zk ;
Antusch:2006vwa ; FernandezMartinez:2007ms ; Goswami:2008mi ; Antusch:2009pm ;
Antusch:2014woa ; Li:2015oal ; Fernandez-Martinez:2016lgt ; Blennow:2016jkn ;
Fong:2016yyh ; Fong:2017gke ; Coutinho:2019aiy . 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.
Parke:2015goa , 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. Parke:2015goa , 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 $\tau$-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 Ellis:2020ehi , 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
$UU^{\dagger}=U^{\dagger}U=\mathbb{I}$ in Section 5.4. 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. In Appendix C, we discuss the Bayesian approach used in many of our
analyses, and the priors that enter this type of analysis. Appendix D includes
the measurement of the phases present in the LMM, a parameterization-dependent
measurement. Lastly, Appendix E 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.

##  2 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.

###  2.1 Unitarity of the Leptonic Mixing Matrix

Neutrino oscillation studies are generally carried out assuming a unitary
$3\times 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
$\mathcal{O}(10^{-13})$ even for an $\mathcal{O}(\text{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. Buchmuller:1991tu ;
Ingelman:1993ve ; Chang:1994hz ; Loinaz:2003gc ; Branco:2019avf ). In
addition, mass models invoking symmetry arguments may also produce large
unitarity violation Wyler:1982dd ; Hewett:1988xc ; Langacker:1988up ;
Nardi:1993ag ; Tommasini:1995ii . It is therefore important to test the
unitarity of the $3\times 3$ lepton mixing matrix with experimental data.

To test unitarity with oscillation data, one could adopt mathematical
assumptions on the $3\times 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 $\nu_{\alpha}(x)$ can be written as a linear
combination of mass eigenstate fields $\nu_{k}(x)$:

$\nu_{\alpha}(x)=\sum_{k}\mathcal{U}_{\alpha k}\nu_{k}(x)\ .$ (2.1)

The flavor index $\alpha\in e,~{}\mu,~{}\tau,~{}\ldots~{}m$ includes the usual
SM flavor fields, along with $m-3$ possible additional right-handed (sterile)
fields, while the mass index $k\in 1,~{}2,~{}3,~{}\ldots~{}n$, allowing for
$n-3$ additional mass eigenstate fields. $\mathcal{U}_{\alpha k}$ is thus an
$m\times 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 Giunti:2004zf ; Antusch:2006vwa
111Reference Giunti:2003qt  provides a thorough explanation for why
$\mathcal{U}_{\alpha k}^{*}$ is the object appearing in Eq. (2.2), contrary to
the definition of the flavor eigenstate field above. The assumption of zero
neutrino masses is reasonable, as the neutrinos are produced ultra-
relativistically. See Ref. Cohen:2008qb  for a discussion of how this
assumption can break down.

$\ket{\nu_{\alpha}}=\frac{1}{\sqrt{\left(\mathcal{U}\mathcal{U}^{\dagger}\right)_{\alpha\alpha}}}\sum_{k}\mathcal{U}_{\alpha
k}^{*}\ket{\nu_{k}}\ ,$ (2.2)

where the sum over $k$ in the normalization factor is implicit. Only active
flavors participate in weak interactions, so the flavor index $\alpha\in
e,~{}\mu,~{}\tau$. 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 Goswami:2008mi ; Li:2015oal ; Blennow:2016jkn .
For pion-decay sources which are used for most experiments, a conservative
cutoff is to include all mass states below 140 MeV.222Heavier sterile
neutrinos may be produced in such sources either off-shell or from heavier
meson decays. However, subtle non-standard oscillation effects from these
neutrinos are not observable because of their reduced fluxes. We therefore
ignore their contribution. 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: $\langle\nu_{\alpha}|\nu_{\alpha}\rangle=1$. Note that if
$\mathcal{U}$ is not unitary, the flavor states are not necessarily
orthogonal: $\langle\nu_{\alpha}|\nu_{\beta}\rangle\neq\delta_{\alpha\beta}$.

In this work, we focus on understanding the structure of the $3\times 3$
(sub-)matrix of the full $m\times n$ mixing matrix $\mathcal{U}$. We will
refer to the $3\times 3$ mixing matrix as the LMM or $U_{\mathrm{LMM}}$, or
simply $U_{\alpha k}$, where now $\alpha=e,~{}\mu,~{}\tau$ 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
Dentler:2018sju ; Diaz:2019fwt ; Boser:2019rta . For heavier steriles that are
still kinematically accessible, see Refs. Fong:2016yyh ; Fong:2017gke  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 Refs. deGouvea:2015euy ; Bolton:2019pcu  for
comprehensive discussions of these experimental searches.

We organize our discussion of oscillations around the following three cases:

  1. 1.

The “standard” case, where $m=n=3$, $U=\mathcal{U}$.

  2. 2.

The “sub-matrix” case, where $m=n>3$. $\mathcal{U}$ is unitary, and $U$ is
not.

  3. 3.

The “agnostic” case, $m\geqslant 3$ and/or $n\geqslant 3$.
$\mathcal{U}_{\alpha 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\times 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 $\mathcal{U}_{\alpha k}$ should be unitary. Additionally, Refs.
Meloni:2009cg ; Blennow:2016jkn  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
$\mathcal{U}_{\alpha 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.

###  2.2 Mixing Matrix Parameterizations

In the standard scenario, $U_{\mathrm{LMM}}$ is a $3\times 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,
$\theta_{12}$, $\theta_{13}$, and $\theta_{23}$, and one complex phase
$\delta_{\text{CP}}$. Often referred to as the PMNS Pontecorvo:1967fh ;
Maki:1962mu  or PDG Tanabashi:2018oca  parameterization, this form of the LMM
is

$U_{\mathrm{LMM}}=U_{\mathrm{PMNS}}\equiv\left(\begin{array}[]{c c
c}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta_{\text{CP}}}\\\
-s_{12}c_{23}-c_{12}s_{13}s_{23}e^{i\delta_{\text{CP}}}&c_{12}c_{23}-s_{12}s_{13}s_{23}e^{i\delta_{\text{CP}}}&s_{23}c_{13}\\\
s_{12}s_{23}-c_{12}c_{23}s_{13}e^{i\delta_{\text{CP}}}&-c_{12}s_{23}-s_{12}c_{23}s_{13}e^{i\delta_{\text{CP}}}&c_{13}c_{23}\end{array}\right),$
(2.3)

where $s_{ij}\equiv\sin\theta_{ij}$ and $c_{ij}\equiv\cos\theta_{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 ($\theta_{12}$),
reactor ($\theta_{13}$), and atmospheric ($\theta_{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
Esteban:2018azc ; deSalas:2017kay ; Capozzi:2020qhw ; Esteban:2020cvm .

More generally, a complex $3\times 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^{\dagger}U=\mathbb{I}$.
This is equivalent to requiring that all columns of the matrix are normalized
to one:

$N_{k}\equiv\left\lvert U_{ek}\right\rvert^{2}+\left\lvert U_{\mu
k}\right\rvert^{2}+\left\lvert U_{\tau
k}\right\rvert^{2}=1\quad(k=1,~{}2,~{}3),$ (2.4)

as well as requiring that the column unitarity triangles close:

$t_{kl}\equiv U_{ek}^{*}U_{el}+U_{\mu k}^{*}U_{\mu l}+U_{\tau k}^{*}U_{\tau
l}=0\quad(k\neq l;\quad k,~{}l=1,~{}2,~{}3).$ (2.5)

Note that these are nine real constraints as $t_{kl}$ can be complex. Because
$U^{\dagger}U$ is Hermitian, the unitarity condition can equivalently be
written as $UU^{\dagger}=\mathbb{I}$, which can be translated to row
normalization conditions:

$N_{\alpha}\equiv\left\lvert U_{\alpha 1}\right\rvert^{2}+\left\lvert
U_{\alpha 2}\right\rvert^{2}+\left\lvert U_{\alpha
3}\right\rvert^{2}=1\quad(\alpha=e,~{}\mu,~{}\tau),$ (2.6)

and the closure of row unitarity triangles:

$t_{\alpha\beta}\equiv U_{\alpha 1}^{*}U_{\beta 1}+U_{\alpha 2}^{*}U_{\beta
2}+U_{\alpha 3}^{*}U_{\beta
3}=0\quad(\alpha\neq\beta;\quad\alpha,~{}\beta=e,~{}\mu,~{}\tau).$ (2.7)

For the general case where $U$ is a non-unitary $3\times 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:

$U_{\mathrm{LMM}}\equiv\left(\begin{array}[]{l l l}\left\lvert
U_{e1}\right\rvert&\left\lvert U_{e2}\right\rvert e^{i\phi_{e2}}&\left\lvert
U_{e3}\right\rvert e^{i\phi_{e3}}\\\ \left\lvert U_{\mu
1}\right\rvert&\left\lvert U_{\mu 2}\right\rvert&\left\lvert U_{\mu
3}\right\rvert\\\ \left\lvert U_{\tau 1}\right\rvert&\left\lvert U_{\tau
2}\right\rvert e^{i\phi_{\tau 2}}&\left\lvert U_{\tau 3}\right\rvert
e^{i\phi_{\tau 3}}\end{array}\right).$ (2.8)

Here, we have nine magnitudes and four CP-violating phases.333The four phases
can be assigned to any $2\times 2$ sub-matrix. 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_{\mathrm{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_{\alpha i}$, which are defined
as

$\varepsilon_{\alpha\beta\gamma}\varepsilon_{ijk}J_{\alpha
i}\equiv\mathrm{Im}\left(U_{\beta j}U_{\gamma k}U_{\beta k}^{*}U_{\gamma
j}^{*}\right),$ (2.9)

where the $\varepsilon$ are Levi-Civita tensors. It is straightforward to see
that

$\displaystyle\sin\phi_{e2}$ $\displaystyle=$
$\displaystyle-\frac{J_{e2}}{\left\lvert U_{e1}\right\rvert\left\lvert
U_{e2}\right\rvert\left\lvert U_{\mu 1}\right\rvert\left\lvert U_{\mu
2}\right\rvert},\quad\sin\phi_{e3}=\frac{J_{e3}}{\left\lvert
U_{e1}\right\rvert\left\lvert U_{e3}\right\rvert\left\lvert U_{\mu
1}\right\rvert\left\lvert U_{\mu 3}\right\rvert},$ (2.10)
$\displaystyle\sin\phi_{\tau 2}$ $\displaystyle=$ $\displaystyle\frac{J_{\tau
2}}{\left\lvert U_{\tau 1}\right\rvert\left\lvert U_{\tau
2}\right\rvert\left\lvert U_{\mu 1}\right\rvert\left\lvert U_{\mu
2}\right\rvert},\quad\sin\phi_{\tau 3}=-\frac{J_{\tau 3}}{\left\lvert U_{\tau
1}\right\rvert\left\lvert U_{\tau 3}\right\rvert\left\lvert U_{\mu
1}\right\rvert\left\lvert U_{\mu 3}\right\rvert}.$ (2.11)

If $U_{\mathrm{LMM}}$ is unitary, all constructible Jarlskog factors must be
equal to each other and equal to the PMNS matrix Jarlskog invariant
Jarlskog:1985ht : 444This can be used to test the unitarity of the LMM. For
details, see Ref. Ellis:2020ehi .

$J_{\rm PMNS}=c_{12}s_{12}c_{23}s_{23}c_{13}^{2}s_{13}\sin\delta_{\rm CP}.$
(2.12)

Enforcing $J_{\tau 2}=J_{\tau 3}=J_{e2}=J_{e3}=J_{\rm PMNS}$ in Eq. (2.11)
allows for the simple derivation of a relation between the phases $\phi_{e2}$,
$\phi_{e3}$, $\phi_{\tau 2}$, and $\phi_{\tau 3}$ and PMNS parameters when
$U_{\rm LMM}=U_{\rm PMNS}$.

Finally, we briefly discuss the sub-matrix case, where $U$ is a $3\times 3$
sub-matrix of a larger unitary matrix $\mathcal{U}$. This introduces two
additional constraints on the structure of $U$: the row and column
normalizations of LMM must not exceed unity:

$N_{\alpha}\leqslant 1,\quad N_{k}\leqslant 1.$ (2.13)

Further, by applying the Cauchy-Schwarz inequality on the vectors
$\big{\\{}\mathcal{U}_{\alpha k}\big{\\}}$, where $\alpha$ or $k$ runs from 4,
5… $n$, we obtain the following inequalities Parke:2015goa  555One can
directly apply Cauchy-Schwarz inequality on the matrix $U$, which leads to the
weaker conditions  $\left|t_{\alpha\beta}\right|^{2}\leqslant
N_{\alpha}N_{\beta},\quad\left|t_{kl}\right|^{2}\leqslant N_{k}N_{j}.$ (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.:

$\displaystyle\left|t_{\alpha\beta}\right|^{2}$ $\displaystyle\leqslant$
$\displaystyle\left(1-N_{\alpha}\right)\left(1-N_{\beta}\right)$
$\displaystyle\left|t_{kl}\right|^{2}$ $\displaystyle\leqslant$
$\displaystyle\left(1-N_{k}\right)\left(1-N_{l}\right).$ (2.15)

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_{\rm 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 Xing:2007zj
; Escrihuela:2015wra ; Blennow:2016jkn ; C:2017scx  is the following:

$U_{\mathrm{LMM}}\equiv NU_{\mathrm{PMNS}}=\left(\begin{array}[]{c c
c}\alpha_{11}&0&0\\\ \alpha_{21}&\alpha_{22}&0\\\
\alpha_{31}&\alpha_{32}&\alpha_{33}\end{array}\right)U_{\mathrm{PMNS}}.$
(2.16)

When there are three active neutrinos and any number of sterile neutrinos, one
can express $\alpha_{kl}$ in terms of mixing angles and phases between active
and sterile neutrino mixing. For the full expressions, see Ref
Escrihuela:2015wra . Here, unitarity is achieved in the limit that
$\alpha_{kl}\to\delta_{kl}$. The off-diagonal $\alpha_{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., $\alpha_{kk}\to 1$ and
$\alpha_{k\neq l}\to 0$, compared to the MP parameterization. However, it is
not straightforward to map between the “$\alpha$” parameterization and the
individual, specific constraints of unitarity – the normalizations and
closures of columns and rows of $U_{\rm LMM}$. As constructed, the map between
the $\alpha_{kl}$ and the normalization of rows and closures between two
different rows is relatively simple:

$\displaystyle N_{e}$ $\displaystyle=$ $\displaystyle\alpha_{11}^{2},$ (2.17)
$\displaystyle N_{\mu}$ $\displaystyle=$
$\displaystyle\alpha_{22}^{2}+\left\lvert\alpha_{21}\right\rvert^{2},$ (2.18)
$\displaystyle N_{\tau}$ $\displaystyle=$
$\displaystyle\alpha_{33}^{2}+\left\lvert\alpha_{32}\right\rvert^{2}+\left\lvert\alpha_{31}\right\rvert^{2},$
(2.19) $\displaystyle t_{e\mu}$ $\displaystyle=$
$\displaystyle\alpha_{11}\alpha_{21}\ ,$ (2.20) $\displaystyle t_{e\tau}$
$\displaystyle=$ $\displaystyle\alpha_{11}\alpha_{32}\ ,$ (2.21)
$\displaystyle t_{\mu\tau}$ $\displaystyle=$
$\displaystyle\alpha_{21}^{*}\alpha_{31}+\alpha_{22}\alpha_{32}\ .$ (2.22)

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 $\alpha_{kl}$ as well as the mixing angles and
$\delta_{\rm CP}$ from the PMNS parameterization.

Since our goal in this work is to determine the current and future constraints
on $N_{\alpha}$, $N_{k}$, $t_{\alpha\beta}$, and $t_{kl}$, we use the MP
parameterization which has a straightforward map between the input parameters
and these quantities. Using the MP paramaterization has a second advantage, in
that the majority of the underlying inputs are parameterization-independent.
Specifically, all of the magnitudes-squared $\left\lvert U_{\alpha
k}\right\rvert^{2}$ are independent of the adopted parameterization, which
makes translation between different experimental results in this context
simpler.

###  2.3 Oscillation Probabilities

We now review the oscillation formulae for different parameterizations. See
Refs. Antusch:2006vwa ; Li:2015oal  for further discussions. In vacuum, the
mass states $|\nu_{k}\rangle$ are eigenstates of the Hamiltonian, such that
they form an orthogonal basis, $\langle\nu_{k}|\nu_{l}\rangle=\delta_{kl}$,
and evolve in time as

$|\nu_{k}(t)\rangle=e^{-iE_{k}t}|\nu_{k}\rangle.$ (2.23)

A flavor state is created as Eq. (2.2). After traveling for time $t$, the
flavor state is evolved to

$|\nu_{\alpha}(t)\rangle=\frac{1}{\sqrt{N_{\alpha}}}\sum_{k=1}^{3}U_{\alpha
k}^{*}e^{-iE_{k}t}|\nu_{k}\rangle.$ (2.24)

The oscillation probability for a neutrino of energy $E$ produced as
$\nu_{\alpha}$ and detected as $\nu_{\beta}$ after propagating $L=ct$ is
therefore

$\displaystyle P(\nu_{\alpha}\rightarrow\nu_{\beta})$ $\displaystyle=$
$\displaystyle|\langle\nu_{\beta}|\nu_{\alpha}(t)\rangle|^{2}=\frac{1}{N_{\alpha}N_{\beta}}\left|\displaystyle\sum_{k=1}^{3}U_{\alpha
k}^{*}U_{\beta k}e^{-iE_{k}t}\right|^{2}.$ (2.25)

Here, $E_{k}=m_{k}^{2}/2E_{\nu}$. In order to better understand the behavior
of the oscillation probability, we separate the discussion here into two
cases, $\alpha=\beta$ (disappearance/survival probability) and
$\alpha\neq\beta$ (appearance probability). Defining $\Delta_{ij}\equiv\Delta
m_{ij}^{2}L/2E_{\nu}$, the disappearance/survival probability may be written
as

$\displaystyle P(\nu_{\alpha}\rightarrow\nu_{\alpha})$ $\displaystyle=$
$\displaystyle 1-\frac{4\left\lvert U_{\alpha
2}\right\rvert^{2}\left(\left\lvert U_{\alpha 1}\right\rvert^{2}+\left\lvert
U_{\alpha
3}\right\rvert^{2}\right)}{N_{\alpha}^{2}}\sin^{2}\left(\frac{\Delta_{21}}{2}\right)-\frac{4\left\lvert
U_{\alpha 3}\right\rvert^{2}\left(\left\lvert U_{\alpha
1}\right\rvert^{2}+\left\lvert U_{\alpha
2}\right\rvert^{2}\right)}{N_{\alpha}^{2}}\sin^{2}\left(\frac{\Delta_{31}}{2}\right)$
(2.26) $\displaystyle+$ $\displaystyle\frac{8\left\lvert U_{\alpha
2}\right\rvert^{2}\left\lvert U_{\alpha
3}\right\rvert^{2}}{N_{\alpha}^{2}}\sin\left(\frac{\Delta_{21}}{2}\right)\sin\left(\frac{\Delta_{31}}{2}\right)\cos\left(\frac{\Delta_{32}}{2}\right).$

For the appearance probability, we define
$\varphi_{k}^{\alpha\beta}\equiv\mathrm{arg}\left(U_{\alpha k}^{*}U_{\beta
k}\right)$ and
$\varphi^{\alpha\beta}\equiv\mathrm{arg}\left(t_{\alpha\beta}\right)$. The
appearance probability may be written as

$\displaystyle P(\nu_{\alpha}\rightarrow\nu_{\beta})$ $\displaystyle=$
$\displaystyle\frac{\left|t_{\alpha\beta}\right\rvert^{2}}{N_{\alpha}N_{\beta}}+\frac{4\left\lvert
U_{\alpha 2}\right\rvert^{2}\left\lvert U_{\beta
2}\right\rvert^{2}}{N_{\alpha}N_{\beta}}\sin^{2}\left(\frac{\Delta_{21}}{2}\right)+\frac{4\left\lvert
U_{\alpha 3}\right\rvert^{2}\left\lvert U_{\beta
3}\right\rvert^{2}}{N_{\alpha}N_{\beta}}\sin^{2}\left(\frac{\Delta_{31}}{2}\right)$
(2.27) $\displaystyle+$ $\displaystyle\frac{8\left\lvert U_{\alpha
2}\right\rvert\left\lvert U_{\beta 2}\right\rvert\left\lvert U_{\alpha
3}\right\rvert\left\lvert U_{\beta
3}\right\rvert}{N_{\alpha}N_{\beta}}\sin\left(\frac{\Delta_{21}}{2}\right)\sin\left(\frac{\Delta_{31}}{2}\right)\cos\left(\frac{\Delta_{32}}{2}+\varphi_{2}^{\alpha\beta}-\varphi_{3}^{\alpha\beta}\right)$
$\displaystyle+$ $\displaystyle\frac{4\left\lvert
t_{\alpha\beta}\right\rvert}{N_{\alpha}N_{\beta}}\left[\left\lvert U_{\alpha
2}\right\rvert\left\lvert U_{\beta
2}\right\rvert\sin\left(\frac{\Delta_{21}}{2}\right)\sin\left(\frac{\Delta_{21}}{2}+\varphi^{\alpha\beta}-\varphi^{\alpha\beta}_{2}\right)\right.$
$\displaystyle+$ $\displaystyle\left.\left\lvert U_{\alpha
3}\right\rvert\left\lvert U_{\beta
3}\right\rvert\sin\left(\frac{\Delta_{31}}{2}\right)\sin\left(\frac{\Delta_{31}}{2}+\varphi^{\alpha\beta}-\varphi^{\alpha\beta}_{3}\right)\right].$

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 $\Delta m_{21}^{2}$
and $\Delta m_{32}^{2}$, 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 $\Delta m_{21}^{2}$ or $\Delta m_{31}^{2}$
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 zero-distance effect Langacker:1988up ; Antusch:2006vwa .
When $L\to 0$ (or in an experiment where $\Delta_{21}$ and $\Delta_{31}\ll
1$), the disappearance probability $P(\nu_{\alpha}\to\nu_{\alpha})\to 1$, but
the appearance probability becomes

$P(\nu_{\alpha}\rightarrow\nu_{\beta},L\sim
0)\simeq\frac{1}{N_{\alpha}N_{\beta}}\left|\sum_{k=1}^{3}U_{\alpha
k}^{*}U_{\beta k}\right|^{2}=\frac{\left\lvert
t_{\alpha\beta}\right\rvert^{2}}{N_{\alpha}N_{\beta}}$ (2.28)

This implies that searches for short-baseline anomalous appearance from one
flavor eigenstate $\alpha$ to another $\beta$ provides a direct constraint on
the closure between the $\alpha$ and $\beta$ 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:

$V_{\alpha\beta}=\sqrt{2}G_{\mathrm{F}}\left(\begin{array}[]{c c
c}n_{e}-n_{n}/2&0&0\\\ 0&-n_{n}/2&0\\\ 0&0&-n_{n}/2\end{array}\right),$ (2.29)

where $n_{e}$ and $n_{n}$ are the electron and neutron density in the medium,
respectively. This potential is rotated by $U_{\rm LMM}$ into the mass basis,
and combined with the (mass-basis-diagonal) energy values
$\Delta_{k1}/(2E_{\nu})$. 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^{\dagger}U=\mathbb{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 non-unitary 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_{\alpha\beta}$ 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 Antusch:2006vwa , and we
adopt it for the remainder of our work.

###  2.4 Normalization Effects

In this subsection we discuss further non-trivial consequences of a non-
unitary LMM. Much of this discussion is adapted from Ref. Antusch:2006vwa .
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 $\alpha$, and the
corresponding CC cross section, may be expressed (relative to their unitary,
SM expectations) as

$\Phi_{\alpha}(E)=N_{\alpha}\Phi^{\rm
SM}_{\alpha}(E),\quad\quad\sigma_{\beta}(E)=N_{\beta}\sigma_{\beta}^{\rm
SM}(E).$ (2.30)

The neutral-current cross section is also modified,

$\displaystyle\sigma_{k}^{\rm NC}(E)$ $\displaystyle=\sigma_{\rm SM}^{\rm
NC}(E)\sum_{l}\Big{|}(U^{\dagger}U)_{kl}\Big{|}^{2}\ .$ (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
$\beta$, $n_{\beta}$, is given by

$n_{\beta}\propto\Phi_{\alpha}P(\nu_{\alpha}\rightarrow\nu_{\beta})\sigma_{\beta}=N_{\alpha}\Phi^{\rm
SM}_{\alpha}P(\nu_{\alpha}\rightarrow\nu_{\beta})N_{\beta}\sigma_{\beta}^{\rm
SM}=\Phi^{\rm
SM}_{\alpha}\hat{P}(\nu_{\alpha}\rightarrow\nu_{\beta})\sigma_{\beta}^{\rm
SM}\ .$ (2.32)

Thus, if an experimental analysis is performed assuming SM predictions of
fluxes and cross-sections as truth, the measurement of $n_{\beta}$ corresponds
to an inferred measurement of $\hat{P}(\nu_{\alpha}\rightarrow\nu_{\beta})$.
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 by

$\hat{P}(\nu_{\alpha}\rightarrow\nu_{\beta})=\left|\sum_{k=1}^{3}U_{\alpha
k}^{*}U_{\beta k}e^{-iE_{k}t}\right|^{2}\ ,$ (2.33)

such that these measurements are not sensitive to the normalization factors
$N_{\alpha}$ and $N_{\beta}$.

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_{\alpha}^{\rm ND}(E)$, where $\alpha$ is the flavor label:

$n_{\alpha}^{\rm
ND}(E)\propto\Phi_{\alpha}(E)\sigma_{\alpha}(E)P_{\alpha\alpha}(E;L=0)\ ,$
(2.34)

At $L=0$, the oscillation probability is $P_{\alpha\alpha}(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
$\Phi_{\alpha}$ and $\sigma_{\alpha}$, so the measured spectrum can be
expressed as

$n_{\alpha}^{\rm ND}(E)\propto\Phi_{\alpha}^{\rm SM}(E)\sigma_{\alpha}^{\rm
SM}(E)N_{\alpha}^{2}.$ (2.35)

Therefore, with a sufficiently precise measurement of the spectrum, as well as
understanding of the SM-expected flux and cross section, a constraint on
$N_{\alpha}^{2}$ can be placed. These results are usually reported in the
context of a limit on disappearance probabilities, i.e.
$P_{\alpha\alpha}(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_{\alpha}$. When considering searches for anomalous short-baseline
appearance of a flavor $\alpha$ to a flavor $\beta$, Eq. (2.35) must be
modified accordingly, and we find that these searches are sensitive to

$n_{\beta}^{\rm ND}(E)\propto\left\lvert t_{\alpha\beta}\right\rvert^{2}\ ,$
(2.36)

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_{\alpha\beta}\equiv
P(\nu_{\alpha}\to\nu_{\beta};E,L)$ by measuring a near detector spectrum
$n_{\alpha}^{\rm ND}$ as above, as well as a far detector spectrum
$n_{\beta}^{\rm FD}$,

$n_{\beta}^{\rm
FD}(E)\propto\Phi_{\alpha}(E)\sigma_{\beta}(E)P_{\alpha\beta}.$ (2.37)

We may express the measured ratio of oscillation probabilities as

$P_{\alpha\beta}=\frac{n_{\beta}^{\rm FD}(E)}{n_{\alpha}^{\rm
ND}(E)}\frac{\sigma_{\alpha}(E)}{\sigma_{\beta}(E)}=\frac{n_{\beta}^{\rm
FD}(E)}{n_{\alpha}^{\rm ND}(E)}\frac{\sigma^{\rm SM}_{\alpha}(E)}{\sigma^{\rm
SM}_{\beta}(E)}\frac{N_{\alpha}}{N_{\beta}}=\hat{P}_{\alpha\beta}\frac{N_{\alpha}}{N_{\beta}}.$
(2.38)

For disappearance measurements where $\alpha=\beta$, the measured
probabilities are the true probabilities, e.g., Eq. (2.27) for the vacuum
case. For appearance measurements where $\alpha\neq\beta$, the measured
probabilities $\hat{P}_{\alpha\beta}$ are the true probabilities
$P_{\alpha\beta}$ multiplied by a factor $N_{\beta}/N_{\alpha}$.

##  3 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 $\Delta
m_{31}^{2}$). While there is a long-standing tension between solar and reactor
experiment measurements of $\Delta m_{21}^{2}$, we find that allowing the
mass-squared splittings to vary does not impact the results of measuring the
elements of the mixing matrix.

Due to the complexity of its simulation (both the computational expense of
simulating the oscillation probabilities and the many relevant systematic
uncertainty parameters), we do not include the atmospheric neutrino results
from Super-Kamiokande in our analysis Abe:2017aap ; Jiang:2019xwn . While
$\chi^{2}$ tables are provided by the collaboration for the results presented
in Ref. Abe:2017aap , these assume three-neutrino mixing and a unitary
leptonic mixing matrix. Since our goal is to keep the data analyzed consistent
whether we are assuming unitarity or not, we do not include this in any of the
following fits.

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”).

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$ | $(\left\lvert U_{e1}\right\rvert^{2}+\left\lvert U_{e2}\right\rvert^{2})N_{2}^{2}+\left\lvert U_{e3}\right\rvert^{2}N_{3}^{2}$
Solar Charged Current | $\sin^{2}\theta_{12}\cos^{4}\theta_{13}+\sin^{4}\theta_{13}$ | $\left\lvert U_{e2}\right\rvert^{2}(\left\lvert U_{e1}\right\rvert^{2}+\left\lvert U_{e2}\right\rvert^{2})+\left\lvert U_{e3}\right\rvert^{4}$
KamLAND | $\cos^{4}\theta_{13}\sin^{2}\left(2\theta_{12}\right)$ | $4\left\lvert U_{e1}\right\rvert^{2}\left\lvert U_{e2}\right\rvert^{2}$
Daya Bay | $\sin^{2}\left(2\theta_{13}\right)$ | $4\left\lvert U_{e3}\right\rvert^{2}(\left\lvert U_{e1}\right\rvert^{2}+\left\lvert U_{e2}\right\rvert^{2})/N_{e}^{2}$
Sterile Neutrino $P_{\alpha\beta}$ ($\alpha\neq\beta$) | $0$ | $\left\lvert t_{\alpha\beta}\right\rvert^{2}$
OPERA | $\cos^{4}\theta_{13}\sin^{2}\left(2\theta_{23}\right)$ | $4\left\lvert U_{\mu 3}\right\rvert^{2}\left\lvert U_{\tau 3}\right\rvert^{2}/N_{\mu}^{2}$
Long-baseline $P_{\mu e}$ | $\sin^{2}\theta_{23}\sin^{2}\left(2\theta_{13}\right)$ | $4\left\lvert U_{e3}\right\rvert^{2}\left\lvert U_{\mu 3}\right\rvert^{2}/N_{\mu}^{2}$
(T2K, NOvA, DUNE, T2HK)
Long-baseline $P_{\mu\mu}$ | $4\cos^{2}\theta_{13}\sin^{2}\theta_{23}\left(1-\cos^{2}\theta_{13}\sin^{2}\theta_{23}\right)$ | $4\left\lvert U_{\mu 3}\right\rvert^{2}(\left\lvert U_{\mu 1}\right\rvert^{2}+\left\lvert U_{\mu 2}\right\rvert^{2})/N_{\mu}^{2}$
(T2K, NOvA, DUNE, T2HK)

###  3.1 Solar Neutrino Measurements

The Sudbury Neutrino Experiment (SNO) Aharmim:2011vm  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 $\ket{\nu_{k}}$ is related to the effective matrix-element-squared
in propagation in matter, which we call $\lvert\widetilde{U}_{ek}\rvert^{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
$\lvert\widetilde{U}_{ek}\rvert^{2}=\left\\{0,(\lvert U_{e1}\rvert^{2}+\lvert
U_{e2}\rvert^{2})/N_{e},\lvert U_{e3}\rvert^{2}/N_{e}\right\\}$.

Following discussion similar to that of Sec. 2.4, we can write the detected
number of neutral current events as

$\displaystyle n_{\rm N}\propto\sum_{k}\Phi_{k}^{\rm exiting\
sun}\sigma_{k}^{\rm NC}$ $\displaystyle=\sum_{k}\left(\Phi_{\rm
SM}N_{e}\lvert\widetilde{U}_{ek}\rvert^{2}\right)\left(\sigma_{\rm SM}^{\rm
NC}\left(\sum_{l}\Big{|}(U^{\dagger}U)_{kl}\Big{|}^{2}\right)\right),$ (3.39)
$\displaystyle=\Phi_{\rm SM}\sigma_{\rm SM}^{\rm
NC}\sum_{k}N_{e}\lvert\widetilde{U}_{ek}\rvert^{2}\left(\sum_{l}\Big{|}(U^{\dagger}U)_{kl}\Big{|}^{2}\right),$
(3.40) $\displaystyle=\Phi_{\rm SM}\sigma_{\rm SM}^{\rm NC}\left[\left(\lvert
U_{e1}\rvert^{2}+\lvert U_{e2}\rvert^{2}\right)\left(N_{2}^{2}+\lvert
t_{12}\rvert^{2}+\lvert t_{23}\rvert^{2}\right)+\lvert
U_{e3}\rvert^{2}\left(N_{3}^{2}+\lvert t_{13}\rvert^{2}+\lvert
t_{23}\rvert^{2}\right)\right],$ (3.41) $\displaystyle=\Phi_{0}\sigma_{\rm
SM}^{\rm NC}\hat{P}_{\rm NC},$ (3.42)

where we define

$\displaystyle\hat{P}_{\mathrm{NC}}$ $\displaystyle\equiv\left(\lvert
U_{e1}\rvert^{2}+\lvert U_{e2}\rvert^{2}\right)\left(N_{2}^{2}+\lvert
t_{12}\rvert^{2}+\lvert t_{23}\rvert^{2}\right)+\lvert
U_{e3}\rvert^{2}\left(N_{3}^{2}+\lvert t_{13}\rvert^{2}+\lvert
t_{23}\rvert^{2}\right),$ (3.43) $\displaystyle=\left(\lvert
U_{e1}\rvert^{2}+\lvert U_{e2}\rvert^{2}\right)N_{2}^{2}+\lvert
U_{e3}\rvert^{2}N_{3}^{2}$ $\displaystyle+\left(\lvert U_{e1}\rvert^{2}+\lvert
U_{e2}\rvert^{2}\right)\left(\lvert t_{12}\rvert^{2}+\lvert
t_{23}\rvert^{2}\right)+\lvert U_{e3}\rvert^{2}\left(\lvert
t_{13}\rvert^{2}+\lvert t_{13}\rvert^{2}\right).$ (3.44)

The terms of the last line in Eq. (3.44) are small relative to those in the
line above it, so we can express the measured oscillation probability at
leading order as $\hat{P}_{\rm NC}\simeq(\lvert U_{e1}\rvert^{2}+\lvert
U_{e2}\rvert^{2})N_{2}^{2}+\lvert U_{e3}\rvert^{2}N_{3}^{2}$. This measurement
is limited by uncertainties associated with the 8B neutrino flux prediction
from the standard solar model Vinyoles:2016djt . We conservatively assume that
this probability is measured at the 25% level, $\hat{P}_{\rm NC}=1\pm 0.25$.

In addition, SNO, and Super-Kamiokande (Super-K) Abe:2016nxk  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}\theta_{12}\cos^{4}\theta_{13}+\sin^{4}\theta_{13}$. We take
the results from a preliminary joint SNO and Super-K analysis SKNu2020 , which
reports $\sin^{2}\theta_{12}=0.306\pm 0.014$ when using a prior of
$\sin^{2}\theta_{13}=0.0219\pm 0.0014$ in their analysis, and interpret it as
a measurement of the survival probability $P_{ee}=0.2932\pm 0.0134$.666The
results in Ref. SKNu2020  result in slightly stronger constraints on $P_{ee}$
than those in the previously reported Ref. Abe:2016nxk . More interestingly,
the results from Ref. SKNu2020  prefer a larger value of $\Delta m_{21}^{2}$
($6.11\times 10^{-5}$ eV2) than those of Ref. Abe:2016nxk  ($4.8\times
10^{-5}$ eV2), more consistent with other observations from KamLAND. If
unitarity is not assumed, the survival probability is

$P_{ee}=\left\lvert U_{e2}\right\rvert^{2}(\left\lvert
U_{e1}\right\rvert^{2}+\left\lvert
U_{e2}\right\rvert^{2})+\left|U_{e3}\right\rvert^{4}.$ (3.45)

We also include the measured value of $\Delta m_{21}^{2}=\left(6.11\pm
1.21\right)\times 10^{-5}$ eV2 from the joint analysis SKNu2020 .

###  3.2 KamLAND

The Kamioka Liquid Scintillator Antineutrino Detector (KamLAND) was an
experiment that observed the oscillation of reactor electron antineutrinos
$\overline{\nu}_{e}$ (with energies between $\simeq$2–10 MeV) at distances of
roughly 180 km. This allows KamLAND to be sensitive to the solar mass-squared
splitting, measuring $\Delta m_{21}^{2}=\left(7.50\pm 0.20\right)\times
10^{-5}$ eV2. We use a slightly older measurement from Ref. Gando:2010aa  to
be consistent with the corresponding measurement of the oscillation
probability. A more recent analysis measures this mass-squared splitting
slightly more precisely Gando:2013nba , but this does not affect our
results777The more recent Ref. Gando:2013nba  reports a slightly smaller
preferred value of $\tan^{2}\theta_{12}=0.436$ than Ref. Gando:2010aa ’s
$\tan^{2}\theta_{12}=0.452$. While this shift is at the ${\sim}0.5\sigma$
level, the more powerful solar neutrino measurements are more important for
the resulting fits than KamLAND..

Appendix B of Ref. Gando:2010aa  gives weighted measurements and uncertainties
on the oscillation probability $P(\overline{\nu}_{e}\to\overline{\nu}_{e})$
for different values of
$x\equiv\langle\sin^{2}{\left(\Delta_{21}/2\right)}\rangle$, 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 $\Delta
m_{31}^{2}$ have averaged out in KamLAND) is written as Gando:2010aa

$P_{ee}=(\cos^{4}\theta_{13}+\sin^{4}\theta_{13})-\cos^{4}\theta_{13}\sin^{2}\left(2\theta_{12}\right)x.$
(3.46)

When not assuming unitarity, this becomes

$P_{ee}=1-\frac{2\left\lvert U_{e3}\right\rvert^{2}(\left\lvert
U_{e1}\right\rvert^{2}+\left\lvert
U_{e2}\right\rvert^{2})}{N_{e}^{2}}-\frac{4\left\lvert
U_{e1}\right\rvert^{2}\left\lvert U_{e2}\right\rvert^{2}}{N_{e}^{2}}x.$ (3.47)

KamLAND does not use a near detector in its analysis, so we need to re-scale
this oscillation probability by $N_{e}^{2}$. The probability measured by
KamLAND is therefore

$\hat{P}_{ee}=|U_{e3}|^{4}+\left(\left\lvert
U_{e1}\right\rvert^{2}+\left\lvert
U_{e2}\right\rvert^{2}\right)^{2}-4\left\lvert
U_{e1}\right\rvert^{2}\left\lvert U_{e2}\right\rvert^{2}x\ .$ (3.48)

The averaging over $x$ performed in Ref. Gando:2010aa  depends on the neutrino
matter effects during propagation, which we discussed in Section 2.3. For
long-baseline $\overline{\nu}_{e}$ oscillations such as those being measured
at KamLAND, the deviations from the (small) unitarity-assumed matter effects
are, to leading order, dependent on quantities such as $N_{2}-N_{1}$ (the
differences between two different column normalizations). Where our global fit
is concerned, these differences are constrained to be greatly sub-dominant
relative to other contributions to the matter-induced mixing angle and mass-
squared splitting. Given this sub-dominance, and the fact that other
experiments (solar neutrino experiments and Daya Bay, specifically) provide
more powerful measurements of the electron row elements of the LMM, we find
this method to be the most complete way of including KamLAND’s results in such
an analysis.

###  3.3 Daya Bay

The Daya Bay experiment observes the disappearance of electron antineutrinos
$P(\overline{\nu}_{e}\to\overline{\nu}_{e})$. Daya Bay operates in the regime
of $\Delta_{21}\ll 1$ and the coefficient of the dominant oscillation term is
given by $\sin^{2}\left(2\theta_{13}\right)$ in the PMNS parameterization.
Given the derivations in Appendix B, we apply Eq. (B.94), which shows that the
disappearance probability (without assuming unitarity, and including near
detector normalization effects) is

$P_{ee}\simeq 1-\frac{4\left\lvert U_{e3}\right\rvert^{2}(\left\lvert
U_{e1}\right\rvert^{2}+\left\lvert
U_{e2}\right\rvert^{2})}{N_{e}^{2}}\sin^{2}\left(\frac{\Delta_{31}}{2}\right)\
.$ (3.49)

The most recent measurement from Daya Bay is $\sin^{2}(2\theta_{13})=0.0856\pm
0.0029$ Adey:2018zwh , which we map on to a measurement of the coefficient in
Eq. (3.49) in our fit. Daya Bay’s measurement of the mass-squared-splitting
$\left\lvert\Delta m_{31}^{2}\right\rvert=\left(2.471\pm 0.070\right)\times
10^{-3}$ eV2 is also included in our analysis Adey:2018zwh .

###  3.4 OPERA

The Oscillation Project with Emulsion-tRacking Apparatus (OPERA) collaboration
has completed data collection and has provided results of searches for
$\nu_{\mu}\to\nu_{\tau}$ oscillations, where the neutrinos have a mean energy
of $17$ GeV and a baseline of $730$ km Agafonova:2018auq . Across four
channels, the experiment observed 10 $\nu_{\tau}$ signal events with an
expectation of $2.0\pm 0.4$ background events, and $6.8\pm 1.4$ signal events
(assuming $\sin^{2}\left(2\theta_{23}\right)=1$ and $\Delta
m_{32}^{2}=2.5\times 10^{-3}$ eV2). This measurement predominantly constrains
the quantity $|t_{\mu\tau}|^{2}$ and the product $|U_{\mu 3}|^{2}|U_{\tau
3}|^{2}$. For our analysis, we assumed that OPERA measures the oscillation
probability $P(\nu_{\mu}\to\nu_{\tau})$ 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_{\mu}N_{\tau}$ 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 $\left\lvert\Delta
m_{32}^{2}\right\rvert=\left(2.7\pm 0.7\right)\times 10^{-3}$ eV2
Agafonova:2018auq .

###  3.5 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. T2KNu2020 ,
using the figures therein to extract the expected signal and background rates
for $\nu_{e}$ and $\overline{\nu}_{e}$ events as a function of oscillation
parameters.888In Ref. Ellis:2020ehi , we had used the detailed results of Ref.
Abe:2018wpn  to perform our simulations. We find that our updated simulation
with the preliminary results of Ref. T2KNu2020  are mostly consistent with the
previous published result, up to the fact that newer data are now included.
Reference Abe:2019vii  also includes newer data than Ref. Abe:2018wpn , and
its results are more-or-less consistent with ours. In Ref. T2KNu2020 ,
expected signal plus background rates are shown for the $\nu_{\mu}\to\nu_{e}$
and $\overline{\nu}_{\mu}\to\overline{\nu}_{e}$ channels for values of
$\delta_{\rm CP}$ of $-\pi/2$, $0$, $\pi/2$, and $\pi$. 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_{\rm T2K}=600$ MeV at a fixed baseline
length $L_{\rm T2K}=295$ km, with a constant matter density of $\rho_{\rm
T2K}=2.6$ g/cm3 Abe:2018wpn . Using the figures provided in Ref. T2KNu2020 ,
we arrive at the following expected event rates:

$\displaystyle N_{\nu_{e}}^{\rm T2K}$ $\displaystyle=21.90+1282.84\times
P(\nu_{\mu}\to\nu_{e},L_{\rm T2K},E_{\rm T2K}),$ (3.50) $\displaystyle
N_{\overline{\nu}_{e}}^{\rm T2K}$ $\displaystyle=10.66+179.59\times
P(\bar{\nu}_{\mu}\to\bar{\nu}_{e},L_{\rm T2K},E_{\rm T2K}).$ (3.51)

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.50)-(3.51), along with
stars that indicate the values given by the figure in Ref. T2KNu2020 . 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. T2KNu2020 ), and its $\pm
1\sigma$ 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_{\mu}$ to account for the use of a near detector for T2K.

Figure 1: Validation of our T2K analysis: each panel displays the number of
expected and observed events for a different data sample at T2K: neutrino mode
$\nu_{e}$ events (left) and antineutrino mode $\overline{\nu}_{e}$ events
(right). All oscillation parameters except $\delta_{\rm CP}$ and $\Delta
m_{31}^{2}$ are fixed to the values listed in the left panel. The green
(purple) lines and stars assume the normal (inverted) mass ordering, with
$\Delta m_{31}^{2}=2.565(-2.46)\times 10^{-3}$ eV2. Stars correspond to the
values given for expected signal and background rates extracted from Ref.
T2KNu2020 , and the lines correspond to our estimates from Eqs. (3.50)-(3.51).
Grey regions indicate the observed number of events in each channel with a
$\pm 1\sigma$ range from Poissonian statistics.

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

$\displaystyle P_{\mu\mu}\simeq 1-\frac{4\left\lvert U_{\mu
3}\right\rvert^{2}(\left\lvert U_{\mu 1}\right\rvert^{2}+\left\lvert U_{\mu
2}\right\rvert^{2})}{N_{\mu}^{2}}\sin^{2}\left(\frac{\Delta_{31}}{2}\right)\
,$ (3.52)

such that the relevant coefficient is

$C_{\mu\mu}^{\mathrm{Dis.}}=\frac{4\left\lvert U_{\mu
3}\right\rvert^{2}\left(\left\lvert U_{\mu 1}\right\rvert^{2}+\left\lvert
U_{\mu 2}\right\rvert^{2}\right)}{N_{\mu}^{2}}.$ (3.53)

We include T2K’s combined $\nu_{\mu}$ and $\overline{\nu}_{\mu}$ disappearance
searches by assuming a measurement of $C_{\mu\mu}=1.0\ \pm\ 0.03$, and find
that it gives us results consistent with Refs. Abe:2018wpn ; Abe:2019vii ;
T2KNu2020 . 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 $\left\lvert\Delta
m_{32}^{2}\right\rvert=\left(2.49\pm 0.082\right)\times 10^{-3}$ eV2 T2KNu2020
.

###  3.6 NOvA

The last long-baseline experiment we include is the NuMI Off-Axis $\nu_{e}$
Appearance (NOvA) experiment. It operates at similar values of $L/E_{\nu}$ 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_{\rm NOvA}=1.9$ GeV, $L_{\rm NOvA}=810$ km, with a constant matter
density along the path of propagation of $\rho_{\rm NOvA}=2.84$ g/cm3
NOvA:2018gge . Using the preliminary results of Ref. NOvANu2020 , 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.999As with T2K, we have
also used the more detailed Ref. Acero:2019ksn  to perform a similar analysis
with less overall data as a cross-check. We find that our simulations match
both Ref. Acero:2019ksn  and Ref. NOvANu2020  well.

Reference NOvANu2020  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,

$\displaystyle N_{\nu,\ {\rm CCQE}}^{\rm NOvA}$
$\displaystyle=29.09+1202.65\times P\left(\nu_{\mu}\to\nu_{e},L_{\rm
NOvA},E_{\rm NOvA}\right),$ (3.54) $\displaystyle N_{\overline{\nu},\ {\rm
CCQE}}^{\rm NOvA}$ $\displaystyle=16.59+438.426\times
P\left(\overline{\nu}_{\mu}\to\overline{\nu}_{e},L_{\rm NOvA},E_{\rm
NOvA}\right).$ (3.55)

Figure 2: Validation of our NOvA analysis: we present the number of expected
signal events in neutrino mode ($x$-axis) and antineutrino mode ($y$-axis),
for $\Delta m_{21}^{2}$, $\sin^{2}\theta_{13}$, and $\sin^{2}\theta_{12}$
values given in the upper-right corner of the figure. The four different
ellipses correspond to different combinations of ($\sin^{2}\theta_{23}$,
$\Delta m_{31}^{2}$) given in the legend where we vary $\delta_{\rm CP}$
between $-\pi$ and $\pi$. 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. NOvANu2020  for $\delta_{\rm CP}=\pm\pi$
(star), $-\pi/2$ (square), $0$ (circle), and $\pi/2$ (triangle).

Instead of presenting the expected event rates for NOvA as a function of
$\delta_{\rm 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
$\Delta m_{21}^{2}$, $\sin^{2}\theta_{12}$, and $\sin^{2}\theta_{13}$ as given
by the values in the top-right of the panel. Each ellipse is generated by
fixing $(\sin^{2}\theta_{23},\Delta m_{31}^{2})$ (given by the legend), and
varying $\delta_{\rm CP}$, while using Eqs. (3.54) and (3.55). The values of
$\sin^{2}\theta_{23}$ and $\Delta m_{31}^{2}$ for the combination correspond
to the upper octant of $\theta_{23}>1/2$, and the normal mass ordering $\Delta
m_{31}^{2}>0$ corresponds to the best-fit according to Ref. NOvANu2020 . Since
corresponding values for NOvA’s preferred values in the inverted mass ordering
and/or lower octant of $\theta_{23}$ are not provided, we simply choose values
such that $\sin^{2}\left(2\theta_{23}\right)$ and $\left\lvert\Delta
m_{32}\right\rvert^{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. NOvANu2020
, displaying how well our results agree with the official collaboration ones.
Finally, the black cross indicates the observed event rate of
$N_{\nu_{e}}^{\rm NOvA}=82$, $N_{\overline{\nu}_{e}}^{\rm NOvA}=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 $\nu_{\mu}$ and $\overline{\nu}_{\mu}$
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.53). Since NOvA measures this
to be slightly away from maximal, we include this as a measurement
$C_{\mu\mu}^{\rm Dis.}=0.99\pm 0.02$, which replicates the results from Ref.
Acero:2019ksn ; NOvANu2020  fairly well. When we analyze NOvA results, either
appearance data alone, disappearance data alone, or combined, we include its
measurement $\left\lvert\Delta m_{32}^{2}\right\rvert=\left(2.41\pm
0.07\right)\times 10^{-3}$ eV2 NOvANu2020 .

###  3.7 Sterile Neutrino Searches

When unitarity is not assumed, there could be additional zero-distance
effects. Given our current knowledge of the mass-squared-splittings $\Delta
m_{21}^{2}$ and $\Delta m_{31}^{2}$ and assuming that there are no additional
neutrinos, any experiment, including near detectors, that operates at
$L/E_{\nu}$ such that

$\frac{\Delta m_{21}^{2}L}{4E_{\nu}}\ll
1,\quad\quad\mathrm{and}\quad\quad\frac{\Delta m_{31}^{2}L}{4E_{\nu}}\ll 1,$
(3.56)

should see no neutrino oscillations, and therefore is sensitive to these zero-
distance effects. In this regime of $L/E_{\nu}$, oscillation probabilities
without assuming unitarity are

$\displaystyle P(\nu_{\alpha}\to\nu_{\alpha})$
$\displaystyle=1,\quad\mathrm{(Disappearance)},$ (3.57) $\displaystyle
P(\nu_{\alpha}\to\nu_{\beta})$ $\displaystyle=\frac{\left\lvert
t_{\alpha\beta}\right\rvert^{2}}{N_{\alpha}N_{\beta}},\quad(\alpha\neq\beta\
\mathrm{Appearance}).$ (3.58)

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 as

$\displaystyle\hat{P}(\nu_{\alpha}\to\nu_{\alpha})$
$\displaystyle=N_{\alpha}^{2},\quad\mathrm{(Disappearance)},$ (3.59)
$\displaystyle\hat{P}(\nu_{\alpha}\to\nu_{\beta})$ $\displaystyle=\left\lvert
t_{\alpha\beta}\right\rvert^{2},\quad(\alpha\neq\beta\ \mathrm{Appearance}).$
(3.60)

Sterile neutrino searches are typically carried out in the following way. If a
fourth neutrino exists with a mass-squared splitting $\Delta
m_{41}^{2}\gg\left\lvert\Delta m_{31}^{2}\right\rvert$, and an experiment
operates in the $L/E_{\nu}$ regime given by Eq. (3.56), then it can search for
oscillations given by the probabilities

$\displaystyle P(\nu_{\alpha}\to\nu_{\alpha})$
$\displaystyle=1-\sin^{2}\left(2\theta_{\alpha\alpha}\right)\sin^{2}\left(\frac{\Delta
m_{41}^{2}L}{4E_{\nu}}\right)=1-4\left\lvert U_{\mu
4}\right\rvert^{2}\left(1-\left\lvert U_{\mu
4}\right\rvert^{2}\right)\sin^{2}\left(\frac{\Delta
m_{41}^{2}L}{4E_{\nu}}\right),\ \mathrm{(Disappearance)}$ (3.61)
$\displaystyle P(\nu_{\alpha}\to\nu_{\beta})$
$\displaystyle=\sin^{2}\left(2\theta_{\alpha\beta}\right)\sin^{2}\left(\frac{\Delta
m_{41}^{2}L}{4E_{\nu}}\right)=4\left\lvert U_{\alpha
4}\right\rvert^{2}\left\lvert U_{\beta
4}\right\rvert^{2}\sin^{2}\left(\frac{\Delta m_{41}^{2}L}{4E_{\nu}}\right),\
(\alpha\neq\beta\ \mathrm{Appearance}).$ (3.62)

Limits or potential observations from these searches are presented in terms of
$\sin^{2}\left(2\theta_{\alpha\beta}\right)$ and $\Delta m_{41}^{2}$ (see, for
example, Refs. Dentler:2018sju ; Diaz:2019fwt ; Boser:2019rta ) For any
particular experiment situated at a baseline $L$ and measuring oscillations
for a specific energy range, the oscillations associated with $\Delta
m_{41}^{2}$ become very rapid (as a function of $E_{\nu}$) as $\Delta
m_{41}^{2}$ becomes larger and larger. Eventually, these oscillations average
out, and the term $\sin^{2}\left(\Delta m_{41}^{2}L/4E_{\nu}\right)\to 1/2$.
In this averaged-out regime, the oscillation probabilities become

$\displaystyle P(\nu_{\alpha}\to\nu_{\alpha})$ $\displaystyle\to
1-\frac{1}{2}\sin^{2}\left(2\theta_{\alpha\alpha}\right)=1-2\left\lvert U_{\mu
4}\right\rvert^{2}\left(1-\left\lvert U_{\mu
4}\right\rvert^{2}\right),\quad\mathrm{(Disappearance)},$ (3.63)
$\displaystyle P(\nu_{\alpha}\to\nu_{\beta})$
$\displaystyle\to\frac{1}{2}\sin^{2}\left(2\theta_{\alpha\beta}\right)=2\left\lvert
U_{\alpha 4}\right\rvert^{2}\left\lvert U_{\beta
4}\right\rvert^{2},\quad(\alpha\neq\beta\ \mathrm{Appearance}).$ (3.64)

These have the same lack of energy-dependence as searches for unitarity
violation. Equating Eq. (3.63) with Eq. (3.59) therefore allows us to map
constraints on $\sin^{2}\left(2\theta_{\alpha\alpha}\right)$ from sterile
neutrino searches in the averaged-out regime onto constraints of
$N_{\alpha}^{2}$, while comparing Eq. (3.64) with Eq. (3.60) allows us to map
$\sin^{2}\left(2\theta_{\alpha\beta}\right)$ onto $\left\lvert
t_{\alpha\beta}\right\rvert^{2}$ for $\alpha\neq\beta$. Table 2 summarizes the
null sterile neutrino searches included in our analysis: KARMEN, NOMAD,
CHORUS, and MINOS/MINOS+.

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
---|---|---|---
KARMEN $\overline{\nu}_{\mu}\to\overline{\nu}_{e}$ Armbruster:2002mp  | $1.8\times 10^{-3}$ | $\sin^{2}\left(2\theta_{\mu e}\right)$ | $\left\lvert t_{\mu e}\right\rvert^{2}$
NOMAD $\nu_{\mu}\to\nu_{e}$ Astier:2003gs  | $1.4\times 10^{-3}$ | $\sin^{2}\left(2\theta_{\mu e}\right)$ | $\left\lvert t_{\mu e}\right\rvert^{2}$
NOMAD $\nu_{e}\to\nu_{\tau}$ Astier:2001yj  | $1.5\times 10^{-2}$ | $\sin^{2}\left(2\theta_{e\tau}\right)$ | $\left\lvert t_{e\tau}\right\rvert^{2}$
NOMAD $\nu_{\mu}\to\nu_{\tau}$ Astier:2001yj  | $3.3\times 10^{-4}$ | $\sin^{2}\left(2\theta_{\mu\tau}\right)$ | $\left\lvert t_{\mu\tau}\right\rvert^{2}$
CHORUS $\nu_{\mu}\to\nu_{\tau}$ Eskut:2007rn  | $4.4\times 10^{-4}$ | $\sin^{2}\left(2\theta_{\mu\tau}\right)$ | $\left\lvert t_{\mu\tau}\right\rvert^{2}$
MINOS/MINOS+ Adamson:2017uda  $\nu_{\mu}\to\nu_{\mu}$ | $2.5\times 10^{-2}$ | $\sin^{2}\left(2\theta_{\mu\mu}\right)$ | $N_{\mu}^{2}$

We do not consider any sterile neutrino searches for
$P(\overline{\nu}_{e}\to\overline{\nu}_{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 Huber:2011wv ; Mueller:2011nm ;
Berryman:2019hme ; Berryman:2020agd , so these experiments do not place
robust, strong limits in the high-$\Delta m_{41}^{2}$ regime.

Both the Liquid Scintillator Neutrino Detector (LSND) Athanassopoulos:1996wc ;
Athanassopoulos:1997pv  and MiniBooNE Aguilar-Arevalo:2018gpe  experiments
have observed an excess of electron-like events in the presence of a beam that
is mostly $\nu_{\mu}$ (or $\overline{\nu}_{\mu}$), which can be interpreted as
a short-baseline oscillation with $P(\nu_{\mu}\to\nu_{e})\approx 2.6\times
10^{-3}$. A combined study of these two experiments favors
$P(\nu_{\mu}\to\nu_{e})\neq 0$ at roughly $6\sigma$. 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 $\left\lvert
t_{\mu e}\right\rvert^{2}\neq 0$ preference from LSND/MiniBooNE in Appendix E.

##  4 Future Experiments and Simulations

In this section, we describe the future experiments that we consider in our
analysis. Specifically, we focus on the IceCube Upgrade Ishihara:2019aao ,
JUNO An:2015jdp ; Abusleme:2020bzt , DUNE Abi:2020evt ; Abi:2020qib , and T2HK
Abe:2018uyc  experiments.

When simulating future data, we assume that the LMM is unitary, and consistent
with the best-fit-point of an analysis of current data with unitarity assumed.
When we analyze all current data assuming unitarity (using the PMNS
parameterization), we obtain the following best-fit point:
$\sin^{2}\theta_{12}=0.308$, $\sin^{2}\theta_{13}=0.02190$,
$\sin^{2}\theta_{23}=0.551$, $\delta_{\rm CP}=-2.78=200.4^{\circ}$, $\Delta
m_{21}^{2}=7.50\times 10^{-5}$ eV2, and $\Delta m_{31}^{2}=2.53\times 10^{-3}$
eV2.101010With the recent update of oscillation data from T2K T2KNu2020  and
NOvA NOvANu2020 , the preference for the normal mass ordering ($\Delta
m_{31}^{2}>0$) over the inverted mass ordering ($\Delta m_{31}^{2}<0$) has
diminished Kelly:2020fkv ; Esteban:2020cvm . 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. We translate
these values in to the values of $\lvert U_{\alpha k}\rvert^{2}$ and
$\phi_{\alpha k}$ and obtain

$\left\lvert U_{\rm LMM}\right\rvert^{2}=\left(\begin{array}[]{c c
c}0.677&0.302&0.022\\\ 0.083&0.378&0.534\\\
0.240&0.320&0.439\end{array}\right);\ \phi_{e2}=3.00,\ \phi_{e3}=-0.47,\
\phi_{\tau 2}=-0.24,\ \phi_{\tau 3}=2.97.$ (4.65)

###  4.1 IceCube Upgrade

The IceCube experiment is capable of detecting atmospheric neutrinos over a
broad range of energies, $1\ \mathrm{GeV}\lesssim E_{\nu}\lesssim 100\
\mathrm{GeV}$. By measuring track-like ($\nu_{\mu}$) and cascade-like
($\nu_{e}$ and $\nu_{\tau}$) events, IceCube is sensitive to effects of the
mass-squared splitting $\Delta m_{31}^{2}$ and oscillations associated with
that splitting. Of importance for constraining leptonic unitarity is IceCube’s
ability to constrain the normalization of $\nu_{\mu}\to\nu_{\tau}$ appearance
in its data sample. Currently, this measurement has a precision of $\sim$40%
Aartsen:2019tjl . However, in the coming years, $\mathcal{O}(10\%)$ precision
will be attainable by considering either eight years of IceCube DeepCore data
or one year of IceCube Upgrade data summer_blot_2020_3959546 . Given the
regime of oscillations that IceCube measures, for $\nu_{\mu}\to\nu_{\tau}$
appearance it is dominantly sensitive to $4\lvert U_{\mu 3}\rvert^{2}\lvert
U_{\tau 3}\rvert^{2}/N_{\mu}^{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.

###  4.2 JUNO

JUNO is an upcoming reactor-based neutrino experiment (scheduled to start
operation in 2022 JUNONu2020 ), where neutrinos from 10 different nuclear
reactors travel for baselines of $\simeq$53 km before reaching a 20 kiloton
detector (with 2 additional reactors at a baseline of $\simeq$200 km),
comprised of liquid scintillator. JUNO is primarily sensitive to the
atmospheric mass-squared splitting $\Delta m_{31}^{2}$, as well as being
sensitive to the solar mass-squared splitting $\Delta m_{21}^{2}$. Its design
goal is to measure $\Delta m_{31}^{2}$ 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 $\Delta_{21}\ll
1$ nor $\Delta_{31}\gg 1$, when considering JUNO we must use the full
oscillation probability in Eq. (B.92). The effects of the matter potential on
the oscillation probability $P(\overline{\nu}_{e}\to\overline{\nu}_{e})$ are
negligible for sensitivity studies on the solar sector mixing An:2015jdp ,
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.92).

We simulate the expected event rate assuming six years of data collection,
corresponding to a total of $1.2\times 10^{5}$ signal events from
$\overline{\nu}_{e}$ inverse beta-decay scattering in the JUNO detector
An:2015jdp ; Abrahao:2015rba . Following Ref. An:2015jdp , 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}\theta_{12}$ in the standard three-flavor oscillation scenario is
0.42%, compared to the official collaboration sensitivity of 0.54%.

JUNO plans on using the Taishan Antineutrino Observatory (TAO)
Abusleme:2020bzt  to constrain the reactor antineutrino spectrum to sub-
percent energy resolution. TAO will be situated 30 meters from the core of the
Taishan Nuclear Power Plant and will be able to observe roughly 2000 reactor
$\overline{\nu}_{e}$ interactions per day. Using TAO and the 20 kt JUNO
detector, the collaboration can perform a near-to-far detector ratio
measurement similar to Daya Bay (see Section 3.3). We also make the
conservative assumption that with a comparison of TAO measurements to reactor
antineutrino flux predictions, the collaboration can constrain the
normalization of its near detector measurement to allow for a $10\%$
measurement of $N_{e}$.111111Without such a near detector measurement, the
near-to-far detector ratio measurement by JUNO-TAO would be subject to an
overall rescaling degeneracy of the electron row elements $\lvert
U_{ek}\rvert^{2}$. Including a 10% measurement from TAO does not impact our
final joint analysis results because $N_{e}$ is already constrained (by
KamLAND, solar neutrino measurements, and Daya Bay) at the few-percent level.

JUNO will in principle be sensitive to each element of the electron row of the
LMM, $\left\lvert U_{e1}\right\rvert^{2}$, $\left\lvert
U_{e2}\right\rvert^{2}$, and $\left\lvert U_{e3}\right\rvert^{2}$. As we will
show in Section 5, JUNO will allow for more precise measurements on
$\left\lvert U_{e1}\right\rvert^{2}$ and $\left\lvert U_{e2}\right\rvert^{2}$
than current experiments (specifically measuring the combination $\left\lvert
U_{e2}\right\rvert^{2}/\left\lvert U_{e1}\right\rvert^{2}$ very precisely),
but will not improve on the precision of $\left\lvert U_{e3}\right\rvert^{2}$
achieved by Daya Bay.

###  4.3 DUNE

DUNE is a future beam-based neutrino experiment, scheduled to begin data
collection in the late 2020s. It consists of an $\mathcal{O}(\text{GeV})$
neutrino beam, consisting primarily of $\nu_{\mu}$ when operating in neutrino
mode and $\overline{\nu}_{\mu}$ when operating in anti-neutrino mode
LauraFluxes , and a 40-kton liquid argon far detector. The baseline length is
1300 km, meaning that DUNE operates near the regime $\Delta_{21}\ll 1$.
Nevertheless, it will be sensitive to the effects of $\Delta m_{21}^{2}$,
allowing for a precise measurement of the CP-violating phase
$\delta_{\text{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(\nu_{\mu}\to\nu_{e})$
(and its CP-conjugate); muon-neutrino disappearance/survival
$P(\nu_{\mu}\to\nu_{\mu})$ (and its CP-conjugate); and tau-neutrino
appearance, sensitive to $P(\nu_{\mu}\to\nu_{\tau})$ (and its CP-conjugate).
To simulate the expected event rates for these different channels, we employ
simulation code developed with Refs. Berryman:2015nua ; deGouvea:2015ndi ;
Berryman:2016szd ; deGouvea:2016pom ; deGouvea:2017yvn  for $\nu_{e}$
appearance and $\nu_{\mu}$ disappearance and Ref. deGouvea:2019ozk  for
$\nu_{\tau}$ 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 $\nu_{e}$ appearance and $\nu_{\mu}$ disappearance
channels, and 25% for $\nu_{\tau}$ 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 $\nu_{\mu}$ disappearance and $\nu_{e}$ appearance, as well as
OPERA for $\nu_{\tau}$ appearance.

Finally, DUNE is capable of improving on existing measurements of the 8B solar
neutrino flux using $\nu_{e}$ CC and elastic electron scattering
Capozzi:2018dat . We simplify the analysis by assuming that DUNE will be able
to measure $P_{ee}\equiv\lvert U_{e2}\rvert^{2}\left(\lvert
U_{e1}\rvert^{2}+\lvert U_{e2}\rvert^{2}\right)+\left\lvert
U_{e3}\right\rvert^{4}$ at the $3\%$ level, consistent with the more complete
analysis of Ref. Capozzi:2018dat .

###  4.4 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_{\nu}$ 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 $\nu_{\mu}\to\nu_{e}$
appearance and $\nu_{\mu}\to\nu_{\mu}$ disappearance, consistent with Refs.
Abe:2015zbg ; Abe:2018uyc , developed from Refs. Kelly:2017kch ;
deGouvea:2017yvn . 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.
Kelly:2017kch ; deGouvea:2017yvn , and signal and background normalization
uncertainties of 5%.

##  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, 5.3, and 5.4. Section 5.2 shows
how subsets of data contribute to the sensitivity in the electron and muon
rows. Section 5.3 presents how well we can currently constrain the LMM matrix
element magnitudes $\left\lvert U_{\alpha k}\right\rvert^{2}$ and how much we
can improve with future data. Section 5.4 demonstrates how well we can
determine the normalization of the rows/columns, and closure of the row/column
triangles of the LMM.

Before presenting our results, we clarify the statistical approaches taken in
our analyses. In Sections 5.1 and 5.2, 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.3 and 5.4, the analyses
depend on 15 parameters, and we use the Bayesian inference tool pyMultiNest
Feroz:2007kg ; Feroz:2008xx ; Feroz:2013hea ; Buchner:2014nha  to construct
credible regions (CR) based on the posterior likelihood density. For further
details, see Appendix C.

###  5.1 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}\theta_{13}$ vs. $\sin^{2}\theta_{12}$ and $\sin^{2}\theta_{13}$ vs.
$\sin^{2}\theta_{23}$.121212For the third combination, $\sin^{2}\theta_{12}$
vs. $\sin^{2}\theta_{23}$, no existing or future measurement is sensitive to
this combination of angles in an interesting and non-trivial way. 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}\theta_{jk}$-$\sin^{2}\theta_{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 $\lvert U_{e3}\rvert^{2}$, or
$\sin^{2}\theta_{13}$ if unitarity is assumed. In what follows, when analyzing
different oscillation coefficients, we will generically allow
$\sin^{2}\theta_{13}$ to be large. This very likely produces results at large
$\sin^{2}\theta_{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}\theta_{13}\gtrsim 0.2$, with dashed lines.

Figure 3: Hypothetical perfect measurements of different oscillation
coefficients as described in the test – the true values of the mixing angles
are assumed to be $\sin^{2}\theta_{12}=0.316$, $\sin^{2}\theta_{13}=0.022$,
and $\sin^{2}\theta_{23}=0.565$. We dash lines above $\sin^{2}\theta_{13}$
where our analytical approximations are likely no longer valid.

In the left panel of Fig. 3, we compare how three types of measurements can
pin down $\sin^{2}\theta_{13}$ vs. $\sin^{2}\theta_{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_{ee}^{\rm SBL}$: short-baseline measurement of
$P(\overline{\nu}_{e}\to\overline{\nu}_{e})$ from a reactor neutrino
experiment, e.g., Daya Bay.

  * •

$P_{ee}^{\rm Solar}$: solar neutrinos measured via CC interactions.

  * •

$P_{ee}^{\rm LBL}$: long-baseline measurement of
$P(\overline{\nu}_{e}\to\overline{\nu}_{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}\theta_{13}$ and $\sin^{2}\theta_{23}$. They are the following:

  * •

$P_{ee}^{\rm SBL}$ as above.

  * •

$P_{\mu\mu}^{\rm LBL}$: long-baseline measurement of
$P(\nu_{\mu}\to\nu_{\mu})$ as performed by the current T2K and NOvA, or future
DUNE and T2HK experiments.

  * •

$P_{\mu e}^{\rm LBL}$: long-baseline measurement of $P(\nu_{\mu}\to\nu_{e})$
as performed by the current T2K/NOvA or future DUNE/T2HK experiments.

  * •

$P_{\mu\tau}^{\rm LBL}$: long-baseline measurement of
$P(\nu_{\mu}\to\nu_{\tau})$ 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}\theta_{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}\theta_{13}$ is needed to break the degeneracy. Also,
we see that an infinitely precise measurement of the $\nu_{\tau}$ 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 $\left\lvert
U_{e3}\right\rvert^{2}$ can be determined fairly robustly. Next, long-baseline
$\nu_{\mu}\to\nu_{e}$ appearance allows us to measure $4\lvert
U_{e3}\rvert^{2}\lvert U_{\mu 3}\rvert^{2}/N_{\mu}^{2}$. Meanwhile, long-
baseline $\nu_{\mu}\to\nu_{\mu}$ disappearance measures $4\lvert U_{\mu
3}\rvert^{2}(\lvert U_{\mu 1}\rvert^{2}+\lvert U_{\mu
2}\rvert^{2})/N_{\mu}^{2}$. If additional information (for instance from
MINOS/MINOS+) on $N_{\mu}^{2}$ is obtained, this combination allows us to
determine $\left\lvert U_{\mu 3}\right\rvert^{2}$ precisely. Finally, long-
baseline $\nu_{\mu}\to\nu_{\tau}$ appearance is sensitive to $4\lvert U_{\mu
3}\rvert^{2}\lvert U_{\tau 3}\rvert^{2}/N_{\mu}^{2}$, from which we can
extract $\lvert U_{\tau 3}\rvert^{2}$. In tandem then, we can measure each of
the elements $\lvert U_{\alpha 3}\rvert^{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 $\nu_{\mu}\to\nu_{\tau}$
appearance, as it is the least well-constrained of these measurements
currently. See the discussion around Fig. 6 for an illustration of this.

Figure 4: Current (left) and projected (right) measurements of the mixing
angles $\sin^{2}\theta_{13}$ and $\sin^{2}\theta_{12}$ at 95% and 99% CL. The
black contours in both panels show the joint-fit region with current data.
Figure 5: Current (left) and projected (right) measurements of the mixing
angles $\sin^{2}\theta_{23}$ and $\sin^{2}\theta_{13}$ at 95% and 99% CL. The
black contours in both panels show the joint-fit region with current data.

We now analyze how well the combination $\sin^{2}\theta_{13}$ vs.
$\sin^{2}\theta_{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}\theta_{12}$. However, it will not measure
$\sin^{2}\theta_{13}$ as precisely as Daya Bay. If we allow
$\sin^{2}\theta_{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}\theta_{13}$ vs. $\sin^{2}\theta_{23}$. In the left panel, because
the T2K and NOvA experiments both measure $\nu_{\mu}\to\nu_{\mu}$
disappearance and $\nu_{\mu}\to\nu_{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 $\delta_{\rm 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\sigma$ level. NOvA and T2K both see modestly higher event rates than
expected in $\nu_{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}\neq 1$. Here, we inject
unitarity violation by making $\lvert U_{\tau 3}\rvert^{2}$ significantly
larger than it should be, making $N_{3}\approx 2$, beyond the edge of the
currently allowed 3$\sigma$ range. We then fit the simulated data using DUNE’s
$\nu_{\tau}$ 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_{\mu\tau}$), orange (DUNE/T2HK $P_{\mu
e}$), and red (DUNE/T2HK $P_{\mu\mu}$) measurements – indicating the need for
a more thorough test of unitarity violation. If such a discrepancy arises,
IceCube’s measurement of $4\lvert U_{\mu 3}\rvert^{2}\lvert U_{\tau
3}\rvert^{2}/N_{\mu}^{2}$ could allow for verification of this arising from
unitarity violation.

Figure 6: Projected measurements of $\sin^{2}\theta_{13}$ vs.
$\sin^{2}\theta_{23}$ when unitarity is violated ($N_{3}\approx 2$). For
DUNE’s long-baseline measurement of $P_{\mu\tau}$ (green), we simulate data
assuming the underlying mixing matrix is non-unitary, and extract the
measurement of these parameters assuming the matrix is unitary.

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.
Ellis:2020ehi , where it was demonstrated how unitarity triangles
$\rho_{xy}+i\eta_{xy}$ can be used. In Ref. Ellis:2020ehi , 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.

###  5.2 Measurements of the Electron and Muon Rows

Before looking at the global fit results of the LMM elements, we show the
results for the absolute-value-squared 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 $\lvert U_{\alpha
k}\rvert^{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_{\tau\tau}$ has never been measured.

Electron Row Only: Figure 7 (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 $\Delta\chi^{2}$ for these CL, after marginalizing over the
third, unseen parameter.

Figure 7: Left: Current measurements coming from experiments that are only
sensitive to electron-row parameters $\left\lvert U_{e1}\right\rvert^{2}$,
$\left\lvert U_{e2}\right\rvert^{2}$, and $\left\lvert
U_{e3}\right\rvert^{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
$\Delta\chi^{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\left\lvert U_{e3}\right\rvert^{2}(\left\lvert
U_{e1}\right\rvert^{2}+\left\lvert U_{e2}\right\rvert^{2})/N_{e}^{2}$, which
is degenerate under $\left\lvert
U_{e1}\right\rvert^{2}\leftrightarrow\left\lvert U_{e2}\right\rvert^{2}$, such
that it appears that Daya Bay places no constraint in the upper-left panel of
Fig. 7 (left), showing constraints in the $\left\lvert
U_{e1}\right\rvert^{2}-\left\lvert U_{e2}\right\rvert^{2}$ plane. However,
Daya Bay constraints on $\left\lvert U_{e3}\right\rvert^{2}$ as a function of
$\left\lvert U_{e1}\right\rvert^{2},~{}\left\lvert U_{e2}\right\rvert^{2}$
combine with solar measurements in the $\left\lvert
U_{e2}\right\rvert^{2}-\left\lvert U_{e3}\right\rvert^{2}$ plane to bound
$\left\lvert U_{e2}\right\rvert^{2}$ and $\left\lvert U_{e1}\right\rvert^{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\leqslant|U_{ek}|^{2}\leqslant 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 upper-right triangular half of the $\lvert U_{e1}\rvert^{2}$
vs. $\lvert U_{e2}\rvert^{2}$ panel would be forbidden, somewhat limiting the
KamLAND, solar CC, and joint fit contours. Compared to Fig. 1 in Ref.
Antusch:2006vwa , the $U_{e1}$-$U_{e2}$ panel is similar, but with some
differences due to how data sets are combined. The $U_{e3}$ panel is
qualitatively different due to the addition of the Daya Bay measurement. The
stars in Fig. 7 represent the best-fit point in this parameter space given by
this combination of datasets: $\left\lvert U_{e1}\right\rvert^{2}=0.670$,
$\left\lvert U_{e2}\right\rvert^{2}=0.306$, and $\left\lvert
U_{e3}\right\rvert^{2}=0.022$.

The upcoming JUNO experiment will measure a combination of these parameters as
well, with the oscillation probability given by Eq. (B.92). JUNO will operate
in the regime of $L/E_{\nu}$ in which the effects of both mass-squared
splittings are relevant, and therefore, with enough statistical power, can be
independently sensitive to each mixing element $|U_{ek}|^{2}$. JUNO’s capacity
to measure each of these three parameters is depicted in Fig. 7, alongside the
combined current measurement from Fig. 7 (left) in black. JUNO will precisely
measure the combination $\left\lvert U_{e2}\right\rvert^{2}/\left\lvert
U_{e1}\right\rvert^{2}$ leading to the very sharp regions in the $\left\lvert
U_{e1}\right\rvert^{2}$ vs. $\left\lvert U_{e2}\right\rvert^{2}$ panel. When
combined with the current fit region, this will lead to impressive
measurements of both $\left\lvert U_{e1}\right\rvert^{2}$ and $\left\lvert
U_{e2}\right\rvert^{2}$.

On its own, JUNO suffers the same degeneracy discussed above for Daya Bay
under the interchange $\lvert U_{e1}\rvert^{2}\leftrightarrow\lvert
U_{e2}\rvert^{2}$, and requires solar neutrino experiments to break the
degeneracy. However, this interchange $\lvert
U_{e1}\rvert^{2}\leftrightarrow\lvert U_{e2}\rvert^{2}$ also requires changing
the neutrino mass ordering (or the sign of $\Delta m_{31}^{2}$). If the mass
ordering can be determined independently of JUNO at high enough significance
(for instance, by DUNE, T2HK) then the solution where $\lvert
U_{e2}\rvert^{2}>\lvert U_{e1}\rvert^{2}$ may be eliminated. The availability
of redundant but independent data to select the right $(\left\lvert
U_{e1}\right\rvert^{2},\left\lvert U_{e2}\right\rvert^{2})$ is a powerful tool
to test new physics scenarios such as possible non-standard interactions of
neutrinos.

In Fig. 7 (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 $\lvert U_{e1}\rvert^{2}$ and $\lvert U_{e2}\rvert^{2}$ will
be dominated by JUNO, whereas the measurements of $\lvert U_{e3}\rvert^{2}$
will be dominated by current experiments, specifically Daya Bay. JUNO will
measure $\lvert U_{e3}\rvert^{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.

Figure 8: Left: Current measurements coming from experiments only sensitive to
the muon-row parameters $\left\lvert U_{\mu 3}\right\rvert^{2}$ and the
combinations $\left\lvert U_{\mu 1}\right\rvert^{2}\pm\left\lvert U_{\mu
2}\right\rvert^{2}$. Here we compare measurements from T2K $\nu_{\mu}$ and
$\overline{\nu}_{\mu}$ disappearance, NOvA $\nu_{\mu}$ and
$\overline{\nu}_{\mu}$ 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 $\nu_{\mu}$ and
$\overline{\nu}_{\mu}$ measurements of these parameters.

Muon Row Only: We perform a similar procedure focusing on the elements
$|U_{\mu k}|^{2}$ in Fig. 8. In the left panel of Fig. 8, we include the
MINOS/MINOS+ sterile neutrino search, NOvA $\nu_{\mu}$ and
$\overline{\nu}_{\mu}$ disappearance, and T2K $\nu_{\mu}$ and
$\overline{\nu}_{\mu}$ disappearance. The combined fit of these data sets is
shown in black. Instead of presenting these measurements in terms of the
elements $\lvert U_{\mu 1}\rvert^{2}$, $\lvert U_{\mu 2}\rvert^{2}$, and
$\lvert U_{\mu 3}\rvert^{2}$, we present the measurements in terms of the
combinations $\lvert U_{\mu 1}\rvert^{2}\pm\lvert U_{\mu 2}\rvert^{2}$ and
$\lvert U_{\mu 3}\rvert^{2}$. This is because the combination $\lvert U_{\mu
1}\rvert^{2}+\lvert U_{\mu 2}\rvert^{2}$ is more precisely measured by
MINOS/MINOS+, NOvA, and T2K, where their difference is not well-constrained by
current data. Again, if we compare against the analogous parameter space in
Ref. Antusch:2006vwa , we find overall consistency, with the stronger
constraints from T2K, NOvA, and MINOS+ contributing to stronger measurements
in the $\lvert U_{\mu 1}\rvert^{2}+\lvert U_{\mu 2}\rvert^{2}$ vs. $\lvert
U_{\mu 3}\rvert^{2}$ space than those in Ref. Antusch:2006vwa .

Again, parameter degeneracies cause the combined fit to appear stronger than
naïve expectations. We see that $\left\lvert U_{\mu 3}\right\rvert^{2}$ is
measured the most precisely of these three parameters, and the combination
$\lvert U_{\mu 1}\rvert^{2}+\lvert U_{\mu 2}\rvert^{2}$ is well measured by
the combination of datasets in the bottom-left panel. The difference $\lvert
U_{\mu 1}\rvert^{2}-\lvert U_{\mu 2}\rvert^{2}$ is not well-measured, and is
predominantly constrained by the requirement that both $\lvert U_{\mu
1}\rvert^{2}$ and $\lvert U_{\mu 2}\rvert^{2}$ are both between $0$ and $1$.
This forces the difference $\lvert U_{\mu 1}\rvert^{2}-\lvert U_{\mu
2}\rvert^{2}$ to be less in magnitude than the sum $\lvert U_{\mu
1}\rvert^{2}+\lvert U_{\mu 2}\rvert^{2}$. The stars in Fig. 8 (left)
correspond to the best-fit points of these elements from the full analysis
discussed in the main text: $\left\lvert U_{\mu 1}\right\rvert^{2}=0.096$,
$\left\lvert U_{\mu 2}\right\rvert^{2}=0.352$, $\left\lvert U_{\mu
3}\right\rvert^{2}=0.552$.

Figure 8 (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) $\nu_{\mu}$ and $\overline{\nu}_{\mu}$
disappearance.131313In order to study how DUNE and T2HK muon-neutrino
disappearance are sensitive to the muon row elements, and only the muon row
elements, we perform this simulation assuming oscillations of $\nu_{\mu}$
occur in vacuum. This allows us to use the expression in Eq. (B.92) for our
calculations. We find this to be a good approximation for the muon
neutrino/antineutrino channels. DUNE and T2HK measure oscillations over a
wider range of $L/E_{\nu}$ 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}\left(\Delta_{31}/2\right)$ in
Eq. (B.92). Indeed, these experiments have some sensitivity to the
interference term, namely the final term of Eq. (B.92). This allows DUNE and
T2HK to be sensitive to $\lvert U_{\mu 1}\rvert^{2}-\lvert U_{\mu
2}\rvert^{2}$ unlike the current data considered in Fig. 8 (left). We see that
these experiments can both demonstrate $\lvert U_{\mu 2}\rvert^{2}>\lvert
U_{\mu 1}\rvert^{2}$ at high significance.

###  5.3 Joint Measurement of All Matrix-Elements-Squared

In this subsection, we present the current constraints on the
parameterization-independent $\left\lvert U_{\alpha k}\right\rvert^{2}$, and
project how well these will be constrained once future data from DUNE/T2HK,
JUNO, and IceCube Upgrade are included. In Appendix D, we show how well the
parameterization-dependent phases $\left\\{\phi_{e2},\phi_{e3},\phi_{\tau
2},\phi_{\tau 3}\right\\}$ are and will be constrained.

Figure 9: Measurements (and projections) of the matrix-elements-squared
$\left\lvert U_{\alpha k}\right\rvert^{2}$ of the LMM. These are one-
dimensional $\Delta\chi^{2}$ measurements of each parameter after
marginalizing a 15 parameter fit down to each individual one. Results are
generated using pyMultiNest Feroz:2007kg ; Feroz:2008xx ; Feroz:2013hea ;
Buchner:2014nha . All results here are obtained under the agnostic assumption.

Figure 9 displays the individual measurements of $\left\lvert U_{\alpha
k}\right\rvert^{2}$ including current data (blue) and current data with future
data (red).141414Note that we present the top-right axis in terms of
$10\times\lvert U_{e3}\rvert^{2}$ rather than $\lvert U_{e3}\rvert^{2}$ for
presentation purposes so that it can share axes with the $\lvert U_{\mu
3}\rvert^{2}$ and $\lvert U_{\tau 3}\rvert^{2}$ panels. Each panel displays
the one-dimensional $\Delta\chi^{2}$ measurement of the element after
marginalizing the 15-parameter fit down to the individual element. Here, we
define $\Delta\chi^{2}=-2\Delta\mathcal{L}$, where $\mathcal{L}$ is the
posterior likelihood obtained in our analysis. The improvement in $\left\lvert
U_{e1}\right\rvert^{2}$ and $\left\lvert U_{e2}\right\rvert^{2}$ is driven
predominantly by JUNO, while the improvement in the muon row is driven by DUNE
and T2HK $\nu_{\mu}\to\nu_{\mu}$ measurements.The tau row improvements are
driven by $\nu_{\mu}\to\nu_{\tau}$ 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\sigma$ 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

$\displaystyle\left\lvert U_{\rm LMM}^{\rm Current}\right\rvert^{2}_{3\sigma}$
$\displaystyle=\left(\begin{array}[]{c c
c}\left[0.606,0.742\right]&\left[0.265,0.337\right]&\left[0.020,0.024\right]\\\
\left[0.051,0.270\right]&\left[0.198,0.484\right]&\left[0.392,0.620\right]\\\
\left[0.028,0.469\right]&\left[0.098,0.685\right]&\left[0.140,0.929\right]\end{array}\right),$
(5.69) $\displaystyle\left\lvert U_{\rm LMM}^{\rm
Future}\right\rvert^{2}_{3\sigma}$ $\displaystyle=\left(\begin{array}[]{c c
c}\left[0.653,0.699\right]&\left[0.291,0.311\right]&\left[0.020,0.024\right]\\\
\left[0.074,0.108\right]&\left[0.355,0.454\right]&\left[0.447,0.561\right]\\\
\left[0.129,0.359\right]&\left[0.212,0.423\right]&\left[0.349,0.595\right]\end{array}\right).$
(5.73)

We can interpret the amount by which each of these measurements will improve,
at the $3\sigma$ level, by computing the reduction in size of the allowed
$3\sigma$ range of each $\lvert U_{\alpha k}\rvert^{2}$, $\Delta\lvert
U_{\alpha k}^{\rm Current}\rvert^{2}/\Delta\lvert U_{\alpha k}^{\rm
Future}\rvert^{2}$:

$\displaystyle\mathrm{Improvement\ Factor:}\ \left(\begin{array}[]{c c
c}3.0&3.6&1.0\\\ 6.4&2.9&2.0\\\ 1.9&2.8&3.2\end{array}\right).$ (5.77)

As is evident in Fig. 9, the improvement is noticeable especially for the
elements $\lvert U_{e1}\rvert^{2}$, $\lvert U_{e2}\rvert^{2}$, $\lvert U_{\mu
1}\rvert^{2}$, and the $\tau$ 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.

###  5.4 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_{\alpha}$ and $N_{k}$ and triangle closures
$t_{\alpha\beta}$ (between two rows) and $t_{kl}$ (between two columns), using
the same analyses as in the previous subsection.

Figure 10: Left: constraints (and projected constraints) on the row
normalizations $N_{e}$, $N_{\mu}$, and $N_{\tau}$ 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.

The left panel of Fig. 10 displays the results of this analysis, projecting
down to two-dimensional CR measuring the row normalizations $N_{e}$,
$N_{\mu}$, and $N_{\tau}$ 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_{\mu}$, some
improvement in $N_{e}$, and significant improvement in $N_{\tau}$.

Similarly, the right panel of Fig. 10 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,
$\left\lvert U_{\tau k}\right\rvert^{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\sigma$ CR. We also show the projected future $3\sigma$ CR, assuming
a true value of $N_{X}=1$, demonstrating the improvement attributable to
future data. Our projected constraint on $N_{e}$ is 1.1%, consistent with the
official JUNO analysis, which reports a 1.2% constraint on $N_{e}$.

Table 3: Summary of current and expected future constraints on the row ($N_{\alpha}$) and column ($N_{k}$) normalizations, under the agnostic assumption.  | Best-fit (current) |  $3\sigma$ (current) |  $3\sigma$ (future)
---|---|---|---
$N_{e}$ | $1.00$ | $[0.94,1.05]$ | $[0.97,1.03]$
$N_{\mu}$ | $0.99$ | $[0.96,1.04]$ | $[0.96,1.03]$
$N_{\tau}$ |  $1.12$. | $[0.32,1.82]$ | $[0.79,1.23]$
$N_{1}$ | $1.01$ | $[0.84,1.22]$ | $[0.89,1.12]$
$N_{2}$ | $1.05$ | $[0.75,1.27]$ | $[0.92,1.10]$
$N_{3}$ | $1.05$ | $[0.67,1.40]$ | $[0.90,1.10]$
Figure 11: Constraints (and projected constraints) on the real ($x$-axes) and
imaginary ($y$-axes) parts of the closures of the six unitarity triangles.
Dashed circles indicate contours corresponding to fixed $\lvert
t_{xy}\rvert^{2}$, with the outer one in each panel as labeled. The inner
dashed circles are an order of magnitude smaller $\lvert t_{xy}\rvert^{2}$
than the outer ones. Here we analyze data under the agnostic assumption.

Figure 11 presents the results on the closures of different triangles
$t_{\alpha\beta}$ and $t_{kl}$. Each panel in this figure presents constraints
on the real and imaginary part of $t_{\alpha\beta}$ (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 $\left\lvert
t_{\alpha\beta}\right\rvert^{2}$ and $\left\lvert t_{kl}\right\rvert^{2}$ as
labeled, where each successive inward circle is an order of magnitude smaller.
When constraining $t_{e\tau}$ and $t_{\mu\tau}$, 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_{\mu 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 $\left\lvert t_{kl}\right\rvert^{2}\lesssim 10^{-1}$
at $95\%$ credibility. With future data, this will improve to roughly
$\left\lvert t_{kl}\right\rvert\lesssim 10^{-2}$ for each of the three
triangles. We summarize the current and future $3\sigma$ credibility upper
limits on the triangle closures in Table 4.

Table 4: Summary of current and expected future constraints on the row closures $\left\lvert t_{\alpha\beta}\right\rvert$ and column closures $\left\lvert t_{kl}\right\rvert$, under the agnostic assumption.  | Current $3\sigma$ Upper Limit | Future $3\sigma$ Upper Limit
---|---|---
$\left\lvert t_{e\mu}\right\rvert$ | $3.2\times 10^{-2}$ | $2.5\times 10^{-2}$
$\left\lvert t_{e\tau}\right\rvert$ | $1.3\times 10^{-1}$ | No Improvement
$\left\lvert t_{\mu\tau}\right\rvert$ | $1.6\times 10^{-2}$ | No Improvement
$\left\lvert t_{12}\right\rvert$ | $2.5\times 10^{-1}$ | $1.0\times 10^{-1}$
$\left\lvert t_{13}\right\rvert$ | $3.2\times 10^{-1}$ | $1.2\times 10^{-1}$
$\left\lvert t_{23}\right\rvert$ | $3.3\times 10^{-1}$ | $1.1\times 10^{-1}$

The analysis yielding Figs. 10 and 11 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. Parke:2015goa , 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.

##  6 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.

###  6.1 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. 11 without short-baseline results.

Figure 12: Constraints (and projected constraints) on the closures of the six
unitarity triangles excluding data from short-baseline sterile neutrino
searches. In each panel, the outer (inner) dashed circle corresponds to
constant $\lvert t_{\alpha\beta}\rvert^{2}=10^{-1}$ ($10^{-2}$) for the
row/column closures.

Figure 12 shows the results. Here we note that the ranges on each of the
panels in Fig. 12 measuring $t_{\alpha\beta}$ and $t_{kl}$ are much larger
than the corresponding ranges in Fig. 11. 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. 11, we draw lines
of constant $\lvert t_{\alpha\beta}\rvert^{2}$ and $\lvert
t_{kl}\rvert^{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.

Table 5 summarizes the numerical results. Comparing 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 Berryman:2015nua ; Kelly:2017kch ; Machado:2019oxb ;
Ghosh:2019zvl  (as discussed in Section 3.7) than those summarized in Table 2.

Comparing the measurements of the individual matrix-elements-squared $\lvert
U_{\alpha k}\rvert^{2}$, as well as the row and column normalization
conditions $N_{\alpha}$ and $N_{k}$, is difficult in this scenario. Short-
baseline sterile neutrino searches, particularly the information from
$\nu_{\tau}$ appearance that $\lvert t_{e\tau}\rvert^{2}$ and $\lvert
t_{\mu\tau}\rvert^{2}$ are small, provide significant information on the
elements $\lvert U_{\tau k}\rvert^{2}$. Additionally, the constraint from
MINOS/MINOS+ that $N_{\mu}\approx 1$ is very important for determining the
muon elements $\lvert U_{\mu k}\rvert^{2}$. If this information is discarded,
every other probe of $\lvert U_{\mu k}\rvert^{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+
$\nu_{\mu}$ disappearance, for precise tests of leptonic unitarity. This
analysis without short-baseline measurements results in a lower limit on
$N_{\mu}$ comparable to the one given in Table 3, however, $N_{\mu}$ can be as
large as $\approx 2$.

Table 5: Summary of current and expected future constraints on the row closures $\left\lvert t_{\alpha\beta}\right\rvert$ and column closures $\left\lvert t_{kl}\right\rvert$, under the agnostic case regarding the LMM, when short-baseline sterile neutrino search results are not included.  | Current $3\sigma$ Upper Limit | Future $3\sigma$ Upper Limit
---|---|---
$\left\lvert t_{e\mu}\right\rvert$ | $8.2\times 10^{-1}$ | $3.0\times 10^{-2}$
$\left\lvert t_{e\tau}\right\rvert$ | $1.0\times 10^{0}$ | $8.1\times 10^{-1}$
$\left\lvert t_{\mu\tau}\right\rvert$ | $6.8\times 10^{-1}$ | $3.2\times 10^{-1}$
$\left\lvert t_{12}\right\rvert$ | $8.2\times 10^{-1}$ | $5.5\times 10^{-1}$
$\left\lvert t_{13}\right\rvert$ | $1.1\times 10^{0}$ | $5.2\times 10^{-1}$
$\left\lvert t_{23}\right\rvert$ | $8.7\times 10^{-1}$ | $2.7\times 10^{-1}$

###  6.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_{\rm LMM}$, specifically requiring

$\displaystyle N_{\alpha}$ $\displaystyle\leqslant 1,\ \quad N_{k}\leqslant
1,$ (6.78) $\displaystyle\lvert t_{\alpha\beta}\rvert^{2}$
$\displaystyle\leqslant\left(1-N_{\alpha}\right)\left(1-N_{\beta}\right),$
(6.79) $\displaystyle\lvert t_{kl}\rvert^{2}$
$\displaystyle\leqslant\left(1-N_{k}\right)\left(1-N_{l}\right).$ (6.80)

For our main analyses in Sections 5.3 and 5.4, 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. 9 for
the current data analyzed, under the two different case assumptions. The
results of this procedure are shown in Fig. 13.

Figure 13: Measurements (and projections) of $\left\lvert U_{\alpha
k}\right\rvert^{2}$ with current experimental data under the agnostic
assumption (blue) and the sub-matrix case (green). Similar to Fig. 9, these
are one-dimensional $\Delta\chi^{2}$ measurements of each parameter after
marginalizing a 15 parameter fit down to each individual one.

We note several effects of the sub-matrix case in Fig. 13. 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, as discussed in Sec. 5.2. We note that the measurement $\lvert
U_{e1}\rvert^{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_{\mu}\leqslant 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 to
the mild excess over expectation of OPERA $\nu_{\tau}$ appearance events, the
current data prefer $N_{\tau}>1$ (at very low significance). When we analyze
the data under the sub-matrix assumption, this solution is forbidden, and the
$\lvert U_{\tau k}\rvert^{2}$ elements are both preferred to be lower in
magnitude, and we end up with smaller measurement regions.

In Section 5.4, 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_{\alpha}$ and $N_{k}$. If we were to repeat this analysis under
the sub-matrix assumption, the resulting figure analogous to Fig. 10 would
look very similar, modulo the regions $N_{\alpha}>1$ and $N_{k}>1$ being
forbidden.

For completeness, Fig. 14 presents the projected future measurements of
$\lvert U_{\alpha k}\rvert^{2}$ under the agnostic (red) and sub-matrix
(purple) assumptions. The red lines are identical to those in Fig. 9. We see
that the sub-matrix hypothesis improves constraints on $\lvert U_{\mu
k}\rvert^{2}$ somewhat, and $\lvert U_{\tau k}\rvert^{2}$ significantly.

Figure 14: Projections of future measurements of $\left\lvert U_{\alpha
k}\right\rvert^{2}$ under the agnostic case assumption (red) and the sub-
matrix case (purple). Similar to Figs. 9 and 13, these are one-dimensional
$\Delta\chi^{2}$ measurements of each parameter after marginalizing a 15
parameter fit down to each individual one.

Finally, we repeat the procedure that generated Fig. 11, 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. 15.

Figure 15: Constraints (and projected constraints) on the closures of the six unitarity triangles with data analyzed under the agnostic (blue) and sub-matrix (green) assumptions. Dashed circles indicate contours corresponding to fixed $\lvert t_{xy}\rvert^{2}$, with the outer one in each panel as labeled. The inner dashed circles are an order of magnitude smaller $\lvert t_{xy}\rvert^{2}$ than the outer ones.  Table 6: Summary of current and expected future constraints on the row closures $\left\lvert t_{\alpha\beta}\right\rvert$ and column closures $\left\lvert t_{kl}\right\rvert$, under the sub-matrix assumption.  | Current $3\sigma$ Upper Limit | Future $3\sigma$ Upper Limit
---|---|---
$\left\lvert t_{e\mu}\right\rvert$ | $2.7\times 10^{-2}$ | $2.1\times 10^{-2}$
$\left\lvert t_{e\tau}\right\rvert$ | $1.0\times 10^{-1}$ | $6.5\times 10^{-2}$
$\left\lvert t_{\mu\tau}\right\rvert$ | $1.6\times 10^{-2}$ | No Improvement
$\left\lvert t_{12}\right\rvert$ | $2.2\times 10^{-1}$ | $8.0\times 10^{-2}$
$\left\lvert t_{13}\right\rvert$ | $2.5\times 10^{-1}$ | $1.0\times 10^{-1}$
$\left\lvert t_{23}\right\rvert$ | $3.3\times 10^{-1}$ | $1.0\times 10^{-1}$

Again, we note several features of this result. First, when looking at the
closure between two rows $\lvert t_{\alpha\beta}\rvert^{2}$ (the top three
panels of Fig. 15), we see that the resulting constraint on $\lvert
t_{\alpha\beta}\rvert^{2}$ is largely unchanged. This is because both analyses
include the searches for short-baseline neutrino appearance discussed in
Section 3.7 which directly constrain the closures $\lvert t_{e\mu}\rvert^{2}$
(from $\nu_{\mu}\to\nu_{e}$ appearance), $\lvert t_{e\tau}\rvert^{2}$ (from
$\nu_{e}\to\nu_{\tau}$ appearance), and $\lvert t_{\mu\tau}\rvert^{2}$ (from
$\nu_{\mu}\to\nu_{\tau}$ appearance). The mild improvement seen in each of
these panels comes from the Cauchy-Schwarz constraint $\lvert
t_{\alpha\beta}\rvert^{2}\leqslant(1-N_{\alpha})(1-N_{\beta})$, where the
normalization constraints assist in these planes. For the bottom panels, the
closure of triangles formed between two columns $\lvert t_{kl}\rvert^{2}$, the
difference between the agnostic and sub-matrix analyses is more drastic. Here,
the Cauchy-Schwarz constraints are of the form $\lvert
t_{kl}\rvert^{2}\leqslant(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.
Parke:2015goa , 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. 16 presents projected future constraints on $\lvert
t_{\alpha\beta}\rvert^{2}$ and $\lvert t_{kl}\rvert^{2}$ under these two
hypotheses, similar to Fig. 15. 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 $\left\lvert t_{\alpha\beta}\right\rvert$ and $\left\lvert
t_{kl}\right\rvert$ 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_{\alpha},N_{k}>1$ in the fits. We do not
reproduce an analogous version of Table 3 here – the same $3\sigma$ lower
limits on $N_{\alpha}$ and $N_{k}$ apply here, where the upper limits are now
simply $1$ from the theoretical constraints of the sub-matrix case.

Figure 16: Projected future constraints on the closures of the six different
unitarity triangles with data analyzed under the agnostic (red) and sub-matrix
(purple) assumptions. Dashed circles indicate contours corresponding to fixed
$\lvert t_{xy}\rvert^{2}$, with the outer one in each panel as labeled. The
inner dashed circles are an order of magnitude smaller $\lvert
t_{xy}\rvert^{2}$ than the outer ones.

##  7 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 mass-squared 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. 1.

The “standard” case, with only 3 neutrino flavors, and the LMM is the PMNS
matrix, and therefore unitary.

  2. 2.

The “sub-matrix” case, where there are $n>3$ neutrino flavors, and the
$3\times 3$ LMM is non-unitary, but the greater $n\times n$ matrix is unitary.

  3. 3.

The “agnostic” case where we assume nothing regarding the unitarity of a
possibly $m\times 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