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

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 isosinglets, the effective $3\times 3$ mixing matrix for the three active neutrinos will be non-unitary. We have analyzed the latest T2K and \nova 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$ C.L.) prefer the non-unitary mixing over unitary mixing. In fact, allowing the non-unitary mixing brings the \nova best-fit point in the $\sin^2\theta_{23}-\delta_{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$ C.L.. 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$ C.L. The combined analysis preferes non-unitarity at $2\, \sigma$ C.L. The stronger tension, which exists between the latest 2020 data of the two experiments, also gets reduced with non-unitary analysis.

As per the latest analysis, the best-fit point for NOνA (T2K) is sin 2 θ 23 = 0.57 (0.528) and δ CP = 0.82π (−1.6π). The tension between the two experiments is even stronger as they exclude each other's allowed region at 1 σ C.L. It has also been shown that although both the experiments individually prefer NH over IH, their combined analysis prefers IH over NH [19].
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ν 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 avg = N obs /N pred = 0.927 ± 0.023 [20]. Recent updates have changed the predictions slightly, giving an average ratio R avg = 0.938 ± 0.023 [21], which is a 2.7 σ 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 ν e disappearance from the SAGE [22] and GALLEX [23] solar neutrino experiments. A combined analysis of the detected and predicted number of neutrino events from the source gives R = 0.86±0.05 [24,25], another 2.7 σ deviation from unity. Both of these deficits of low energyν e events can be explained by an oscillation at ∆m 2 ≥ 1 eV 2 over very short baseline.
2. LSND and MiniBooNE anomalies: The LSND experiment [26] at the Los Alamos National laboratory was designed to observeν µ →ν e oscillations over a baseline of 30 m. After 5 years of data taking, it observed 89.7 ± 22.4 ± 6.0ν e candidate events over background, providing a 3.8 σ evidence ofν µ →ν e oscillations at ∆m 2 = 1 eV 2 region. Therefore, this result cannot be accommodated in the three flavor scenario.
The MiniBooNE experiment [27] at the Fermilab was designed to observe ν µ → ν e and ν µ →ν 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 σ signal excess of ν e candidate and 2.8 σ signal excess ofν 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. [28] for a comprehensive review. The minimal model consists of 3+1 neutrino mixing, dominated by ν e , ν µ and ν τ , Recent results from the IceCube experiment constrain the sterile neutrino mass and mixing using atmospheric neutrino fluxes [29]. In 2018, however, MiniBooNE has again confirmed a 4.7 σ excess of combined ν e andν e events [30]. The present significance of the excess from a combined analysis of MiniBooNE and LSND is 6 σ.
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 × 3 non-unitary mixing matrix [39]. NHL would induce charged lepton flavor violation processes [39,40]. If Majorana type, these NHLs will modify the rate of neutrinoless double beta decay [41,42]. 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. Ref. [43] 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. [44]. In ref. [45], physics potential of long baseline experiments with non-unitary mixing has been studied, while ref. [46] has made a theoretical study of the non-unitary oscillation at the probability level. In ref. [47], the CP violation measurement potential of short baseline (SBL) experiments in the presence of nonunitary ν µ → ν τ 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νA experiments as we have done in this work. However, in ref. [48,49] an effort has been made to resolve the tension between NOνA and T2K with non-standard NC interaction during propagation and in ref. [50] 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. [51].
In this paper, we have explored whether the 2019 and 2020 T2K and NOνA data can exclude the non-unitary 3 × 3 mixing, and if not, whether the non-unitary 3 × 3 mixing hypothesis can lead to better agreement between the T2K and NOνA data. A combined analysis with the non-unitary hypothesis has also been done. In Section 2, we have discussed the theory of neutrino oscillation probability in the presence of a non-unitary 3 × 3 mixing matrix. Details of the simulation method have been discussed in Section 3. In Sections 4 and 5, we have presented our results for the 2019 and 2020 data respectively, and the conclusion has been drawn in Section 6.

II. NON-UNITARY OSCILLATION PROBABILITIES
In the standard case of 3 active neutrinos, the flavor basis ν f (f = e, µ, τ ) is related to the mass basis ν m (m = 1, 2, 3) by the relation ν f = U ν m , where U is the unitary PMNS matrix. The Schrödinger equation for neutrino propagation in matter can be written as where H is the Hamiltonian and A is the matter potential. The solution to the above equation, after propagation over a distance L in matter, is is the S-matrix in the mass basis, with W being the transformation matrix between the mass basis ν m and the mass basis in matterν m such that ν m = Wν m . The energiesẼ = diag(Ẽ 1 ,Ẽ 2 ,Ẽ 3 ), whereẼ i are the eigenvalues in theν 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 Flavor change in terms of the S-matrix is ν f (L) = S f (L)ν f . Correspondingly, the oscillation probability from flavor α to flavor β can be written as One can extend this formalism to n×n unitary mixing matrix U n×n with 3 active neutrinos and (n − 3) heavy singlet neutrinos [51]. In that case, where N is a 3×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) × (n − 3) matrix in the heavy neutrino sector. The interaction potential matrix A can be written as where, for the usual matter potential, the term A 3×3 is the same as in eq. (1) and A (n−3)×(n−3) = 0. An explicit formalism for the potential matrix has been discussed in details in refs. [44]. The Hamiltonian in vacuum in this case can be written as where H 3×3 is the Hamiltonian in vacuum from eq. (1). Similarly, the neutrino flavor vector can be rewritten as where we split the flavor vector in active ν f 3×1 and heavy ν f (n−3)×1 parts. The W matrix similar to that in eq. (2) can be written as and theẼ matrix from eq. (2) can be written as The Hamiltonian in the presence of matter potential, following eq. (1), is Similarly, the S-matrix in the mass basis, following eq. (2), is 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×(n−3) and W (n−3)×3 satisfy conditions that they are 1, following the typical see-saw mechanisms [52,53]. Therefore, we can neglect these terms and write the S-matrix as The S-matrix in the flavor basis is written, following eq. (3), as In case of 3 active neutrinos the flavor changes, following eq. (8), as Here S m 3×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 [51,53]. 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 The corresponding oscillation probability can be calculated using eq. (4).
We use the following parametrization of the non-unitary mixing matrix N [51] 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 α parameters vary at least at the percent scale [39,51,54]. There exists severe constraint on the parameter |α 10 | < 10 −5 from non-observation of the µ → eγ decay [54].
Anti-neutrino oscillation probability Pμē can be obtained by changing the sign of A and δ CP in eq. 19. The oscillation probability mainly depends on hierarchy (sign of ∆ 31 ), octant of θ 23 and δ CP .
From eq. (19), we can see that for the standard unitary case, the oscillation probability P µe (Pμē) gets a double boost (double suppression) when the hierarchy is NH and δ CP is in the lower half plane (LHP). Similarly P µe (Pμē) gets a double suppression (double boost) when the hierarchy is IH and δ CP is in the upper half plane (UHP). Therefore P µe (Pμē) is maximum (minimum) for NH and δ CP in the LHP and minimum (maximum) for IH and δ CP in the UHP.
In the case of non-unitary mixing, an analytic expression for ν µ → ν 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., [51]. NOν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 δ CP < 0 with smaller values of L/E or δ 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νA and T2K, along with DUNE, have discrimination capability between the non-unitary and unitary cases for a very small range of δ 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.

