Searching for non-unitary neutrino oscillations in the present T2K and NO$\nu$A data

Abstract

The mixing of three active neutrino flavors is parameterized by the unitary PMNS matrix. If there are more than three neutrino flavors and if the extra generations are heavy iso-singlets, the effective $3\times 3$ mixing matrix for the three active neutrinos will be non-unitary. We have analyzed the latest T2K and NO$\nu$A data with the hypothesis of non-unitary mixing of the active neutrinos. We found that the 2019 NO$\nu$A data slightly (at $\sim 1\, \sigma$ CL) prefer the non-unitary mixing over unitary mixing. In fact, allowing the non-unitary mixing brings the NO$\nu$A best-fit point in the $\sin ^2{\theta _{23}}-\delta _{\mathrm {CP}}$ plane closer to the T2K best-fit point. The 2019 T2K data, on the other hand, cannot rule out any of the two mixing schemes. A combined analysis of the NO$\nu$A and T2K 2019 data prefers the non-unitary mixing at $1\, \sigma$ CL. We derive constraints on the non-unitary mixing parameters using the best-fit to the combined NO$\nu$A and T2K data. These constraints are weaker than previously found. The latest 2020 data from both the experiments prefer non-unitarity over unitary mixing at $1\, \sigma$ CL. The combined analysis prefers non-unitarity at $2\, \sigma$ CL. The stronger tension, which exists between the latest 2020 data of the two experiments, also gets reduced with non-unitary analysis.

1 Introduction

Neutrino oscillation phenomena has been well established by the solar [1, 2], atmospheric [3] and reactor [4,5,6] neutrino experiments. The oscillations between 3-neutrino flavors depend on three mixing angles ${\theta _{12}}$, ${\theta _{13}}$ and ${\theta _{23}}$; two independent mass-squared differences ${\Delta _{21}}=m_{2}^{2}-m_{1}^{2}$ and ${\Delta _{31}}=m_{3}^{2}-m_{1}^{2}$; and a CP violating phase $\delta _{\mathrm {CP}}$. The angle ${\theta _{12}}$ and mass-squared difference ${\Delta _{21}}$ were measured in the solar neutrino experiments [1, 2], $|{\Delta _{31}}|$ was measured in the accelerator neutrino experiment MINOS [7], and ${\theta _{13}}$ was measured in the reactor neutrino experiments [4,5,6]. However, the sign of ${\Delta _{31}}$, the octant of ${\theta _{23}}$ and the value of $\delta _{\mathrm {CP}}$ are still unknown. Currently operating long baseline accelerator neutrino experiments, namely T2K [8] and NO$\nu$A [9], are expected to measure these unknown quantities. Both the experiments published their data in 2018 and 2019 [10,11,12,13]. According to an analysis of these data sets, the T2K best-fit point is at $\sin ^2 {\theta _{23}}\!=\!0.53^{+0.03}_{-0.04}$ for both hierarchies, and $\delta _{\mathrm {CP}}/\pi \!=\!-1.89^{+0.70}_{-0.58}$ ($-1.38^{+0.48}_{-0.54}$) for normal (inverted) hierarchy [13]. On the other hand the NO$\nu$A best-fit point is at $\sin ^2 {\theta _{23}}\!=\! 0.56^{+0.04}_{-0.03}$, and $\delta _{\mathrm {CP}}/\pi =0^{+1.3}_{-0.4}$ for normal hierarchy. Therefore, there is a visible difference between the $\delta _{\mathrm {CP}}$ values at the NO$\nu$A and T2K best-fit points. A previous result from the NO$\nu$A collaboration [14,15,16] had some mild tension with the T2K data, which has been discussed in Ref. [17]. In the 2019 analysis with the NO$\nu$A data, the best-fit point is now closer to the T2K best-fit point, but differences still exist as NO$\nu$A disfavors the T2K best-fit point at $1\, \sigma$ CL and vice versa. Recently NO$\nu$A [18] and T2K [19] have published their latest results in the Neutrino, 2020 conference. As per the latest analysis, the best-fit point for NO$\nu$A (T2K) is $\sin ^2{\theta _{23}}= 0.57$ (0.528) and $\delta _{\mathrm {CP}}=0.82 \pi$ ($-1.6\pi$). The tension between the two experiments is even stronger as they exclude each other’s allowed region at $1\, \sigma$ CL. It has also been shown that although both the experiments individually prefer NH over IH, their combined analysis prefers IH over NH [20].

Apart from the unknowns in the three flavor neutrino mixing, there are anomalies from the short baseline experiments, which cannot be accounted for in the three flavor oscillation formalism. These anomalies are namely

• 1. Reactor anomalies: It implies a deficit of observed $\bar{\nu }_e$ event numbers in different detectors situated at a few meters away from the reactor sources, compared to the predicted number. In particular, the average ratio is $R_{\mathrm{avg}}= N_\mathrm{obs}/N_{\mathrm{pred}}=0.927\pm 0.023$ [21]. Recent updates have changed the predictions slightly, giving an average ratio $R_{\mathrm{avg}}=0.938\pm 0.023$ [22], which is a $2.7\, \sigma$ deviation from unity. However, there is a lack of knowledge of the reactor neutrino fluxes and a detailed study of the forbidden transition in the reactor neutrino spectra may increase the systematic uncertainties to a few percentage. Moreover, there are similar indications of $\nu _e$ disappearance from the SAGE [23] and GALLEX [24] solar neutrino experiments. A combined analysis of the detected and predicted number of neutrino events from the source gives $R=0.86\pm 0.05$ [25, 26], another $2.7\, \sigma$ deviation from unity. Both of these deficits of low energy $\bar{\nu }_e$ events can be explained by an oscillation at $\Delta m^2 \ge 1\, \mathrm{eV}^2$ over very short baseline.

• 2. LSND and MiniBooNE anomalies: The LSND experiment [27] at the Los Alamos National Laboratory was designed to observe $\bar{\nu }_{\mu }\rightarrow \bar{\nu }_e$ oscillations over a baseline of 30 m. After 5 years of data taking, it observed $89.7 \pm 22.4 \pm 6.0$ $\bar{\nu }_e$ candidate events over background, providing a $3.8\, \sigma$ evidence of $\bar{\nu }_{\mu }\rightarrow \bar{\nu }_e$ oscillations at $\Delta m^2 = 1\, \mathrm{eV}^2$ region. Therefore, this result cannot be accommodated in the three flavor scenario. The MiniBooNE experiment [28] at the Fermilab was designed to observe $\nu _{\mu } \rightarrow \nu _e$ and $\bar{\nu }_{\mu }\rightarrow \bar{\nu }_e$ oscillations over 540 m baseline using the Booster Neutrino Beam (BNB), a predominantly muon-neutrino beam, peaking at 700 MeV. It observed a $3.4\, \sigma$ signal excess of $\nu _e$ candidate and $2.8\, \sigma$ signal excess of $\bar{\nu }_e$ candidates.

The most common explanation of these anomalies is the existence of one or more “sterile” neutrino states with mass at or below a few eV range, see Ref. [29] for a comprehensive review. The minimal model consists of $3+1$ neutrino mixing, dominated by $\nu _e$, $\nu _{\mu }$ and $\nu _\tau$, with very small perturbative contribution from the new sterile flavor $\nu _x$. The $\nu _x$ mainly consists of a very heavy eigenstate $\nu _4$ with mass $m_4$, such that $m_1,\, m_2,\, m_3\ll m_4$ and $\Delta _{41} = m_{4}^{2}-m_{1}^{2}=[0.1-10]\, \mathrm{eV}^2$. Recent results from the IceCube experiment constrain the sterile neutrino mass and mixing using atmospheric neutrino fluxes [30]. In 2018, however, MiniBooNE has again confirmed a $4.7\, \sigma$ excess of combined $\nu _e$ and $\bar{\nu }_e$ events [31]. The present significance of the excess from a combined analysis of MiniBooNE and LSND is $6\, \sigma$. Constraints on the existence of sterile neutrino have been discussed in Refs. [32,33,34], while Refs. [35,36,37,38,39] discuss about the effects of light sterile neutrino on present and future long baseline experiments.

