Global oscillation data analysis on the $3\nu$ mixing without unitarity

We present results of a combined analysis in neutrino oscillations without unitarity assumption in the $3\nu$ mixing picture. Constraints on neutrino mixing matrix elements are based on recent data from the reactor, solar and long-baseline accelerator neutrino oscillation experiments. The current data are consistent with the standard $3\nu$ scheme. The precision on different matrix elements can be as good as a few percent at $3\sigma$ CL, and is mainly limited by the experimental statistical uncertainty. The $\nu_e$ related elements are the most precisely measured among all sectors with the uncertainties $<20\%$. The measured leptonic CP violation is very close to the one assuming the standard $3\nu$ mixing. The deviations on normalization and the unitarity triangle closure are confined within $\mathcal{O}(10^{-3})$, $\mathcal{O}(10^{-2})$ and $\mathcal{O}(10^{-1})$, for $\nu_e$, $\nu_{\mu}$ and $\nu_{\tau}$ sectors, respectively. We look forward to the next-generation neutrino oscillation experiments \textit{such as} DUNE, T2HK, and JUNO, especially the precise measurements on $\nu_\tau$ oscillations, to significantly improve the precision of unitarity test on the $3\nu$ mixing matrix.


Introduction
The observation of neutrino oscillation phenomena [1][2][3], which can be described by neutrino masses and the leptonic mixing, points to new physics beyond the Standard Model (SM). The leptonic mixing can be incorporated into the SM by involving the neutrino mixing matrix U so that neutrino flavour eigenstates are identified by mass eigenstates with the superposition |ν α = U αi |ν i , where U αi is an element of U . So far, ν e , ν µ and ν τ are three observed neutrino flavour eigenstates that participate in standard weak interactions. They are called 'active neutrinos'. Without any experimental evidence of the new neutrino states in the neutrino mixing, the matrix U is generally assumed to be a unitary 3×3 matrix U 3ν , and is commonly parameterized by the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [4][5][6] with three mixing angles θ 12 , θ 13 and θ 23 , a Dirac phase δ CP , and two Majorana phases. In principle, two Majorana phases never affect the neutrino oscillation formalism. Therefore, we will only talk about the Dirac CP phase in the following discussion.
In the standard 3ν mixing framework, the recent data from solar, atmospheric, reactor and long-baseline accelerator neutrino oscillation experiments can be combined to constrain these neutrino mixing parameters [7]. Two out of four parameters in the PMNS matrix have already been measured in a relatively high precision, rather than the other two. One is the parameter sin 2 θ 12 primarily determined by solar neutrino experiments like SNO and the long-baseline reactor experiment like KamLAND. The other is sin 2 θ 13 , reaching a precision of ∼ 3%, predominantly derived from the reactor neutrino experiments such as Daya Bay, Double Chooz and RENO, with baselines at ∼ 1 km. Apart from the current achievement, there are still puzzles in neutrino physics to be resolved: θ 23 -octant problem (θ 23 >, = or < 45 • ), the neutrino mass ordering problem (if ∆m 2 31 > or < 0), whether there is any CP violation in the leptonic sector and the size of δ CP if so. According to the current result, it seems that, by the combining atmospheric and accelerator neutrino oscillation data, we will have sin 2 θ 23 in the higher octant and the CP violation scenario δ CP = 0, π. Recent updates from running accelerator neutrino experiments, NOvA and T2K [8,9] have significantly improved the precision of these two parameters. Given the standard 3ν mixing, the current global fit results prefer the solution with sin 2 θ 12 = 0.307 +0.013 −0.012 , sin 2 θ 23 = 0.542 +0.019 −0.022 (normal ordering, upper octant) and sin 2 θ 13 = 0.0218 ± 0.0007. In addition, the remaining parameter δ CP shows an inclination to be non-zero, though the uncertainty ∆δ CP still needs to be reduced [10]. Entering a new era of precise measurements in neutrino oscillation physics, whether the 3ν mixing framework conserves unitarity deserves further scrutiny. Moreover, our understanding of neutrino mixing does not necessitate the foundation of unitary 3ν mixing assumption. There are theoretical extensions in the neutrino sector in the SM, which allow the non-unitarity of 3ν mixing, e.g. the sterile-active neutrino mixing [11][12][13][14][15][16][17][18][19]. The motivations for this new mixing are (1) to explain experimental anomalies [20][21][22], (2) to explain the origin of neutrino masses e.g., seesaw I [23][24][25][26] and seesaw III [27,28], and (3) to explain warm dark matter [29][30][31]. In addition, unknown couplings involving neutrinos may result in effective the 3ν non-unitarity [32][33][34]. As the precision gets better and better, one might be always questioning how much we know about the neutrino mixing if taking away the unitarity assumptions, and whether the unitarity assumption is valid. In the quark sector, similar issues have been discussed for decades, such as Refs. [35][36][37][38][39][40]. Some of the neutrino phenomenological studies answer these questions [41][42][43][44][45][46]. Along with Ref. [44], testing the unitarity hypothesis in a general manner, authors of current works [45,46] have updated the results with the current data, and their predicted results emphasis the prospect of future experiments: DUNE, Hyper, and JUNO, as well as the τ neutrino measurements.
In this paper, we present an analysis beyond the standard 3ν mixing scheme, using current neutrino oscillation data. In presence of sterile neutrinos, the neutrino mixing matrix is enlarged to a unitary (3 + N ) × (3 + N ) mixing matrix, where there are N types of sterile neutrinos. This results in the non-unitarity of the 3ν mixing matrix. The nonunitary 3ν mixing assumption is compared against data from accelerator, reactor and solar neutrino measurements. Our understanding of the 3ν mixing matrix is so far mostly limited to the ν e and ν µ sectors, by reactor and solar neutrino experiments and accelerator neutrino experiments, respectively. In our analysis, we will try to answer how the neutrino mixing looks like beyond the standard 3ν mixing scheme and how much the current data prefer the 3ν unitarity assumption. This paper is organized as follows. In section 2 we introduce the non-unitary 3ν mixing framework. Section 3 describes the oscillation data used in this analysis, as well as how we analyse the data. Section 4 presents the resultant non-unitary 3ν mixing parameters, their correlations and CP-violation, as well as the results of unitarity test. We present allowed parameter regions of the e − µ triangle in Fig. 4 without assuming unitarity. Unitarity conditions are tested and shown in Fig. 5. Finally, we summarize and conclude this work in section 5.

