Model-Independent Constraints on Non-Unitary Neutrino Mixing from High-Precision Long-Baseline Experiments

Our knowledge on the active 3$\nu$ mixing angles ($\theta_{12}$, $\theta_{13}$, and $\theta_{23}$) and the CP phase $\delta_{\mathrm{CP}}$ is becoming accurate day-by-day enabling us to test the unitarity of the leptonic mixing matrix with utmost precision. Future high-precision long-baseline experiments are going to play an important role in this direction. In this work, we study the impact of possible non-unitary neutrino mixing (NUNM) in the context of next-generation long-baseline experiments DUNE and T2HKK/JD+KD having one detector in Japan (T2HK/JD) and a second detector in Korea (KD). We estimate the sensitivities of these setups to place direct, model-independent, and competitive constraints on various NUNM parameters. We demonstrate the possible correlations between the NUNM parameters, $\theta_{23}$, and $\delta_{\mathrm{CP}}$. Our numerical results obtained using only far detector data and supported by simple approximate analytical expressions of the oscillation probabilities in matter, reveal that JD+KD has better sensitivities for $|\alpha_{21}|$ and $\alpha_{22}$ as compared to DUNE, due to its larger statistics in the appearance channel and less systematic uncertainties in the disappearance channel, respectively. For $|\alpha_{31}|$, $|\alpha_{32}|$, and $\alpha_{33}$, DUNE gives better constraints as compared to JD+KD, due to its larger matter effect and wider neutrino energy spectrum. For $\alpha_{11}$, both DUNE and JD+KD give similar bounds. We also show how much the bounds on the NUNM parameters can be improved by combining the prospective data from DUNE and JD+KD setups. We find that due to zero-distance effects, the near detectors alone can also constrain $\alpha_{11}$, $|\alpha_{21}|$, and $\alpha_{22}$ in both these setups. Finally, we observe that the $\nu_\tau$ appearance sample in DUNE can improve the constraints on $|\alpha_{32}|$ and $\alpha_{33}$.


Introduction and motivation
The Standard Model (SM) of particle physics is the most successful gauge theory which has been rigorously tested in several laboratory experiments, including those which are being performed at the Large Hadron Collider at CERN in Geneva [1]. In spite of this huge success, the observed input parameters of the SM, such as the Higgs mass, the fermion masses and mixings, and the QCD theta term are quite offbeat and certainly require further explanations [2,3]. These pressing issues point towards the existence of a microscopic theory beyond the Standard Model which should address the electroweak hierarchy problem [4,5], the flavor puzzle, and the strong CP problem [6][7][8]. Various next-generation high-precision experiments in the energy, intensity, and cosmic frontiers are expected to provide crucial information on these issues. In the intensity frontier, several high-precision neutrino oscillation experiments are currently running and also in the pipeline to measure the mass-mixing parameters with unprecedented precision. Marvelous data from several pioneering neutrino experiments like Super-K [9][10][11], IceCube-DeepCore [12], ANTARES [13], Daya Bay [14], RENO [15], MINOS [16], Tokai-to-Kamioka (T2K) [17], and NuMI Off-axis ν e Appearance (NOνA) [18] have been decisive for our understanding of neutrino flavor mixing. Global fit analyses of the world neutrino data [19][20][21][22] have already been able to determine the values of the neutrino oscillation parameters with reasonable accuracy. At present, the relative 1σ precision on the mixing angles θ 23 , θ 13 , and θ 12 lies in the range of O(3 − 7)%. For the mass-squared differences, the achieved relative 1σ precision is around O(1 − 3)%. Another important result that is being emerged from these global fit studies is that now we have an overall 1.6σ hint for leptonic CP-violation (sin δ CP < 0) [20,22].
The upcoming medium-baseline reactor experiment Jiangmen Underground Neutrino Observatory (JUNO) [23,24] is going to address the issue of neutrino mass hierarchy at high confidence level. JUNO will also improve our knowledge on the mixing angle θ 12 and mass-squared differences ∆m 2 21 and ∆m 2 31 . The Fermilab-based Deep Underground Neutrino Experiment (DUNE) [25][26][27][28] is expected to resolve the issue of neutrino mass hierarchy and leptonic CP violation with a confidence level never achieved before. DUNE will measure the value of CP phase δ CP and atmospheric oscillation parameters θ 23 and ∆m 2 31 with high precision using its wide band neutrino beam, which provides information on oscillation parameters at several L/E values. Another major candidate in the race of next-generation high-precision long-baseline neutrino oscillation experiments is Tokai to Hyper-Kamiokande (T2HK) with one detector in Japan, which we refer as JD in the present work [29,30]. The possibility of having a second detector in Korea (KD) is also being explored actively. The combination of these two detectors (one in Japan and other in Korea) exposed to a common high-intensity, off-axis (2.5 • ), narrow-band beam from J-PARC, is widely known as T2HKK [31], which we refer as JD+KD setup in the present study. The JD setup with a relatively shorter baseline as compared to DUNE and high statistics is going to address the issue of leptonic CP violation with unmatched sensitivity and measure the value of δ CP quite precisely without facing the issue of fake CP violation due to matter effect. On the other hand, the KD setup having a baseline roughly around four times that of JD contains some information on Earth matter effect and measure δ CP around the second oscillation maximum with reasonable precision.
With the excellent foreseen precision on the neutrino mass-mixing parameters in the coming years, it will be interesting to explore whether tiny new physics effects beyond the standard three-flavor framework are present in the scenario and what are their impact on the neutrino oscillation probabilities [32]. In the present paper, in the context of highprecision long-baseline experiments, we study the possible consequences of the fact that the 3 × 3 active neutrino mixing matrix N no longer respects the unitarity condition N † N = I. Theoretically, there are various neutrino mass models in the literature [33][34][35][36], which allows the possibility of non-unitary neutrino mixing (NUNM). The most appealing ones are the so-called see-saw models, in which new heavy neutral leptons and/or scalars are introduced in the basic Standard Model to explain tiny neutrino masses. In these models, the standard 3×3 active neutrino mixing matrix becomes a non-unitary sub-matrix of the larger mixing matrix (see discussion in Appendix of Ref. [37]). If the heavy lepton is very massive in the see-saw model [38][39][40], then it can cause a departure from the unitarity of the order of 10 −3 . The amount of non-unitarity can be even larger in the low-scale see-saw models [41,42].
In this work, our aim is to study the sensitivity of various long-baseline experiments to probe the non-unitarity of the standard 3 × 3 active neutrino mixing matrix via oscillation, which allows us to constrain the NUNM parameters considering one at a time or all of them at the same time. Here, we adopt a completely model-independent approach to invoke this possible deviation from unitarity without relying on the underlying mechanism.
In the present work, we study the potential of two next-generation long-baseline experiments DUNE and T2HKK (JD+KD), to constrain the NUNM parameters in a complete model-independent approach. We present the expected bounds, as obtained from these experiments individually as well as from their combination, assuming no priors in the fit procedure. Compared to other similar studies, we differ in several aspects; first of all, to carry on our numerical simulations, we use the most recent configuration of the DUNE [27] and T2HKK [31] experiments, also investigating the role of near detectors (ND) and τ neutrino detection in the case of DUNE. Then, for a better understanding of the obtained bounds as well as of the observed correlations between the standard oscillation parameters θ 23 and δ CP with the various NUNM parameters, we provide simple and useful analytical expressions of the relevant transition probabilities in the regime of small matter effects, which is a correct approximation for the neutrino facilities under investigation.
The manuscript is organized as follows. In Sec. 2, we discuss the formalism and the parameterization of non-unitary neutrino mixing used in this work, and derive simple approximate analytical expressions of the oscillation probability in ν µ → ν e and ν µ → ν µ oscillation channels in the presence of standard matter interaction. In Sec. 3, we give a brief description of the future long-baseline experiments discussed in this work. Sec. 4 provides the details of the numerical simulations performed in our analysis and some results on the expected signal events. The main results of our study are presented in Sec. 5, where we illustrate the correlation of the standard oscillation parameters θ 23 and δ CP with various NUNM parameters and give the expected constraints achievable by DUNE and T2HKK separately as well as their combination. In Sec. 6, we show the possible improvements in the bounds due to the presence of near detectors. Section 7 shows the improved constraints on the NUNM parameters when we add ν τ events sample in the DUNE simulation. Our concluding remarks are discussed in Sec. 8. Two appendices complete our work: in Appendix A, we compare our analytical expressions of the oscillation probabilities and those computed numerically in the presence of non-unitarity. We also provide the analytical expressions of the oscillation probabilities for all the remaining channels. In Appendix B, we discuss in detail the impact on the bounds discussed in the main text of α 11 α 22 α 33 |α 21 | |α 31 | |α 32 | < 0.031 < 0.005 < 0.110 < 0.013 < 0.033 < 0.009 Table 1: 90% confidence level (C.L.) limits on the NUNM parameters using data from various shortbaseline and long-baseline experiments, as obtained from the Ref. [53].
the simultaneous marginalization over all NUNM parameters.