If extra neutrino generations exist as iso-singlet neutral heavy leptons (NHL), in the minimal extension of the standard model, they would not take part in neutrino oscillations, however. The admixture of such leptons in the charged current weak interactions would affect the neutrino oscillation, and the neutrino oscillation would be described by an effective $3\times 3$ non-unitary mixing matrix [40]. NHL would induce charged lepton flavor violation processes [40, 41]. If Majorana type, these NHLs will modify the rate of neutrinoless double beta decay [42,43,44]. The theory of neutrino oscillation in the presence of non-unitarity in the 3-generation scheme and its effect on long baseline accelerator neutrino, have been studied in several references. Reference [45] has studied the CP violation measurement potential of T2K and future experiment T2HK in presence of non-unitary oscillation, whereas a similar study in the context of DUNE has been done in Ref. [46]. In Ref. [47], physics potential of long baseline experiments with non-unitary mixing has been studied, while Ref. [48] has made a theoretical study of the non-unitary oscillation at the probability level. In Ref. [49], the CP violation measurement potential of short baseline (SBL) experiments in the presence of non-unitary $\nu _{\mu }\rightarrow \nu _\tau$ oscillation has been studied. All these works were done with simulations for the future experiments, but they did not analyse already available data, e.g., from the T2K and NO$\nu$A experiments as we have done in this work. However, in Refs. [50, 51] an effort has been made to resolve the tension between NO$\nu$A and T2K with non-standard NC interaction during propagation and in Ref. [52] the same has been done with CPT conserving Lorentz invariance violation. A global analysis assuming non-unitary hypothesis has been done and limits on non-unitary parameters have been given in Ref. [53].

In this paper, we have explored whether the 2019 and 2020 T2K and NO$\nu$A data can exclude the non-unitary $3\times 3$ mixing, and if not, whether the non-unitary $3\times 3$ mixing hypothesis can lead to better agreement between the T2K and NO$\nu$A data. A combined analysis with the non-unitary hypothesis has also been done. In Sect. 2, we have discussed the theory of neutrino oscillation probability in the presence of a non-unitary $3\times 3$ mixing matrix. Details of the simulation method have been discussed in Sect. 3. In Sects. 4 and 5, we have presented our results for the 2019 and 2020 data respectively, and the conclusion has been drawn in Sect. 6.

2 Non-unitary oscillation probabilities

In the standard case of 3 active neutrinos, the flavor basis $\nu _f$ ($f=e, \mu , \tau$) is related to the mass basis $\nu _m$ ($m = 1, 2, 3$) by the relation $\nu _f = U\nu _m$, where U is the unitary PMNS matrix. The Schrödinger equation for neutrino propagation in matter can be written as

$$\begin{aligned} i\frac{d}{dt}\nu _m = (H + U^\dagger {\mathcal {A}} U)\nu _m, \end{aligned}$$

(1)

where H is the Hamiltonian and ${\mathcal {A}}$ is the matter potential. The solution to the above equation, after propagation over a distance L in matter, is $\nu _m (L) = S_m(L)\nu _m$, where

$$\begin{aligned} S_m(L) = W e^{-i L {\tilde{E}}} W^\dagger = e^{-i LW {\tilde{E}} W^\dagger } \end{aligned}$$

(2)

is the S-matrix in the mass basis, with W being the transformation matrix between the mass basis $\nu _m$ and the mass basis in matter ${\tilde{\nu }}_m$ such that $\nu _m = W{\tilde{\nu }}_m$. The energies ${\tilde{E}} = diag(\tilde{E_1}, \tilde{E_2}, \tilde{E_3})$, where $\tilde{E_i}$ are the eigenvalues in the ${\tilde{\nu }}_m$ basis. The S-matrix in the flavor basis can be found from the mass basis in Eq. (2) by using the unitary property of the PMNS matrix as

$$\begin{aligned} S_f(L) = US_m(L)U^\dagger = e^{-iLUW{\tilde{E}}W^\dagger U^\dagger }. \end{aligned}$$

(3)

Flavor change in terms of the S-matrix is $\nu _f (L) = S_f(L)\nu _f$. Correspondingly, the oscillation probability from flavor $\alpha$ to flavor $\beta$ can be written as

$$\begin{aligned} P_{\alpha \beta }=|(S_f(L))_{\beta \alpha }|^2 . \end{aligned}$$

(4)

Fig. 1

Comparison between unitary and non-unitary $3\times 3$ mixing of $\nu _{\mu } \rightarrow \nu _e$ ($\nu _{\mu } \rightarrow \nu _{\mu }$) probabilities in the upper (lower) panel for the NO$\nu$A experiment (810 km baseline). Left (right) panel is for normal hierarchy (inverted hierarchy). The unitary CP-violating phase $\delta _{\mathrm {CP}}$ has been varied in its total range $[-180^\circ :180^\circ ]$. Non-unitary complex phase $\phi _{10}=0$. Non-unitary parameters are fixed at their boundary values in Eq. (18) taken from Ref. [46]

Fig. 2

Same as Fig. 1 but for the T2K experiment (295 km baseline)

One can extend this formalism to $n\times n$ unitary mixing matrix $U_{n\times n}$ with 3 active neutrinos and $(n-3)$ heavy singlet neutrinos [53]. In that case,

$$\begin{aligned} U_{n\times n}= \left[ {\begin{array}{cc} N_{3\times 3} &{} Q_{3\times (n-3)} \\ V_{(n-3)\times 3} &{} T_{(n-3)\times (n-3)} \end{array} } \right] , \end{aligned}$$

(5)

where N is a $3\times 3$ matrix in the light neutrino sector; Q and V depict the coupling parameters of the extra iso-singlet state, expected to be heavy; and T is a $(n-3)\times (n-3)$ matrix in the heavy neutrino sector. The interaction potential matrix ${\mathcal {A}}$ can be written as

$$\begin{aligned} {\mathcal {A}}_{n\times n}= \left[ {\begin{array}{cc} {\mathcal {A}}_{3\times 3} &{} 0 \\ 0 &{} {\mathcal {A}}_{(n-3)\times (n-3)} \end{array} } \right] , \end{aligned}$$

(6)

where, for the usual matter potential, the term ${\mathcal {A}}_{3\times 3}$ is the same as in Eq. (1) and ${\mathcal {A}}_{(n-3)\times (n-3)} = 0$. An explicit formalism for the potential matrix has been discussed in details in Ref. [46]. The Hamiltonian in vacuum in this case can be written as

$$\begin{aligned} H_{n\times n}= \left[ {\begin{array}{cc} H_{3\times 3} &{} 0 \\ 0 &{} H_{(n-3)\times (n-3)} \end{array} } \right] , \end{aligned}$$

(7)

where $H_{3\times 3}$ is the Hamiltonian in vacuum from Eq. (1). Similarly, the neutrino flavor vector can be rewritten as

$$\begin{aligned} \nu _{f_{n}}= \left[ {\begin{array}{c} \nu _{f_{3}} \\ \nu _{f_{(n-3)}} \end{array} } \right] , \end{aligned}$$

(8)

where we split the flavor vector in active $\nu _{f_{3\times 1}}$ and heavy $\nu _{f_{(n-3)\times 1}}$ parts. The W matrix similar to that in Eq. (2) can be written as

$$\begin{aligned} W_{n\times n}= \left[ {\begin{array}{cc} W_{3\times 3} &{} W_{3\times (n-3)} \\ W_{(n-3)\times 3} &{} W_{(n-3)\times (n-3)} \end{array} } \right] , \end{aligned}$$

(9)

and the ${\tilde{E}}$ matrix from Eq. (2) can be written as

$$\begin{aligned} {\tilde{E}}_{n\times n}= \left[ {\begin{array}{cc} {\tilde{E}}_{3\times 3} &{} 0 \\ 0 &{} {\tilde{E}}_{(n-3)\times (n-3)} \end{array} } \right] . \end{aligned}$$

(10)

The Hamiltonian in the presence of matter potential, following Eq. (1), is