Non-unitarity and neutrino oscillation
The neutrino oscillation is a quantum coherent phenomenon described by the neutrino mixing and the mass-squared differences. The presence of sterile neutrinos in neutrino oscillations requires extra neutrino states and the extension of the mixing matrix from 3 × 3 to the larger. Assuming the system with three active neutrinos and N sterile neutrinos, the whole 3 + N mixing matrix as shown in Eq. (2.1) is a complete and unitary mixing matrix. As a result, the left-up 3 × 3 subset that couples e, µ, τ and ν 1 , ν 2 , ν 3 is non-unitary. The effect of non-unitarity is detectable, for example, deficits in disappearance measurements as a part of active neutrinos change into invisible sterile neutrinos.
In this section, we will introduce the standard 3ν scheme and the basic conditions for its unitarity assumption. Then a general 3ν mixing matrix relaxing the unitarity conditions will be presented. Additional conditions in the non-unitarity assumption will be also discussed. After presenting the way to discuss the CP violation, we will come up to discuss neutrino oscillations under the assumption of non-unitarity.

The 3 × 3 unitary mixing matrix
In the standard three-neutrino-mixing scheme, the mixing matrix U 3ν between the flavour and mass eigenstates is defined as Neglecting the two Majorana CP phases, which have no effect on neutrino oscillations, U PMNS is given by three-dimensional rotation matrix described with three mixing angles and one Dirac CP phase, 3) where the θ ij is the mixing angle for each rotation, and δ CP is the Dirac CP phase. This parametrization satisfies the unitarity conditions, We note that the first two conditions Eq. (2.4) and (2.5), are the normalization conditions, while the other two Eq. (2.6) and (2.7) are for the unitarity triangle closure. These conditions break in the (3 + N )-neutrino scenario. For example, in the 3 + 1 case, As a result, by testing the validity of these conditions, we can exam the standard 3ν scheme.

The 3 × 3 non-unitary mixing matrix
A general 3 × 3 mixing matrix consists of nine complex matrix elements. Each of them has real and imaginary parts. In total, there are 18 parameters in this matrix. Among them, five phases can be removed by redefining the phase in eigenstates, leaving nine absolute values and four phases, and the physics is unchanged. One possible parametrization is given as The assignment of the four phases (or image parts) can be arbitrary, while this makes no difference in neutrino oscillations. We use such parametrization Eq. (2.8) in the following analysis.
Once the neutrino mixing goes beyond the standard 3ν scheme, equalities in Eqs. (2.4)-(2.7) are not preserved. Given U N U is a 3 × 3 sub-matrix of a (3 + N ) × (3 + N ) unitary matrix, the sum of squares of each row or column in U N U shall never exceeds 1. Thus the 3 × 3 matrix U N U satisfies: (2.9) Considering any non-unitarity in U N U induced by the active-sterile mixing, additional constraints can be applied, to further bound the parameter space. As discussed in Ref. [43,44], the active-sterile mixing matrix elements in the (3 + N ) × (3 + N ) matrix satisfy Cauchy-Schwarz inequalities for α, β = e, µ, τ, α = β, (2.10) where s k denotes the kth sterile state. By noticing the unitarity of the (3 + N ) × (3 + N ) matrix, Eqs.(2.10) and (2.11) can be rewritten as |U N U βi | 2 , for α, β = e, µ, τ, α = β, (2.13) We note here that the inequalities Eqs. (2.12) and (2.13) need to be satisfied as requiring at least three active neutrinos join in the neutrino mixing. As a result, we will impose these inequality conditions in our later analysis. We need to keep in mind that the constraint on each matrix element of U N U is expected to be worse than those for U 3ν . This is not only because the degrees of freedom are extended for U N U , but also the restrictions on these 9 elements are weaker (compare Eqs. (2.12,2.13) to Eqs. (2.4-2.7)). In addition, it is also interesting to ask if more degeneracy solutions will appear in U N U , or if the existing degeneracy solutions in U 3ν is harder to be separated from the current best fit. These will be included in our later analysis.