Oscillation probabilities in matter with non-unitary neutrino mixing
The non-unitarity of the PMNS matrix can be parameterized in different ways [49,56,[59][60][61][62]. One possibility, which turned out to be very useful in oscillation analysis, consists of factorizing the deviation from unitarity into a matrix α in such a way that the non-unitary neutrino mixing matrix N is expressed as 1 : In order to compare our numerical results with those already presented in the literature, we adhere here at the widely used lower triangular structure of the matrix α, containing nine free parameters organized as follows: This parameterization simplifies the oscillation probabilities and let the parameter α ij to be the main source of non-unitarity for the oscillation channel ν i → ν j (i, j = e, µ, τ ).
Bounds on the α ij parameters have been recently computed, among others, in Ref. [53] and reported for convenience in Table 1. These results have been obtained using data from the short-baseline experiments NOMAD and NuTeV, and the long-baseline experiments MINOS/MINOS+, T2K, and NOνA. For the off-diagonal NUNM parameters, the authors also used the triangular inequalities 2 [40], which appear due to the assumption that the standard 3×3 active neutrino mixing matrix is a non-unitary sub-matrix of a larger unitary mixing matrix. Note that in the present study, we do not take into account these inequalities in order to study the capability of long-baseline experiments alone to put bounds on these NUNM parameters in a modelindependent fashion. 1 We choose the convention in which the matrix α, which invokes non-unitarity, is added to the identity matrix. Note that other phenomenological studies adopt the relation N = (1 − α) UP M N S [39]. Our results can be compared to the others just changing the sign in front of the diagonal elements. 2 Note that in this model, the diagonal parameters can only be negative.  Taking into account the expression in Eq. 2.1, the complete effective neutrino propagation Hamiltonian in the mass-eigenstate is: As usual, the matter potential parameters are given by a e = 2 √ 2E ν G F N e and a n = − √ 2E ν G F N n where, N e and N n are the electron and the neutron number densities, respectively. Note that in this framework, the neutral current (NC) matter potential is necessary since the non-unitarity of the matrix N does not allow the subtraction of an identity matrix proportional to a n . From the Schroedinger equation, the transition probability at a given baseline L is obtained from the following expression: The relevant probabilities for long-baseline experiments are the ν µ → ν e appearance and ν µ → ν µ disappearance channels. We will discuss these probabilities here in details, while for the sake of completeness, we also quote the τ appearance and electron disappearance into the Appendix A.
In order to get approximate analytical expressions for the transition probabilities, we observe that the vacuum approximation cannot be sufficiently precise in experiments like DUNE, since the matter effects can modify the appearance probability up to about 10%. For this reason, we derive approximate analytical expressions in the presence of matter. We use perturbation theory in the small expansion parameters (r, s, and a) defined as follows: where, r, s, and a represent the deviation from the tri-bimaximal mixing values of the neutrino mixing parameters, namely, sin θ 13 = 0, sin θ 23 = 1/ √ 2, and sin θ 12 = 1/ √ 3 [63,64]. Given the recent global fit of neutrino oscillation data, we can assume that r, s, a ∼ O(0.1) and we can further expand them up to the second order [19][20][21][22]. To further simplify the notation, we also introduce ∆ 31 = ∆m 2 31 L/4E ν , ∆ e = a e L/4E ν and ∆ n = a n L/4E ν ; at the DUNE peak energy, namely, E ν = 2.5 GeV, ∆ e ∼ 0.36 and ∆ n ∼ 0.18, we can further expand in the small matter potentials up to the first order. Note that for the other experimental facilities discussed in this paper, this approximation is even better; in fact, for the Tokai to Hyper-Kamiokande (T2HK) setup with a far detector in Japan (JD), at beam energy of E ν = 0.6 GeV, we have ∆ e ∼ 0.08 and ∆ n ∼ 0.04, while with a second in Korea (KD) at a distance of 1100 km from the source with same energy, we get ∆ e ∼ 0.30 and ∆ n ∼ 0.15. Also, we use one mass scale dominance (OMSD) approximation (∆ 31 >> ∆ 21 , where ∆ 21 = ∆m 2 21 L/4E ν ) in our derivation; which is a valid approximation in the atmospheric regime.
In the case of the ν µ → ν e appearance probability, we thus obtain: From the above expression, it is clear that the ν µ → ν e appearance probability strongly depends on |α 21 | and |α 31 |. The parameter |α 21 | survives in the vacuum case while |α 31 | always appears with the matter potential ∆ n . This essentially means that an experiment in which the matter effect is not negligible is able to put strong bounds also to |α 31 | which, otherwise, would not be accessible by P µe . Note that due to the loss of unitarity property of the neutrino mixing, some terms remain non-zero in ν µ → ν e appearance probability expression even when the neutrino propagation length L is zero, which is known as zerodistance effect: So, it is clear that even the near detectors (ND) of long-baseline experiments could contribute to the bounds of non-unitarity parameters, as it will be discussed in Sec. 6. Finally, we point out that the vacuum limit of Eq. 2.6 assumes a particularly simple expression: which, in the limit of vanishing solar mass difference, agrees with the results presented in Ref. [37]. In order to assess the modifications in the ν µ → ν e appearance probabilities caused by the presence of non-unitary neutrino mixing, in Fig. 1, we show the exact ν µ → ν e appearance probabilities as a function of energy, obtained numerically using the General long-baseline Experiment Simulator (GLoBES) software [65,66]. In the extreme left column, we consider a baseline of 1300 km for DUNE. In the middle column, we show the results for the JD baseline of 295 km. In the extreme right column, we deal with the KD setup having a baseline of 1100 km. Note that for both JD and KD, we consider the neutrino energy range of 0 to 1.5 GeV having a peak around 0.6 GeV. In every rows, we switch-on one off-diagonal NUNM parameters at a time, while maintaining others to zero. The effect of |α 21 | is shown in the top panels, |α 31 | in the middle panels, and |α 32 | in the bottom panels. In each panel, the solid black curves correspond to the probabilities in the unitary neutrino mixing (UNM) case, while the colored curves correspond to the NUNM cases with four benchmark values of the phases associated with each off-diagonal NUNM parameters, as reported in the legend. As we can see, the impact of the NUNM parameter |α 21 | (top panels) is comparatively larger than the other two off-diagonal NUNM parameters |α 31 | and |α 32 | even though the strength of the |α 21 | is much smaller than the other two. This feature is more clear for JD in the middle columns. In fact, from the approxi-  Table 2.
mated analytical expression of ν µ → ν e appearance probability in Eqs. 2.6 and 2.8, we see that only terms containing |α 21 | survive in vacuum, while the effect of |α 31 | and |α 32 | is linked to ∆ e and ∆ n which are very small at the considered baseline. This is not the case for DUNE (L 1300 km) where, given the largest baseline under consideration, ∆ e and ∆ n are no longer negligible and the impact of |α 31 | and |α 32 | on P µe is of the same order as |α 21 |. For the ν µ → ν µ disappearance channel, we get: from which we learn that all the NUNM parameters α ij but α 11 enter into the probability expression. Also, in this case, we get the zero-distance expression of the ν µ → ν µ survival probability, given by: 10) and the vacuum approximation (also in agreement with Ref. [37] in the limit of vanishing ∆m 2 21 ): P vac µµ = cos 2 ∆ 31 1 + 2|α 21 | 2 + 4α 22 + 6α 2 22 + 4a 2 sin 2 ∆ 31 . (2.11) In Fig. 2, we show the exact ν µ → ν µ oscillation probabilities as a function of energy for the baseline lengths corresponding to DUNE (left panel), JD (middle panel), and KD (right panel) setups. In each panel, the black solid curves correspond to the UNM case, while the red and blue curves show the presence of α 22 and α 32 , respectively with strength reported in the legend 3 , one at a time. The impact of these two NUNM parameters can be understood from our approximated analytical expressions in Eq. 2.9. When matter effects are negligible (for example, in the middle panel of Fig. 2), we expect that the parameter α 22 dominates the deviation from UNM since it appears already at first order in P µµ . This remains true when matter parameters are switched-on; the relevant difference compared to the vacuum case relies on the fact that also |α 32 | enter at first order, although suppressed by ∆ n . Thus, we expect that for DUNE and KD, one can see deviation from the UNM predictions, as visible in Fig. 2. Note that the impact of α 32 is amplified by the larger benchmark value compared to the choice for α 22 .