$$\begin{aligned} H_{m_{n\times n}}= & {} H_{n\times n} + U_{n\times n}^\dagger {\mathcal {A}}_{n\times n} U_{n\times n} \nonumber \\= & {} \left[ {\begin{array}{cc} H_{3\times 3}+N_{3\times 3}^\dagger {\mathcal {A}}_{3\times 3} N_{3\times 3} \,+ &{} N_{3\times 3}^\dagger {\mathcal {A}}_{3\times 3} Q_{3\times (n-3)} \,+ \\ V_{(n-3)\times 3}^\dagger {\mathcal {A}}_{(n-3)\times (n-3)}V_{(n-3)\times 3} &{} \, V_{(n-3)\times 3}^\dagger {\mathcal {A}}_{(n-3)\times (n-3)}T_{(n-3)\times (n-3)} \\ \\ Q_{3\times (n-3)}^\dagger {\mathcal {A}}_{3\times 3} N_{3\times 3} \, + &{} H_{(n-3)\times (n-3)}+Q_{3\times (n-3)}^\dagger {\mathcal {A}}_{3\times 3} Q_{3\times (n-3)}\, + \\ T_{(n-3)\times (n-3)}^\dagger {\mathcal {A}}_{(n-3)\times (n-3)}V_{(n-3)\times 3} &{} T_{(n-3)\times (n-3)}^\dagger {\mathcal {A}}_{(n-3)\times (n-3)}T_{(n-3)\times (n-3)} \end{array} } \right] .\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \end{aligned}$$

(11)

Similarly, the S-matrix in the mass basis, following Eq. (2), is

$$\begin{aligned} S_{m_{n\times n}}(L)= & {} W_{n\times n}e^{-iL{\tilde{E}}}W_{n\times n}^\dagger \nonumber \\= & {} \left[ {\begin{array}{cc} W_{3\times 3}e^{-iL{\tilde{E}}_{3\times 3}}W_{3\times 3}^\dagger \, + &{} W_{3\times 3}e^{-iL{\tilde{E}}_{3\times 3}}W_{(n-3)\times 3}^\dagger \, + \\ W_{3\times (n-3)} e^{-iL{\tilde{E}}_{(n-3)\times (n-3)}} W_{3\times (n-3)}^\dagger &{} W_{3\times (n-3)} e^{-iL{\tilde{E}}_{(n-3)\times (n-3)}} W_{(n-3)\times (n-3)}^\dagger \\ \\ W_{(n-3)\times 3}e^{-iL{\tilde{E}}_{3\times 3}}W_{3\times 3}^\dagger \, + &{} W_{(n-3)\times 3}e^{-iL{\tilde{E}}_{3\times 3}}W_{(n-3)\times 3}^\dagger \, + \\ W_{(n-3)\times (n-3)}e^{-iL{\tilde{E}}_{(n-3)\times (n-3)}} W_{3\times (n-3)}^\dagger &{} W_{(n-3)\times (n-3)}e^{-iL{\tilde{E}}_{(n-3)\times (n-3)}} W_{(n-3)\times (n-3)}^\dagger \end{array} } \right] . \end{aligned}$$

(12)

The masses of the heavy right handed neutrinos are expected to be several orders of magnitude larger than the masses of the active neutrinos and the potential terms. Under this condition the elements $W_{3 \times (n-3)}$ and $W_{(n-3)\times 3}$ satisfy conditions that they are $\ll 1$, following the typical see-saw mechanisms [54, 55]. Therefore, we can neglect these terms and write the S-matrix as

$$\begin{aligned} S_{m_{n\times n}}(L)\approx & {} \left[ {\begin{array}{cc} W_{3\times 3}e^{-iL{\tilde{E}}_{3\times 3}}W_{3\times 3}^\dagger &{} 0\\ 0 &{} W_{(n-3)\times (n-3)} e^{-iL{\tilde{E}}_{(n-3)\times (n-3)}} W_{(n-3)\times (n-3)}^\dagger \end{array} } \right] \nonumber \\\equiv & {} \left[ {\begin{array}{cc} S_{m_{3\times 3}} &{} 0\\ 0 &{} S_{m_{(n-3)\times (n-3)}} \end{array} } \right] . \end{aligned}$$

(13)

The S-matrix in the flavor basis is written, following Eq. (3), as

$$\begin{aligned} S_{f_{n\times n}}(L)= & {} U_{n\times n} S_{m_{n\times n}}(L)U_{n\times n}^\dagger \nonumber \\= & {} \left[ {\begin{array}{cc} N S_{m_{3\times 3}} N^\dagger + Q S_{m_{(n-3)\times (n-3)}} Q^\dagger &{} N S_{m_{3\times 3}} V^\dagger + Q S_{m_{(n-3)\times (n-3)}} T^\dagger \\ V S_{m_{3\times 3}} N^\dagger + T S_{m_{(n-3)\times (n-3)}} Q^\dagger &{} V S_{m_{3\times 3}} V^\dagger + T S_{m_{(n-3)\times (n-3)}} T^\dagger \end{array} } \right] . \end{aligned}$$

(14)

In case of 3 active neutrinos the flavor changes, following Eq. (8), as

$$\begin{aligned} \nu _{f_{3}}(L)= & {} (N S_{m_{3\times 3}} N^\dagger + Q S_{m_{(n-3)\times (n-3)}} Q^\dagger )\nu _{f_3} \nonumber \\&+ (N S_{m_{3\times 3}} V^\dagger + Q S_{m_{(n-3)\times (n-3)}} T^\dagger ) \nu _{f_{(n-3)}}. \end{aligned}$$

(15)

Here $S_{m_{3\times 3}}$ is the same as in Eq. (2). With a similar see-saw mechanism argument due to the connection between the mass and flavor bases, the elements of Q and V are much lower compared to the terms containing N and T matrices [53, 55]. As a good approximation, we can therefore eliminate terms with Q and V matrices and write the non-unitary solution for 3-active neutrino flavor states as

$$\begin{aligned} \nu _{f_{3}}(L) \approx N e^{-iLW_{3\times 3}{\tilde{E}}_{3\times 3}W_{3\times 3}^{\dagger }} N^\dagger \nu _{f_3} \equiv S_f(L)\nu _{f_3}. \end{aligned}$$

(16)

The corresponding oscillation probability can be calculated using Eq. (4).

We use the following parametrization of the non-unitary mixing matrix N [53]

$$\begin{aligned} N=N_{NP}U_{3\times 3}= \left[ {\begin{array}{ccc} \alpha _{00} &{} 0 &{} 0 \\ \alpha _{10} &{} \alpha _{11} &{} 0 \\ \alpha _{20} &{} \alpha _{21} &{} \alpha _{22} \end{array} } \right] U_{\mathrm{PMNS}} , \end{aligned}$$

(17)

where the diagonal term(s) must deviate from unity and/or the off-diagonal term(s) deviate from zero to allow non-unitarity effect. The oscillation experiments can probe non-unitarity only if the $\alpha$ parameters vary at least at the percent scale [40, 53, 56]. There exists severe constraint on the parameter $|\alpha _{10}| < 10^{-5}$ from non-observation of the $\mu \rightarrow e\gamma$ decay [56]. However, this bound can be relaxed in certain neutrino mass-generation models involving inverse or linear see-saw mechanism [40]. In this paper, to calculate the probabilities, we have kept the $\alpha$ parameters fixed at their $3\,\sigma$ boundary values [46]:

$$\begin{aligned}&\alpha _{00}> 0.93;\quad \alpha _{11}> 0.95;\quad \alpha _{22} > 0.61 \nonumber \\&|\alpha _{10}|< 3.6\times 10^{-2};\quad |\alpha _{20}|< 1.3\times 10^{-1};\nonumber \\&\qquad \qquad \qquad \qquad \qquad \qquad |\alpha _{21}| < 2.1\times 10^{-2}. \end{aligned}$$

(18)

The standard unitary parameter values have been fixed at their best-fit values [57].

For the standard unitary case, the $\nu _{\mu } \rightarrow \nu _e$ oscillation probability can be written as [58]

$$\begin{aligned} P_{\mu e}\simeq & {} \sin ^2 2 {\theta _{13}}\sin ^2 {\theta _{23}}\frac{\sin ^2{\hat{\Delta }}(1-{\hat{A}})}{(1-{\hat{A}})^2}\nonumber \\&+ \alpha \cos {\theta _{13}}\sin 2{\theta _{12}}\sin 2{\theta _{13}}\sin 2{\theta _{23}}\cos ({\hat{\Delta }}+\delta _{\mathrm {CP}})\nonumber \\&\times \frac{\sin {\hat{\Delta }}{\hat{A}}}{{\hat{A}}}\frac{\sin {\hat{\Delta }}(1-{\hat{A}})}{1-{\hat{A}}}, \end{aligned}$$