CP violation
The CP violation can happen in the active-neutrino sector, and is defined as where U 3×3 can be U 3ν or U N U . The size of CP violation for neutrino oscillations is described by the Dirac CP phase δ CP in the standard 3ν scheme. As PMNS parametrization is not adoptable in U N U , we should use a more general way to scale this size in the case of 3ν non-unitarity.
Under the 3ν unitarity assumption, a rephasing-invariant quantity, so-called Jarlskog factor, describes the CP violation in the neutrino sector, and is the imaginary part of the four-element combination, This is obvious that once the Jarlskog factor is non-zero, the CP is violated as Eq. 2.14.
The Jarlskog factor is an invariant in the standard 3ν scheme, and the factor J can be expressed, J = sin θ 13 cos 2 θ 13 sin θ 12 cos θ 12 sin θ 23 cos θ 23 sin δ. (2.16) As our analysis is beyond the standard 3ν scheme, Jarlskog factors do not need to be the same for different combinations. In other words, the second equality in Eq. (2.15) does not hold. In the ν e − ν µ sector we have three Jarlskog factors: ], which will be discussed in the following analysis.

Oscillation probabilities
With the non-unitary 3ν mixing matrix U N U , the probability of a neutrino of flavour α to be detected as a neutrino of flavour β is given by where ∆m 2 ji is the mass-squared difference between ν j and ν i . L is the baseline distance and E ν is the neutrino energy. The second term is a flavour-changing term, which is a function of ∆m 2 ji L/4E ν , and is independent of CP violation. The last term, which is a combination of the imaginary part of the quartic product (U αi U βj U * αj U * βi ) and a periodic oscillation, is the CP violating term. It vanishes when CP is conserved. For antineutrinos, the oscillation probability is the same as Eq. (2.17), but the matrix elements need to be replaced by their complex conjugate partners (U αi → U * αi ). Eq. (2.17) also shows that for the purpose of measuring a certain matrix element we need to choose a proper L/E ν setup, and even the flavour of neutrino is also an important factor in the measurement. Therefore, we will include data from different experiments in the analysis, trying to constrain all elements in U N U , including real and imaginary parts. These experiments will be introduced in the next section.

Simulation details and analysis setup
To analyse each of the matrix elements in U N U , we include data from a variety of neutrino oscillation experiments, such as the Daya Bay [47], Double Chooz [48], KamLAND [49], NOvA [8], OPERA [50], RENO [51], SNO [1,52,53], and T2K experiments [9]. Extra results from sterile neutrino searches are also taken into account. In the following of this section, we will introduce how the neutrino mixing parameters are determined by different types of experiments, and present simulation details for the reactor, accelerator, and solar experiments. We will finally show how we combine data from all of the experiments in our statistical analysis, and how we include the extra information from sterile-neutrinosearching experiments. Table 1. All relevant categories, experiments, and the corresponding measurements in this work. The abbreviations MBL and LBL denote medium and long baseline, respectively.

Types
Exps Measurements The purpose of combining different experiments is to try to measure all matrix elements of U N U independently as much as possible. These experiments are sensitive to different U N U elements, because of measuring different oscillation channels, focusing on neutrino energy ranges or using different oscillation baselines. We classify them according to their measurements. Moreover, to better understand how many matrix elements are measured independently, we introduce these experiments in this classification method. We summarise all experiments in Table 3.1. More details are introduced in the following.
Medium-baseline (MBL) reactor neutrino experiments: Medium-baseline reactor neutrino experiments play crucial roles in examining the unitarity condition in the ν e sector. We include latest data from Daya Bay, Double Chooz, and RENO. With the baseline of ∼ 1 km and 0.8 − 12 MeV neutrino energy, these experiments are sensitive to theν e →ν e oscillation driven by mass squared splittings |∆m 2 31 | and |∆m 2 32 |. Thus these data can be used to constrain |U e3 |, and |U e1 | 2 + |U e2 | 2 .
Long-baseline (LBL) reactor neutrino experiments: KamLAND is the only longbaseline reactor neutrino experiment by far. KamLAND studies the reactor anti-neutrino oscillation, utilizing more than 50 nuclear power reactors. The flux-weighted average baseline to the reactors is ∼ 180 km. The leading term of the KamLAND measurement is Solar neutrino experiments: Solar neutrino experiments study ν e flavour conversions from the sun. In this analysis, we consider the observations of solar 8 B neutrinos. Generated at the center of the sun, solar 8 B neutrinos travel through the dense matter to the sun surface. The high matter density in the solar core enables solar 8 B neutrinos with energy of ∼6 MeV to undergo MSN resonant transitions. Also, in such environment with high matter density, the superposition of ν e in the effective mass states is directly ν 2 . Then, some of these neutrinos travel through the vacuum to the earth. The distance between the sun and the earth is so long that the oscillation probability average out at the surface of the earth. Thus, the fraction of solar 8 B neutrinos that remain ν e at the earth is the effectively the same as the fraction of solar 8 B neutrinos being ν e at the surface of the sun. Measurements of solar 8 B neutrino charged-current flux verify this fraction, directly probe |U e2 | 2 , and therefore disentangle the degeneracy between |U e1 | and |U e2 |.
Long-baseline accelerator neutrino experiments (NOvA and T2K): The non- are mainly constrained by accelerator experiments measuring ν µ (ν µ ) disappearance channels and ν e (ν e ) appearance channels. We consider here the NOvA and T2K experiments, which are designed to study accelerator neutrino oscillations driven by |∆m 2 31 | and |∆m 2 32 |. Neglecting sub-leading terms and matter effects, their disappearance measurements and appearance measurements constrain , respectively. An ambiguity presents in measuring |U µ3 | and |U µ1 | 2 + |U µ2 | 2 , as known as "octant degeneracy" in the context of 3ν unitarity [54].