Key features of the experiments
Long-baseline experiments provide a clear environment to study the phenomenology of neutrino oscillations. Indeed, the possibility to fix the baseline and have a focused neutrino beam allow us to have different features of the experiments under control. The nextgeneration long-baseline experiments will increase the neutrino events statistics with the help of high precision detectors. Two experiments currently under construction are DUNE and T2HK. For the latter is being taken into account the possibility of a second far detector placed approximately at the second oscillation maximum. These two complementary experiments, which will be described in more details in the following subsections, are going to improve the measurements of the oscillation parameters and the searches of new physics effects in a significant way. In the non-unitarity framework, the unprecedented statistics collected by the detectors of the two experiments, as well as the possibility to look for ν µ disappearance, ν e appearance but also ν τ appearance in DUNE, will provide the chance to bound at a very good level all the non-standard parameters.

DUNE
DUNE (Deep Underground Neutrino Experiment) is a next-generation long-baseline experiment which will play an important role in solving the existing puzzles in neutrino oscillation as well as in increasing the precision on the measured neutrino oscillation parameters. It will use an on-axis, high-intensity, wide-band neutrino beam traveling a distance of 1284 km from Fermilab to South Dakota. For this baseline, they consider an average matter density of 2.85 g/cm 3 . The far detector (FD) is a liquid argon time projection chamber (LArTPC) of 40 kt fiducial mass placed underground at the Homestake mine. According to the recent Technical Design Report (TDR) [27], DUNE will use a proton beam having 1.2 MW power which will deliver 1.1 × 10 21 proton on target (P.O.T.) per year. DUNE considers a total run-time of 10 years with 5 years each in neutrino and antineutrino modes.
However, following most of the existing studies related to DUNE, we present our results using a total run-time of 7 years with 3.5 years each in neutrino and antineutrino modes. Beside the far detector, the possibility to have three modules of near detectors (ND) has been proposed [67]. The first one would be a liquid Argon Time Projection chamber (TPC) situated at a distance of 574 m from the source. It will measure the flux and cross-section of the neutrinos. The second one would be a multi purpose detector (MPD) equipped with a magnetic spectrometer with one ton High-Pressure Gaseous Argon TPC, which will be useful to study possible new physics signals. The third one would be the System for On-Axis Neutrino Detection (SAND). It will be made up of a former KLOE magnet and a calorimeter which will track the outgoing particles.