(19)

where $\alpha =\frac{{\Delta _{21}}}{{\Delta _{31}}}$, ${\hat{\Delta }}=\frac{{\Delta _{31}}L}{4E}$ and ${\hat{A}}=\frac{A}{{\Delta _{31}}}$. A is the Wolfenstein matter term [59], given by $A=2\sqrt{2}G_FN_eE$, where E is the neutrino beam energy and L is the length of the baseline. Anti-neutrino oscillation probability $P_{\bar{\mu } {\bar{e}}}$ can be obtained by changing the sign of A and $\delta _{\mathrm {CP}}$ in Eq. (19). The oscillation probability mainly depends on hierarchy (sign of ${\Delta _{31}}$), octant of ${\theta _{23}}$ and $\delta _{\mathrm {CP}}$.

From Eq. (19), we can see that for the standard unitary case, the oscillation probability $P_{\mu e}$ ($P_{\bar{\mu }{\bar{e}}}$) gets a double boost (double suppression) when the hierarchy is NH and $\delta _{\mathrm {CP}}$ is in the lower half plane (LHP). Similarly $P_{\mu e}$ ($P_{\bar{\mu }{\bar{e}}}$) gets a double suppression (double boost) when the hierarchy is IH and $\delta _{\mathrm {CP}}$ is in the upper half plane (UHP). Therefore $P_{\mu e}$ ($P_{\bar{\mu }{\bar{e}}}$) is maximum (minimum) for NH and $\delta _{\mathrm {CP}}$ in the LHP and minimum (maximum) for IH and $\delta _{\mathrm {CP}}$ in the UHP.

In the case of non-unitary mixing, an analytic expression for $\nu _{\mu } \rightarrow \nu _e$ oscillation with matter effect is not available in literature and it is difficult to calculate them. However, approximate analytic oscillation probability in vacuum with non-unitary mixing is given in different articles, see, e.g., [53]. From the expression of vacuum oscillation probability with non-unitary mixing, we can see that the parameters ij’s, where $i,\, j=0,\, 1, \, 2$, are multiplied with different unitary parameters. Since $\alpha _{ij}<1$, we can expect them to reduce the oscillation probability in case of non-unitary mixing. This feature is visible in Figs. 1, 2, 3 and 4.

In Fig. 1, we have shown the comparison of the $\nu _{\mu } \rightarrow \nu _e$ oscillation probabilities and $\nu _{\mu } \rightarrow \nu _{\mu }$ survival probabilities between the unitary (labeled u) and non-unitary (labeled n-u) $3\times 3$ mixing, for both the normal hierarchy (NH) and inverted hierarchy (IH) in the case of the NO$\nu$A experiment baseline. The unitary CP violating phase $\delta _{\mathrm {CP}}$ has been kept as a floating parameter. The phases associated with the non-diagonal elements of the $\alpha$ matrix are zero. Similar comparison has been done for the T2K experiment as well in Fig. 2. In the case of unitary mixing, the upper (lower) boundary for $P_{\mu e}$ denotes $\delta _{\mathrm {CP}}=-90^\circ$ ($90^\circ$). For the anti-neutrino oscillation, this characteristic is just the opposite. The probability for all other $\delta _{\mathrm {CP}}$ values fall in between. In both figures, the disappearance plots in lower panel look similar for both the hierarchies. This is because the $\nu _{\mu }$ survival probability does not get affected by matter and hence is not sensitive to hierarchy.

From both Figs. 1 and 2 bottom panels, it is obvious that for $\nu _{\mu }$ disappearance channels the survival probability for unitary $3\times 3$ mixing matrix can be differentiated from the non-unitary mixing at energies away from the oscillation maxima and flux peak energy of NO$\nu$A (2.0 GeV flux peak energy) and T2K (0.7 GeV flux peak energy) experiments. It is also observable that there are overlaps between the unitary and non-unitary probabilities for both $\nu _e$ and $\bar{\nu }_e$ appearances, and therefore, the appearance channels reduce sensitivity to differentiate between the unitary and non-unitary mixing in these experiments.

In Figs. 3 and 4, we have compared between the oscillation probabilities of the unitary and non-unitary case for $\phi _{10}=\pi /2$ and $-\pi /2$, respectively, where $\phi _{10}$ is the phase associated with $\alpha _{10}=|\alpha _{10}|e^{i\phi _{10}}$. For both NH and IH, the oscillation probabilities in the neutrino channel, can be mimicked by the non-unitary oscillation probability with $\phi _{10}=\pi /2$ and it holds for both NO$\nu$A and T2K. For this value of $\phi _{10}$, the discrimination between the unitary and non-unitary cases is better in the anti-neutrino channel. However, for $\phi _{10}=-\pi /2$ and for both the experiments, the unitary oscillation probabilities in the anti-neutrino channel can be mimicked by the non-unitary probability. The neutrino channel has better discrimination capability between the unitary and non-unitary cases for both the experiments in this case. Moreover, for T2K, the unitary oscillation probabilities in the neutrino (anti-neutrino) channel, with NH and almost for the total range of $\delta _{\mathrm {CP}}$, can be mimicked by the non-unitary probabilities with IH and $\phi _{10}=\pi /2$ ($-\pi /2$). It is clear from this discussion that data from the NO$\nu$A and T2K experiments may not constrain the non-unitarity significantly.

Fig. 3

Comparison between the unitary and non-unitary $\nu _{\mu } \rightarrow \nu _e$ probabilities for $\phi _{10}=\pi /2$ ($-\pi /2$) in the upper (lower) panel for the NO$\nu$A experiment. Left (right) panel is for Normal hierarchy (inverted hierarchy). The unitary CP violating phase $\delta _{\mathrm {CP}}$ has been varied in its total range $[-180^\circ :180^\circ ]$. The non-unitary parameters are fixed at their boundary values in Eq. (18) taken from Ref. [46]

Fig. 4

Same as Fig. 3 but for the T2K experiment

In Figs. 5, 6, 7 and 8, we have shown the ratio of the non-unitary to unitary $\nu _{\mu } \rightarrow \nu _e$ oscillation probabilities as a function of L/E and $\delta _{\mathrm {CP}}$, where L is the baseline of an experiment and E is the energy of the neutrino beam. To do this, we fixed E to four values, namely 0.7, 2, 2.5 and 5 GeV and for each E we varied L from 100 km to no more than 1500 km to obtain different L/E values. The non-unitary parameters are fixed at their boundary values as before and $\phi _{10}=0$. The first three energy values are respectively for the NO$\nu$A, T2K and DUNE flux peak points. The farther away the ratio is from 1, the better the potential to differentiate between the two mixing scenarios. In general the sensitivity to non-unitarity is higher at the two opposing corners of the plots, e.g., for $\delta _{\mathrm {CP}}< 0$ with smaller values of L/E or $\delta _{\mathrm {CP}}> 0$ with larger values of L/E for the neutrino channel in NH. The opposite is true for the antineutrino channels and so on. Hence, it is clear from the figures that NO$\nu$A and T2K, along with DUNE, have discrimination capability between the non-unitary and unitary cases for a very small range of $\delta _{\mathrm {CP}}$. If a future experiment can be built with smaller baseline and larger flux peak energy, it will be possible to differentiate the non-unitary mixing from the unitary one better.

Fig. 5

Ratio of the non-unitary to unitary oscillation probabilities as a function of L/E and $\delta _{\mathrm {CP}}$. The reference energy $E=0.7$ GeV has been fixed to the T2K flux peak energy. For this peak energy and T2K baseline, $L/E = 421$. The upper (lower) panel shows the ratio for NH (IH), the left (right) panel shows it for neutrino (anti-neutrino). The non-unitary parameters are fixed to their boundary values in Eq. (18) taken from Ref. [46] and we have set $\phi _{10}=0$

Fig. 6

Same as Fig. 5 but the reference energy $E=2.0$ GeV has been fixed to the NO$\nu$A flux peak energy. For this peak energy and NO$\nu$A baseline, $L/E = 405$