Long-baseline accelerator neutrino experiments (OPERA):
To get the independent measurement in the ν τ sector, we further include the observation of 5 ν τ events from a ν µ beam at OPERA. With a baseline of 730 km and a beam energy peaked at ∼ 20 GeV, the experiment can only measure the tail of the first cycle of the oscillation driven by |∆m 2 31 | and |∆m 2 32 |. This measurement can be used to constrain the combination

Medium-baseline reactor neutrino experiments
The Daya Bay, Double Chooz, and RENO experiments adopted relative near-far measurements, allowing inter-detector ratio fits, so that uncertainties and potential biases from neutrino sources are minimized. We take the far-to-near ratios R obs i from Double Chooz [48] and RENO [51], and the measuredν e survival probability versus L ef f / E ν from Daya Bay [55], along with the corresponding error σ obs i . 1 L ef f is the effective propagation distance at Daya Bay, and E ν is the average trueν e energy. Background has been removed from these data points.
The data points from Double Chooz and RENO are ratios of far spectra to expected nooscillation spectra, where the expected no-oscillation spectra are obtained by weighting the observed spectra in the near detectors with no-oscillation assumptions. Noting correlated terms cancel out for a relative measurement, our prediction is given by: Eq. (2.17), given distance L jk and neutrino energy E i . Φ j is the number of neutrino generated in the jth reactor, and L jk denotes the distance from the jth reactor to the kth detector 2 .
The Daya Bay data points are measured survival probabilities, and can be expressed as 2) where N obs i is the observed number of event in bin i. N no−osc i is the expected number of event with no oscillation, derived from near-site measurements. We redefine ω i ≡ (L/E) i and use the following equation for Daya Bay prediction .

(3.3)
A χ 2 quantity is constructed for a medium-baseline reactor neutrino measurement, comparing our predictions with the observed ratios, where the subscript M R denotes different medium-baseline reactor experiments. All the three measurements share the same formula. Eventually χ 2 s for each measurements will be summed up,

Long-baseline reactor neutrino experiment
KamLAND experiment uses 1 kton of liquid scintillator to monitorν e flux from more than 50 nuclear power reactors at long baselines. We take the information of 23 major reactors from Ref. [57], including their thermal power and distances to KamLAND detector. Our analysis is using the ratio of data to no-oscillation expectation (Fig. 5 in Ref. [49]) from KamLAND. As there is only one detector, correlated terms do not cancel out. To analyse KamLAND data, we use a function similar to Eq. 3.1 to predict the ratio, with an additional normalization factor θ KLD which deals with the correlated uncertainties, Then we construct the following χ 2 LR function, with the uncertainty of θ KLD being σ KLD = 5%. The notation σ obs i is the published uncertainty of the corresponding observed data. (3.7)

Solar neutrino experiment
For solar neutrinos, we include the experimental results from the SNO experiment. The SNO 8 B solar neutrino flux measured with charged-current (CC) interactions, Φ obs CC , are analysed 3 . We take Φ obs CC from all three phases of the SNO experiments in Ref. [1,52,53], alongside the published errors σ CC . Assuming adiabatic evolution, the predicted 8 B solar neutrino CC flux Φ pred CC is expressed as where | U ei | is an effective mixing matrix element where 8 B neutrinos are produced in the sun, and |U ei | is the mixing matrix element in vacuum. The effective mixing matrix elements depend on the charged-current potential at the neutrino source. The NC potential is not considered, as in this work active-sterile mixing angles (θ A−S ) are assumed to be tiny compared to the mixing angle in the active section (θ ij ), i.e. θ A−S θ ij . We consider 8 B neutrinos generated at a single point where the matter density is 93.11 g/cm 3 and with the energy 6.44 MeV [58]. Φ8 B is the predicted 8 B solar neutrino flux from solar model BS05(AGS,OP) [58], which is 4.51 × 10 6 cm −2 s −1 . Simultaneously fitting results from all three SNO operational phases, the χ 2 solar is defined as (3.9)