T2HK (JD) and T2HKK (JD+KD)
T2HK (Tokai to Hyper-Kamiokande) is another promising next-generation long-baseline experiment which will play a very important role in δ CP measurement as well as in the study of various BSM physics [68][69][70]. In the T2HK setup [29,71,72], an intense neutrino beam from the J-PARC proton synchroton facility will be detected at the Hyper-Kamiokande (HK) detector situated at a distance of 295 km from the source. The power of the proton beam at the source is 1.3 MW which will produce 27 × 10 21 P.O.T. in its total 10 years run-time. The detector is a water Cherenkov (WC) detector with 187 kt of fiducial mass.
To have an equal contribution from the neutrino and the antineutrino signal events, the proposed running time of the experiment is 2.5 years (25% exposure) in neutrino mode and 7.5 years (75% exposure) in antineutrino mode. The HK detector will be placed at 2.5 • off-axis angle from the neutrino beam-line to receive a narrow band beam with energy peaked at around the first oscillation maximum (∼ 0.6 GeV).
There is also a proposal [73] to have another identical detector which will be situated in Korea, 1100 km far from J-PARC. We assume that this detector will also be placed at an off-axis angle of 2.5 • from the neutrino beamline (similar to JD) and operate around the second oscillation maximum with a baseline of 1100 km and an average neutrino energy of 0.6 GeV. Combination of the T2HK setup along with the detector in Korea is known as T2HKK. However, in our work, we call the Hyper-Kamiokande as the Japanese detector (JD) and the detector in Korea as the Korean detector (KD). For both JD and KD baselines, we assume an average matter density of 2.8 g/cm 3 . In the following, we will show our numerical results on the sensitivity to the various NUNM parameters separately for each detectors as well as for their combination.
In addition to these two detectors, two more near detectors, namely, the ND280 detector and Intermediate Water Cherenkov detector (IWCD), have been proposed. ND280 will be located at 280 m from the source with the same off-axis angle as the far detectors. It will be helpful to measure the flux of the unoscillated neutrino beam and study the neutrino cross-sections [74]. IWCD is a water Cherenkov detector of mass of 1 kt and possibly located at a distance of 1 km from the source [75,76]. One advantage of this detector is that it can be moved vertically to take data at different off-axis angles. In Table 3, we summarize the relevant information of all the experimental setups discussed in this section. 3.5 yrs + 3.5 yrs 2.5 yrs + 7.5 yrs Table 3: Essential features of DUNE [27] and JD/KD [31] experiments used in our simulation.

Simulation details and discussion at the event level
In order to perform our numerical simulations, we use the General long-baseline Experiment Simulator (GLoBES) package [65,66] along with the plug-in MonteCUBES [39]. For the standard oscillation parameters, we use the value as given in Table 2. For the simulation of the DUNE experiment, we use the most recent GLoBES configuration files as given in Ref. [27]. The simulation details of T2HK and T2HKK setups are taken from Ref. [31].
In Table 3, we summarize the relevant information of the two experiments considered in our simulation. Note that while showing the sensitivity of the FDs in constraining various NUNM parameters, we do not take into account the NDs in our simulation. The only exception is in Sec. 6, where we discuss the possible improvement in the bounds due to presence of near detectors without considering any possible correlations among FDs and NDs (see Sec. 6 for further details). Both DUNE and T2HKK are expected to have access to two oscillation channels, namely the ν µ (ν µ ) disappearance and the ν e (ν e ) appearance channels. In particular: • for DUNE, backgrounds to the appearance channel are the ν e beam contamination and the misidentified ν µ , ν τ , and NC events. The signal systematic normalization error have been chosen to be 2%. For the disappearance channel, background to the signal are misidentified ν τ and NC events. The signal error is 5%. Efficiency functions as well as smearing matrices have been provided by the DUNE collaboration in Ref. [27]; • for the two T2HKK far detectors JD and KD, background in the appearance channels are the ν e from the beam contamination and misidentified ν µ and NC events. Signal systematic uncertainties are 5% normalization and 5% calibration errors. In the disappearance channel, backgrounds to the signal are misidentified ν e and NC events. Systematic uncertainties are 3.5% normalization and 5% calibration errors. Efficiencies and energy resolutions are taken from the Ref. [31].
The number of expected signal events for both channels simulated here are summarized in Tables 4 and 5 Table 3. The values of the standard oscillation parameters used to calculate event rate are quoted in Table 2. The phases associated with the off-diagonal NUNM parameters are considered to be zero.
is fully in agreement with our analytical discussions. First of all, as shown in Eq. 2.6, the appearance channel is mainly influenced by |α 21 |, even in vacuum. This reflects in an enhancement of the number of events by roughly 5% in both JD and DUNE. On the other hand, the NUNM parameters |α 31 | and α 33 are also relevant but they are coupled to the matter potentials, so we expect them to be relevant primarily for DUNE, where matter effects are more important: in fact, |α 31 | causes an increase in the number of events up to 10%, while α 33 provokes a small but visible reduction of the order of 4%. Finally, some impact on the number of signal events is also given by α 11 , even though it only appears at higher orders in our perturbative expansion and has not been displayed (but it present in the vacuum probabilities reported in Ref. [37]). Note that the number of ν e andν e events in KD is only slightly influenced by the NUNM parameters due to fact that the experiment works close to the second oscillation maximum of the atmospheric oscillation (ν µ → ν τ ), where the ν µ → ν e appearance probability approached one of its minima and the effects of new physics are suppressed.
For the disappearance channel, the parameter α 22 , which enters at the first perturbative order in Eq. 2.9, produces a reduction of about 6% in the number of events for all three detectors. This can be roughly understood from the fact that the standard disappearance probability is multiplied by 4α 22 = 0.06, which causes a reduction by a similar factor in the number of events. The other relevant NUNM parameter is |α 32 | which, being coupled to matter potential in Eq. 2.9, can cause a ∼ 6% increase of events especially in DUNE. The other parameters at their benchmark values only have a negligible impact on the number of disappearance events.
In Fig. 3, we show the relation between the appearance events in neutrino and antineutrino modes at the three detectors discussed in this paper, namely DUNE (left panel),  The plots show that DUNE (left panels) is in principle expected to be able to distinguish the mass hierarchy even in presence of NUNM. On the other hand, the two standard model curves and the two shadowed regions overlap in the case of JD (center panels) and KD (right panels). For these, the amplitude of the shadowed regions when φ 31 is varied is much smaller than the DUNE one since, as it is clear from the analytical formulas, |α 31 | is always coupled to matter effects, which are smaller for the two J-PARC based experiments.
In our numerical simulations on the sensitivity to the NUNM parameters, the true values of the standard oscillation parameters have been chosen as in Table 2, while the true values of the NUNM parameters are set to zero. Since the presence of non-unitarity has basically the same effect on the number of events in the case of NH and IH, as shown in Fig. 3, we consider only the NH case. The computation of the ∆χ 2 is based on the pull method [77][78][79] implemented in the GLoBES software. We study the NUNM parameters by fixing the mixing angles θ 12 , θ 13 , and two mass-squared differences ∆m 2 21 , and ∆m 2 31 both in data and theory at their best fit values as given in Table 2. We check that the marginalization over the atmospheric mass-squared difference ∆m 2 31 does not have any significant effect on our analysis. On the other hand, the only notable effect of the marginalization over the reactor mixing angle θ 13 (which has a very small experimental  Table 2. uncertainty of 3%), is the worsening of the α 11 bound at the level of 15%. This is due to the fact that there is a correlation between these two parameters, which appear in a term proportional to α 2 11 sin 2 2θ 13 in the ν µ → ν e transition probability as shown in Ref. [37]. Finally, we marginalized θ 23 in its current 3σ allowed range [22], which is approximately [40 • , 50 • ] and the CP phase δ CP in its entire possible range [−180 • , 180 • ]. We keep both these parameters with true values as in Table 2. Moreover, we consider one NUNM parameter at a time, i.e., when a parameter is taken into account the others are considered to be zero. In Appendix B, we also discuss in detail how our results would be affected by the marginalization over multiple NUNM parameters. Finally, the phases associated with each off-diagonal NUNM parameters are marginalized over the entire possible range from −180 • to 180 • .