Fig. 7

Same as Fig. 5 but the reference energy $E=2.5$ GeV has been fixed to the DUNE flux peak energy. For this peak energy and DUNE baseline, $L/E = 520$

Fig. 8

Same as Fig. 5 but the reference energy $E=5.0$ GeV has been fixed to the flux peak energy of an arbitrary future experiment

3 Simulation details

The T2K experiment [8] uses the $\nu _{\mu }$ beam from the J-PARC accelerator at Tokai and the water Cerenkov detector at Super-Kamiokande, which is 295 km away from the source. The detector is situated $2.5^\circ$ off-axis. The flux peaks at 0.7 GeV, which is also close to the first oscillation maximum. T2K started taking data in 2009 and until 2019 release of results, these [10, 11, 13] correspond to $14.9 \times 10^{20}$ $(16.4\times 10^{20})$ protons on target (POT) in neutrino (anti-neutrino) mode. The NO$\nu$A detector [9] is a 14 kt totally active scintillator detector (TASD), placed 810 km away from the neutrino source at the Fermilab and it is situated at $0.8^\circ$ off-axis of the NuMI beam. The flux peaks at 2 GeV, close to the oscillation maxima at 1.4 GeV for NH and at 1.8 GeV for IH. NO$\nu$A started taking data in 2014 and took data until 2019 release [12], these correspond to $8.85 \times 10^{20}$ $(6.9 \times 10^{20})$ POTs, for neutrino (anti-neutrino) mode.

To analyse the T2K and NO$\nu$A data, we have taken the solar neutrino parameters ${\Delta _{21}}$ and $\sin ^2 {\theta _{12}}$ to be fixed at $7.50 \times 10^{-5}\, \mathrm{eV}^2$ and 0.30, respectively. For the reactor neutrino angle, $\sin ^2 2{\theta _{13}}$ has been varied in its $3~\sigma$ range around the central value of 0.084 with $3.5 \%$ uncertainty [60]. For the atmospheric mixing angle, $\sin ^2 {\theta _{23}}$ has been varied in the $3\, \sigma$ range [0.40 : 0.63] [57, 61]. The atmospheric effective mass squared difference $\Delta m^2_{\mathrm{eff}}$ has been varied in the MINOS $3~\sigma$ range around the best-fit value of $2.32 \times 10^{-3}\, \mathrm{eV}^2$ [62]. The effective mass-squared difference $\Delta m^{2}_{\mathrm{eff}}$ is related with ${\Delta _{31}}$ by the following relation [63]:

$$\begin{aligned} \Delta m^{2}_{\mathrm{eff}}= & {} \sin ^2 {\theta _{23}}{\Delta _{31}}+ \cos ^2 {\theta _{12}}\Delta _{32} \nonumber \\&+\cos \delta _{\mathrm {CP}}\sin 2{\theta _{12}}\sin {\theta _{13}}\tan {\theta _{12}}{\Delta _{21}}. \end{aligned}$$

(20)

Among the non-unitary parameters, $\alpha _{00},\, \alpha _{11}$, $|\alpha _{10}|$ and $\phi _{10}$ have been varied, while the other non-unitary parameters have been kept constant at their boundary values. This is because only these non-unitary parameters affect appreciably the $\nu _{\mu } \rightarrow \nu _e$ oscillation probability and the $\nu _{\mu } \rightarrow \nu _{\mu }$ survival probability. The choice of the values of $\alpha _{20}$, $\alpha _{21}$ and $\alpha _{22}$ is justified in Fig. 9. In this plot, we have shown the comparison between probabilities for NO$\nu$A when all non-unitary parameters are fixed at their boundary values (denoted by NOvA1) and when $\alpha _{00}$, $\alpha _{11}$, and $\alpha _{10}$ are at their boundary values, but other non-unitary parameters are fixed at their unitary values (denoted by NOvA2). We can see that the difference between these two probabilities are negligible. The unitary parameters are fixed at their present best-fit values taken from Refs. [57, 61] and the hierarchy is NH. We further calculated the $\chi ^2$ between these two parameter sets for neutrino and anti-neutrino run in NO$\nu$A with the latest POT. The $\chi ^2$ value between these two sets of parameter with the latest POT in NO$\nu$A happens to be 0.2. Therefore, it is safe to say that the choice of $\alpha _{20}$, $\alpha _{21}$ and $\alpha _{22}$ do not have any significant effect on the present accelerator neutrino data and we can fix them at their boundary values without any controversy.

Fig. 9

Comparison of probabilities for NO$\nu$A with NH and different sets of values for $\alpha _{20}$, $\alpha _{21}$ and $\alpha _{22}$. For NOvA1 (NOvA2), $\alpha _{00}=0.93$; $\alpha _{11}=0.95$; $|\alpha _{10}=3.6\times 10^{-2}|$; $\alpha _{20}=1.3\times 10^{-1}$ (0); $\alpha _{21}=2.1\times 10^{-2}$ (0); $\alpha _{22}=0.61$ (1). The standard oscillation parameters are fixed at their NH best-fit values, taken from Ref. [57, 61]. The left (right) panel shows the probabilities for neutrino (anti-neutrino) channel and the upper (lower) panel shows probabilities for appearance (disappearance)

We have used GLoBES [64, 65] to calculate the binned theoretical event rates as a function of the test values of the oscillation parameters. The energy dependent efficiencies of the detector for both signal and background events have been fixed according to the expected event rates plots as a function of energy, given in [10,11,12,13]. For $\nu _e$ appearance data, we considered backgrounds from $\nu _{\mu }$ CC interactions, beam contamination and NC interactions. For $\nu _{\mu }$ disappearance data, the backgrounds were from $\nu _e$ CC interactions, beam contamination and NC interactions. Automatic bin based energy smearing for generated theoretical events has been implemented in the same way as described in the GLoBES manual [64, 65]. For this purpose, we used a Gaussian smearing function

$$\begin{aligned} R^c (E,E^{\prime })=\frac{1}{\sqrt{2\pi }}e^{-\frac{(E-E^{\prime })^2}{2\sigma ^2(E)}}, \end{aligned}$$

(21)

where $E^{\prime }$ is the reconstructed energy. The energy resolution function is given by

$$\begin{aligned} \sigma (E)=\alpha E+\beta \sqrt{E}+\gamma , \end{aligned}$$

(22)

where $\alpha =0$, $\beta =0.075$, $\gamma =0.05$ for T2K. For NO$\nu$A, however, we used $\alpha =0$, $\beta =0.085$ (0.06), and $\gamma =0$ for $\nu _e$ ($\nu _{\mu }$) events. For NC events, in NO$\nu$A, we used migration matrices as discussed in [66]. The similar energy smearing techniques have been used in Refs. [17, 67, 68].

Fig. 10

Allowed regions of the $\alpha$ parameters in a triangle plot for NO$\nu$A (T2K) in the upper (lower) panel. The left (right) panel is for NH (IH). The parameter values at the best fit point are listed in Table 1 (2) for NO$\nu$A (T2K). These results are for 2019 data

The experimental event rates have been taken from NO$\nu$A [12] and T2K [10, 11, 13] collaboration papers. The $\chi ^2$ between the theory and experiments have been calculated for the appearance and disappearance channels for both the neutrino and anti-neutrino runs of both the experiments. To generate $\chi ^2$, we have used 30 million data points in the parameter ranges stated above. We have used the Poissonian $\chi ^2$ formula:

$$\begin{aligned} \chi ^2= & {} 2\sum _i \left\{ (1+z) N_i^{\mathrm{th}} - N_i^{\mathrm{exp}} + N_i^{\mathrm{exp}} \ln \left[ \frac{N_i^{\mathrm{exp}}}{(1+z) N_i^{\mathrm{th}}} \right] \right\} \nonumber \\&+ 2 \sum _j (1+z) N_j^{\mathrm{th}} + z^2 \end{aligned}$$

(23)

where i stands for the bins for which $N_i^{\mathrm{exp}}\ne 0$ and j stands for the bins for which $N_j^{\mathrm{exp}} = 0$. The parameter z defines the additional systematic uncertainties. For each of the two experiments, we have included systematic uncertainties of $10\%$, using the pull method. We have varied the pull parameters in their $3\sigma$ range and have marginalized over it to calculate the $\chi ^{2}_{m}$ as a function of the test values of the oscillation parameters and mass hierarchies. For a particular experiment, the total $\chi ^2$ is calculated by