Long-baseline accelerator experiments
Latest data from NOvA [8] and T2K [9,59] from ν µ (ν µ ) disappearance and ν e (ν e ) appearance channels, for both neutrino and anti-neutrino modes are included in this analysis. The NOvA data consist of an exposure of 0.89 × 10 21 POT neutrino and 1.23 × 10 21 POT anti-neutrino beams. The T2K data include measurements from a neutrino beam mode with an exposure of 1.49 × 10 21 POT, and an anti-neutrino beam with an exposure of 1.64 × 10 21 POT. Each experiment collects several data samples, e.g. the T2K experiment has ν µ (ν µ )-enriched samples, ν e (ν e )-enriched samples, and ν e CC1π + sample. We take the spectra sampled by far detectors.
(ν e +ν e )CCQE NC, Beam (ν e +ν e ) ν e appearance (ν e +ν e )CCQE NC, Beam (ν e +ν e ) ν e CC1π + appearance ν e CC1π + NC, Beam (ν e +ν e ) NOvA (ν e +ν e )CCQE NC, Cosmic ν e appearance (ν e +ν e )CCQE NC, Cosmic To analyses the published results, we refer to Ref. [60] and predict the spectra according to available information, including neutrino flux spectra, cross-sections, energy responses, and backgrounds. We summarise the data samples used in this analysis in Table. 2, along with the signal and background components we analysed. For a ν α -enriched data sample α, the estimated number of event N i,α in the ith energy bin is given by where N α bkg,i is the number of background events in the energy bin i, which we have extracted from NOvA [8] and T2K [61,62] predicted spectra. We do not include any oscillations for the background components. E ν and E rec are true and reconstructed neutrino energy, respectively. η (E ν ), as a function of E ν , is the detection efficiency for the event η. [E i , E i+1 ] are bin edges. dΦ dEν is the incident ν µ (ν µ ) flux. We extract the neutrino flux spectra for T2K from Ref. [63], and for NOvA from Ref. [64]. The flux spectra are rescaled by the number of protons on target and the fiducial number of nucleons in the detector. R η (E rec , E ν ) is the energy response function given by: and η denotes different interactions such as CCQE, CC1π + , or CC-nonQE. The NOvA far detector energy resolution δE ν /E ν is 9% for ν µ (ν µ ) CCQE events, and 11% for ν e (ν e ) CCQE events [65]. The energy resolution for the T2K far detector is described in the formula, which fits to 14% at 0.7 GeV, 8.4% at 1 GeV, and 5% at 2 GeV taken from Refs. [66,67]. Concerning the energy loss as hadron energy may not be properly measured in T2K far detector, we assume a −0.4 GeV energy shift for CC non-QE components. σ η (E ν ) is the cross-section for interaction η. We use the default cross section generated by GENIE [68] for NOvA, and extracted cross-section for T2K from Ref. [63]. P νµ→να is the ν µ → ν α oscillation probability, taking into account matter effects for neutrino passing through the Earth's crust. However, only charged-current potential is taken into account. Similar to solar neutrino experiments, we assume that these NC matter effects from the active-sterile mixing are negligible, as these mixing angles are assumed to be tiny in this work. Including NC matter potentials, the neutrino oscillation patter will depend on the number of sterile neutrinos, which is set to be unknown in this analysis. This assumption can be realised, and we leave it to the future work.
We also take into consideration 5 ν τ events observed in a ν µ beam in OPERA. A rateonly fit to the observed number of events is performed. The prediction for the ν τ event number in OPERA is expanded where the no-oscillation prediction dN ντ no−osc (θ best ) dEν is obtained from the best-fit prediction dN ντ osc (θ best ) dEν [69], divided by the best-fit oscillation probability P 3ν νµ→ντ (θ best , E ν ). The estimated number of background N bkg [50] is 0.25 and is assumed to be unchanged.
Eventually, the estimated spectra are compared to the observed spectra N obs i (total number of ν τ events in the case of OPERA) by constructing a Poissonian χ 2 acc , where σ j are uncertainties [60] of nuisance parameters θ x and their values are summarized in Table.3. The normalization factors for NOvA ν µ disappearance channel andν µ disappearance channel are fully correlated. In other words, these two channels share a common nuisance parameter.
It is worthy to note that in analyses for long-baseline accelerator experiments, some of the inputs are derived assuming unitarity. For example, cross-sections for the T2K experi- NOvA ν µ (ν µ ) disappearance 6% ν e appearance 5% ν e appearance 6% ments were tuned according to T2K near detector data. In presence of sterile neutrinos, the true oscillation probability deviates from that assuming unitarity. This can be a potential improvement into our future work.