Numerical results
In this section, we provide a detailed discussion of the numerical results obtained with the GLoBES simulations. We study the correlations between the NUNM parameters and the two less constrained standard oscillation parameters θ 23 and δ CP . Then, we determine the expected bounds on all the non-standard parameters at 90% confidence level (C.L.) when the expected data from DUNE, JD, KD, and their various combinations are taken into account.   Table 2).

Correlations in test
In Fig. 4, we show the correlation between various NUNM parameters and the atmospheric mixing angle θ 23   planes, the allowed regions for JD+KD are smaller than the ones corresponding to DUNE, which suggest that the J-PARC based experiments will have a better sensitivity on these parameters (see the following subsection). For the other three NUNM parameters, DUNE shows much better sensitivity than JD+KD, due to the larger matter effects which couples to the NUNM parameters α 33 , |α 31 |, and |α 32 |. In particular, for maximal value of θ 23 , there are no lower limits on α 33 from JD+KD. One interesting point to be noted is that, in the presence of non-unitary mixing, sin 2 θ 23 can be constrained better for non-maximal true values of θ 23 .
In Fig. 5, we show the correlation of various NUNM parameters with the CP phase δ CP . In each panel, the black dots correspond to the true values of δ CP = 0 • , 90 • , 180 • , and − 90 • and the benchmark value chosen for the α ij under consideration. The red contours refer to the sensitivity at 1σ (2 d.o.f.) C.L. expected to be obtained by DUNE setup and the blue curve shows the same for JD+KD combination. The plots follow a similar behavior already seen in Fig. 4. In particular, the allowed regions for the two experiments are almost the same in the case of α 11 , smaller for JD+KD in the case of α 22 and |α 21 | and, finally, much smaller for DUNE in the case of the NUNM parameters α 33 , |α 31 |, and |α 32 |. Note that when true δ CP = 90 • , 180 • , and − 90 • , the J-PARC based detectors will not be able to set any lower limit on α 33 .

Constraints on non-unitary neutrino mixing parameters
In this subsection, we present our numerical results showing the expected constraints on the six NUNM parameters (α ij ) that DUNE, JD, KD, and JD+KD setups can place. We derive bounds on α ij following the simulation method as discussed in Sec. 4. The statistical significance with which we can constrain the NUNM parameters (α ij ) in a given experiment is defined as where, χ 2 (α ij = 0) and χ 2 (α ij = 0) are calculated by fitting the prospective data assuming NUNM (α ij = 0) and UNM (α ij = 0). Note that χ 2 (α ij = 0) ≈ 0 because the statistical fluctuations are suppressed to obtain the median sensitivity of a given experiment in the frequentist approach [80]. While estimating the constraints, we marginalize over the most uncertain oscillation parameters (θ 23 , δ CP ) and the phases associated with the off-diagonal NUNM parameters (φ ij ) in the fit. We also minimize over the systematic pulls on signal (λ 1 ) and background (λ 2 ). In Fig. 6, we plot the ∆χ 2 function for the six NUNM parameters analyzed in our paper considering only the far detectors in a given setup. Upper (Lower) panels show the results for the diagonal (off-diagonal) NUNM parameters. The red curves in each panel refer to the sensitivities obtained for the DUNE setup considering a total of 336 kt-MW-yrs exposure, corresponding to a total 7 years of data collection with equal run-time in neutrino and antineutrino modes. The green curves show the results for JD for which we consider a total exposure of 2431 kt-MW-yrs with 10 years of total run-time (2.5 years in neutrino mode and 7.5 years in antineutrino mode). The magenta curves correspond to KD assuming the same exposure. We also estimate the results for JD+KD as shown by the blue curves. From the upper left panel, we observe that DUNE and JD+KD place similar constraints on α 11 . The sensitivity to this parameter comes from two contributions: disappearance of intrinsic ν e beam and the ν e appearance. For both of these contributing channels, DUNE has better systematics and JD+KD has more statistics. As a result, limits on α 11 is found to be almost the same for the two setups. However, for α 22 (upper middle panel), JD+KD has significantly better sensitivity compared to DUNE setup. This is because α 22 is mainly constrained by the disappearance channels, which due to the large statistics, is primarily limited by the systematic uncertainties. Since the normalization error for this channel is 3.5% (5%) for JD+KD (DUNE), it is clear that JD+KD can put better limit than DUNE. We have checked that if we consider the same amount of systematic uncertainties for both the setups, DUNE shows slightly better sensitivity than T2HKK.
For |α 21 |, DUNE and JD+KD have comparable sensitivities (see lower left panel). The slightly better limit on |α 21 | achieved in the case of JD+KD as compared to DUNE, is due to the fact that JD+KD has larger statistics in the appearance channels. For the other three parameters |α 31 |, |α 32 | and α 33 , which enter the ν µ → ν e appearance channel through matter parameters ∆ e and ∆ n (see Eq. 2.6), DUNE outperforms JD+KD setups because of its large matter effects.
We summarize our results in Table 6, where we give the bounds on the six NUNM parameters at 90% C.L. for the various long-baseline experimental setups discussed in   this paper. As clear from our previous discussion, the expected constraints on NUNM parameters from DUNE is better than the other two experiments JD and KD (and their combination) except for the parameters α 22 and |α 21 |, where JD has better sensitivity than DUNE. Finally, in the sixth column of the Table, we give the final constraints on the NUNM parameters by combining the expected results from DUNE and JD+KD setups. As we have anticipated, the bounds experience a general improvements by ∼ 20%, with the precise magnitude depending on the parameter under consideration.
For a comparison with the ongoing long-baseline experiments, we also add the expected constraints from the combination of the T2K and NOνA setup in the last column. For T2K, we consider a total exposure of 84.4 kt-MW-yrs with 22.5 kt detector mass, 750 kW proton beam power, 5 years run-time (2.5 years each for neutrino and antineutrino mode). For NOνA, the considered exposure is 58.8 kt-MW-yrs with 14 kt detector mass, 700 kW proton beam power, and 6 years for total run-time (3 years each for neutrino and antineutrino modes). Due to the limited statistics, we observe that the expected constraints from T2K+NOνA setup is worse than the DUNE or JD+KD setup.
Note that if some information coming from the near detector (for an example, measurement of the initial neutrino flux) are used to analyze the far detector data then the constraints on α 11 and α 22 may be modified (see Sec. 6 for a detailed discussion). However, in this section, we adopt a different strategy, where we simulate the far detector data alone to set limits on the NUNM parameters. In principle, this approach is valid if the initial flux can be predicted by some theoretical calculation or measured at an experiment which is insensitive to neutrino oscillation phenomena. In fact, the MINOS/MINOS+ experiment adopted this strategy where the oscillation data at the far detector was analyzed using information from the MINERvA flux predictions [81]. In future, if somehow we can apply this approach for DUNE and T2HKK, then in principle, one can estimate the limits on the NUNM parameters using only the far detector data.
At the same time, we understand that the assumptions that we take in our paper for the systematic uncertainties at the far detector may be too optimistic if we do not use the near detector to measure the initial flux. To have a better understanding on this issue, we perform some study and find that limits on α 11 and α 22 are mainly governed by the systematic uncertainties. In other words, the expected constraints on these two NUNM parameters are proportional to the systematic uncertainties that we consider in our simulation. Therefore, it may be possible to predict what would be the limits on α 11 and α 22 for a given assumptions on systematic uncertainties.
We compare our results summarized in Table 6 with the bounds reported in Table 1 4 . We observe that the bound we achieve from the DUNE+JD+KD (DUNE+T2HKK) setup for the diagonal α 11 is ∼ 80% better than the bound quoted in Ref. [53]. In the α 22 case, the two results are comparable, with a slightly better limit when the global neutrino data analysis is performed. For the remaining diagonal parameter α 33 , NC data from MINOS/MINOS+ give a 60% stronger bound [53] compared to the one expected from the DUNE+JD+KD setup. As for α 21 , the authors of the Ref. [53] make use of the triangular inequality as well as the data from the short-baseline experiments; this allows to constrain the mentioned parameter very tightly. However, due to the large statistics and good systematics of DUNE and JD+KD setups, we can achieve an almost similar bound without using any external hypothesis on the relations between the α ij . On the other hand, constraints on α 31 and α 32 in Table 1 are substantially better than the ones we obtain from DUNE+JD+KD setup. Also, in these cases, the triangular inequalities which link them to the diagonal NUNM parameters play an important role, together with the short-baseline experiments limits on the ν τ appearance. However, it is important to stress that all our results are obtained in a complete model independent fashion, relying only on the expected data from DUNE and T2HKK. We check that for our best setup, namely DUNE+T2HKK, the only parameter whose limit gets improved when we consider these inequalities is |α 32 | because of the stringent bound on α 22 . Since |α 21 | is already tightly constrained, we do not see any improvement in its limit using these inequalities. As far as the bound on |α 31 | is concerned, since the limit on the diagonal parameter α 33 is very poor, we also do not see any improvement. 4 In order to get their results, the authors of Ref. [40] left free the standard oscillation parameters θ23, δCP, and ∆m 2 31 and the NUNM parameters α11, α21, α22. Conversely, in our work we marginalize over δCP and θ23 only, but we have checked that the marginalization over ∆m 2 31 does not have any significant impact. In Appendix B, we also show that the marginalization over α11, α21, and α22 does not cause remarkable change in the results.
Recently, the DUNE collaboration [27] exploited the possibility of increasing the exposure of the experiment from 336 kt-MW-yrs to 480 kt-MW-yrs (corresponding to an increase of the data taking time from 7 years to 10 years with 5 years in neutrino mode and 5 years in antineutrino mode). In Table 7, we compare our previous constraints from the DUNE experiment, Table 6, with those obtained in the (5 + 5) years configuration. We observe that the constraints on all six NUNM parameters improve by small amount except for |α 21 |, which shows a significant improvement. This happens because the higher runtime increases statistics of the ν µ → ν e appearance channel, which is the one driving the α 21 sensitivity. On the other hand, the ν µ → ν µ disappearance channel is almost already saturated by systematics after 3.5 years + 3.5 years of running. This leads to only small improvements on the other NUNM parameters sensitivities.