$$\begin{aligned} \chi ^2(\mathrm{tot})= & {} \chi ^{2}_{m}(\nu \, \mathrm{app})+ \chi ^{2}_{m}({{\bar{\nu }}}\, \mathrm{app})+ \chi ^{2}_{m}(\nu \, \mathrm{disapp})\nonumber \\&+ \chi ^{2}_{m}({{\bar{\nu }}}\, \mathrm{disapp})+\chi ^2(\mathrm{prior}). \end{aligned}$$

(24)

During the calculation of $\chi ^2(\mathrm{tot})$, we have to keep in mind that the test values of the oscillation parameters are same for all the individual $\chi _{m}^{2}$s. The $\chi ^2(\mathrm{tot})$ is a function of the test values of the oscillation parameters and hierarchies. The definitions of the $\chi ^2(\mathrm{prior})$ and its significance have been discussed in details in Ref. [69]. In our analysis, we have used priors to $\sin ^2 2{\theta _{13}}$, $\sin ^2 {\theta _{23}}$, and $|\Delta m^{2}_{\mathrm{eff}}|$. Then, we found out the minimum $\chi ^{2}(\mathrm{tot})$ and subtracted it from the $\chi ^{2}(\mathrm{tot})$ values to calculate the $\Delta \chi ^2$ as a function of the oscillation parameters.

To do a combined analysis of the NO$\nu$A and T2K data, we define the total $\chi ^2$ as:

$$\begin{aligned} \chi ^2(\mathrm{tot})&=\chi ^{2}_{m}(\mathrm{NO}\nu \mathrm{A}\, \nu \, \mathrm{app})+ \chi ^{2}_{m}(\mathrm{NO}\nu \mathrm{A}\, {{\bar{\nu }}}\, \mathrm{app}) \nonumber \\&\quad +\chi ^{2}_{m}(\mathrm{T2K}\, \nu \, \mathrm{app})+ \chi ^{2}_{m}(\mathrm{T2K}\, {{\bar{\nu }}}\, \mathrm{app}) \nonumber \\&\quad + \chi ^{2}_{m}(\mathrm{NO}\nu \mathrm{A}\, \nu \, \mathrm{disapp}) + \chi ^{2}_{m}(\mathrm{NO}\nu \mathrm{A}\, {{\bar{\nu }}}\, \mathrm{disapp}) \nonumber \\&\quad +\chi ^{2}_{m}(\mathrm{T2K}\, \nu \, \mathrm{disapp})+ \chi ^{2}_{m}(\mathrm{T2K}\, {{\bar{\nu }}}\, \mathrm{disapp})+\chi ^2(\mathrm{prior}). \end{aligned}$$

(25)

Just like the separate analysis, priors have been added for the $\sin ^2 2{\theta _{13}}$, $\sin ^2 {\theta _{23}}$ and $|\Delta m^{2}_{\mathrm{eff}}|$. The $\Delta \chi ^2$ has been calculated as before.

In this work, we did not simulate the near detector in detail. The effect of the near detector is included in both T2K and NO$\nu$A as errors in the systematic uncertainties. This procedure is the common approach taken in the literature. This approximation is valid in some specific regimes in the presence of extra heavy neutrinos, see Ref. [56] for the discussion. We also remark that in order for the bounds on non-unitarity in any regime from near detectors measurements to be competitive, it is necessary to have a very good knowledge of the flux arriving at the near detectors (see a more detailed discussion on the impact of the flux uncertainties in non-unitary and near detectors in [70]). Because the uncertainties in the near detector flux for both T2K and NO$\nu$A are much larger than the present bounds on those parameters, our sensitivity comes entirely from the non-unitary effects of the propagation of the three (active) neutrinos. This implies that our bounds are valid for heavy neutrinos whose oscillations are averaged out or not produced because their masses are heavier than the experimental energy. Therefore, the systematic-like near detector effect approximation that we have used gives conservative bounds.

4 Analysis of 2019 data

In this section, we discuss the analysis of the 2019 T2K and NO$\nu$A data individually with the hypothesis of non-unitary $3\times 3$ mixing matrix. We have also done a combined analysis of the T2K and NO$\nu$A 2019 data with the non-unitary hypothesis. To do so, we have varied the parameters $\alpha _{00}$ and $\alpha _{11}$ from 0.7 to 1. The parameter $|\alpha _{10}|$ has been varied from 0 to 0.2. The phase $\phi _{10}$, associated with $\alpha _{10}$, has been varied in its total range of $[-180^\circ :180^\circ ]$. These are the only non-unitary parameters that matters in the $\nu _{\mu } \rightarrow \nu _e$ oscillation or $\nu _{\mu } \rightarrow \nu _{\mu }$ survival probabilities [46]. Therefore, all other $\alpha$ parameters have been kept constant at their boundary values. Moreover, only those values of the parameters were chosen, for which $|\alpha _{10}|\le \sqrt{(1-\alpha _{00}^{2})(1-\alpha _{11}^{2})}$ bound is obeyed [46, 71].

4.1 Individual analyses

For the non-unitary case, the minimum $\chi ^2$ for the NO$\nu$A and T2K data are 44.32 and 121.37 for 50 (46 d.o.f.) and 104 (100 d.o.f.) bins respectively and they both occur at the NH. Since we have done similar analyses for both the experiments separately with the unitary $3\times 3$ mixing hypothesis and got the minimum $\chi ^2$ for the NO$\nu$A and T2K data as 47.92 (42 d.o.f.) and 123.71 (96 d.o.f.), respectively, both for NH. We have noted down the values of the unitary and non-unitary parameters at the best fit point for NO$\nu$A and T2K in Tables 1 and 2 respectively. It can be seen that the T2K data cannot rule out either of the two mixing schemes at $1\, \sigma$ CL. For the IH, the minimum $\Delta \chi ^2$ for the NO$\nu$A (T2K) data in the case of non-unitary mixing is 1.36 (2.49), and that for the unitary mixing is 2.7 (6.46). In Fig. 10, we have plotted the allowed region of the $\alpha$ parameters in triangular plots, showing the $1\,\sigma$ and $3\,\sigma$ contours. The grey triangle, labeled $\Delta \chi _{\mathrm{min}}$, corresponds to the best-fit point in the non-unitary case. The cyan triangle, labeled STD, corresponds to the best-fit point in the standard $3\times 3$ unitary case.

Fig. 11

Allowed regions in the $\sin ^2{\theta _{23}}-\delta _{\mathrm {CP}}$ plane for NO$\nu$A (T2K) in the upper (lower) panel. The left (right) panel is for NH (IH). The red (blue) lines indicate $1\, \sigma$ ($3\, \sigma$) CL. The minimum $\chi ^2$ for NO$\nu$A (T2K) in the case of unitary and non-unitary mixing are 47.92 and 44.38 (123.71 and 121.37), respectively. The parameter values at the best-fit point have been mentioned in Table 1 (2) for NO$\nu$A (T2K). These results are for 2019 data

In Fig. 11, we have shown the analysis of individual NO$\nu$A (upper panels) and T2K (lower panels) data analysis in the $\delta _{\mathrm {CP}}-\sin ^2{\theta _{23}}$ plane for both unitary and non-unitary hypothesis. For the standard unitary case of NO$\nu$A, we have got the exact same best-fit value as published by the collaboration [12]. The allowed region is also qualitatively same as the collaboration, though we got a smaller allowed region for the CP conserving $\delta _{\mathrm {CP}}$ values. For T2K, our best fit point is close to the collaboration best fit point [13] and we have got exact same allowed region as had been found by the collaboration.