Limits to sterile neutrinos
Results from sterile neutrino searches are used in this analysis to provide extra information of the sterile sector. Global fits with data from various experiments [70,71] report exclusion limits on the sterile neutrino mixing. We treat the exclusion limits as constraints to the nonunitarity. The constraints are considered to be Gaussian, with central values being zero. We note for large mass-squared differences, the sterile neutrino mixing driven oscillations are averaged out, and the exclusion limits on sterile neutrino mixing parameters will be independent from mass-squared differences. Proper choices are made with exclusion limits for ∆m 2 41 ≥ 0.1eV 2 [43]. Constraint provided by Ref. [70,71] are for |U α4 | 2 and 4|U e4 | 2 |U µ4 | 2 , which cannot be used directly in our analysis. Therefore, we need a translation. Limits on |U α4 | 2 are interpreted as limits on the variation from unity of the normalization of row α, i.e. 1 − 3 i=1 |U αi | 2 (α = e, µ, τ ). Limits on 4|U e4 | 2 |U µ4 | 2 are interpreted as limits on To constrain each of the non-unitarities, we define the χ 2 sterile as where σ α are the exclusion limits at 1σ confidence level. Shown in Table.4 are the exclusion limits and the non-unitarities they constrain. Exclusion limits reported at certain confidence levels are recast to limits at 1σ confidence level, assuming Gaussian distribution. Non-unitarity Data Limit (1σ) CDHS+MNS+NOvA 0.0659 +MB+SK+DC+IC+SNO

χ 2 function
To combine all data from above experiments to constrain matrix elements of U N U , we sum up all χ 2 values from Eqs. (3.5, 3.7, 3.9, 3.14, and 3.15) and define a total χ 2 value, total is a function of U N U , two mass-squared differences ∆m 2 21 and ∆m 2 31 and a vector θ. In more details, inside U N U there are 13 mixing parameters, including 9 absolute values and 4 phases. The vector of nuisance parameters θ collects all normalization factors applied to the prediction for KamLAND Eq. (3.6) and accelerator neutrino experiments Eq. (3.10). For the ease of calculation, we fixed the mass-squared differences as ∆m 2 21 = 2.51×10 −3 eV 2 and ∆m 2 32 = 7.53×10 −5 eV 2 , assuming normal mass ordering. Except for the constrained parameters and two mass-squared differences, the χ 2 total is minimized over all parameters by the TMinuit.Migrad minimizer [72]. Finally, in our simulation, we impose the input hypothesis satisfying constraints such as Eq. (2.9) and Cauchy-Schwarz inequalities Eq. (2.12)-(2.13).

Results
After combining all data from medium-baseline reactor, long-baseline reactor, solar, longbaseline accelerator experiments that are introduced in Sec. 3, we will firstly present the resultant constraints on U N U with Eq. (3.18) in this section. Both 3ν unitarity and nonunitarity assumptions are considered. We will also visit the CP violation in the case of 3ν non-unitarity, before showing the goodness of fit of the current data to the 3ν unitarity assumption.
3 τ Figure 1. The constraints on the moduli of the mixing matrix elements with (black dashed) and without (red solid) unitarity. The x-axes are the moduli of matrix elements, ranging from 0 to 1. The y-axes are the ∆χ 2 , ranging from 0 to 9. A fit setting χ 2 sterile = 0, i.e. without using limits on sterile neutrinos, is also shown (blue dotted).
We present the ∆χ 2 values against the absolute values of each matrix element in Fig. 1 for different scenarios. We perform several groups of fittings with different assumptions: 'N U ' (blue-dotted) and 'N U ⊕ ν s ' (red) are for the fit assuming non-unitarity with and without sterile search results, respectively; '3ν' (black-dashed) is for the fit assuming 3ν unitarity.
At the first glance of Fig. 1, the 'N U ⊕ ν s ' fit is in good agreement with the '3ν' fit. In addition, the ∆χ 2 behaves the same in all scenarios in the ν e sector. This indicates that the current reactor and solar data seem to well match the unitarity of the 3ν mixing matrix. From the narrowness of the 1-D ∆χ 2 curves, one sees that elements |U ei | have been precisely measured, with the 3σ errors being 0.065, 0.103, and 0.013 for |U e1 |, |U e2 |, and |U e3 |, respectively. The precision measurements in the ν e sector are passed to the ν µ sector via the Cauchy-Schwarz inequalities Eqs. (2.12) and (2.13).
Let us draw readers attention to the spike-like structure in the 1-D ∆χ 2 curves at |U µ1 | ∼ 0.5 and the dip structure around |U µ2 | ∼ 0.65 and |U µ3 | ∼ 0.7, in Fig. 1. In this analysis, ν µ → ν µ disappearance and ν µ → ν e appearance measurements at NOvA and T2K provide predominant constraints to the ν µ sector. The leading terms of the oscillation probabilities in NOvA and T2K are easily derived from Eq. (2.17), read (4.1) With the 3ν unitarity assumption, the leading term for ν µ disappearance channel is which is symmetric with respect to |U µ3 | = √ 2 2 , known as the octant degeneracy in the standard 3ν scheme. This degeneracy can be disentangled by combining the ν µ disappearance measurements with the ν e appearance measurements. Our unitary fit '3ν' shows a preference for the upper octant, which is consistent with the published results from NOvA and T2K, with the best-fit value |U µ3 | = 0.733. Without the 3ν unitarity assumption, the same global best-fit is also found, while for 'lower octant' 6 |U µ3 | = 0.689(∆χ 2 = 0.84), and is consistent with |U µ3 | = √ 2 2 within 1σ. The spike-like structure around |U µ1 | = 0.5 and the dip structure around |U µ2 | = 0.6 are due the conversion of |U µ3 | from 'upper octant' to 'lower octant', left-to-right. More details will be introduced in Fig. 2.
We summarise our result as follows. The best-fit points are the same for both nonunitary and unitary cases, We intend to investigate the correlations between the 3ν mixing matrix elements when the conditions of 3ν unitarity are not imposed. In Fig. 2, we plot 1 to 3σ allowed regions on the |U µi | − |U µj | plane, with (black) and without unitarity (red). On the |U µ1 | − |U µ2 | 6 Here "octant" means |Uµ3| = in the case of unitarity, to show the maximal mixing. All of the contours are depressed around the maximal mixing. We also plot the local minimum (light blue and orange curves) with fixed U µi value, and find the discontinuity in Fig. 1, happens at the place where the local minimum switches the octant.