III. SIMULATION DETAILS
The T2K experiment [8] uses the ν µ beam from the J-PARC accelerator at Tokai and the water Cerenkov detector at Super-Kamiokande, which is 295 km away from the source. taken from ref. [55,59] and the hierarchy is NH. We further calculated the χ 2 between these two parameter sets for neutrino and anti-neutrino run in NOνA with the latest POT. The χ 2 value between these two sets of parameter with the latest POT in NOνA happens to be 0.2. Therefore, it is safe to say that the choice of α 20 , α 21 and α 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.
We have used GLoBES [62,63] 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 ν e appearance data, we considered backgrounds from ν µ CC interactions, beam contamination and NC interactions. For ν µ disappearance data, the backgrounds were from ν 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 [62,63]. For this purpose, we used a Gaussian smearing function where E is the reconstructed energy. The energy resolution function is given by where α = 0, β = 0.075, γ = 0.05 for T2K. For NOνA, however, we used α = 0, β = 0.085 (0.06), and γ = 0 for ν e (ν µ ) events. For NC events, in NOνA, we used migration matrices as discussed in [64]. The similar energy smearing techniques have been used in refs. [16,65,66].
The experimental event rates have been taken from NOνA [12] and T2K [10,11,13] collaboration papers. The χ 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 χ 2 , we have used 30 million data points in the parameter ranges stated above. We have used the Poissonian χ 2 formula: where i stands for the bins for which N exp the χ 2 m as a function of the test values of the oscillation parameters and mass hierarchies. For a particular experiment, the total χ 2 is calculated by During the calculation of χ 2 (tot), we have to keep in mind that the test values of the oscillation parameters are same for all the individual χ 2 m s. The χ 2 (tot) is a function of the test values of the oscillation parameters and hierarchies. The definitions of the χ 2 (prior) and its significance have been discussed in details in ref. [67]. In our analysis, we have used priors to sin 2 2θ 13 , sin 2 θ 23 , and |∆m 2 eff |. Then, we found out the minimum χ 2 (tot) and subtracted it from the χ 2 (tot) values to calculate the ∆χ 2 as a function of the oscillation parameters.
To do a combined analysis of the NOνA and T2K data, we define the total χ 2 as: Just like the separate analysis, priors have been added for the sin 2 2θ 13 , sin 2 θ 23 and |∆m 2 eff |. The ∆χ 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ν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. [54] 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 [68]). Because the uncertainties in the near detector flux for both T2K and NOν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.