It is clear from the plots in Fig. 11 that the agreement between the two experiments in the $\sin ^2{\theta _{23}}-\delta _{\mathrm {CP}}$ plane is better, though not perfect, in the case of non-unitarity. In fact, the NO$\nu$A data can put better constraints on the $\delta _{\mathrm {CP}}$ values in the case of non-unitarity. While the unitary hypothesis allows the whole $\delta _{\mathrm {CP}}$ range, the non-unitarity can exclude almost all of the UHP in the $\delta _{\mathrm {CP}}$ range. The best-fit value of the $\sin ^2 {\theta _{23}}$ increases a bit in the case of non-unitarity. Unlike the unitary mixing, the T2K best fit point is allowed by the NO$\nu$A data, in the case of non-unitary mixing. However, the unitary hypothesis can rule out the allowed region for the IH at $1\, \sigma$ CL, whereas in the case of non-unitary, a small region of the IH with $0.46<\sin ^2 {\theta _{23}}<0.52$ and $-140^\circ<\delta _{\mathrm {CP}}<-30^\circ$ is allowed with a minimum $\Delta \chi ^2$ of 1.4.

In the case of T2K (Fig. 11, bottom panels), the non-unitarity does not have any significant effect on the best-fit point or the allowed region, as the best-fit points are similar to the best-fit points of the unitary case. The non-unitarity just makes slightly larger significance region in the $\sin ^2{\theta _{23}}-\delta _{\mathrm {CP}}$ plane to be allowed at $3\, \sigma$ CL for IH. Because of that, the T2K data continue to exclude (include) the NO$\nu$A best-fit point for the NH (IH) at $1\, \sigma$ ($3\, \sigma$) CL. Unlike NO$\nu$A, both the unitary and non-unitary mixing can exclude the IH allowed region at $1\, \sigma$ CL for the T2K data. Therefore, T2K has a better hierarchy sensitivity than NO$\nu$A in the case of non-unitary hypothesis.

In Fig. 12, we have shown the allowed region in the $\delta _{\mathrm {CP}}-\phi _{10}$ plane for both the NO$\nu$A and T2K data. It can be seen that for the NH, both NO$\nu$A and T2K cannot exclude any value of $\phi _{10}$ at $1\, \sigma$ CL. For the IH, at $1\, \sigma$ CL, data from NO$\nu$A allow a very small region of $\delta _{\mathrm {CP}}$ in the LHP and $\phi _{10}$ in the UHP.

Fig. 12

Allowed regions in the $\delta _{\mathrm {CP}}-\phi _{10}$ plane for NO$\nu$A (T2K) in the upper (lower) panel. The left (right) panel is for NH (IH). The minimum $\chi ^2$ for NO$\nu$A (T2K) in the case of unitary and non-unitary mixing are 47.92 and 44.38 (123.71 and 121.37), respectively. The parameter values at the best-fit point have been mentioned in Table 1 (2) for NO$\nu$A (T2K). These results are for 2019 data

4.2 Combined analysis

Figure 13 shows the allowed region in the $\delta _{\mathrm {CP}}-\sin ^2 {\theta _{23}}$ and $\delta _{\mathrm {CP}}-\phi _{10}$ planes for the combined analysis of the T2K and NO$\nu$A data. The minimum $\chi ^2$ for the combined analysis with non-unitary mixing is 170.90 for 146 d.o.f. Same analysis has been done for the unitary $3\times 3$ mixing, and the minimum $\chi ^2$ for the unitary case has been found out to be 173.40 for 150 d.o.f. The combined analysis, just like NO$\nu$A, prefers non-unitary mixing at $1\, \sigma$ CL. The parameter values at the best-fit point for this combined analysis has been given in Table 3. The minimum $\Delta \chi ^2$ for IH is 8.16 (4.89) for (non-) unitary case. Therefore in Fig. 13 (top right panel), there is no allowed region at $1\, \sigma$ CL for IH. In Fig. 14, we have shown the allowed region of $\alpha$ parameters in a triangular plot, similar to Fig. 10.

Fig. 13

Allowed regions in the $\sin ^2{\theta _{23}}-\delta _{\mathrm {CP}}$ ($\delta _{\mathrm {CP}}-\phi _{10}$) plane for the combined analysis in the upper (lower) panel. The left (right) panel is for NH (IH). The red (blue) lines indicate $1\, \sigma$ ($3\, \sigma$) CL. The minimum $\chi ^2$ for the unitary (non-unitary) case is 173.40 (170.90). The parameter values at the best-fit point have been mentioned in Table 3. These results are for 2019 data

Fig. 14

Allowed regions of the $\alpha$ parameters in a triangular plot for the combined analysis. The left (right) panel is for NH (IH). The parameter values at the best-fit point have been mentioned in Table 3. These results are for 2019 data

We have shown $\Delta \chi ^2$ as a function of $\delta _{\mathrm {CP}}$, in the case of non-unitarity, for the individual T2K and NO$\nu$A, and the combined analysis of the two experiments in Fig. 15. It is obvious that the individual T2K analysis can exclude $60\%$ of the $\delta _{\mathrm {CP}}$ plane at $2\, \sigma$ for the NH. It can also exclude the IH for $90\%$ of the $\delta _{\mathrm {CP}}$ plane at $2\, \sigma$ CL. But, NO$\nu$A can exclude the IH only for $50\%$ of the $\delta _{\mathrm {CP}}$ plane at $2\, \sigma$, and for the NH, it cannot disfavor any value of $\delta _{\mathrm {CP}}$ at $2\, \sigma$ CL. The combined analysis can exclude the IH at $2\, \sigma$ CL. for every value of $\delta _{\mathrm {CP}}$, and for the NH, it can exclude the UHP of $\delta _{\mathrm {CP}}$ at $2\, \sigma$ CL.

Fig. 15

$\Delta \chi ^2$ as a function of $\delta _{\mathrm {CP}}$ for the individual T2K and NO$\nu$A, and the combined analysis. These results are for 2019 data. The limiting values of $\Delta \chi ^2$ for $1\, \sigma$ and $2\, \sigma$ CL are 1 and 4 respectively

Finally, in Fig. 16, we have plotted $\Delta \chi ^2$ as a function of the individual non-unitary $\alpha$ parameters, so that one can have an idea of the bounds put on these parameters by the individual T2K and NO$\nu$A analyses, and their combined analysis. It is clear from Fig. 16 that analyses of T2K and the combined data from both the experiments prefer rather tiny deviation from unitarity at the best-fit point. However, large deviation from unitarity is allowed at $1\, \sigma$ CL. Both NO$\nu$A and the combined analysis also rules out the unitary values of $\alpha _{00}$, $|\alpha _{10}|$ and $\alpha _{11}$ at $1\, \sigma$ CL.

Fig. 16

$\Delta \chi ^2$ as a function of the $\alpha$ parameters for the individual T2K and NO$\nu$A, and the combined analysis. The x-axis of $|\alpha _{10}|$ plot is in log scale. These results are for 2019 data

We list the best-fit parameter values with $1\, \sigma$ CL intervals for NO$\nu$A in Table 1. The best fit parameter values are reported with $1\, \sigma$ CL intervals, and the $90\%$ and $3\, \sigma$ CL limit for $\alpha$ parameters from the analyses of T2K and the combined data have been mentioned in Tables 2 and 3, respectively.

5 Analysis of the 2020 data

In June, 2020, NO$\nu$A [18] and T2K [19] have published their new data in the Neutrino 2020 conference. So far, NO$\nu$A data have been analysed for $1.36 \times 10^{21}$ ($1.25 \times 10^{21}$) POT in $\nu$ ($\bar{\nu }$) mode. T2K data have been analysed for $1.97 \times 10^{21}$ ($1.63 \times 10^{21}$) POT in $\nu$ ($\bar{\nu }$) mode. According to the present data, the tension between the two experiments are even stronger. The best-fit point for NO$\nu$A (T2K) is $\sin ^2{\theta _{23}}= 0.57$ (0.528) and $\delta _{\mathrm {CP}}=0.82 \pi$ ($-1.6\pi$). Moreover, there is no overlap between the $1\, \sigma$ allowed regions of the two experiments. In such a scenario, it is even more important to test new physics hypotheses with the new data from T2K and NO$\nu$A experiments.

As before, in the GLoBES software we have tuned the signal and background efficiencies according to the Monte-Carlo simulations given by the collaborations. This time, GLoBES in its latest update has included data analysis facility and we have used GLoBES completely to analyse the data.

To analyse the 2020 data, in Eq. (22), for NO$\nu$A, we have used $\alpha =0.11\, (0.09)$, $\beta =\gamma =0$ for electron (muon) like events. For T2K, we have used $\alpha =0$, $\beta =0.075$, $\gamma =0.05$ for both electron and muon like events.