Benefits of having near detectors
Near detectors (ND) are a fundamental component for long-baseline neutrino experiments. Indeed, a detector placed very close to the beam source (from hundreds of meters to a few kilometers) is able to monitor the neutrino beam, measuring the flavor composition, and the total number of neutrinos emitted from the source. Near detectors are not expected to improve any of the standard oscillation parameter measurements, since at such short distances, oscillations do not develop for neutrinos with energies in the GeV range. However, in some new physics scenarios, in which, oscillation probabilities contain zerodistance terms, near detectors can be used to constrain non-standard parameters in a very straightforward way. This is the case of the non-unitarity framework under discussion where, as already mentioned in Sec. 2, at vanishing distances we have zero-distance terms in case ν µ → ν e appearance channel: P L=0 µe ∼ |α 21 | 2 , and ν µ → ν µ disappearance channel: P L=0 µµ ∼ 1 + 2|α 21 | 2 + 6α 2 22 + 4α 22 . Thus, we can expect that T2HKK and DUNE near detectors would be able to constrain two parameters |α 21 | and α 22 from ν µ → ν e appearance and ν µ → ν µ disappearance channel, respectively, but also α 11 considering the ν e beam contamination (see Eq. A1). So, in this section, we analytically infer the order of magnitude of bounds implied by ND measurements. Let us consider the total number of events of a given channel as [82]: where, the normalization factor N 0 includes all the detector properties. For an oscillation channel ν α → ν β , N 0 can be defined as: where, σ β denotes the production cross-section of the β lepton, ε β represents the detector efficiency, and φ α stands for the initial neutrino flux of flavor α. If we want to put bounds on new physics parameters, we can use a simple χ 2 test with a gaussian χ 2 defined as   where, σ represents the uncertainty on the number of events; in this case, neglecting the backgrounds, we get: For the ν µ → ν µ disappearance channel, the leading term of the probability is P L=0 µµ = 1 + 4α 22 . Therefore, the χ 2 assumes the form: At a chosen confidence level, represented by a cut at χ 2 0 , it is possible to exclude the region satisfying: Since the disappearance channel is expected to produce a huge number of events at the near detector, one can consider the uncertainty to be dominated by systematic errors σ sys . Thus, it is possible to approximate σ ∼ N 0 σ sys , where N 0 represents the number of events in absence of zero-distance effects, being the disappearance probability in that case equal to 1. This allows to simplify Eq. 6.6 as follows: which tells us that, neglecting backgrounds effects, the near detector limits would be of the order of the chosen systematic uncertainty. A similar approach can be used for the ν e → ν e oscillation channel (see Eq. A1), which arises from the ν e beam contamination, obtaining an inequality for |α 11 | of the similar form: where σ νe sys refers to the systematic uncertainty on the ν e → ν e transition. For the appearance channel, the zero-distance probability reads P L=0 µe = |α 21 | 2 ; the χ 2 function can therefore be written as: and the excluded region is expected to be determined by the following relation: Since the number of events at the near detector is in principle very small (being only caused by new physics), the uncertainty is dominated by statistics. Thus, given a certain number of observed events, σ ∼ √ N obs and the excluded values of |α 21 | reduced to:  suggesting that the bounds are very sensitive to the number of events and to the running time of the experiment. Both DUNE and T2HKK will have near detectors [55,67,83] which may play a crucial role to probe various new physics scenarios including the possibility of non-unitarity of the PMNS matrix which is the main thrust of this work. In our analysis, for DUNE, we consider a 67 tons LArTPC near detector placed at a baseline of 574 meters from Fermilab [67]. For JD+KD, we consider a 1 kt water Cherenkov near detector located at a baseline of 1 km from J-PARC which is known as IWCD [75,76]. In order to simulate their responses, we scale the far detector fluxes for ND baselines and take into account their fiducial masses. We follow a very conservative approach as far as the systematic uncertainties at the near detectors are concerned. We multiply the FD systematic uncertainties by a factor of three and consider them as inputs for the ND. In DUNE near detector, we expect O(10 7 ) ν µ andν µ events, which provide bounds on α 22 . DUNE can place stringent constraints on α 11 and α 21 using O(10 6 ) ν e andν e events at ND, which stem from both intrinsic ν e (ν e ) beam contamination and via ν µ → ν e (ν µ →ν e ) appearance caused due to zero-distance effect. For the NDs, we consider their appropriate baselines, fiducial masses, and systematic uncertainties which we assume to be larger than the systematic uncertainties considered for the FDs.
Before discussing the limits that the near detectors would be able to set using their own data, we want to study the effect of the ND flux measurements on the far detector constraints. Indeed, if the initial neutrino flux is measured at the near detector and then extrapolated to the far detector, the probability which could be inferred at the far detector is the effective probability defined as The P L=0 αα term that appears in the denominator is the survival probability initial neutrino flavor at the source or the zero-distance term which act as a normalization factor. If we normalize the ν µ → ν µ survival probability at the far detector using the zero-distance term in Eq. 2.10, it is observed that the contribution from α 22 gets canceled at the leading order. As a result, sensitivity to the parameter α 22 is worsened for a given setup. The same happens for α 11 , whose contribution in the effective ν e → ν e disappearance probability is canceled at the leading order (see Eq. A1 and A4). Since the sensitivity to this parameter arises partially due to the intrinsic ν e that we have in the beam to begin with, the near detector normalization causes a deterioration of α 11 limits. In Table 8, we show how the constraints on α 11 and α 22 would be modified when taking into account the FD and ND correlation for the three setups namely, DUNE, JD+KD, and DUNE+JD+KD. We observe that the bound on α 11 is increased by a factor of almost two when we consider the correlation between the FD and ND. For α 22 , the bound is deteriorated at least three times compared to the FD case only. We have checked that no other NUNM parameter is affected significantly if we consider the FD and ND correlation. Indeed, the non-diagonal parameters and α 33 can be constrained using appearance channels (for which we do not have full cancellations in the effective probabilities) or using the interplay with matter effects, which are not developed at the near site.
The bounds obtained using the above-mentioned near detectors are shown in Fig. 7 for DUNE and T2HKK together with the results we got using the FD data with effective probabilities. For α 11 , NDs of the two setups can put bounds better than the one set by the FDs considering the ND normalization due to the very high statistics and the strong α 11 dependence of the zero-distance probability. The improvement is roughly a factor of two for DUNE and 60% for T2HKK (see Table 9). Note that the obtained limits are in agreement with the predictions deduced from Eq. 6.8, once we insert a normalization uncertainty of 6% (15%) for DUNE (JD+KD).
For the second diagonal parameter α 22 we also observe a similar situation, in which the ND alone can put more stringent bounds than the far detector when the normalization is considered, despite of the increased systematics. The improvement can be quantified as roughly 25% in DUNE and a factor of two in T2HKK. Once again, the analytical predictions from Eq. 6.8 are sufficiently recovered by the numerical simulations, considering that the near detectors normalization systematics are 15% for DUNE and 10.5% for T2HKK. Note that the bounds from the far detector data alone would be considerably better than the near detector ones.
Finally, for the NUNM parameter |α 21 | the near detectors bounds are considerably better than the far detector ones, due to the zero-distance effect outlined in Eq. 2.7. In particular, the limits are ∼ 3 times smaller than the one set by the far detector in the DUNE facility and ∼ 70% smaller in the case of T2HKK. Considering a number of observed events of O(10), and taking into account that we expect N 0 ∼ 10 6 per year [67], our analytic estimate for |α 21 | is comparable with the numerical results.
7 Improvement due to the ν τ sample in DUNE The ν τ production at accelerator experiments is very challenging since the charged current interactions of such particles with nuclei have an energy threshold of 3.1 GeV. Thus, many proposed long-baseline experiments are not able to detect such neutrinos 5 . However, the DUNE neutrino spectra will have peak at around 2.5 GeV (differently from the T2HKK where E ν 0.6 GeV) and the most energetic neutrinos of the beam will have enough energy to produce τ leptons. Recently, different studies [51,[85][86][87] take into account the possibility of including the ν τ (ν τ ) sample in the DUNE analysis.  The recognition of such events could be in principle possible due to the imaging capabilities of LArTPC detectors. Because of a relatively small number of neutrinos with an energy above the production threshold and of the short lifetime of the τ leptons which could make the recognition of the τ interaction and the decay vertices a difficult task, the ν τ appearance channel is not really useful to constrain the standard oscillation parameters. However, when new physics affects the oscillation probabilities, it has been shown in Refs. [51,86,88,89] that the NUNM parameters α 33 and |α 32 | constraints can be improved even by the small number of ν τ events. In case of non-unitary neutrino mixing, the ν µ → ν τ oscillation probability mainly depends on the three parameters |α 32 |, α 33 and α 22 (see Eq. A3 in Appendix A). While we expect that the sensitivity to α 22 will not be improved by the small number of ν τ events, the other two poorly constrained NUNM parameters |α 32 | and α 33 could take advantage of ν µ → ν τ oscillation channel. We include τ events in our analysis in the following fashion.
• For the hadronic decays of τ events having a branching ratio of 65%, a 30% signal efficiency has been assumed and 10% of the NC events are considered as background [51].
• For τ decaying to electron with a branching ratio of 17.4%, we assume a signal efficiency of 30% considering ν e events as possible background. We consider the signal to background ratio of 2.45 in our analysis [86].
The muonic decays have not been taken into account since the discrimination of the number of background events would be too large compared to the signal events [51,86]. The total number of ν τ (ν τ ) events in DUNE is expected to be roughly 72 (37) per year. The normalization error for the signal is taken to be 20%. The results of our analysis for α 33 and |α 32 | are shown in Fig. 8 and Table 10. The allowed range of α 33 , which appears only at the second order in the probability, is reduced of ∼ 13% by the inclusion of the new oscillation channel, and the new limits are set into the range [-0.16, 0.15] (see Table 10).
On the other hand, the sensitivity to |α 32 |, which impacts linearly the ν τ appearance probability, is significantly improved: in this case, the new upper bound is roughly 60% smaller than the one set by the standard oscillation channels, namely |α 32 | < 0.19.

Summary and conclusions
Neutrino oscillation has been one of the most important discoveries over the last few decades and a large number of pioneering neutrino experiments have confirmed this phenomena of mass-induced flavor transition. Excellent data from various neutrino oscillation experiments having different baselines and energies have enabled us to achieve remarkable precision on almost all the three-flavor oscillation parameters. However, some of the features of neutrino mass-mixing parameters, like, neutrino mass hierarchy, the value of the CPviolating phase δ CP , and octant of θ 23 are still poorly known and next-generation neutrino experiments will play an important role to address these issues with high confidence level. Given the magnificent precision on neutrino mixing angles and rapidly increasing knowledge on δ CP , it seems quite natural to ask if there is any violation of the unitary property of 3 × 3 PMNS mixing matrix, which may be related to the existence of new mass eigenstates of neutrinos. In this context, there exists one well-known parameterization in the literature, which takes into account the non-unitary neutrino mixing (NUNM) by introducing a lower triangular matrix with three real and three complex parameters, α ij .
To study the impact of these NUNM parameters on various oscillation channels, for the first time, we derive simple approximate analytical expressions for the oscillation probabilities in matter in the atmospheric regime (∆ 31 >> ∆ 21 ). Our perturbative expansions are valid up to second order in the small deviations from the mixing angles by their tribimaximal values, second order in the NUNM parameters, and first order in the matter potential. We have shown that the ν µ → ν e appearance probability (see Eq. 2.6) mainly depends on |α 21 |, but when matter potential is large, the impact of |α 31 | and |α 32 | can also be significant. On the other hand, the ν µ → ν µ disappearance probability (see Eq. 2.9) mainly relies on α 22 , but sub-leading dependencies on |α 21 | in vacuum and α 33 , |α 31 |, and |α 32 | in matter are also present. The only parameter that does not appear in our formulae is α 11 , which is only relevant at the higher-orders in our perturbative expansions as shown in Ref. [37] for the vacuum case.
In this work, we analyze in detail, the impact of possible NUNM in the context of longbaseline experiments DUNE and T2HK having one detector in Japan (JD) and a second detector in Korea (KD), and the combination of these two detectors, popularly known as T2HKK or JD+KD. First, we show how the various NUNM parameters (α ij ) are correlated with the oscillation parameters θ 23 and δ CP for these setups. Then, we estimate in detail the sensitivities of these experiments to place direct, model-independent, competitive constraints on the six NUNM parameters at 90% C.L. (see Fig. 6 and Table 6). The wide-band neutrino beam in DUNE encompassing both first and second oscillation maxima allows us to measure the NUNM parameters at several L/E values in the presence of significant matter effect due to its 1300 km baseline. Indeed, DUNE shows better sensitivity than JD+KD in constraining the NUNM parameters, which are influenced by matter effects, namely, α 33 , |α 31 |, and |α 32 |. We observe that JD+KD will provide a stringent constraint on α 22 as compared to DUNE because it has less systematic uncertainties in the disappearance channel. Also, due to the larger statistics in the appearance channel, JD+KD will give significantly better limit on the NUNM parameter |α 21 | in comparison to DUNE. Lastly, α 11 is expected to be constrained almost in the same way by these two experiments. We show how the limits on the NUNM parameters get improved in case of DUNE, if the total exposure is increased from 336 kt-MW-yrs to 480 kt-MW-yrs (corresponding to an increase in the total run-time from 7 years to 10 years), as proposed in the recent TDR [27]. We also estimate how much the sensitivities can be improved by adding the prospective data from DUNE and JD+KD. Finally, we compare our results with the constraints that can be achieved using the full exposure of currently running experiments T2K and NOνA.
Due to the so-called zero-distance effects which are induced by the non-unitary neutrino mixing in neutrino oscillation probabilities, the prospective data from near detectors in both DUNE and JD+KD experiments could be in principle used to bound the three NUNM parameters |α 21 |, α 11 and α 22 . However, the zero-distance effect should also be taken into account if the near detectors data will be used to measure the initial neutrino flux for both experiments. This would lead to a substantial deterioration of the limits that the far detectors could set on α 11 and α 22 , as summarized in Table 8.
Moreover, in DUNE, the expected limits on |α 32 | and α 33 get improved by O(20)% when we also add the ν τ appearance sample in our analysis.
Long-baseline experiments have contributed significantly in our journey to establish the three-neutrino paradigm on a strong footing. Future high-precision long-baseline experiments such as DUNE and T2HKK (JD+KD) will continue this legacy and rigorously test the unitarity of the PMNS matrix by measuring the mixing angles and CP phase δ CP with remarkable precision. We hope our present work will provide further boost along this direction. In this appendix, we provide a check of the goodness of our expanded probabilities, Eqs. 2.6 to 2.9. To this aim, in the upper panels of Fig. A1, we show the ν µ → ν e appearance probabilities as a function of neutrino energy considering the NUNM parameters |α 21 | and |α 31 | one at a time. We portray the same in the lower panels for ν µ → ν µ disappearance channel in presence of the diagonal NUNM parameter α 22 . We consider three different baselines, corresponding to DUNE (left panels), JD (middle panels), and KD (right panels). In each panel, solid curves are obtained numerically from the full (not expanded) Hamiltonian, while dashed curves correspond to the oscillation probabilities obtained analytically using Eq. 2.6 and Eq. 2.9. The value of the CP phase δ CP and phases φ 21 , φ 31 are fixed to zero, while the values of the other standard oscillation parameters are taken from Table 2. As we can see, the agreement between the full solid and the analytical dashed curves is generally very good for both probabilities.
Perturbative expansion of P ee , P eτ , and P µτ For the sake of completeness, we give here the perturbative expressions of the probabilities P ee , P eτ and P µτ , using the same approximations discussed in Sec. 2. We start with the ν e → ν e disappearance channel:  The ν e → ν τ transition is governed by the following expression: Eventually, the ν τ appearance from a ν µ beam is regulated by the following probability expression: It is easy to check that, when all α parameters are set to zero, we recover the unitary relation P µe + P µµ + P µτ = 1.

B Impact of marginalization over the other NUNM parameters while showing the constraints on α ij
In this appendix, we derive the constraints on the NUNM parameters obtained, differently from the procedure adopted in the main text, marginalizing over the other NUNM parameters α ij . Our numerical results are shown in Fig. B1, where we show the impact of the marginalization over various NUNM parameters on the sensitivity of the DUNE+JD+KD  Table B1. We observe that the marginalization over all other NUNM parameters worsen the bounds on all of them but α 22 . This is due to the fact that, as it can be seen in Eq. 2.9, α 22 appears at the leading order, with no correlations to the other α ij . The other diagonal parameter α 11 shows only a marginal deterioration at the level of 15% since, as before, its correlation with the other NUNM parameters are mild in ν e appearance probability as shown in Ref. [37]. We see from Table B1 that there is a considerable deterioration in the sensitivities for α 33 , |α 21 |, |α 31 |, and |α 32 | when we marginalize over the un-displayed NUNM parameters in the fit. This happens because all these four parameters are strongly correlated among them and with other two NUNM parameters α 11 Table 2. We marginalize over θ23 and δCP in the fit (see text for details). and 2.9). Note that |α 21 | does not have any correlation with the other NUNM parameters in ν µ → ν e appearance channel but it is strongly correlated with α 22 , |α 31 |, and |α 32 | in ν µ → ν µ disappearance channel which reduces its sensitivity by roughly 40%, when we keep the other NUNM parameter free in the fit. Similar correlations among the NUNM parameters are also responsible for worsening the bounds on |α 31 | by 60%, |α 32 | by 45%, and by a factor of two for α 33 when we marginalize over the other NUNM parameters in the fit.