IV. ANALYSIS OF 2019 DATA
In this section, we discuss the analysis of the 2019 T2K and NOνA data individually with the hypothesis of non-unitary 3 × 3 mixing matrix. We have also done a combined analysis of the T2K and NOνA 2019 data with the non-unitary hypothesis. To do so, we have varied the parameters α 00 and α 11 from 0.7 to 1. The parameter |α 10 | has been varied from 0 to 0.2. The phase φ 10 , associated with α 10 , has been varied in its total range of [−180 • : 180 • ]. These are the only non-unitary parameters that matters in the ν µ → ν e oscillation or ν µ → ν µ survival probabilities [44]. Therefore, all other α parameters have been kept constant at their boundary values. Moreover, only those values of the parameters were chosen, for which |α 10 | ≤ (1 − α 2 00 )(1 − α 2 11 ) bound is obeyed [44,69]. in Tables I and II  In fig. 11, we have shown the analysis of individual NOνA (upper panels) and T2K

A. Individual analyses
(lower panels) data analysis in the δ CP − sin 2 θ 23 plane for both unitary and non-unitary hypothesis. For the standard unitary case of NOν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 δ CP values.    Table I (II) for NOνA (T2K). These results are for 2019 data.
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 θ 23 − δ CP plane is better, though not perfect, in the case of non-unitarity. In fact, the 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 θ 23 − δ CP plane to be allowed at 3 σ C.L. for IH. Because of that, the T2K data continue to exclude (include) the NOνA best-fit point for the NH (IH) at 1 σ (3 σ) C.L. Unlike NOνA, both the unitary and non-unitary mixing can exclude the IH allowed region at 1 σ C.L. for the T2K data. Therefore, T2K has a better hierarchy sensitivity than NOνA in the case of non-unitary hypothesis.
In fig. 12, we have shown the allowed region in the δ CP −φ 10 plane for both the NOνA and T2K data. It can be seen that for the NH, both NOνA and T2K cannot exclude any value of φ 10 at 1 σ C.L. For the IH, at 1 σ C.L., data from NOνA allow a very small region of δ CP in the LHP and φ 10 in the UHP.   Table III. These results are for 2019 data.
We list the best-fit parameter values with 1 σ C.L. intervals for NOνA in Table I. The best fit parameter values are reported with 1 σ C.L. intervals, and the 90% and 3 σ C.L. limit for α parameters from the analyses of T2K and the combined data have been mentioned in tables II and III, respectively.  Table III. These results are for 2019 data.
According to the present data, the tension between the two experiments are even stronger.
Moreover, there is no overlap between the 1 σ 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νA experiments. For both NOν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 µ-like events.
Implementing systematic uncertainties has been discussed in details in GLoBES manual [62,63]. Here also we have used more than 30 million test data points in the parameter ranges discussed in section III.     2 σ C.L. 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 nonunitary hypothesis, as there are overlaps between the allowed regions in the sin 2 θ 23 − δ CP plane at 1 σ. However, The δ CP best-fit points are still far away from each other and there is a new, mild tension between the θ 23 octant at best-fit points between the two experiments,    It should be noted that the tension between NOνA and T2K 2019 data is reduced when both the experiments are analysed with non-unitary hypothesis. The 90% confidence level limit on the α parameters, from these two long baseline accelerator based neutrino experiments is weaker compared to the constraints given from the global analysis in ref. [44] for the neutrinos only.
For the latest 2020 data set, both the experiments individually prefer non-unitary mixing over the unitary one at 1 σ C.L., 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 σ C.L. The tension between the two experiments is also reduced when analysed with non-unitary hypothesis. Both the experiments lose θ 23 octant sensitivity when analysed with non-unitarity. The constraints on the nonunitary 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