CP violation
The CP violation in the case of non-unitarity is currently measured in the form of Jarlskog factors J αβij . For U 3ν , an invariant Jarlskog factor characterizes CP violation, as expressed in Eq. (2.15). Without the 3ν unitarity assumption, the Jarlskog factor is no longer an invariant. Taking different combinations of αβij, there can be at most nine different Jarlskog factors. As we have little direct measurements of the ν τ sector, six of the Jarlskog factors related to the ν τ sector are not analysed. We show in Fig. 3 three Jarlskog factors associated to the ν e − ν µ sector: , along with the 3ν Jarlskog invariant.
] (green). A minus sign is appended to J eµ13 . Also shown is the 3ν Jarlskog invariant (black dashed).
The deviation from the CP conservation scenario J αβij = 0 is found in three results. The best-fit values are -0.0137, -0.0128, and -0.0125, for J eµ12 , J eµ13 , and J eµ23 , respectively. These values are consistent with the best-fit value in the standard 3ν scheme (J = −0.0148). In addition, the CP conservation hypothesis is excluded with ∆χ 2 = 1.2(0.8) with(without) unitarity.

Conditions for 3ν mixing unitarity
The closure of each unitarity triangles are direct tests of the 3ν unitarity hypothesis. In Fig. 4, we show the the e − µ unitarity triangle on the left panel, and a zoom-in of the comparison between the best-fits with and without unitarity on the right panel. We present here only the e − µ unitarity triangle for the same reason as the Jarlskog factor analyses Fig. 3. The e − µ unitarity triangle corresponds to the unitarity condition, We parameterize the three sides of the e − µ unitarity triangle as where each side consists of a combination of U ei U * µi and is normalized so that side C is unity. If U is unitary, the three sides form a close triangle with A + B + C ≡ A + B + 1 = 0. However, this is not always true when the unitarity conditions are not imposed. There could be three different triangles, depending on which side to be normalized to unity. Here side C is chosen for the ease of calculation, because we did not put phases on U e3 and U µ3 , as shown in Eq. (2.8). We plot the allowed regions for the complex numbers In the following we would like to understand the shapes of the contours. Before discussing the features, we would like to remark on the e − µ unitarity triangle. The vertex of the e − µ unitarity triangle is given by sin θ 12 cos θ 12 cos θ 23 sin θ 13 sin θ 23 e iδ CP . (4.7) One finds the vertex goes along a circle, with the center at (sin 2 θ 12 , 0) in the complex plane, and the radii being (sin θ 12 cos θ 12 cos θ 23 )/(sin θ 13 sin θ 23 ). The vertex points towards (−∞, 0) when δ CP = 0 o , and rotates around the center counter-clockwise as δ CP increases. Without the 3ν unitarity assumption, the endpoints of side A and B seem to have the same behaviours as in the case of unitarity. The contours, i.e. the allowed regions of the endpoints, are roughly ring shaped. The ambiguity in the radii of both rings are caused by the octant degeneracy when measuring |U µ3 |, which we have discussed in Sec. 4.1. Currently only ν e appearance measurements provide knowledge of CP phases, and therefore result in preferred directions of the e − µ triangle. In our case, the e − µ triangle prefers the top-right direction, and the bottom-left parts of the rings are unfavored. More data from running experiments like NOvA and T2K, and future experiments like DUNE and T2HK, help shrinking the contours to pin down the CP violation, and test the unitarity. Alternative methods have been discussed, which proposed to determine the combinations of U ei U * separately [73].
We also present tests of the unitarity by verifying the following quantities: and (4.9) The normalization deviations from the 3ν unity condition for each row and column are defined in Eq. (4.8), and the unitarity triangle closure deviations from zero are defined in Eq. (4.9). If the 3ν unitarity holds, these quantities are zero. In Fig. 5, we plot 1-D ∆χ 2 for δ α (δ i ) on the left panel, and 1-D ∆χ 2 for ζ αβ (ζ ij ) on the right panel. No signal of non-unitarity is found. We find the violation from unitarity in the ν e sector (δ e ) and in the ν µ sector δ µ well constrained to be less than 0.003 and 0.02, at 3σ CL. These have benefited from the precise measurements to reactorν e disappearance oscillations, accelerator ν µ (ν µ ) disappearance oscillations, and the solar 8 B ν e flavour conversions. δ τ is less known with an upper bound of ∼ 0.2, for we lack of knowledge to the ν τ sector. Two out of three elements in a column, |U ei | and |U µi |, can be measured. Therefore, δ i for each column are constrained better than δ τ , being less than 0.06-0.2 at 3σ CL.
The non-closure of the e−µ unitarity triangle ζ eµ is constrained to be less than ∼ 0.005. With little(no) experimental data of the ν µ → ν τ (ν e → ν τ ) channel, ζ µτ and ζ eτ are known to be 0.05 and < 0.02, worse than ζ eµ . The precision measurements in the ν e and ν µ sectors provide dominant constraints to ζ eτ and ζ µτ via Cauchy-Schwarz inequalities. The remaining three column-wise unitarity triangles, bounded only by the Cauchy-Schwarz inequalities, are known to be < 0.07 − 0.08, at 3σ CL. One can roughly estimate these constraints from δ 1 , δ 2 , and δ 3 , with ζ ij ≈ δ i δ j .