For both NO$\nu$A and T2K, we have used

• $5\%$ normalisation and $5\%$ energy calibration systematic uncertainties for the e-like events, and

• $5\%$ normalisation and $0.01\%$ energy calibration systematic uncertainties for the $\mu$-like events.

Implementing systematic uncertainties has been discussed in details in GLoBES manual [64, 65]. Here also we have used more than 30 million test data points in the parameter ranges discussed in Sect. 3.

Table 1 Parameter values at the best-fit points for NO$\nu$A.The $1\, \sigma$ error bars have been mentioned where possible. The $90\%$ limits for 1 d.o.f have also been mentioned. In each box, the result with 2019 (2020) data has been mentioned at the top (bottom) of the box

Table 2 Parameter values at the best-fit points for T2K.The $1\, \sigma$ error bars have been mentioned where possible. The $90\%$ limits for 1 d.o.f have also been mentioned. In each box, the result with 2019 (2020) data has been mentioned at the top (bottom) of the box

Table 3 Parameter values at the best-fit points for the combined analysis of NO$\nu$A and T2K.The $1\, \sigma$ error bars have been mentioned where possible. The $90\%$ limits for 1 d.o.f have also been mentioned. In each box, the result with 2019 (2020) data has been mentioned at the top (bottom) of the box

At first we have analysed the data with unitary mixing hypothesis and the results have been shown in Fig. 17. The minimum $\chi ^2$ for NO$\nu$A (T2K) with 46 (84) d.o.f. is 48.65 (95.85) and it occurs at NH. For the combined analysis, minimum $\chi ^2$ for 138 bins (134 d.o.f.) is 147.14 and it occurs at IH. From Fig. 17, it can be seen that the tension between the two experiments continue. Each experiment excludes the other’s allowed region at $1\, \sigma$ CL. The combined analysis prefers IH over NH. It allows NH for a very tiny CP conserving region at $1\, \sigma$ CL.

At the next step, we have analysed the data with non-unitary hypothesis. The minimum $\chi ^2$ for NO$\nu$A (T2K) is 45.88 (93.36) 42 (80) d.o.f. and it is at NH. The minimum $\chi ^2$ for the combined analysis, however is at IH and its value is 142.72 for 130 d.o.f. The IH best fit points for individual analysis and the NH best fit point for the combined analysis are mentioned in Fig. 18.

Fig. 17

Allowed region in the $\sin ^2 {\theta _{23}}-\delta _{\mathrm {CP}}$ plane after analysing NO$\nu$A and T2K complete 2020 data sets with unitary mixing hypothesis. The left (right) panel represents the test hierarchy as NH (IH). The red (blue) lines indicate the results for NO$\nu$A (T2K) and the black lines indicate the combined analysis of both. The solid (dashed) lines indicate the $1\, \sigma$ ($3\, \sigma$) allowed regions. The minimum $\chi ^2$ for NO$\nu$A (T2K) with 46 (84) d.o.f. is 48.65 (95.85) and it occurs at NH. For the combined analysis, the minimum $\chi ^2$ with 134 d.o.f. is 147.14 and it prefers IH

Fig. 18

Allowed region in the $\sin ^2 {\theta _{23}}-\delta _{\mathrm {CP}}$ plane after analysing NO$\nu$A and T2K complete 2020 data set with non-unitary mixing hypothesis. The left (right) panel represents test hierarchy NH (IH). The red (blue) lines indicate the results for NO$\nu$A (T2K) and the black line indicates the combined analysis of both. The solid (dashed) lines indicate the $1\, \sigma$ ($3\, \sigma$) allowed regions. The minimum $\chi ^2$ for NO$\nu$A (T2K) with 42 (80) d.o.f. is 45.88 (93.36) and it occurs at NH. For the combined analysis, the minimum $\chi ^2$ with 130 d.o.f. is 142.72 and it prefers IH

In Fig. 19, we have shown $\Delta \chi ^2$ as a function of individual non-unitary parameters. To do this, we have marginalised $\Delta \chi ^2$ over all parameters, except the one against which we have plotted.

It is quite certain that with the new data, both T2K and NO$\nu$A prefer non-unitary mixing over unitary mixing at $1\, \sigma$ CL. The combined analysis excludes unitary mixing at more than $2\, \sigma$ CL. In all three cases, the preference for non-unitarity is stronger compared to the 2019 data set. The tension between the two experiments is also reduced when analysed with non-unitary hypothesis, as there are overlaps between the allowed regions in the $\sin ^2{\theta _{23}}-\delta _{\mathrm {CP}}$ plane at $1\, \sigma$. However, The $\delta _{\mathrm {CP}}$ best-fit points are still far away from each other and there is a new, mild tension between the ${\theta _{23}}$ octant at best-fit points between the two experiments, T2K prefers the lower octant while NO$\nu$A prefers the higher octant. But NO$\nu$A (T2K) cannot rule out lower (higher) octant at $1\,\sigma$ CL. Although, NH is preferred over IH by both the experiments in separate analyses, there is an almost degenerate IH best-fit point with $\Delta \chi ^2=0.66$ (0.14) for NO$\nu$A (T2K). The combined analysis prefers IH over NH at $1\, \sigma$ CL. It should also be noted that the two experiments have strong agreement for IH best-fit point with $\delta _{\mathrm {CP}}$ in the LHP and ${\theta _{23}}$ in the lower octant.

The best-fit points along with the $1\,\sigma$ error bars and the $90\%$ confidence level constraints on the new physics parameters have been mentioned in Tables 1, 2 and 3. It can be seen the constraints on the non-unitary parameters are even weaker for the 2020 data compared to the 2019 data.

6 Conclusions

With the 2019 data, the NO$\nu$A experiment disfavors the unitary mixing at $~1\, \sigma$ CL in favor of the non-unitary mixing. The T2K experiment, however, cannot exclude any of the two hypotheses at $1\, \sigma$ CL. With the non-unitary hypothesis, NO$\nu$A includes the T2K best-fit point at $1\,\sigma$ CL, but T2K still continues to disfavor the NO$\nu$A best-fit point at $1\, \sigma$ CL. Unitary hypothesis can exclude the IH at $1.5\, \sigma$ CL for NO$\nu$A and at $2\, \sigma$ CL for T2K. With the non-unitary hypothesis, hierarchy can be determined only at $1\, \sigma$ CL for the NO$\nu$A data. However, T2K can determine hierarchy at $2\, \sigma$ for $90\%$ of the $\delta _{\mathrm {CP}}$ plane with the non-unitary hypothesis. The combined analysis prefers non-unitary mixing over unitary mixing at $1\, \sigma$ CL.

It should be noted that the tension between NO$\nu$A and T2K 2019 data is reduced when both the experiments are analysed with non-unitary hypothesis. The $90\%$ confidence level limit on the $\alpha$ parameters, from these two long baseline accelerator based neutrino experiments is weaker compared to the constraints given from the global analysis in Ref. [46] for the neutrinos only.

For the latest 2020 data set, both the experiments individually prefer non-unitary mixing over the unitary one at $1\, \sigma$ CL, preferring NH in both cases. The combined analysis, however, prefers IH both for unitary and non-unitary mixing but the non-unitary mixing is favored over unitary mixing at more than $2\, \sigma$ CL. The tension between the two experiments is also reduced when analysed with non-unitary hypothesis. Both the experiments lose ${\theta _{23}}$ octant sensitivity when analysed with non-unitarity.

Fig. 19

$\Delta \chi ^2$ as a function of individual non-unitary parameters for 2020 long baseline data

The constraints on the non-unitary parameters are even weaker with the 2020 data as compared to the 2019 data. As a consequence, the preference for non-unitarity is stronger with the 2020 data.

It can be commented that the present long baseline data prefer non-unitary mixing giving a hint of the possibility of new physics. It is important that the future analysis of NO$\nu$A and T2K data are done with the non-unitary $3\times 3$ mixing hypothesis, besides the standard unitary mixing hypothesis in order to find new physics signatures. If the two experiments continue to have better agreement with the non-unitary hypothesis, that can be a strong hint of the presence of new physics in the neutrino sector.