Conclusions
The development in neutrino oscillations in the past decades allows us to conduct precision measurements of the neutrino mixing in the active sector (U N U ). Entering the new precision era, we are able to explore other possibilities, e.g. the 3ν unitarity-violated neutrino mixing hypothesis.
An analysis of neutrino oscillations was performed without unitarity assumption in the 3ν picture. We have combined the medium and long-baseline reactor, solar, long-baseline accelerator neutrino data to constrain the mixing matrix in the active sector U N U . We have found that elements U ei are measured to be the best among all sectors (3σ uncertainty < 20%). At the same confidence level, the uncertainties > 20% were obtained in |U µi | in the current global analysis. Though currently data for the ν τ sector are limited, via Cauchy-Scharz ineqauilities Eqs. (2.10) and (2.11) the constraints for this sector can be passed from that in the µ sector. And, therefore the size of uncertainties for the ν µ and ν τ sectors are similar. Our result prefers |U µ3 | 2 > 1/2, which corresponds to the upper-octant solution in the standard 3ν scheme. A negative correlation was noticed between |U µ1 | and |U µ2 |, as the ν µ disappearance measurements determine the combination |U µ1 | 2 + |U µ2 | 2 . The ν e appearance measurements help with distinguishing |U µ1 | and |U µ2 | at 1σ C.L., as well as the "octant-like degeneracy".
Other properties of U N U were further discussed. The CP-violation in the e − µ sector was investigated by measuring three Jarlskog factors J eµ12 , J eµ13 , and J eµ23 . The result is similar to that in the case of 3ν unitarity, though the uncertainties are slightly worse. Furthermore, the e − µ unitary triangle coincides with the unitary 3ν mixing scheme Fig. 4. This encouraged us to test the assumption of 3ν unitarity Eq. (4.8) and (4.9), as given in Fig. 5. Normalizations: δ e and δ µ (Eq. (4.8)) are relatively well-constrained. Closure: we know ζ eµ (Eq. (4.9)) the best. No significant violations to these conditions are observed.
Our results also bring some prospects of the future experiments on the 3ν unitarity hypothesis testing, while we summarise the current status in Table 3.1. The precision of |U e1 | and |U e2 | can be further improved by theν e detection in the JUNO experiment. Future accelerator experiments DUNE and T2HK can provide more data relevant to the ν µ sector, and reduce the statistical uncertainties. Further, the large matter density and long baseline of DUNE can be used to measure the NC matter effects due to the active-sterile mixing. The lack of ν τ events brings the large uncertainties of |U τ i |, which might be improved by the possibility of ν τ detection in DUNE. On the phenomenological analysis itself, we need more information from experiments that is not based on the 3ν assumption, such as systematics uncertainties.
We close up this work with a conclusion that we have not found any significant evidence for the 3ν non-unitarity, though the uncertainties (mainly from statistical uncertainties) still needs to be improved. We see the prospect of future DUNE, T2HK, and JUNO data, and are looking forward to any experimental proposals for ν τ oscillation channels. At the mean time, we will keep sharpening our simulation and analysis tools.