# A review of the tension between the T2K and NO$\nu$A appearance data and
hints to new physics

Ushak Rahaman<EMAIL_ADDRESS>Centre for Astro-Particle Physics (CAPP) and
Department of Physics, University of Johannesburg, PO Box 524, Auckland Park
2006, South Africa  Soebur Razzaque<EMAIL_ADDRESS>Centre for Astro-
Particle Physics (CAPP) and Department of Physics, University of Johannesburg,
PO Box 524, Auckland Park 2006, South Africa  S. Uma Sankar
<EMAIL_ADDRESS>Dept. of Physics, IIT Bombay, Powai, Mumbai-400076, India

###### Abstract

In this article, we review the status of the tension between the long-baseline
accelerator neutrino experiments T2K and NO$\nu$A. The tension arises mostly
due to the mismatch in the appearance data of the two experiments. We explain
how this tension arises based on $\nu_{\mu}\to\nu_{e}$ and
$\bar{\nu}_{\mu}\to\bar{\nu}_{e}$ oscillation probabilities. We define the
reference point of vacuum oscillation, maximal ${\theta_{23}}$ and
$\delta_{\mathrm{CP}}=0$ and compute the $\nu_{e}/\bar{\nu}_{e}$ appearance
events for each experiment. We then study the effects of deviating the unknown
parameters from the reference point and the compatibility of any given set of
values of unknown parameters with the data from T2K and NO$\nu$A. T2K observes
a large excess in the $\nu_{e}$ appearance event sample compared to the
expected $\nu_{e}$ events at the reference point, whereas NO$\nu$A observes a
moderate excess. The large excess in T2K dictates that $\delta_{\mathrm{CP}}$
be anchored at $-90^{\circ}$ and that ${\theta_{23}}>\pi/4$ with a preference
for normal hierarchy. The moderate excess at NO$\nu$A leads to two degenerate
solutions: A) NH, $0<\delta_{\mathrm{CP}}<180^{\circ}$, and
${\theta_{23}}>\pi/4$; B) IH, $-180^{\circ}<\delta_{\mathrm{CP}}<0$, and
${\theta_{23}}>\pi/4$. This is the main cause of the tension between the two
experiments. We have reviewed the status of three beyond standard model (BSM)
physics scenarios, (a) non-unitary mixing, (b) Lorentz invariance violation
and (c) non-standard neutrino interactions, to resolve the tension.

Neutrino Mass Hierarchy, Long-Baseline Experiments

###### pacs:

14.60.Pq,14.60.Lm,13.15.+g

##  I Introduction

Neutrino oscillations have provided the first signal for physics beyond the
standard model (SM). They were first proposed to explain the deficit in the
solar neutrino flux observed by the pioneering Homestake experiment Cleveland
et al. (1998). Oscillations between two neutrino flavours require them to mix
and form two mass eigenstates. The survival probability of a neutrino with
energy $E$ and given flavour $\alpha$, after propagation over a distance $L$
in vacuum, is given by

$P(\nu_{\alpha}\to\nu_{\alpha})=1-\sin^{2}2\theta\sin^{2}\left(1.27\frac{\Delta
m^{2}L}{E}\right),$ (1)

where $\Delta m^{2}$ is the difference between the squares of the neutrino
masses and $\theta$ is the mixing angle. In eq. (1), the units are chosen such
that $\Delta m^{2}$ should be specified in eV2, $L$ in meters and $E$ in MeV.
The solar neutrino deficit was confirmed by the water Cerenkov detector
Kamiokande Hirata et al. (1991), which detected the solar neutrinos in real
time. Radio-chemical Gallium experiments, GALLEX Hampel et al. (1999), SAGE
Abdurashitov et al. (2009) and GNO Altmann et al. (2005), which were mostly
sensitive to the low energy $pp$ solar neutrinos, also observed a deficit. The
high statistics water Cerenkov detector Super-Kamiokande Abe et al. (2016a)
and the heavy water Cerenkov detector SNO Aharmim et al. (2013) and Borexino
Agostini et al. (2019) also have made detailed spectral measurements of the
solar neutrino fluxes. Analysing the solar neutrino data in a two flavour
oscillation framework gives the oscillation parameters

$\Delta m^{2}_{\rm sol}\sim 10^{-4}~{}{\rm eV}^{2}~{}~{}~{}~{}{\rm
and}~{}~{}~{}~{}\sin^{2}\theta_{\rm sol}\sim 0.33.$ (2)

Observation of proton decay is one of the main physics motivations for the
construction of the water Cerenkov detectors, IMB Casper et al. (1991);
Becker-Szendy et al. (1992) and Kamiokande Hirata et al. (1992); Fukuda et al.
(1994). The interactions of atmospheric neutrinos in the detector, especially
those of $\nu_{e}$ and $\bar{\nu}_{e}$, could mimic the proton decay signal.
Hence, these experiments made a detailed study of the atmospheric neutrino
interactions. They did not find any signal for proton decay but instead
observed a deficit of up-going atmospheric $(\nu_{\mu}+\bar{\nu}_{\mu})$ flux,
relative to the down-going flux. It was proposed that the up-going neutrinos,
which travel thousands of km inside the Earth, oscillate into another flavour
whereas the down-going neutrinos, which travel tens of km, do not. Super-
Kamiokande experiment Fukuda et al. (1998) observed a zenith angle dependence
of the deficit, which is expected from neutrino oscillations. An analysis of
the atmospheric neutrino data in a two-flavour oscillation framework gives the
oscillation parameters

$|\Delta m^{2}_{\rm atm}|\sim 3\times 10^{-3}~{}{\rm eV}^{2}~{}~{}~{}~{}{\rm
and}~{}~{}~{}~{}\sin^{2}\theta_{\rm atm}\sim 0.5.$ (3)

We note that $\Delta m^{2}_{\rm sol}\ll\Delta m^{2}_{\rm atm}$.

It is known that there are three flavours of neutrinos, $\nu_{e}$, $\nu_{\mu}$
and $\nu_{\tau}$ Steinberger (1991). They mix to form three mass eigenstates,
$\nu_{1}$, $\nu_{2}$ and $\nu_{3}$, with mass eigenvalues $m_{1}$, $m_{2}$ and
$m_{3}$. The $3\times 3$ unitary matrix $U$, connecting the flavour basis to
the mass basis,

$\left[\begin{array}[]{c}\nu_{e}\\\ \nu_{\mu}\\\
\nu_{\tau}\end{array}\right]=\left[\begin{array}[]{ccc}U_{e1}&U_{e2}&U_{e3}\\\
U_{\mu 1}&U_{\mu 2}&U_{\mu 3}\\\ U_{\tau 1}&U_{\tau 2}&U_{\tau
3}\end{array}\right]\left[\begin{array}[]{c}\nu_{1}\\\ \nu_{2}\\\
\nu_{3}\end{array}\right],$ (4)

is called the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix Maki et al.
(1962); Bilenky and Pontecorvo (1978). Naively, it seems desirable to label
the lightest mass $m_{1}$, the middle mass $m_{2}$ and the heaviest mass
$m_{3}$. However, no method exists at present to directly measure these
masses. What can be measured in oscillation experiments are the mixing matrix
elements $U_{\alpha j}$, where $\alpha$ is a flavour index and $j$ is a mass
index. In particular, the three elements of the first row $U_{e1}$, $U_{e2}$
and $U_{e3}$ are well measured. The labels $1$, $2$ and $3$ are chosen such
that $|U_{e3}|<|U_{e2}|<|U_{e1}|$.

Given the three masses $m_{1}$, $m_{2}$ and $m_{3}$, it is possible to define
two independent mass-squared differences, $\Delta_{31}=m_{3}^{2}-m_{1}^{2}$
and $\Delta_{21}=m_{2}^{2}-m_{1}^{2}$. The third mass-squared difference is
then $\Delta_{32}=m_{3}^{2}-m_{2}^{2}=\Delta_{31}-\Delta_{21}$. Without loss
of generality, we can choose $\Delta_{21}=\Delta m^{2}_{\rm sol}$ and
$\Delta_{31}=\Delta m^{2}_{\rm atm}$. Since $\Delta m^{2}_{\rm sol}\ll\Delta
m^{2}_{\rm atm}$, we find that $\Delta_{32}\approx\Delta m^{2}_{\rm atm}$.
Because of the way the labels $1$, $2$ and $3$ are chosen, these mass-squared
differences, in principle, can be either positive or negative. Their signs
have to be determined by experiments. The solar neutrinos, produced at the
core of the sun, undergo forward elastic scattering as they travel through the
solar matter. This scattering leads to matter effect Wolfenstein (1978, 1979),
which modifies the solar electron neutrino survival probability $P_{\rm ee}$.
Super-Kamiokande Abe et al. (2018a) and SNO Aharmim et al. (2013) have
measured $P_{\rm ee}$ as a function of neutrino energy for $E>5$ MeV and found
it to be of a constant value $\simeq 0.3$. SNO has also measured Aharmim et
al. (2013) the neutral current interaction rate of solar neutrinos to be
consistent with predictions of the standard solar model Bahcall and
Pinsonneault (2004). The measurements of the Gallium experiments imply that
$P_{\rm ee}>0.5$ for neutrino energies $E<0.5$ MeV. This increase in $P_{\rm
ee}$ at lower solar neutrino energies can be explained only if $\Delta
m^{2}_{\rm sol}=\Delta_{21}$ is positive. At present, there is no definite
experimental evidence for either positive or negative sign of $\Delta
m^{2}_{\rm atm}=\Delta_{31}\approx\Delta_{32}$.

The PMNS matrix is similar to the quark mixing matrix introduced by Kobayashi
and Maskawa Kobayashi and Maskawa (1973). It can be parameterized in terms of
three mixing angles, $\theta_{12}$, $\theta_{13}$ and $\theta_{23}$, and one
CP-violating phase $\delta_{\mathrm{CP}}$. The following parameterization of
the PMNS matrix is found to be the most convenient to analyze the neutrino
data

$U=\left(\begin{array}[]{ccc}1&0&0\\\ 0&c_{23}&s_{23}\\\
0&-s_{23}&c_{23}\end{array}\right)\left(\begin{array}[]{ccc}c_{13}&0&s_{13}e^{-i\delta_{\mathrm{CP}}}\\\
0&1&0\\\
-s_{13}e^{i\delta_{\mathrm{CP}}}&0&c_{13}\end{array}\right)\left(\begin{array}[]{ccc}c_{12}&s_{12}&0\\\
-s_{12}&c_{12}&0\\\ 0&0&1\end{array}\right),$ (5)

where $c_{ij}=\cos\theta_{ij}$ and $s_{ij}=\sin\theta_{ij}$. Among these
mixing angles, CHOOZ experiment set a strong upper bound on the middle angle
$\theta_{13}$ Apollonio et al. (1998)

$\sin^{2}2\theta_{13}\leq 0.1.$ (6)

By combining this limit with the solar and atmospheric neutrino data, it can
be shown that $\theta_{13}\ll 1$ Narayan et al. (1998).

The mixing angles $\theta_{ij}$ and the mass-squared differences $\Delta
m^{2}_{ij}$ have been determined in a series of precision experiments with
man-made neutrino sources, which we briefly describe below.

  * •

The long baseline reactor neutrino experiment KamLAND Araki et al. (2005) has
$L\simeq 180$ km. At this long distance, it can observe the oscillations due
to the small mass-squared difference $\Delta_{21}$. In the limit of neglecting
$\theta_{13}$ in the three-flavour oscillation, the expression for the anti-
neutrino survival probability $P(\bar{\nu}_{e}\to\bar{\nu}_{e})$ reduces to an
effective two flavour expression

$P(\bar{\nu}_{e}\to\bar{\nu}_{e})=1-\sin^{2}2\theta_{12}\sin^{2}\left(1.27\frac{\Delta_{21}L}{E}\right).$
(7)

KamLAND measured the spectral distortion $P(\bar{\nu}_{e}\to\bar{\nu}_{e})$
precisely. A combined analysis of KamLAND and solar neutrino data yields the
results

$\Delta_{21}=(7.9\pm 0.5)\times 10^{-5}~{}{\rm eV}^{2}~{}~{}~{}~{}{\rm
and}~{}~{}~{}~{}\tan^{2}\theta_{12}=0.4^{+0.10}_{-0.07}$ (8)

  * •

The short baseline reactor neutrino experiments, Daya Bay Adey et al. (2018),
RENO Bak et al. (2018) and DoubleCHOOZ Abe et al. (2016b), have baselines of
the order of 1 km. At this distance, the oscillating term in
$P(\bar{\nu}_{e}\to\bar{\nu}_{e})$ containing $\Delta_{21}$ is negligibly
small. In this approximation, $P(\bar{\nu}_{e}\to\bar{\nu}_{e})$ again reduces
to an effective two-flavour expression

$P(\bar{\nu}_{e}\to\bar{\nu}_{e})=1-\sin^{2}2\theta_{13}\sin^{2}\left(1.27\frac{\Delta_{31}L}{E}\right).$
(9)

High statistics measurement from Daya Bay gives the measurement

$|\Delta_{31}|=(2.47\pm 0.07)\times 10^{-3}~{}{\rm eV}^{2}~{}~{}~{}~{}{\rm
and}~{}~{}~{}~{}\sin^{2}2\theta_{13}=(0.0856\pm 0.0029).$ (10)

  * •

The long baseline accelerator experiment MINOS Adamson et al. (2014) has a
baseline of 730 km and it measured the survival probability of the accelerator
$\nu_{\mu}$ beam. For this baseline and for accelerator neutrino energies, the
oscillating term in $P(\nu_{\mu}\to\nu_{\mu})$ due to $\Delta_{21}$ is
negligibly small. Setting this term and $\theta_{13}$ to be zero in the three
flavour expression for $P(\nu_{\mu}\to\nu_{\mu})$, once again we obtain an
effective two flavour expression

$P(\nu_{\mu}\to\nu_{\mu})=1-\sin^{2}2\theta_{23}\sin^{2}\left(1.27\frac{\Delta_{31}L}{E}\right).$
(11)

MINOS data gives the results

$|\Delta_{31}|=(2.3-2.5)\times 10^{-3}~{}{\rm eV}^{2}~{}~{}~{}~{}{\rm
and}~{}~{}~{}~{}\sin^{2}\theta_{23}=(0.35-0.65).$ (12)

Note that both short baseline reactor neutrino experiments and long baseline
accelerator experiments as well as atmospheric neutrino data determine only
the magnitude of $\Delta_{31}$ but not its sign. Hence, we must consider both
positive and negative sign possibilities in the data analysis. The case of
positive $\Delta_{31}$ is called normal hierarchy (NH) and that of negative
$\Delta_{31}$ is called inverted hierarchy (IH). Both atmospheric neutrino
data and accelerator neutrino data are functions of $\sin^{2}2\theta_{23}$ and
they prefer $\sin^{2}2\theta_{23}\simeq 1$. For such values, there are two
possibilities: $\theta_{23}<45^{\circ}$ called lower octant (LO) and
$\theta_{23}>45^{\circ}$ called higher octant (HO). At present, the data is
not able to make a distinction between these two cases.

A number of groups have done global analysis of neutrino oscillation data from
all the available sources: solar, atmospheric, reactor and accelerator Esteban
et al. (2020); de Salas et al. (2021). In Table 1, we present the latest
results obtained by the nu-fit collaboration Esteban et al. (2020).

| Normal Hierarchy | Inverted Hierarchy
---|---|---
$\sin^{2}\theta_{12}$ | $0.304^{+0.012}_{-0.012}$ | $0.304^{+0.013}_{-0.012}$
$\sin^{2}\theta_{13}$ | $0.02246\pm 0.00062$ | $0.02241^{+0.00074}_{-0.00062}$
$\sin^{2}\theta_{23}$ | $0.45^{+0.019}_{-0.016}$ | $0.570^{+0.016}_{-0.022}$
$\frac{\Delta_{21}}{10^{-5}~{}{\rm eV^{2}}}$ | $7.42\pm 0.2$ | $7.42\pm 0.2$
$\frac{\Delta_{31}}{10^{-3}~{}{\rm eV^{2}}}$ | $2.51\pm 0.027$ | $-2.49^{+0.026}_{-0.028}$
Table 1: Neutrino oscillation parameters determine by nu-fit collaboration
using global neutrino oscillation data Esteban et al. (2020).

Among the neutrino oscillation parameters, the mass-squared differences and
the mixing angles (except for $\theta_{23}$) are measured to a precision of a
few percent. On the other hand, the CP-violating phase $\delta_{\mathrm{CP}}$
still eludes measurement. In addition, we also need to resolve the issues of
the sign of $\Delta_{31}$ (also called the problem of neutrino mass hierarchy)
and the octant of $\theta_{23}$. Thus, there are currently three main unknowns
in the three-flavour neutrino oscillation paradigm.

It can be shown that the survival probability of neutrinos,
$P(\nu_{\alpha}\to\nu_{\alpha})$ is equal to that of the anti-neutrinos
$P(\bar{\nu}_{\alpha}\to\bar{\nu}_{\alpha})$ due to CPT invariance. But the
oscillation probabilities $P(\nu_{\alpha}\to\nu_{\beta})$ and
$P(\bar{\nu}_{\alpha}\to\bar{\nu}_{\beta})$ are not equal if there is CP
violation. A measurement of the difference between these two probabilities
will establish CP-violation in neutrino oscillations and will also determine
$\delta_{\mathrm{CP}}$. While considering oscillation probabilities, in
principle, we can enumerate six possibilities, with $\alpha=e/\mu/\tau$ and
$\beta$ taking two possible values other than $\alpha$. It can be shown that,
in all six cases, the difference

$\Delta P^{\alpha\beta}_{\rm
CP}=P(\nu_{\alpha}\to\nu_{\beta})-P(\bar{\nu}_{\alpha}\to\bar{\nu}_{\beta})$
(13)

is proportional to

$\frac{\Delta_{21}}{\Delta_{31}}\cos\theta_{13}\sin 2\theta_{12}\sin
2\theta_{13}\sin
2\theta_{23}\sin\delta_{\mathrm{CP}}=\frac{\Delta_{21}}{\Delta_{31}}J,$ (14)

where $J$ is the Jarlskog invariant of the PMNS matrix Jarlskog (1985).

Practically speaking, $\alpha=e/\tau$ at production are not possible because
there are no intense sources of $\nu_{e}/\nu_{\tau}/\bar{\nu}_{\tau}$. Nuclear
reactors do produce $\bar{\nu}_{e}$ copiously but their energies are in the
range of a few MeV. When such $\bar{\nu}_{e}$ oscillate into
$\bar{\nu}_{\mu}$, the resultant anti-neutrinos are not energetic enough to
produce $\mu^{+}$ by their interactions in the detector. Hence,
$\alpha=e/\tau$ are not practical choices. The neutrino beams produced by
accelerators yield intense fluxes of $\nu_{\mu}/\bar{\nu}_{\mu}$. In
principle, it is possible to search for CP-violation by contrasting
$P(\nu_{\mu}\to\nu_{e})$ with $P(\bar{\nu}_{\mu}\to\bar{\nu}_{e})$ or by
contrasting $P(\nu_{\mu}\to\nu_{\tau})$ with
$P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})$. However, the second option is much
more difficult compared to the first for the reasons listed below.

  * •

To produce $\tau^{\mp}$ in the detector, due to the interactions of
$\nu_{\tau}/\bar{\nu}_{\tau}$, the neutrinos should have energies of tens of
GeV. At these energies, the oscillation probabilities are quite small.

  * •

Even when we have energetic-enough beams to produce $\tau^{\mp}$ in the
detector, the efficiency of reconstructing these particles is very poor. Thus
the event numbers will be very limited.

  * •

From eq. (14), we see that the CP violating asymmetry

$A_{CP}^{\mu\tau}=\frac{P(\nu_{\mu}\to\nu_{\tau})-P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})}{P(\nu_{\mu}\to\nu_{\tau})+P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})}$

will be very small because the numerator is a product of two small quantities
$(\Delta_{21}/\Delta_{31})$ and $\sin 2\theta_{13}$, whereas the denominator
is close to $1$. Hence a measurement of this CP-asymmetry requires very high
statistics.

Thus the most feasible method to establish CP-violation in neutrino
oscillations and to determine $\delta_{\mathrm{CP}}$ is to measure the
difference between $P(\nu_{\mu}\to\nu_{e})$ and
$P(\bar{\nu}_{\mu}\to\bar{\nu}_{e})$. The neutrinos do not require large
energies to produce electrons/positrons on interacting in the detector. So the
neutrino beam energy can be tuned to the oscillation maximum. The produced
electrons and positrons are relatively easy to identify in the detector. The
dominant term in the expression for $P(\nu_{\mu}\to\nu_{e})$ is proportional
to $\sin^{2}\theta_{23}\sin^{2}2\theta_{13}$ Cervera et al. (2000). The
expression for the CP-asymmetry in $\nu_{\mu}\to\nu_{e}$ oscillations has the
form

$A_{\rm CP}^{\mu
e}\sim\frac{\Delta_{21}}{\Delta_{31}}\frac{J}{\sin^{2}2\theta_{13}}\sim\frac{\Delta_{21}}{\Delta_{31}}\frac{1}{\sin
2\theta_{13}},$ (15)

which is much larger than $A_{\rm CP}^{\mu\tau}$. Thus, CP-violation in these
oscillations can be established with moderate statistics.

There are, however, some other difficulties to overcome before the goal of
establishing CP-violation in neutrino oscillations can be achieved. The matter
effect, which modifies the solar neutrino oscillation probabilities, modifies
$P(\nu_{\mu}\to\nu_{e})$ and $P(\bar{\nu}_{\mu}\to\bar{\nu}_{e})$ Lipari
(2000); Narayan and Sankar (2000). These modifications depend on the sign of
$\Delta_{31}$. Since the dominant terms in these oscillation probabilities are
proportional to $\sin^{2}\theta_{23}$, they are also subject to the octant
ambiguity of $\theta_{23}$. That is, the two oscillation probabilities,
$P(\nu_{\mu}\to\nu_{e})$ and $P(\bar{\nu}_{\mu}\to\bar{\nu}_{e})$, depend on
all the three unknowns of the three flavour neutrino oscillation parameters.
In such a situation, the change in the probabilities induced by changing one
of the unknowns can be compensated by changing another unknown. This leads to
degenerate solutions which can explain a given set of measurements.
Unravelling these degeneracies and making a distinction between the degenerate
solutions requires a number of careful measurements and moderately high
statistics.

In this review article, we have discussed the theory of oscillation
probability and parameter degeneracy in section II. Details of the $\chi^{2}$
analysis have been discussed in section III. In section IV, we have discussed
the chronology of NO$\nu$A and T2K data. In the same section, we have also
tried to explain the results of the analysis of data from NO$\nu$A and T2K in
the past and the present ones as well with the help of parameter degeneracy
and explain the cause for the tension between the data of the two experiments.
The resolution of the tension in terms of BSM physics has been discussed in
section V. A summary of the article has been drawn in section VI.

##  II Oscillation probability and parameter degeneracy

Two accelerator experiments, T2K Itow et al. (2001) and NO$\nu$A Ayres et al.
(2007), are taking data with the aim of establishing CP-violation as well as
determining neutrino mass hierarchy and the octant of $\theta_{23}$. Both
experiments share the following common features.

  * •

They aim a beam of $\nu_{\mu}/\bar{\nu}_{\mu}$ to a far detector a few hundred
km away, which is at an off-axis location.

  * •

The off-axis location leads to a sharp peak in neutrino spectrum McDonald,
Kirk T (2001), which is crucial to suppress the $\pi^{0}$ events, produced in
via the neutral current reaction $\nu_{\mu}N\to\nu_{\mu}\pi^{0}N$, that form a
large background for the $\nu_{\mu}\to\nu_{e}$ oscillation signal.

  * •

They have a near detector, a few hundred meters from the accelerator, which
measures the neutrino flux accurately.

  * •

The energy of the neutrino beam is tuned to be close to the oscillation
maximum.

  * •

They measure the two survival probabilities, $P(\nu_{\mu}\to\nu_{\mu}$ and
$P(\bar{\nu}_{\mu}\to\bar{\nu}_{\mu})$, and the oscillation probabilities
$P(\nu_{\mu}\to\nu_{e}$ and $P(\bar{\nu}_{\mu}\to\bar{\nu}_{e})$.

The survival probabilities lead to further improvement in the precision of
$|\Delta_{31}|$ and $\sin 2\theta_{23}$. A careful analysis of the oscillation
probabilities can lead to the determination of three unknowns of the neutrino
oscillation parameters. The crucial parameters of T2K and NO$\nu$A experiments
are summarised in table 2. Note that the integrated flux of accelerator
neutrinos is specified in units of protons on target (POT).

| T2K | NO$\nu$A
---|---|---
Baseline | 295 km | 810 km
Energy of Peak Flux | 0.7 GeV | 2.0 GeV
Detector Type | Water Cerenkov | Liquid Scintillator
Detector Mass | 22.5 ktons | 14 ktons
Total Flux in $\nu$ mode |  $13.6\times 10^{20}$ POT |  $16.3\times 10^{20}$ POT
Total Flux in $\bar{\nu}$ mode |  $19.7\times 10^{20}$ POT |  $12.5\times 10^{20}$ POT
Table 2: Summary of important information of T2K and NO$\nu$A experiments.

We first start with a discussion of the oscillation probabilities,
$P(\nu_{\mu}\to\nu_{e})$ and $P(\bar{\nu}_{\mu}\to\bar{\nu}_{e})$ and describe
how they vary with each of the unknown neutrino oscillation parameters. The
three flavour $\nu_{\mu}\to\nu_{e}$ oscillation probability in the presence of
matter effect with constant matter density can be written as Cervera et al.
(2000)

$\displaystyle P_{\mu e}$ $\displaystyle\simeq$
$\displaystyle\sin^{2}2{\theta_{13}}\sin^{2}{\theta_{23}}\frac{\sin^{2}\hat{\Delta}(1-\hat{A})}{(1-\hat{A})^{2}}$
(16) $\displaystyle+$ $\displaystyle\alpha\cos{\theta_{13}}\sin
2{\theta_{12}}\sin 2{\theta_{13}}\sin
2{\theta_{23}}\cos(\hat{\Delta}+\delta_{\mathrm{CP}})\frac{\sin\hat{\Delta}\hat{A}}{\hat{A}}\frac{\sin\hat{\Delta}(1-\hat{A})}{1-\hat{A}}$
$\displaystyle+$
$\displaystyle\alpha^{2}\sin^{2}2\theta_{12}\cos^{2}\theta_{13}\cos^{2}\theta_{23}\frac{\sin^{2}\hat{\Delta}\hat{A}}{\hat{A}^{2}},$

where $\alpha=\frac{{\Delta_{21}}}{{\Delta_{31}}}$,
$\hat{\Delta}=\frac{{\Delta_{31}}L}{4E}$ and
$\hat{A}=\frac{A}{{\Delta_{31}}}$, with $E$ being the energy of the neutrino
and $L$ being the length of the baseline. The parameter $A$ is the Wolfenstein
matter term Wolfenstein (1978), given by $A=2\sqrt{2}G_{F}N_{e}E$, where
$G_{F}$ is the Fermi coupling constant and $N_{e}$ is the number density of
the electrons in the matter. Anti-neutrino oscillation probability
$P_{\bar{\mu}\bar{e}}$ can be obtained by changing the sign of $A$ and
$\delta_{\mathrm{CP}}$ in eq. (16). The oscillation probability mainly depends
on hierarchy (sign of ${\Delta_{31}}$), octant of ${\theta_{23}}$ and
$\delta_{\mathrm{CP}}$, and precision in the value of ${\theta_{13}}$. $P_{\mu
e}$ is enhanced if $\delta_{\mathrm{CP}}$ is in the lower half plane (LHP,
$-180^{\circ}<\delta_{\mathrm{CP}}<0$), and it is suppressed if
$\delta_{\mathrm{CP}}$ is in the upper half plane (UHP,
$0<\delta_{\mathrm{CP}}<180^{\circ}$), compared to the CP conserving
$\delta_{\mathrm{CP}}$ values. In the following paragraph, for the sake of
discussion, we will treat $\delta_{\mathrm{CP}}$ as a binary variable that
either increases or decreases oscillation probability.

The dominant term in $P_{\mu e}$ is proportional to $\sin^{2}2{\theta_{13}}$.
Therefore, the oscillation probability is rather small. It can be enhanced
(suppressed), by $8\%$ for T2K and $22\%$ for NO$\nu$A, due to the matter
effect if ${\Delta_{31}}$ is positive (negative). This dominant term is also
proportional to $\sin^{2}{\theta_{23}}$. If $\sin^{2}2{\theta_{23}}<1$, there
can be two possible cases: i) $\sin^{2}{\theta_{23}}<0.5$ which will suppress
$P_{\mu e}$, and ii) $\sin^{2}{\theta_{23}}>0.5$ which will enhance $P_{\mu
e}$ relative to the maximal ${\theta_{23}}$. Since each of the unknowns can
take 2 possible values, there are 8 different combinations of three unknowns.
Any given value of $P_{\mu e}$ can be reproduced by any of these eight
combinations of the three unknowns with the appropriate choice of
${\theta_{13}}$ value. Thus, if the value of ${\theta_{13}}$ is not known
precisely, it will lead to eight fold degeneracy in $P_{\mu e}$. Given that
$\theta_{13}$ has been measured quite precisely, this degeneracy becomes less
severe.

####

a. Hierarchy-$\delta_{\mathrm{CP}}$ degeneracy: To start with, we assume that
${\theta_{23}}$ is maximal and the values of ${\theta_{13}}$ and
${\theta_{12}}$ are precisely known. With these assumptions, the only two
unknowns are hierarchy and $\delta_{\mathrm{CP}}$. From table 1, we see that,
according to the current measurements, $\sin 2{\theta_{13}}\approx 0.3\pm
0.005$ whereas $|\alpha|\approx 0.03\pm 0.001$. Therefore, the first term in
the expression of $P_{\mu e}$ (and in $P_{\bar{\mu}\bar{e}}$) in eq. (16) has
the maximum matter effect contribution. This term is much larger than second
term and the third term is extremely small. We will neglect the third term in
all further discussions.

For NH (IH), the first term in $P_{\mu e}$ becomes larger (smaller). For
$P_{\bar{\mu}\bar{e}}$, the situation is reversed. These changes in $P_{\mu
e}$ and in $P_{\bar{\mu}\bar{e}}$ can be amplified or canceled by the second
term, depending on the value of $\delta_{\mathrm{CP}}$. Because of the
dependence on $\hat{A}$ term, $P_{\mu e}$ ($P_{\bar{\mu}\bar{e}}$) for NH is
always larger (smaller) than that for IH. At the oscillation maxima,
$\hat{\Delta}\simeq 90^{\circ}$. Thus for $P_{\mu e}$
($P_{\bar{\mu}\bar{e}}$), the term $\cos(\hat{\Delta}+\delta_{\mathrm{CP}})$
is maximum (minimum) for $\delta_{\mathrm{CP}}=-90^{\circ}$ and it is minimum
(maximum) for $\delta_{\mathrm{CP}}=90^{\circ}$. Therefore, for NH and
$\delta_{\mathrm{CP}}=-90^{\circ}$, $P_{\mu e}$ ($P_{\bar{\mu}\bar{e}}$) is
maximum (minimum) and for IH and $\delta_{\mathrm{CP}}=90^{\circ}$, it is
minimum (maximum). These two hierarchy-$\delta_{\mathrm{CP}}$ combinations,
for both $P_{\mu e}$ and $P_{\bar{\mu}\bar{e}}$, are well separated from each
other. It can be shown that oscillation probability, for the NH and
$\delta_{\mathrm{CP}}$ in the LHP, is well separated from that for the IH and
$\delta_{\mathrm{CP}}$ in the UHP, for both neutrino and anti-neutrino. But
for the other two hierarchy-$\delta_{\mathrm{CP}}$ combinations, NH and
$\delta_{\mathrm{CP}}$ in the UHP, and IH and $\delta_{\mathrm{CP}}$ in the
LHP, $P_{\mu e}$ and $P_{\bar{\mu}\bar{e}}$ are quite close to each other,
leading to hierarchy-$\delta_{\mathrm{CP}}$ degeneracy. This is illustrated in
Fig. 1, where $P_{\mu e}$ and $P_{\bar{\mu}\bar{e}}$ are plotted for the
NO$\nu$A experiment baseline. For these plots, we have used maximal
${\theta_{23}}$, i.e., $\sin^{2}{\theta_{23}}=0.5$. The other mixing angle
values are $\sin^{2}{\theta_{12}}=0.32$ and $\sin^{2}2{\theta_{13}}=0.089$.
For the mass-squared differences, we have used ${\Delta_{21}}=7.50\times
10^{-5}\ {\rm eV}^{2}$ and $|\Delta_{\mu\mu}|=2.40\times 10^{-3}\ {\rm
eV}^{2}$. $\Delta m^{2}_{\mu\mu}$ is related with ${\Delta_{31}}$ by the
following equation Nunokawa et al. (2005):

$\Delta_{\mu\mu}=\sin^{2}{\theta_{23}}{\Delta_{31}}+\cos^{2}{\theta_{12}}\Delta_{32}+\cos\delta_{\mathrm{CP}}\sin
2{\theta_{12}}\sin{\theta_{13}}\tan{\theta_{12}}{\Delta_{21}}.$ (17)

$\Delta_{\mu\mu}$ is positive (negative) for NH (IH). For the NH and
$\delta_{\mathrm{CP}}$ in the LHP, the values of $P_{\mu e}$
($P_{\bar{\mu}\bar{e}}$) are reasonably greater (lower) than the values of
$P_{\mu e}$ ($P_{\bar{\mu}\bar{e}}$) for the IH and any value of
$\delta_{\mathrm{CP}}$. Similarly, for the IH and $\delta_{\mathrm{CP}}$ in
the UHP the values of $P_{\mu e}$ ($P_{\bar{\mu}\bar{e}}$) are reasonably
lower (greater) than the values of $P_{\mu e}$ ($P_{\bar{\mu}\bar{e}}$) for
the NH and any value of $\delta_{\mathrm{CP}}$. Hence, for these favourable
combinations, NO$\nu$A is capable of determining the hierarchy at a confidence
level (C.L.) of $2\sigma$ or better, with 3 years run each for $\nu$ and
${\bar{\nu}}$. However, as mentioned above, the change in the first term in
eq. (16) can be canceled by the second term for unfavourable values of
$\delta_{\mathrm{CP}}$. This leads to hierarchy-$\delta_{\mathrm{CP}}$
degeneracy Barger et al. (2002); Minakata and Sugiyama (2004); Mena and Parke
(2004). From Fig. (1), it can be seen that $P_{\mu e}$ and
$P_{\bar{\mu}\bar{e}}$ for NH and $\delta_{\mathrm{CP}}$ in the UHP are very
close to or degenerate with those of IH and $\delta_{\mathrm{CP}}$ in the LHP.
For these unfavourable combinations, NO$\nu$A has no hierarchy sensitivity
Prakash et al. (2012). Since the neutrino energy of T2K is only one third of
the energy of NO$\nu$A, the matter effect of T2K is correspondingly smaller.
Therefore, T2K has very little hierarchy sensitivity. The cancellation of
change due to matter effect occurs for different values of
$\delta_{\mathrm{CP}}$ in the case of NO$\nu$A and T2K. Therefore, combining
the data of NO$\nu$A and T2K leads to a small hierarchy discrimination
capability for unfavourable hierarchy-$\delta_{\mathrm{CP}}$ combinations Mena
and Parke (2004); Prakash et al. (2012); Agarwalla et al. (2012).

Figure 1: $P_{\mu e}$ (left panel) and $P_{\bar{\mu}\bar{e}}$ (right panel)
vs. neutrino energy for the NO$\nu$A baseline. Variation of
$\delta_{\mathrm{CP}}$ leads to the blue (red) bands for NH (IH). The plots
are drawn for maximal ${\theta_{23}}$ and other neutrino parameters given in
the main text. Figure 2: Illustration of degenerate $P_{\mu e}$ and non-
degenerate $P_{\bar{\mu}\bar{e}}$ for the following two cases. Left: (LO-NH,
$\delta_{\mathrm{CP}}=-45^{\circ}$) and (HO-IH,
$\delta_{\mathrm{CP}}=-45^{\circ}$), Right: (LO-NH,
$\delta_{\mathrm{CP}}=-90^{\circ}$) and (HO-IH,
$\delta_{\mathrm{CP}}=-45^{\circ}$).

####

b. Octant-hierarchy degeneracy: Even though the atmospheric neutrino
experiments prefer maximal ${\theta_{23}}$ ($\sin^{2}2{\theta_{23}}=1$), MINOS
experiment prefers non-maximal values, $\sin^{2}2{\theta_{23}}=0.96$ Nichol
(2012a). The global fits, before the NO$\nu$A and T2K experiment started
taking data, also favour a non-maximal value of ${\theta_{23}}$ Forero et al.
(2012); Fogli et al. (2012); Gonzalez-Garcia et al. (2012), leading to two
degenerate solutions: ${\theta_{23}}$ in the lower octant (LO)
($\sin^{2}{\theta_{23}}=0.41$) and ${\theta_{23}}$ in the higher octant (HO)
($\sin^{2}{\theta_{23}}=0.59$). Given the two hierarchy and two octant
possibilities, there are four possible octant-hierarchy combinations: LO-NH,
HO-NH, LO-IH and HO-IH. The first term of $P_{\mu e}$ in eq. (16) becomes
larger (smaller) for NH (IH). The same term also becomes smaller (larger) for
LO (HO). If HO-NH (LO-IH) is the true octant-hierarchy combination, then the
values of $P_{\mu e}$ are significantly higher (smaller) than those for IH
(NH) and for any octant. For these two cases, only $\nu$ data has good
hierarchy determination capability. But the situation is very different for
the two cases LO-NH and HO-IH. The increase (decrease) in the first term of
$P_{\mu e}$ due to NH (IH) is canceled (compensated) by the decrease
(increase) for LO (HO). Therefore the two octant-hierarchy combinations, LO-NH
and HO-IH, have degenerate values for $P_{\mu e}$. However, this degeneracy is
not present in $P_{\bar{\mu}\bar{e}}$, which receives a double boost
(suppression) for the case of HO-IH (LO-NH). Thus the octant-hierarchy
degeneracy in $P_{\mu e}$ is broken by $P_{\bar{\mu}\bar{e}}$ (and vice-
verse). Therefore $\nu$-only data has no hierarchy sensitivity if the cases
LO-NH or HO-IH are true, but a combination of $\nu$ and ${\bar{\nu}}$ data
will have a good sensitivity.

This has been illustrated in Figure 2. From the figure, we can see that
$P_{\mu e}$ has degeneracy for the octant-hierarchy combinations LO-NH and HO-
IH. This degeneracy does not exist in the case of $P_{\bar{\mu}\bar{e}}$
Prakash et al. (2014).

####

c. Octant-$\delta_{\mathrm{CP}}$ degeneracy: The possibility of two octants of
${\theta_{23}}$ also leads to octant-$\delta_{\mathrm{CP}}$ degeneracy. To
highlight this degeneracy, we rewrite the expression for $P_{\mu e}$ in eq.
(16) as Agarwalla et al. (2013)

$P_{\mu
e}=\beta_{1}\sin^{2}{\theta_{23}}+\beta_{2}\cos(\hat{\Delta}+\delta_{\mathrm{CP}})+\beta_{3}\cos^{2}{\theta_{23}},$
(18)

where

$\displaystyle\beta_{1}=\sin^{2}2{\theta_{13}}_{13}\frac{\sin^{2}\hat{\Delta}(1-\hat{A})}{(1-\hat{A})^{2}},$
$\displaystyle\beta_{2}=\alpha\cos{\theta_{13}}\sin 2{\theta_{12}}\sin
2{\theta_{13}}\sin 2{\theta_{23}}\frac{\sin(\hat{\Delta}\hat{A})}{\hat{A}}$
$\displaystyle\beta_{3}=\alpha^{2}\sin^{2}2{\theta_{12}}\cos^{2}{\theta_{13}}\frac{\sin^{2}(\hat{\Delta}\hat{A})}{\hat{A}^{2}}.$
(19)

In Figure 3 we have plotted $P_{\mu e}$ ($P_{\bar{\mu}\bar{e}}$) for the
NO$\nu$A experiment as a function of neutrino energy $E_{\nu}$, for normal
hierarchy and for different values of $\delta_{\mathrm{CP}}$. In our
calculation, $\sin^{2}{\theta_{23}}=0.41$, when ${\theta_{23}}$ is in the LO
and $\sin^{2}{\theta_{23}}=0.59$, when ${\theta_{23}}$ is in the HO. The
$\sin^{2}2{\theta_{13}}$ has been taken equal to 0.089.

Figure 3: $P_{\mu e}$ as a function of neutrino energy for the NO$\nu$A
baseline. The left (right) panel is for $\nu$ (${\bar{\nu}}$). The plots have
been drawn for different possible $\delta_{\mathrm{CP}}$ values between
$-180^{\circ}$ and $180^{\circ}$. The $\sin^{2}2{\theta_{13}}$ value is 0.089.
The value of $\sin^{2}{\theta_{23}}$ is 0.41 (0.59) for LO (HO).

From the left panel of the figure, we can see that when ${\theta_{23}}$ is in
the LO and $\delta_{\mathrm{CP}}$ is $90^{\circ}$, $P_{\mu e}$ is quite
distinctive from probability values with other ${\theta_{23}}$ and
$\delta_{\mathrm{CP}}$ combinations. Similar arguments hold for $P_{\mu e}$
with ${\theta_{23}}$ in the HO and $\delta_{\mathrm{CP}}=-90^{\circ}$.
Therefore, ${\theta_{23}}$ in the LO (HO) and
$\delta_{\mathrm{CP}}=90^{\circ}$ ($-90^{\circ}$) is a favourable
octant-$\delta_{\mathrm{CP}}$ combination to determine the octant of
${\theta_{23}}$. However $P_{\mu e}$ for ${\theta_{23}}$ in LO and
$\delta_{\mathrm{CP}}=-90^{\circ}$ overlaps with that for ${\theta_{23}}$ in
HO and $\delta_{\mathrm{CP}}=90^{\circ}$. Therefore, ${\theta_{23}}$ in the LO
(HO) and $\delta_{\mathrm{CP}}=-90^{\circ}$ ($90^{\circ}$) is an unfavourable
combination to determine the octant.

But these unfavourable combinations, become favourable for
$P_{\bar{\mu}\bar{e}}$ and vice-versa, as can be seen from the right panel of
Figure 3. Thus the octant-$\delta_{\mathrm{CP}}$ degeneracy, present in
neutrino data, can be removed by anti-neutrino data and vice-versa. This
aspect is different from the hierarchy-$\delta_{\mathrm{CP}}$ degeneracy.

The above features can also be understood by following algebraic analysis of
eq. (18). From that equation, we see that $P_{\mu e}$ increases with the
increase in ${\theta_{23}}$. But $P_{\mu e}$ can increase or decrease with
change in $\delta_{\mathrm{CP}}$. In the case of octant-$\delta_{\mathrm{CP}}$
degeneracy, for different $\delta_{\mathrm{CP}}^{\rm LO}$ and
$\delta_{\mathrm{CP}}^{\rm HO}$, we can have $P_{\mu e}(\rm
LO,\delta_{\mathrm{CP}}^{\rm LO})=P_{\mu e}(\rm HO,\delta_{\mathrm{CP}}^{\rm
HO})$. It leads to

$\cos(\hat{\Delta}+\delta_{\mathrm{CP}}^{\rm
LO})-\cos(\hat{\Delta}+\delta_{\mathrm{CP}}^{\rm
HO})=\frac{\beta_{1}-\beta_{3}}{\beta_{2}}(\sin^{2}{\theta_{23}}^{\rm
HO}-\sin^{2}{\theta_{23}}^{\rm LO}).$ (20)

The NO$\nu$A experiment has a baseline of 810 km and the flux peaks at an
energy 2 GeV. Now for the NH and neutrino,

$\cos(\hat{\Delta}+\delta_{\mathrm{CP}}^{\rm
LO})-\cos(\hat{\Delta}+\delta_{\mathrm{CP}}^{\rm HO})=1.7.$ (21)

This equation can have solutions only if

$\displaystyle 0.7\leq\cos(\hat{\Delta}+\delta_{\mathrm{CP}}^{\rm LO})\leq
1.0,$ $\displaystyle-1.0\leq\cos(\hat{\Delta}+\delta_{\mathrm{CP}}^{\rm
HO})\leq-0.7.$ (22)

From the above equation, we have the ranges of $\delta_{\mathrm{CP}}$ as

$\displaystyle-116^{\circ}\leq\delta_{\mathrm{CP}}^{\rm LO}\leq-26^{\circ}$
$\displaystyle 64^{\circ}\leq\delta_{\mathrm{CP}}^{\rm HO}\leq 161^{\circ}.$
(23)

Therefore, for the NO$\nu$A experiment, for NH and neutrino, $P_{\mu e}(\rm
LO,-116^{\circ}\leq\delta_{\mathrm{CP}}\leq-26^{\circ})$ is close to $P_{\mu
e}(\rm HO,64^{\circ}\leq\delta_{\mathrm{CP}}\leq 161^{\circ})$ Agarwalla et
al. (2013). Similar equations for $P_{\bar{\mu}\bar{e}}$ show that this
degeneracy can be removed by anti-neutrino run.

##  III Details of data analysis

In this section, we describe the methodology by which we have done data
analyses of T2K and NO$\nu$A data in Sections IV and V. We have computed
$\chi^{2}$ and $\Delta\chi^{2}$ between the data and a given theoretical
model. In section IV, we have discussed the evolution of NO$\nu$A and T2K data
with time. To do so, we have presented the analysis of the latest data from
both the experiments in standard 3-flavour oscillation scenario. In section V,
we have discussed how different BSM scenarios alleviate the tension between
the two experiments. This is done by analyzing the latest data from NO$\nu$A
and T2K in each of the different BSM frameworks. In both Sections IV and V,
the results have been presented in the form of $\Delta\chi^{2}$.

The observed and the expected number of events in the $i$-th energy bin of a
given experiment are denoted by $N_{i}^{\rm obs}$ and $N_{i}^{\rm th}$
respectively. The $\chi^{2}$ between these two distributions is calculated as

$\displaystyle\chi^{2}$ $\displaystyle=$ $\displaystyle
2\sum_{i}\left\\{(1+z)N_{i}^{\rm th}-N_{i}^{\rm exp}+N_{i}^{\rm
exp}\ln\left[\frac{N_{i}^{\rm exp}}{(1+z)N_{i}^{\rm
th}}\right]\right\\}+2\sum_{j}(1+z)N_{j}^{\rm th}+z^{2},\,\,,$

where $i$ stands for the bins for which the observed event numbers are non-
zero, $j$ stands for bins for which the observed event numbers are zero and
$z$ is the parameter defining systematic uncertainties.

The theoretical expected event numbers for each energy bin and the $\chi^{2}$
between theory and experiment have been calculated using the software GLoBES
Huber et al. (2005, 2007). To do so, we fixed the signal and background
efficiencies of each energy bin according to the Monte-Carlo simulations
provided by the experimental collaborations Himmel (2020); Acero et al.
(2021); Dunne (2020); Abe et al. (2021). We kept $\sin^{2}{\theta_{12}}$ and
${\Delta_{21}}$ at their best-fit values $0.304$ and $7.42\times 10^{-5}\,{\rm
eV}^{2}$, respectively. We varied $\sin^{2}2{\theta_{13}}$ in its $3\,\sigma$
range around its central value $0.084$ with $3.5\%$ uncertainty Dohnal (2021).
$\sin^{2}{\theta_{23}}$ has been varied in its $3\,\sigma$ range $[0.41:0.62]$
(with $2\%$ uncertainty on $\sin^{2}2{\theta_{23}}$ Esteban et al. (2019)). We
varied $|\Delta_{\mu\mu}|$ in its $3\,\sigma$ range around the MINOS best-fit
value $2.32\times 10^{-3}\,{\rm eV}^{2}$ with $3\%$ uncertainty Nichol
(2012b). The CP-violating phase $\delta_{\mathrm{CP}}$ has been varied in its
complete range $[-180^{\circ}:180^{\circ}]$. In case of BSM physics, we
modified the software to include new physics. The ranges of different new
parameters for each of the BSM scenarios have been discussed in section V.

Automatic bin-based energy smearing for the generated theoretical events has
been implemented within GLoBES Huber et al. (2005, 2007) using a Gaussian
smearing function

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

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

$\sigma(E)=\alpha E+\beta\sqrt{E}+\gamma.$ (26)

For NO$\nu$A, we used $\alpha=0.11\,(0.09)$, $\beta=\gamma=0$ for electron
(muon) like events Acero et al. (2018, 2019). For T2K, we used $\alpha=0$,
$\beta=0.075$, $\gamma=0.05$ for both electron and muon like events. For T2K,
the relevant systematic uncertainties are

  * •

$5\%$ normalisation and $5\%$ energy calibration systematics uncertainty for
$e$-like events, and

  * •

$5\%$ normalisation and $0.01\%$ energy calibration systematics uncertainty
for $\mu$-like events.

For NO$\nu$A, we used $8.5\%$ normalization and $5\%$ energy calibration
systematic uncertainties for both the $e$ like and $\mu$ like events Acero et
al. (2018). Details of systematic uncertainties have been discussed in the
GLoBES manual Huber et al. (2005, 2007).

During the calculation of $\chi^{2}$ we added (for the older, pre-2020 data)
priors on $\sin^{2}2{\theta_{13}}$, $\sin^{2}2{\theta_{23}}$, and
$|\Delta_{\mu\mu}|$, in cases where we have not included muon disappearance
data. In all other cases, priors have been added on $\sin^{2}2{\theta_{13}}$
only (to account for electron disappearance data from reactor neutrino
experiments).

We calculated the $\chi^{2}$ for both the hierarchies. Once the $\chi^{2}$s
had been calculated, we subtracted the minimum $\chi^{2}$ from them to
calculate the $\Delta\chi^{2}$. The parameter values and hierarchy, for which
the $\Delta\chi^{2}=0$, is called the best-fit point.

##  IV Evolution of the tension between NO$\nu$A and T2K data

In this section, we consider the appearance data of T2K and NO$\nu$A in both
$\nu/\bar{\nu}$ modes and discuss how they give rise to degenerate solutions.
We will highlight the growing tension between these appearance data when they
are interpreted within the three flavour oscillation paradigm. In the next
section, we will consider new physics scenarios which can reduce that tension.

###  IV.1 Evolution of the NO$\nu$A data

In 2017, NO$\nu$A published first result with combined analysis of $\nu_{e}$
appearance and $\nu_{\mu}$ disappearance data, corresponding to $6.05\times
10^{20}$ POT Adamson et al. (2017). This analysis gave three almost degenerate
solutions

  * •

NH, $\sin^{2}{\theta_{23}}=0.4$, $\delta_{\mathrm{CP}}=-90^{\circ}$ (NH, LO,
$-90^{\circ}$),

  * •

NH, $\sin^{2}{\theta_{23}}=0.62$, $\delta_{\mathrm{CP}}=135^{\circ}$ (NH, HO,
$135^{\circ}$), and

  * •

IH, $\sin^{2}{\theta_{23}}=0.62$, $\delta_{\mathrm{CP}}=-90^{\circ}$ (NH, HO,
$-90^{\circ}$).

In ref. Bharti et al. (2018), an effort was made to understand this degeneracy
in the NO$\nu$A data from 2017. To do so, the authors of ref. Bharti et al.
(2018) first calculated the expected $\nu_{e}$ appearance event numbers for
vacuum oscillation, maximal ${\theta_{23}}=45^{\circ}$, and
$\delta_{\mathrm{CP}}=0$ for $6.05\times 10^{20}$ POT. This case was labelled
as ${}^{\prime}000^{\prime}$. Then they considered changes in this number due
to a) matter effect, b) octant of non-maximal ${\theta_{23}}$, and c) large
value of $\delta_{\mathrm{CP}}$. The parameter value for which $P_{\mu e}$ is
increased (decreased) was labelled as ${}^{\prime}+^{\prime}$
(${}^{\prime}-^{\prime}$). First, one change at a time was introduced in the
following manner:

  * •

NH which increases $P_{\mu e}$ (labelled as ${}^{\prime}+^{\prime}$) or IH
which decreases it (labelled as ${}^{\prime}-^{\prime}$),

  * •

HO which increases $P_{\mu e}$ (labelled as ${}^{\prime}+^{\prime}$) or LO
which decreases it (labelled as ${}^{\prime}-^{\prime}$),

  * •

$\delta_{\mathrm{CP}}=-90^{\circ}$ which increases $P_{\mu e}$ (labelled as
${}^{\prime}+^{\prime}$) or $\delta_{\mathrm{CP}}=+90^{\circ}$ which decreases
it (labelled as ${}^{\prime}-^{\prime}$).

The event numbers for $6.05\times 10^{20}$ POT were calculated using the
software GLoBES Huber et al. (2005, 2007). Other parameters were fixed at
following constant values: ${\Delta_{21}}=7.5\times 10^{-5}\,{\rm eV}^{2}$,
$\sin^{2}{\theta_{12}}=0.306$, ${\Delta_{31}}(\rm NH)=2.74\times 10^{-3}\,{\rm
eV}^{2}$, ${\Delta_{31}}(\rm IH)=-2.65\times 10^{-3}\,{\rm eV}^{2}$, and
$\sin^{2}2{\theta_{13}}=0.085$. The values of ${\Delta_{31}}$(NH) and
${\Delta_{31}}$(IH) were taken from the analysis of NO$\nu$A $\nu_{\mu}$
disappearance data. The results are presented in table 3. From this table, it
is obvious that the increase (decrease) in $P_{\mu e}$ for any single
${}^{\prime}+^{\prime}$ (${}^{\prime}-^{\prime}$) change in the unknown
parameters is essentially same.

Hierarchy-$\sin^{2}{\theta_{23}}$-$\delta_{\mathrm{CP}}$ | Label |  $\nu_{e}$ Appearance events
---|---|---
Vacuum-$0.5$-$0$ | $000$ | $26.49$
NH-$0.5$-$0$ | $+00$ | $31.28$
IH-$0.5$-$0$ | $-00$ | $21.08$
Vac-$0.5$-$-90^{\circ}$ | $00+$ | $31.04$
Vac-$0.5$-$+90^{\circ}$ | $00-$ | $21.18$
Vac-$0.62$-$0$ | $0+0$ | $32.88$
Vac-$0.4$-$0$ | $0-0$ | $24.59$
Table 3: Expected $\nu_{e}$ appearance events of NO$\nu$A for $6.05\times
10^{20}$ POT. They are listed for the reference point and for change of one
unknown parameter at a time. NO$\nu$A observed $33$ $\nu_{e}$ appearance
events in 2017.

Next all eight possible combinations in the changes of the three unknown
parameters were considered. All three unknown parameters can shift in a way
that each change leads to an increase in $P_{\mu e}$. This can be lebelled as
${}^{\prime}+++^{\prime}$. In this case, one gets the maximum number of
$\nu_{e}$ appearance events. Another case is that two of the unknown
parameters change to increase $P_{\mu e}$, whereas the third one decreases it.
This can happen in three possible ways, labelled as ${}^{\prime}+-+^{\prime}$,
${}^{\prime}++-^{\prime}$, and ${}^{\prime}-++^{\prime}$. These three
combinations lead to a moderate increase in $\nu_{e}$ appearance events
compared to the reference ${}^{\prime}000^{\prime}$ case. Similarly a moderate
decrease in $\nu_{e}$ appearance event numbers compared to the reference case
can occur due to increase in $P_{\mu e}$ by one unknown parameter and decrease
by the other two. This can occur in three possible ways:
${}^{\prime}+--^{\prime}$, ${}^{\prime}-+-^{\prime}$ and
${}^{\prime}--+^{\prime}$. Finally there is a possible case where each of the
three changes lowers $P_{\mu e}$ and we get the minimum number of $\nu_{e}$
appearance events. This can be labelled as ${}^{\prime}---^{\prime}$. In table
4, the number of $\nu_{e}$ appearance events for $6.05\times 10^{20}$ POT and
for all the eight combinations, mentioned above, have been listed. From the
table, we can see that the number of events for ${}^{\prime}++-^{\prime}$,
${}^{\prime}+-+^{\prime}$ and ${}^{\prime}-++^{\prime}$ are nearly the same. A
similar statement can be made about ${}^{\prime}+--^{\prime}$,
${}^{\prime}-++^{\prime}$ and ${}^{\prime}--+^{\prime}$. The predicted number
of events for ${}^{\prime}+++^{\prime}$, and ${}^{\prime}---^{\prime}$ are
totally unique. Thus, the eight-fold degeneracy, present before the
${\theta_{13}}$ was measured, breaks itself down to $1+3+3+1$ pattern after
the precise measurement of ${\theta_{13}}$. The 2017 NO$\nu$A data saw a
moderate increase in $\nu_{e}$ appearance events, compared to the expected
event numbers for the reference case ${}^{\prime}000^{\prime}$. Hence, there
was a three-fold degeneracy in the NO$\nu$A solutions. In table 5, the
expected number of $\nu_{e}$ appearance events for $6.05\times 10^{20}$ POT at
each of the NO$\nu$A solutions have been listed. The predictions for two NH
solutions matched the experimentally observed event numbers $33$. The
prediction of the IH solution (which was $0.5\,\sigma$ away from the NH
solutions) was lower by $3$ (half the statistical uncertainty in the expected
number). The occurrence of three-fold degeneracy in the NO$\nu$A data, based
on the inherent degeneracy in $P_{\mu e}$ was also discussed in ref. Lindner
et al. (2018).

Hierarchy-$\sin^{2}{\theta_{23}}$-$\delta_{\mathrm{CP}}$ | Label |  $\nu_{e}$ Appearance events
---|---|---
NH-$0.62$-$-90^{\circ}$ | $+++$ | $43.67$
NH-$0.4$-$-90^{\circ}$ | $+-+$ | $33.54$
NH-$0.62$-$+90^{\circ}$ | $++-$ | $33.04$
IH-$0.62$-$-90^{\circ}$ | $-++$ | $30.94$
NH-$0.4$-$+90^{\circ}$ | $+--$ | $22.79$
IH-$0.4$-$-90^{\circ}$ | $--+$ | $24.47$
IH-$0.62$-$+90^{\circ}$ | $-+-$ | $22.63$
IH-$0.4$-$+90^{\circ}$ | $---$ | $16.07$
Table 4: Expected $\nu_{e}$ appearance events of NO$\nu$A for $6.05\times 10^{20}$ POT, and for eight different combinations of unknown parameters. Hierarchy-$\sin^{2}{\theta_{23}}$-$\delta_{\mathrm{CP}}$ | Label |  $\nu_{e}$ Appearance events
---|---|---
NH-$0.404$-$-86^{\circ}$ | $+-+$ | $33.55$
NH-$0.62$-$+135^{\circ}$ | $++-$ | $34.36$
IH-$0.62$-$-90^{\circ}$ | $-++$ | $30.94$
Table 5: Expected $\nu_{e}$ appearance events of NO$\nu$A for $6.05\times
10^{20}$ POT, and for the three solutions in ref. Adamson et al. (2017)

In ref. Bharti et al. (2018) the authors also considered the possible
resolution of the three-fold degeneracy with the anti-neutrino run. One can
obtain the anti-neutrino oscillation probability $P_{\bar{\mu}\bar{e}}$ by
reversing the signs of matter term $A$, and $\delta_{\mathrm{CP}}$ in eq.
(16). $P_{\bar{\mu}\bar{e}}$ decreases (increases) for NH (IH). Similarly,
$\delta_{\mathrm{CP}}$ in the UHP (LHP) increases (decreases)
$P_{\bar{\mu}\bar{e}}$. But, we will continue to label the NH (IH) as
${}^{\prime}+^{\prime}$ (${}^{\prime}-^{\prime}$). In the same way,
$\delta_{\mathrm{CP}}$ in the LHP (UHP) will be labelled as
${}^{\prime}+^{\prime}$ (${}^{\prime}-^{\prime}$). However, it should be noted
that ${}^{\prime}+^{\prime}$ (${}^{\prime}-^{\prime}$) sign in hierarchy and
$\delta_{\mathrm{CP}}$ leads to an decrease (increase) in
$P_{\bar{\mu}\bar{e}}$. But for the octant of ${\theta_{23}}$,
${}^{\prime}+^{\prime}$ (${}^{\prime}-^{\prime}$) sign leads to an increase
(decrease) in $P_{\bar{\mu}\bar{e}}$. Now, for the ${}^{\prime}++-^{\prime}$
solution, $P_{\bar{\mu}\bar{e}}$ decreases due to the hierarchy, and increases
due to the octant of ${\theta_{23}}$ and $\delta_{\mathrm{CP}}$. Similarly,
for ${}^{\prime}-++^{\prime}$ solution, $P_{\bar{\mu}\bar{e}}$ increases due
to the hierarchy and the octant of ${\theta_{23}}$, and decreases due to
$\delta_{\mathrm{CP}}$. Hence these two solutions are degenerate for anti-
neutrino data also, and the NO$\nu$A $\bar{\nu}_{e}$ appearance data would not
be able to distinguish between them. However, for the third solution
${}^{\prime}+-+^{\prime}$, $P_{\bar{\mu}\bar{e}}$ decreases due to all three
unknown parameters. The expected $\bar{\nu}_{e}$ appearance events for this
case would be the smallest. In principle, the NO$\nu$A $\bar{\nu}_{e}$
appearance data for $6.05\times 10^{20}$ POT would be able to distinguish this
solution for the other two. Since the expected number of events for this
particular scale is rather small, as can be seen from table 6, the statistical
uncertainty would be large, and NO$\nu$A $\bar{\nu}_{e}$ appearance data with
$6.05\times 10^{20}$ POT would not be able to distinguish this solution from
the other two at the $3\,\sigma$ level.

Hierarchy-$\sin^{2}{\theta_{23}}$-$\delta_{\mathrm{CP}}$ | Label |  $\bar{\nu}_{e}$ Appearance events
---|---|---
NH-$0.404$-$-86^{\circ}$ | $+-+$ | $6.85$
NH-$0.62$-$+135^{\circ}$ | $++-$ | $11.78$
IH-$0.62$-$-90^{\circ}$ | $-++$ | $12.86$
Table 6: Expected $\bar{\nu}_{e}$ appearance events of NO$\nu$A for
$6.05\times 10^{20}$ POT, and for the three solutions in ref. Adamson et al.
(2017)

In the Neutrino 2018 conference, NO$\nu$A published results after analysing
data corresponding to $8.85\times 10^{20}$ ($6.9\times 10^{20}$) POTs in the
neutrino (anti-neutrino) mode Sanchez (2018); Acero et al. (2018). They found
out the best-fit point at $\delta_{\mathrm{CP}}=30.6^{\circ}$
($-95.4^{\circ}$), $\sin^{2}{\theta_{23}}=0.58\pm 0.03$ ($0.58\pm 0.04$),
$\Delta_{32}=2.51^{+0.12}_{-0.08}\times 10^{-3}\,{\rm eV}^{2}$ ($-2.56\times
10^{-3}\,{\rm eV}^{2}$) for NH (IH). NH was preferred over IH at $1.8\,\sigma$
C.L. Also $\delta_{\mathrm{CP}}=90^{\circ}$ in the IH was excluded at more
than $3\,\sigma$ C.L. The neutrino disappearance data were consistent with the
maximal ${\theta_{23}}$, whereas the $\bar{\nu}_{\mu}$ disappearance data
preferred a non-maximal mixing Nizam et al. (2019). Therefore, there was a
mild tension between the two different data sets of the same experiment.
However, since the statistics from the anti-neutrino disappearance data were
very low, this tension was not statistically significant. In the case of
appearance data, the expected number of $\nu_{e}$ ($\bar{\nu}_{e}$) appearance
events at the reference point ${}^{\prime}000^{\prime}$ was $39$ ($15.5$)
Nizam et al. (2019). The observed $\nu_{e}$ ($\bar{\nu}_{e}$) appearance
events were $58$ ($18$). Hence, there was a moderate excess in both these
channels. As we have already seen, the moderate excess in both $\nu_{e}$ and
$\bar{\nu}_{e}$ appearance channels is possible when the $20\%$ change due to
the hierarchy and $\delta_{\mathrm{CP}}$ cancels each other, and there is an
increase in both channels induced by ${\theta_{23}}$ in the HO. We have
already labelled these possibilities by A) ${}^{\prime}++-^{\prime}$ and B)
$-++$. These are two of the three degenerate solutions of the previous 2017
data. The other degenerate solution of 2017 data, namely
${}^{\prime}+-+^{\prime}$ was ruled out by 2018 data because although this
solution leads to moderate excess in the $\nu_{e}$ appearance events, it
causes minimum number of $\bar{\nu}_{e}$ appearance events. In ref. Nizam et
al. (2019), an analysis of the $\nu_{e}$ and $\bar{\nu}_{e}$ appearance data
from 2018 was performed. It was found out that there were two degenerate best-
fit solutions: i) NH, $\sin^{2}{\theta_{23}}=0.65$,
$\delta_{\mathrm{CP}}=120^{\circ}$, and ii) IH, $\sin^{2}{\theta_{23}}=0.67$,
$\delta_{\mathrm{CP}}=-50^{\circ}$. The first solution is in the form of A)
and the second solution is in the form of B). The best-fit points given by the
NO$\nu$A collaboration were also in these two forms. However, because of the
inclusion of the disappearance data in the analysis of the NO$\nu$A
collaboration, a smaller value, compared to those mentioned above, of
$\sin^{2}{\theta_{23}}$ was obtained.

In 2020, NO$\nu$A published the analysis of their data corresponding to
$1.36\times 10^{21}$ ($1.25\times 10^{21}$) POT in the $\nu$ ($\bar{\nu}$)
mode Himmel (2020); Acero et al. (2021). The best-fit points were
$\Delta_{32}=+(2.41\pm 0.07)\times 10^{-3}\,{\rm eV}^{2}$ ($-2.45\times
10^{-3}\,{\rm eV}^{2}$), $\sin^{2}{\theta_{23}}=0.57^{+0.03}_{-0.04}$
($0.56$), $\delta_{\mathrm{CP}}/\pi=0.82^{+0.27}_{-0.87}$ ($1.52$). In ref.
Rahaman and Raut (2021), a detailed analysis of the NO$\nu$A data has been
done. The expected $\nu_{e}$ and $\bar{\nu}_{e}$ appearance event numbers at
the reference point ${}^{\prime}000^{\prime}$ are $76.14$ and $32.93,$
respectively. The observed event numbers in these channels are $82$ and $33$,
respectively. Thus, there is a moderate excess in the $\nu_{e}$ appearance
channel, and this excess can be explained by, as explained before, three
possible solutions: A) ${}^{\prime}+-+^{\prime}$, B) ${}^{\prime}++-^{\prime}$
and C) ${}^{\prime}-++^{\prime}$. As for the $\bar{\nu}_{e}$ appearance data,
the observed event number $33$ is consistent with the
${}^{\prime}000^{\prime}$ solution. However, due to the lack of statistics in
the $\bar{\nu}_{e}$ appearance channel, other solutions are also allowed at
the $1\,\sigma$ C.L.. Exceptions are the cases labelled as
${}^{\prime}-+-^{\prime}$, and ${}^{\prime}+-+^{\prime}$, since these cases
lead to the maximum and minimum number of expected $\bar{\nu}_{e}$ appearance
events, respectively. Thus the solution in the form of B) is excluded when
$\nu_{e}$ and $\bar{\nu}_{e}$ appearance data are analysed together. The
result is shown in Fig. 4. Best-fit points are of the forms A) and C). The
best-fit points obtained by the NO$\nu$A collaboration after analysing
appearance and disappearance data together are of the same forms as well.

Figure 4: Expected allowed regions in the
$\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane from the NO$\nu$A
appearance data, as given in ref. Himmel (2020); Acero et al. (2021), for both
the $\nu_{e}$ and $\bar{\nu}_{e}$ channels. In the left (right) panel, the
hierarchy is assumed to be NH (IH). The best-fit point is at NH with a minimum
$\chi^{2}=6.77$ for 12 energy bins. IH has a minimum $\Delta\chi^{2}=0.2$.

###  IV.2 Evolution of the T2K data

In 2013 T2K published their first analysis of $\nu_{e}$ appearance data
corresponding to $6.57\times 10^{20}$ POT Abe et al. (2014a, b). They found
out their best-fit point at NH and $\delta_{\mathrm{CP}}=-90^{\circ}$. Both
the hierarchies with $\delta_{\mathrm{CP}}$ in the LHP were allowed at
$2\,\sigma$ C.L., whereas $\delta_{\mathrm{CP}}$ values in the UHP were
disfavoured at $2\,\sigma$ C.L. for both the hierarchies. From the $\nu_{\mu}$
disappearance data published in 2014 Abe et al. (2014b), ${\theta_{23}}$ was
found to be close to the maximal mixing. In ref. Bharti et al. (2016), a
detailed analysis of the physics potential of NO$\nu$A in the presence of the
information from T2K data were made. It was shown that if the hierarchy and
$\delta_{\mathrm{CP}}$ are in favourable combinations, T2K data have no effect
on the hierarchy determination potential of NO$\nu$A. Among the unfavourable
hierarchy-$\delta_{\mathrm{CP}}$ combinations, T2K data picked out the correct
(incorrect) solution from the two degenerate solutions allowed by the NO$\nu$A
data for the hierarchy being IH (NH) and $\delta_{\mathrm{CP}}$ in the LHP
(UHP). Therefore it was concluded that if the combination of NO$\nu$A and T2K
prefer IH and $\delta_{\mathrm{CP}}$ in the LHP as the correct solution, one
needs to be careful, because the actual combination might be NH and
$\delta_{\mathrm{CP}}$ in the LHP. We will find out that this prediction in
ref. Bharti et al. (2016) was quite accurate in the context of the latest
published data from NO$\nu$A and T2K.

In 2018, T2K published the analysis of their data with $14.7\times 10^{20}$
($7.6\times 10^{20}$) POT in the neutrino (anti-neutrino) mode Abe et al.
(2018b). In ref. Nizam et al. (2019) a detailed analysis of T2K disappearance
and appearance data was done separately. It was found out that the analysis of
T2K disappearance data gave the best-fit point at
$\sin^{2}{\theta_{23}}=0.51$. The $3\,\sigma$ C.L. on $\sin^{2}{\theta_{23}}$
was constrained in the range $[0.43:0.6]$. The constraint on
$\sin^{2}{\theta_{23}}$ was valid for all values of $\delta_{\mathrm{CP}}$,
since the disappearance data do not have $\delta_{\mathrm{CP}}$ sensitivity.
As for the appearance data it was found out that the expected $\nu_{e}$
appearance event number at the reference point ${}^{\prime}000^{\prime}$ for
the given neutrino run was found out to be $60$. Inclusion of matter effect
changed the number by $4$, and inclusion of maximal CP violation changed the
number by $11$. Therefore, for NH and $\delta_{\mathrm{CP}}=-90^{\circ}$, the
expected number was increased to 80 Nizam et al. (2019). T2K observed $89$
$\nu_{e}$ events. Hence, the $\nu_{e}$ appearance data of T2K pulled the
$\sin^{2}{\theta_{23}}$ to a value larger than $0.5$. An analysis of T2K
$\nu_{e}$ appearance data was done in ref. Nizam et al. (2019). The
$\bar{\nu}_{e}$ data were not included because the observed number of events
in this channel was too small to have any statistical significance. The best-
fit point was found to be at $\sin^{2}{\theta_{23}}=0.63$, although
$\sin^{2}{\theta_{23}}=0.5$ was allowed at $1\,\sigma$ C.L. Therefore there
was a mild tension between the T2K appearance and disappearance data. The T2K
collaboration, after analysing the appearance and disappearance data together
got the best-fit point at $\sin^{2}{\theta_{23}}=0.53$ Abe et al. (2018b).
Because of the larger statistical weight of disappearance data, the final
value of $\sin^{2}{\theta_{23}}$ was determined by the disappearance data and
the final value of $\sin^{2}{\theta_{23}}$ was found to be close to the
maximal value. Because of the large excess in the observed $\nu_{e}$
appearance events in T2K, the $\delta_{\mathrm{CP}}$ was found to be in the
vicinity of $-90^{\circ}$. For $\delta_{\mathrm{CP}}$ in the UHP, the expected
number of $\nu_{e}$ events was smaller than that at the reference point. Thus
$\delta_{\mathrm{CP}}$ in the UHP was highly disfavoured. This data also
disfavoured IH, because the corresponding matter effect reduced the number of
expected $\nu_{e}$ events. IH with $\delta_{\mathrm{CP}}=-90^{\circ}$ was
barely allowed at $2\,\sigma$ C.L. Abe et al. (2018b).

In 2019, T2K published data corresponding to $14.9\times 10^{20}$ ($16.4\times
10^{20}$) POT in the neutrino (anti-neutrino) mode. The best-fit point was at
$\delta_{\mathrm{CP}}/\pi=-1.89^{+0.70}_{-0.58}$, ($-1.38^{+0.48}_{-0.54}$)
$\sin^{2}{\theta_{23}}=0.53^{+0.03}_{-0.04}$ for NH (IH) Abe et al. (2020).
They also found out that $\delta_{\mathrm{CP}}=0$ was excluded at $99\%$ C.L.

In 2020, T2K published their latest data with $1.97\times 10^{21}$
($1.63\times 10^{21}$) POT in the neutrino (anti-neutrino) mode Dunne (2020);
Abe et al. (2021). The best-fit point is at
$\delta_{\mathrm{CP}}/\pi=-2.14^{+0.90}_{-0.69}$ ($-1.26^{+0.61}_{-0.69}$),
$\sin^{2}{\theta_{23}}=0.512^{+0.045}_{-0.042}$ for NH (IH). This result can
be explained with the change in event number due to the change in unknown
parameter values from their reference point values Rahaman and Raut (2021). At
the reference point ${}^{\prime}000^{\prime}$, the expected number of
$\nu_{e}$ ($\bar{\nu}_{e}$) appearance event is $78$ ($19$). T2K observed
$113$ ($15$) $\nu_{e}$ ($\bar{\nu}_{e}$) appearance events. The large excess
in the $\nu_{e}$ appearance channel observed by T2K can only be explained by
making the choice of unknowns to be ${}^{\prime}+++^{\prime}$. Choosing NH
leads to only an $8\%$ boost and we need ${\theta_{23}}$ to have a large value
in HO to explain the large excess $\nu_{e}$ events. But the disappearance data
limits $\sin^{2}{\theta_{23}}\leq 0.59$. Given that only about $20\%$ boost is
possible from the hierarchy and octant, $\delta_{\mathrm{CP}}$ has to be
firmly anchored around $\delta_{\mathrm{CP}}=-90^{\circ}$ to accommodate the
large excess in the $\nu_{e}$ appearance channel. The $\bar{\nu}_{e}$
appearance data see a reduction in the observed events. This reduction is
consistent with the event numbers expected from ${}^{\prime}+++^{\prime}$
choice of the unknowns. However, the number of events in this channel is too
small to have any statistical significance. In ref. Rahaman and Raut (2021),
an analysis of the T2K $\nu_{e}$ and $\bar{\nu}_{e}$ appearance data has been
done and the result is presented in Fig. 5. It is obvious that the large
excess in the T2K $\nu_{e}$ appearance channel is responsible for
$\delta_{\mathrm{CP}}$ being close to $-90^{\circ}$, and this is what leads to
the tension between the NO$\nu$A and T2K data.

Figure 5: Expected allowed regions in the
$\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane from the appearance data
from T2K, as given in ref. Dunne (2020); Abe et al. (2021), in both the
$\nu_{e}$ and $\bar{\nu}_{e}$ channels. In the left (right) panel, the
hierarchy is assumed to be NH (IH). The best-fit point is at NH with a minimum
$\chi^{2}=20.23$ for 18 energy bins.

###  IV.3 Combined analysis of NO$\nu$A and T2K data

Joint analyses have been planned between the NO$\nu$A and T2K collaborations
with the aim of obtaining better constraints on the oscillation parameters due
to resolved degeneracy and to understand the non-trivial systematic
correlations between them Berns (2021). However, independent joint analyses of
the two experiments have been done in ref. Kelly et al. (2021), and it has
been found out that the combined analysis prefers IH over NH. The result has
been shown in Fig. 6. Note that there is no overlap between the $1\,\sigma$
individual allowed regions of T2K and NO$\nu$A. This also leads to an
extremely tiny $1\,\sigma$ allowed region for the combined analysis.

Figure 6: Allowed region in the $\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$
plane after analysing NO$\nu$A and T2K complete data sets. The left (right)
panel represents the test hierarchy to be 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
50 (88) bins is 48.65 (95.85) and it occurs at NH. For the combined analysis,
the minimum $\chi^{2}$ with 138 bins is 147.14, which occurs for IH.

##  V Resolution of the NO$\nu$A -T2K tension with new physics

One of the possible reasons for the tension between NO$\nu$A and T2K is the
existence of new physics in the neutrino oscillation. Effect of new physics on
the determination of unknown oscillation parameters in the long-baseline
accelerator neutrino experiments have been studied in details in the
literature Argüelles et al. (2020). In recent times, efforts have been made to
resolve the tension between NO$\nu$A and T2K experiments with the help of non-
standard physics Miranda et al. (2021); Rahaman (2021); Chatterjee and Palazzo
(2021); Denton et al. (2021). In this section, we will discuss the present
status of different non-standard physics to resolve this tension.

###  V.1 Non-unitary mixing

The anomalies observed in LSND Aguilar-Arevalo et al. (2001) and MiniBooNE
Aguilar-Arevalo et al. (2007) experiments can be explained with the existence
of one or more “sterile” neutrino states with mass at or below a few eV
Abazajian et al. (2012)111A recent results from the MicroBooNE experiment
Abratenko et al. (2021) rules out any excess electron-like events at $94.8\%$
C.L. However, they do not rule out the complete parameter space suggested by
the MiniBooNE experiment and other data, nor do they probe the $\nu_{e}$
interpretation of MiniBooNE result in a model independent way Argüelles et al.
(2021).. The effect of sterile neutrino on the long-baseline accelerator
neutrino experiments have been discussed in detail in Ref. Gupta et al.
(2018); Chatla et al. (2018); Choubey et al. (2018, 2017); Berryman et al.
(2015). If the sterile neutrinos exist as iso-singlet neutral heavy leptons
(NHL), then in the minimum extension of the standard model, they do not take
part in neutrino oscillation. However, their admixture in the charged current
weak interactions affects neutrino oscillation. In such scenarios, neutrino
oscillation will be governed by a non-unitary mixing matrix. The mixing
matrix, in this case, can be parameterised as Escrihuela et al. (2015)

$N=N_{\rm 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_{\rm PMNS}\,.$ (27)

Here $U_{\rm PMNS}$ is the standard unitary PMNS mixing matrix. The diagonal
(off-diagonal) terms of $N_{\rm NP}$ matrix are real (complex). To allow the
effect of non-unitarity, the diagonal terms of $N_{\rm NP}$ matrix must
deviate from unity, and/or the off-diagonal terms must deviate from $0$. The
present constraints on non-unitary parameters at $3\,\sigma$ C.L. are
Escrihuela et al. (2017):

$\displaystyle\alpha_{00}>0.93\,;\,\alpha_{11}>0.95\,;\,\alpha_{22}>0.61$
$\displaystyle|\alpha_{10}|<3.6\times 10^{-2}\,;\,|\alpha_{20}|<1.3\times
10^{-1}\,;\,|\alpha_{21}|<2.1\times 10^{-2}\,.$ (28)

The calculation of oscillation probability in the presence of matter effect in
the case of non-unitary mixing has been discussed in Ref. Escrihuela et al.
(2015); Miranda et al. (2021). The effects of non-unitary mixing on the
determination of the unknown oscillation parameters in the present and future
long-baseline accelerator neutrino experiments have been discussed in
literature Ge et al. (2017); Escrihuela et al. (2015); Soumya and Rukmani
(2018); Fong et al. (2019); Verma and Bhardwaj (2018).

In the context of non-unitary mixing, the mixing between flavour states and
mass eigenstates can be written as

$|\nu_{\beta}\rangle=\sum_{i=1}^{3}N_{\beta i}|\nu_{i}\rangle,$ (29)

where $\beta$ denotes the flavour states and $i$ denotes the mass eigenstates.
The evolution of neutrino mass eigenstates during the propagation, can be
written as

$i\frac{d}{dt}|\nu_{i}\rangle=H_{\rm vac}|\nu_{i}\rangle,$ (30)

where $H_{\rm vac}$ is the Hamiltonian in vacuum and it is defined as

$H_{\rm vac}=\frac{1}{2E}\left[{\begin{array}[]{ccc}0&0&0\\\
0&{\Delta_{21}}&0\\\ 0&0&{\Delta_{31}}\end{array}}\right]$ (31)

The non-unitary neutrino oscillation probability in vacuum can be written as

$\displaystyle P_{\mu e}^{\rm NU}({\rm vac})$ $\displaystyle=$
$\displaystyle\sum_{i,j}^{3}N_{\mu i}^{*}N_{ei}N_{\mu
j}N_{ej}^{*}-4\sum_{j>i}^{3}Re\left[N_{\mu j}N_{ej}N_{\mu
i}N_{ei}^{*}\right]\sin^{2}\left(\frac{\Delta_{ji}L}{4E}\right)$ (32)
$\displaystyle+2\sum_{j>i}^{3}Im\left[N_{\mu j}^{*}N_{ej}N_{\mu
i}N_{ei}\right]\sin\left(\frac{\Delta_{ji}L}{2E}\right).$

If written explicitly, neglecting the cubic products of $\alpha_{10}$,
$\sin{\theta_{13}}$, and ${\Delta_{21}}$, the oscillation probability takes
the form Escrihuela et al. (2015)

$P_{\mu e}^{\rm NU}({\rm vac})=(\alpha_{00}\alpha_{11})^{2}P_{\mu e}^{\rm
SO}+\alpha_{00}^{2}\alpha_{11}|\alpha_{10}|P_{\mu e}^{\rm
I}+\alpha_{00}^{2}|\alpha_{10}|^{2}.$ (33)

$P_{\mu e}^{\rm SO}$ is the standard three-flavour unitary neutrino
oscillation probability in vacuum and can be written as

$\displaystyle P_{\mu e}^{\rm SO}$ $\displaystyle=$ $\displaystyle
4\cos^{2}{\theta_{12}}\cos^{2}{\theta_{23}}\sin^{2}{\theta_{12}}\sin^{2}\left(\frac{{\Delta_{21}}L}{4E}\right)$
(34)
$\displaystyle+4\cos^{2}{\theta_{13}}\sin^{2}{\theta_{13}}\sin^{2}{\theta_{23}}\sin^{2}\left(\frac{{\Delta_{31}}L}{4E}\right)$
$\displaystyle+\sin(2{\theta_{12}})\sin{\theta_{13}}\sin(2{\theta_{23}})\sin\left(\frac{{\Delta_{21}}L}{2E}\right)\sin\left(\frac{{\Delta_{31}}L}{4E}\right)\cos\left(\frac{{\Delta_{31}}L}{4E}-I_{123}\right),$

and

$\displaystyle P_{\mu e}^{\rm I}$ $\displaystyle=$
$\displaystyle-2\left[\sin(2{\theta_{13}})\sin{\theta_{23}}\sin\left(\frac{{\Delta_{31}}L}{4E}\right)\sin\left(\frac{{\Delta_{31}}L}{4E}+I_{\rm
NP}-I_{012}\right)\right]$ (35)
$\displaystyle-\cos{\theta_{13}}\cos{\theta_{23}}\sin(2{\theta_{12}})\sin\left(\frac{{\Delta_{21}}L}{2E}\right)\sin(I_{\rm
NP}),$

where $I_{012}=-\delta_{\mathrm{CP}}=\phi_{10}-\phi_{20}+\phi_{21}$,
$I_{123}=\phi_{21}-\phi_{31}+\phi_{32}$ and $I_{\rm NP}=\phi_{10}-{\rm
Arg}(\alpha_{10})$. $\phi_{ij}$’s are the phases associated with the off-
diagonal terms $\alpha_{ij}=|\alpha_{ij}|e^{i\phi_{ij}}$. It should be noted
that non-unitary parameters $\alpha_{00}$, $\alpha_{11}$, and $\alpha_{10}$
have the most significant effects on $P_{\mu e}^{\rm NU}({\rm vac})$.

While propagating through matter, the neutrinos undergo forward scattering and
the neutrino oscillation probability gets modified due to interaction
potential between neutrino and matter. In case of non-unitary mixing, the CC
and NC terms in the interaction Lagrangian is given as

$\mathcal{L}=V_{\rm
CC}\sum_{i,j}N_{ei}^{*}N_{ej}\bar{\nu}_{i}\gamma^{0}\nu_{j}+V_{\rm
NC}\sum_{\alpha,i,j}N_{\alpha i}^{*}N_{\alpha
j}\bar{\nu}_{i}\gamma^{0}\nu_{j},$ (36)

where $V_{\rm CC}=\sqrt{2}G_{F}N_{e}$, and $V_{\rm NC}=-G_{F}N_{n}/\sqrt{2}$
are the potentials for CC and NC interactions, respectively. Therefore, the
effective Hamiltonian becomes

$H_{\rm matter}^{\rm NU}=\frac{1}{2E}\left[{\begin{array}[]{ccc}0&0&0\\\
0&{\Delta_{21}}&0\\\
0&0&{\Delta_{31}}\end{array}}\right]+\frac{1}{2E}N^{\dagger}\left[{\begin{array}[]{ccc}V_{\rm
CC}+V_{\rm NC}&0&0\\\ 0&V_{\rm NC}&0\\\ 0&0&V_{\rm NC}\end{array}}\right]N\,.$
(37)

The non-unitary neutrino oscillation probability, after neutrinos travel
through a distance $L$, is given as

$P_{\alpha\beta}(E,L)=|\langle\nu_{\beta}|\nu_{\alpha}(L)\rangle|^{2}=\left|\left(Ne^{-iH_{\rm
matter}^{\rm NU}L}N^{\dagger}\right)_{\beta\alpha}\right|^{2}\,.$ (38)

A detailed description of non-unitary neutrino oscillation probability in the
presence of matter effect has been discussed in Ref. Miranda et al. (2021),
where an effort has been made to resolve the tension between the two
experiments with non-unitary mixing. To do so, the authors first analysed the
present data from NO$\nu$A and T2K with standard unitary oscillation
hypothesis. We have already presented the result in Fig. 6 using the 2020
data. As mentioned in the caption of Fig. 6, the minimum $\chi^{2}$ for
NO$\nu$A (T2K) is $48.65$ ($95.85$) for $50$ ($88$) energy bins. For the
combined analysis, the minimum $\chi^{2}$ was $147.14$ for $138$ energy bins.
In the next step, data were analysed with non-unitary mixing hypothesis. To do
so, only the effects of $\alpha_{00}$, $|\alpha_{10}|$, and $\alpha_{11}$ were
considered, as these three parameters have the most significant effects on the
oscillation probabilities $P_{\mu e}$, and $P_{\bar{\mu}\bar{e}}$. Only those
values were chosen for which the condition
$|\alpha_{10}|\leq\sqrt{(1-\alpha_{00}^{2})(1-\alpha_{11}^{2})}$ Antusch and
Fischer (2014); Escrihuela et al. (2017) is satisfied. All other non-unitary
parameters have been kept fixed at their unitary values. The software GLoBES
Huber et al. (2005, 2007) was used to analyse the data for both standard and
non-unitary oscillations. In the later case, the software was modified to
include non-unitary mixing. We have already talked about the standard unitary
parameter values in section III. After analysing the data with non-unitary
mixing hypothesis, the minimum $\chi^{2}$ was found to be $45.88$ ($93.36$)
for NO$\nu$A (T2K). For the combined analysis, the minimum $\chi^{2}$ was
$142.72$. Thus it can be said that each of the two experiments individually
prefer non-unitary mixing over unitary mixing by $1\,\sigma$ C.L., and the
combined analysis rules out unitary mixing at $2\,\sigma$ C.L. The result has
been presented in Fig. 7. It can be observed from the figure that there is a
large overlap between the $1\,\sigma$ allowed regions of the two experiments
for both the hierarchies. Both the experiments lose their hierarchy
sensitivity when the data are analysed with non-unitary mixing. For the NH
best-fit point, NO$\nu$A (T2K) prefers ${\theta_{23}}$ to be in LO (HO), but a
nearly degenerate best-fit point occurs at HO (LO). Thus the experiments lose
their octant determination sensitivity as well. Best-fit points at IH coincide
for both the experiments.

Figure 7: Allowed region in the $\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$
plane after analysing NO$\nu$A and T2K complete data set with non-unitary
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 50 (88) bins is 45.88 (93.36) and it occurs at NH. For the combined
analysis, the minimum $\chi^{2}$ with 138 bins is 142.72.

In Fig. 8, $\Delta\chi^{2}$ as a function of individual non-unitary parameters
has been represented. It is obvious that each experiment rules out unitary
mixing at $1\,\sigma$ C.L., and the combined analysis does the same at
$2\,\sigma$ C.L.

Figure 8: $\Delta\chi^{2}$ as a function of individual non-unitary parameters
after analysing 2020 data.

In Fig. 9, we have represented the bi-event plots with $x$-axis ($y$-axis)
denoting the $\nu_{e}$ ($\bar{\nu}_{e}$) appearance events. These plots are
helpful in understanding the origin of the tension and their resolution with
the help of new physics. As we have already seen, in NO$\nu$A and T2K, all the
information from appearance channels can be represented by total number of
events, because of the limited statistics, the information extracted from the
shape of the energy spectrum is limited. The ellipses have been generated by
fitting the combined data from NO$\nu$A and T2K. The black ellipses,
corresponding to the standard unitary oscillation, have been generated for the
combined best-fit values of $|{\Delta_{31}}|$, $\sin^{2}2{\theta_{13}}$, and
$\sin^{2}{\theta_{23}}$. The red ellipses, corresponding to the non-unitary
oscillation, have been generated for the combined best-fit values of
$|{\Delta_{31}}|$, $\sin^{2}2{\theta_{13}}$, and $\sin^{2}{\theta_{23}}$,
$\alpha_{00}$, $|\alpha_{10}|$, $\alpha_{11}$, and $\phi_{10}$. In both cases,
only the $\delta_{\mathrm{CP}}$ has been varied in its complete range
$[-180^{\circ}:180^{\circ}]$ while keeping all other parameters fixed at their
best-fit values from the combined analysis. For the black (red) ellipse, the
square (circle) represent the best-fit event numbers for
$\delta_{\mathrm{CP}}=-170^{\circ}$ ($-134.56^{\circ}$) for NH, and
$\delta_{\mathrm{CP}}=-90^{\circ}$ ($-102.39^{\circ}$) for IH. These best-fit
values come as a compromise between T2K and NO$\nu$A . However, for the
standard case, T2K appearance events strongly prefer NH and
$\delta_{\mathrm{CP}}=-80^{\circ}$, and NO$\nu$A prefers
NH-$\delta_{\mathrm{CP}}=150^{\circ}$ or IH-$\delta_{\mathrm{CP}}=-90^{\circ}$
(Table II of Ref. Rahaman and Raut (2021)). Thus we see that for NH, neither
of the two experiments can give a good fit to the data at the combined best-
fit point. Non-unitary mixing hypothesis, brings the expected appearance event
numbers closer to the observed event numbers. For IH, however, T2K alone
cannot give a good fit to the data at the combined best-fit point for standard
oscillation, and non-unitary hypothesis brings the expected event number
closer to the observed event number. But NO$\nu$A can give a good fit to the
data at the combined best-fit point for standard oscillation, and including
non-unitary mixing does not affect it in any significant way. Hence, the
tension between the two experiments gets reduced by the non-unitary mixing.

Figure 9: Bi-event plot for T2K (NO$\nu$A) in left (right) panel. The upper
(lower) panel is for NH (IH). To generate the ellipses, $\delta_{\mathrm{CP}}$
has been varied in the range $[-180^{\circ}:180^{\circ}]$ while keeping all
other parameters at their combined best-fit values. The black (red) ellipses
represent the SM (non-unitary) case with the best-fit points indicated by
black square (red circle). The ellipses and the best-fit points have been
determined by fitting the combined data from NO$\nu$A and T2K. The black
circle with error bars represent the experimental data.

###  V.2 Lorentz invariance violation

Neutrino oscillation is the first experimental signature towards BSM physics,
as it requires neutrinos to be massive, albeit extremely light. Without loss
of any generality, SM can be considered as the low-energy effective theory
derived from a more general theory – governed by the Planck mass ($M_{P}\simeq
10^{19}\,{\rm GeV}$) – which unify the gravitational interactions along with
the weak, strong, and electromagnetic interactions. There are models which
include spontaneous Lorentz invariance violation (LIV) and CPT violations in
that more complete framework at the Planck scale Kostelecky and Samuel (1989a,
b); Kostelecky and Potting (1991, 1995, 1996). At the observable low energy,
these violations can give rise to minimal extension of SM through perturbative
terms suppressed by $M_{P}$. CPT invariance imposes that particles and anti-
particles have the same mass and lifetime. Observation of difference between
masses and lifetimes of particles and anti-particles would be a hint of CPT
violation. The present upper limit on CPT violation from the kaon system is
$|m_{K^{0}}-m_{\bar{K}_{0}}|/m_{K}<6\times 10^{-18}$ Tanabashi et al. (2018).
Since kaons are bosons and the natural mass term appearing in the Lagrangian
is mass squared term, the above constraints can be rewritten as
$|m_{K^{0}}^{2}-m_{\bar{K}_{0}}^{2}|<0.25\,{\rm eV}^{2}$. Current neutrino
oscillation data provide the bounds
$|{\Delta_{21}}-\bar{\Delta}_{21}|<5.9\times 10^{-5}\,{\rm eV}^{2}$ and
$|{\Delta_{31}}-\bar{\Delta}_{31}|<1.1\times 10^{-3}\,{\rm eV}^{2}$ Ohlsson
and Zhou (2015). The non-zero differences are manifestations of some kind of
CPT violation and this can change neutrino oscillation probability Kostelecky
and Mewes (2004a); Diaz and Kostelecky (2012); Kostelecky and Mewes (2004b);
Katori et al. (2006).

Several studies have been done about the LIV/CPT violation with neutrinos
Dighe and Ray (2008); Barenboim and Lykken (2009); Rebel and Mufson (2013); de
Gouvêa and Kelly (2017); Barenboim et al. (2018, 2019); Majhi et al. (2020);
Giunti and Laveder (2010); Datta et al. (2004); Chatterjee et al. (2014);
Singh Koranga and Khurana (2014); Diaz and Schwetz (2016); Hooper et al.
(2005); Tomar et al. (2015); Liao and Marfatia (2018); Kumar Agarwalla and
Masud (2020). Different neutrino oscillation experiments have looked for
LIV/CPT violations and put on constraints on the LIV/CPT violating parameters
Auerbach et al. (2005); Adamson et al. (2008, 2012); Aguilar-Arevalo et al.
(2018); Abe et al. (2012, 2015); Aartsen et al. (2018); Abe et al. (2017).
Constraints on all the relevant LIV/CPT violating parameters have been listed
in Ref. Kostelecky and Russell (2008). In Ref. Rahaman (2021), an effort has
been made to resolve the tension between NO$\nu$A and T2K by considering the
changes in neutrino oscillation probability due to CPT-violating LIV.

The effective Lagrangian for the Lorentz invariance violating neutrinos and
antineutrinos can be written as Kostelecky and Mewes (2004a); Kostelecky and
Mewes (2012)

$\mathcal{L}=\bar{\Psi}_{A}\left(i\gamma_{\mu}\partial_{\mu}\delta_{AB}-M_{AB}+\hat{\mathcal{Q}}_{AB}\right)\Psi_{B}+{\rm
h.c.}.$ (39)

$\Psi_{A(B)}$ is a $2N$ dimensional spinor containing $\psi_{\alpha(\beta)}$,
which is a spinor field with $\alpha(\beta)$ ranging over $N$ spinor flavours,
and their charge conjugates given by
$\psi_{\alpha(\beta)}^{C}=C\bar{\psi}_{\alpha(\beta)}^{T}$. Therefore,
$\Psi_{A(B)}$ can be expressed as

$\Psi_{A(B)}=\left(\psi_{\alpha(\beta)},\psi_{\alpha(\beta)}^{C}\right)^{T}.$
(40)

$\hat{\mathcal{Q}}$ in eq. (39) is a generic Lorentz invariance violating
operator. The first term in the right side of eq. (39) is the kinetic term,
the second term is the mass term involving the mass matrix $M$ and the third
term gives rise to the LIV effect. $\hat{\mathcal{Q}}$ is small and
perturbative in nature.

Restricting ourselves only to the renormalizable Dirac couplings in the theory
(terms only with mass dimension $\leq 4$ will be incorporated), the LIV
Lagrangian in the flavour basis can be written as Kostelecky and Mewes (2004a)

$\mathcal{L}_{\rm
LIV}=-\frac{1}{2}\left[a^{\mu}_{\alpha\beta}\bar{\psi}_{\alpha}\gamma_{\mu}\psi_{\beta}+b^{\mu}_{\alpha\beta}\bar{\psi}_{\alpha}\gamma_{5}\gamma_{\mu}\psi_{\beta}-ic_{\alpha\beta}^{\mu\nu}\bar{\psi}_{\alpha}\gamma_{\mu}\partial_{\nu}\psi_{\beta}-id_{\alpha\beta}^{\mu\nu}\bar{\psi}_{\alpha}\gamma_{\mu}\gamma_{5}\partial_{\nu}\psi_{\beta}\right],$
(41)

where $a^{\mu}_{\alpha\beta}$, $b^{\mu}_{\alpha\beta}$,
$c^{\mu\nu}_{\alpha\beta}$ and $d^{\mu\nu}_{\alpha\beta}$ are Lorentz
invariance violating parameters. Considering that only left handed neutrinos
are present in the SM, these terms can be written as

$\left(a_{L}\right)^{\mu}_{\alpha\beta}=\left(a+b\right)^{\mu}_{\alpha\beta},\left(c_{L}\right)^{\mu\nu}_{\alpha\beta}=\left(c+d\right)^{\mu\nu}_{\alpha\beta}.$
(42)

$\left(a_{L}\right)^{\mu}_{\alpha\beta}$, and
$\left(c_{L}\right)^{\mu}_{\alpha\beta}$ are constant Hermitian matrices which
can modify the standard Hamiltonian in vacuum. In Ref. Rahaman (2021), only
direction-independent isotropic terms were considered where $\mu=\nu=0$. From
now on, for simplicity, we will call $a_{\alpha\beta}^{0}$ terms as
$a_{\alpha\beta}$ and $c_{\alpha\beta}^{00}$ term as $c_{\alpha\beta}$.
$a_{\alpha\beta}$ involves CPT violating terms, and $c_{\alpha\beta}$ involves
CPT conserving Lorentz invariance violating terms. Taking into account only
these isotropic LIV terms, the neutrino Hamiltonian with LIV effect becomes:

$H=H_{\rm vac}+H_{\rm mat}+H_{\rm LIV},$ (43)

where

$H_{\rm vac}=\frac{1}{2E}U\left[\begin{array}[]{ccc}m_{1}^{2}&0&0\\\
0&m_{2}^{2}&0\\\ 0&0&m_{3}^{2}\\\ \end{array}\right]U^{\dagger};H_{\rm
mat}=\sqrt{2}G_{F}N_{e}\left[\begin{array}[]{ccc}1&0&0\\\ 0&0&0\\\ 0&0&0\\\
\end{array}\right];$ (44) $H_{\rm
LIV}=\left[\begin{array}[]{ccc}a_{ee}&a_{e\mu}&a_{e\tau}\\\
a_{e\mu}^{*}&a_{\mu\mu}&a_{\mu\tau}\\\
a_{e\tau}^{*}&a_{\mu\tau}^{*}&a_{\tau\tau}\\\
\end{array}\right]-\frac{4}{3}E\left[\begin{array}[]{ccc}c_{ee}&c_{e\mu}&c_{e\tau}\\\
c_{e\mu}^{*}&c_{\mu\mu}&c_{\mu\tau}\\\
c_{e\tau}^{*}&c_{\mu\tau}^{*}&c_{\tau\tau}\\\ \end{array}\right].$ (45)

Here $G_{F}$ is the Fermi coupling constant and $N_{e}$ is the electron
density along the neutrino path. The $-4/3$ in front of the second term arises
due to non observability of the Minkowski trace of the CPT conserving LIV term
$c_{L}$ which relates the $xx$, $yy$, and $zz$ component to the $00$ component
Kostelecky and Mewes (2004a). The effects of $a_{\alpha\beta}$ are
proportional to the baseline $L$ and the effects of $c_{\alpha\beta}$ are
proportional to the product of energy and baseline $LE$. In Ref. only the CPT
violating LIV was considered. More specifically, the authors restricted
themselves to the effects of $a_{e\mu}=|a_{e\mu}|e^{i\phi_{e\mu}}$, and
$a_{e\tau}=|a_{e_{\tau}}|e^{i\phi_{e\tau}}$, because these two terms have the
maximum effects on the $\nu_{\mu}\to\nu_{e}$ oscillation probability Kumar
Agarwalla and Masud (2020). The current constraint on these parameters from
Super-kamiokande experiment at $95\%$ C.L. is Abe et al. (2015)

$|a_{e\mu}|<2.5\times 10^{-23}\,{\rm GeV};\,|a_{e\tau}|<5\times 10^{-23}\,{\rm
GeV}$ (46)

The $\nu_{\mu}\to\nu_{e}$ oscillation probability in matter after inclusion of
LIV can be approximately written as Kumar Agarwalla and Masud (2020)

$P_{\mu e}^{\rm SM+LIV}\simeq P_{\mu e}({\rm SM})+P_{\mu e}(a_{e\mu})+P_{\mu
e}(a_{e\tau}).$ (47)

The $P_{\mu e}(\rm SM)$ term in eq. (47) has been given in eq. (16). The other
two terms can be written as Kumar Agarwalla and Masud (2020)

$\displaystyle P_{\mu e}(a_{e\beta})$ $\displaystyle=$
$\displaystyle\frac{4|a_{e\beta}|\hat{A}\hat{\Delta}\sin\theta_{13}\sin
2\theta_{23}\sin\hat{\Delta}}{\sqrt{2}G_{F}N_{e}}\left[Z_{e\beta}\sin(\delta_{\mathrm{CP}}+\phi_{e\beta})+W_{e\beta}\cos(\delta_{\mathrm{CP}}+\phi_{e\beta})\right]$
$\displaystyle=$ $\displaystyle\frac{4|a_{e\beta}|L\sin\theta_{13}\sin
2\theta_{23}\sin\hat{\Delta}}{2}\left[Z_{e\beta}\sin(\delta_{\mathrm{CP}}+\phi_{e\beta})+W_{e\beta}\cos(\delta_{\mathrm{CP}}+\phi_{e\beta})\right]$

where $\beta=\mu,\,\tau$;

$\displaystyle Z_{e\beta}$ $\displaystyle=$
$\displaystyle-\cos\theta_{23}\sin\hat{\Delta},\,{\rm if}\,\beta=\mu$ (49)
$\displaystyle=$ $\displaystyle\sin\theta_{23}\sin\hat{\Delta},\,{\rm
if}\,\beta=\tau$

and

$\displaystyle W_{e\beta}$ $\displaystyle=$
$\displaystyle\cos\theta_{23}\left(\frac{\sin^{2}{\theta_{23}}\sin\hat{\Delta}}{\cos^{2}{\theta_{12}}\hat{\Delta}}+\cos\hat{\Delta}\right),\,{\rm
if}\,\beta=\mu$ (50) $\displaystyle=$
$\displaystyle\sin\theta_{23}\left(\frac{\sin\hat{\Delta}}{\hat{\Delta}}-\cos\hat{\Delta}\right),\,{\rm
if}\,\beta=\tau.$

From, eq. (LABEL:pmue-a), it can be concluded that the LIV effects considered
in this paper are matter independent. The oscillation probability
$P(\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e})$ for antineutrino can be
calculated from eqs. (16) and (LABEL:pmue-a) by substituting $A\to-A$,
$\delta_{\mathrm{CP}}\to-\delta_{\mathrm{CP}}$, $|a_{e\beta}|\to-|a_{e\beta}|$
and $\phi_{e\beta}\to-\phi_{e\beta}$, where $\beta=\mu,\,\tau$.

Since the effects of the CPT-violating LIV terms are proportional to $L$,
NO$\nu$A will be more sensitive to CPT-violating LIV than T2K. Thus, NO$\nu$A
might be sensitive to CPT-violating LIV which T2K is insensitive to. Hence,
LIV can be a possible explanation of the disagreement between the two
experiments.

In Ref. Rahaman (2021), the individual data from NO$\nu$A and T2K, as well as
the combination of data from both the experiments have been analysed with LIV.
As before, GLoBES Huber et al. (2005, 2007) software was used to analyse the
data. The software was modified to include LIV, and the oscillation
probabilities and the event numbers were calculated in case of LIV without the
approximations required to calculate oscillation probabilities given in eqs.
(47)-(50). As we have mentioned before, among the LIV parameters, only
$|a_{e\mu}|$ ($[0:20\times 10^{-23}]$ GeV), $|a_{e\tau}|$ ($[0:20\times
10^{-23}]$ GeV), $\phi_{e\mu}$ ($-180^{\circ}:180^{\circ}$), and
$\phi_{e\tau}$ ($-180^{\circ}:180^{\circ}$) have been varied. All other LIV
parameters have been kept fixed to zero. The standarad parameter values are
same as described in the subsection III. We have represented the result on the
$\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane in Fig. 10. The minimum
$\chi^{2}$ for NO$\nu$A (T2K) is $47.71$ ($93.14$). For the combined analysis,
the minimum $\chi^{2}$ is $145.09$. Both the experiments individually prefer
NH as the best-fit hierarchy, although there are degenerate best-fit points at
IH for both of them. The experiments lose their hierarchy sensitivity when
analysed with LIV. The best-fit values of $\delta_{\mathrm{CP}}$ at NH are
close to each other. Moreover, there is a large overlap between the
$1\,\sigma$ allowed regions of the two experiments. Thus, it can be concluded
that the tension between the two experiments for NH has been reduced.
Although, there is a new mild tension between the best-fit values of
$\sin^{2}{\theta_{23}}$ from NO$\nu$A and T2K as the former prefers
${\theta_{23}}$ to be in the HO while the later prefers it to be in the LO,
both of them have a degenerate best-fit point at the other octant (LO for
NO$\nu$A, and HO for T2K) as well. The combined experiment prefers IH over NH.
However, just like the individual data, the combined data lose hierarchy (as
well as octant) sensitivity and has a degenerate best-fit point at NH.

Figure 10: Allowed region in the $\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$
plane after analysing NO$\nu$A and T2K complete data set with LIV hypothesis.
The left (right) panel represents test hierarchy to be 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
50 (88) bins is 47.71 (93.14) and it occurs at NH. For the combined analysis,
the minimum $\chi^{2}$ with 138 bins is 145.09.

In Fig. 11, we have expressed $\Delta\chi^{2}$ as a function of the LIV
parameters. It can be seen from the figure that the present T2K data
disfavours standard oscillation at $1\,\sigma$ C.L., whereas the NO$\nu$A data
do not have any preference. The combined analysis disfavours standard
oscillation at $1\,\sigma$ C.L. as well.

Figure 11: $\Delta\chi^{2}$ as a function of individual LIV parameters.

To emphasize the result, we have presented a similar bi-event plot like Fig. 9
for the non-unitary oscillation analysis in Fig. 12 for the LIV analysis. It
is obvious that inclusion of LIV in the theory brings the expected $\nu_{e}$,
and $\bar{\nu}_{e}$ event numbers, of both the experiments, at the combined
best-fit point at NH closer to the observed event numbers, and thus reduces
the tension between the two experiments. For IH, inclusion of LIV brings the
observed $\nu_{e}$, and $\bar{\nu}_{e}$ event numbers of T2K at the combined
best-fit point closer to the observed event numbers. Although, LIV takes the
observed $\nu_{e}$, and $\bar{\nu}_{e}$ event numbers of NO$\nu$A at the
combined IH best-fit point farther away from the observed event numbers, the
change is minuscule.

Figure 12: Bi-event plot for T2K (NO$\nu$A) in left (right) panel. The upper
(lower) panel is for NH (IH). To generate the ellipses, $\delta_{\mathrm{CP}}$
has been varied in the range $[-180^{\circ}:180^{\circ}]$ while keeping all
other parameters fixed at the combined best-fit values. The black (red)
ellipses represent the SM (LIV) case with the best-fit points indicated by
black square (red circle). The ellipses and the best-fit points have been
determined by fitting the combined data from NO$\nu$A and T2K. The black
circle with error bars represent the experimental data.

###  V.3 Non-standard interaction (NSI)

Non-standard interactions can arise as a low-energy manifestation of new heavy
states of a more complete model at high energy Farzan and Tortola (2018);
Biggio et al. (2009); Ohlsson (2013); Miranda and Nunokawa (2015); Pro (2019)
or it can arise due to light mediators Farzan (2015); Farzan and Shoemaker
(2016). NSI can modify the neutrino and antineutrino flavour conversion in
matter Wolfenstein (1978); Mikheev and Smirnov (1985, 1986). The effect of NSI
on the present and future long-baseline accelerator neutrino experiments have
been discussed in details in literature Friedland and Shoemaker (2012); Rahman
et al. (2015); Coelho et al. (2012); Masud et al. (2019); Deepthi et al.
(2018); Blennow et al. (2016). Ref. Chatterjee and Palazzo (2021); Denton et
al. (2021) have tried to invoke NSI to resolve the tension between NO$\nu$A
and T2K data.

Neutral current NSI during neutrino propagation can be represented by a
dimension 6 operator Wolfenstein (1978):

$\mathcal{L}_{\rm NC-
NSI}=-2\sqrt{2}G_{F}\epsilon_{\alpha\beta}^{fC}\left(\bar{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\beta}\right)\left(\bar{f}\gamma_{\mu}P_{C}f\right),$
(51)

where $\alpha,\,\beta=e,\,\mu,\,\tau$ denote the neutrino flavour,
$f=e,\,\mu,\,\tau$ denotes the fermions inside matter, $P$ is the projection
operator with the superscript $C$ referring to the $L$ or $R$ chirality of the
$ff$ current, and $\epsilon_{\alpha\beta}^{fC}$ denotes the strength of the
NSI. From the hermiticity of the interaction,

$\epsilon_{\beta\alpha}^{fC}=\left(\epsilon_{\alpha\beta}^{fC}\right)^{*}\,.$
(52)

For neutrino propagation through earth matter, the relevant expression is

$\epsilon_{\alpha\beta}\equiv\sum_{f=e,u,d}\epsilon_{\alpha\beta}^{f}\frac{N_{f}}{N_{e}}\equiv\sum_{f=e,u,d}\left(\epsilon_{\alpha\beta}^{fL}+\epsilon_{\alpha\beta}^{fR}\right)\frac{N_{f}}{N_{e}},$
(53)

where $N_{f}$ is the density of $f$ fermion. If we consider earth matter to be
neutral and isoscalar, then $N_{n}\simeq N_{p}=N_{e}$. Thus,

$\epsilon_{\alpha\beta}\simeq\epsilon_{\alpha\beta}^{e}+3\epsilon_{\alpha\beta}^{u}+3\epsilon_{\alpha\beta}^{d}\,.$
(54)

The effective Hamiltonian for neutrino propagation in matter in presence of
NSI can be written in the flavour basis as

$H=H_{\rm vac}+H_{\rm mat}+H_{\rm NSI},$ (55)

where

$H_{\rm vac}=\frac{1}{2E}U\left[\begin{array}[]{ccc}m_{1}^{2}&0&0\\\
0&m_{2}^{2}&0\\\ 0&0&m_{3}^{2}\\\ \end{array}\right]U^{\dagger};H_{\rm
mat}=\sqrt{2}G_{F}N_{e}\left[\begin{array}[]{ccc}1&0&0\\\ 0&0&0\\\ 0&0&0\\\
\end{array}\right];$ (56) $H_{\rm
NSI}=\sqrt{2}G_{F}N_{e}\left[\begin{array}[]{ccc}\epsilon_{ee}&\epsilon_{e\mu}&\epsilon_{e\tau}\\\
\epsilon_{e\mu}^{*}&\epsilon_{\mu\mu}&\epsilon_{\mu\tau}\\\
\epsilon_{e\tau}^{*}&\epsilon_{\mu\tau}^{*}&\epsilon_{\tau\tau}\\\
\end{array}\right].$ (57)

It is important to note that the NC NSI Hamiltonian during propagation
presented in eqs. (55)-(57) are analogous to the CPT violating LIV Hamiltonian
given in eqs. (43)-(45). A relationship between CPT-violating LIV and NSI can
be found by the following relation Diaz (2015):

$\epsilon_{\alpha\beta}=\frac{a_{\alpha\beta}}{\sqrt{2}G_{F}N_{e}}.$ (58)

The $\nu_{\mu}\to\nu_{e}$ oscillation probability with matter effect in
presence of NSI during propagation can be written in the similar way as in eq.
(47) Kikuchi et al. (2009); Agarwalla et al. (2016); Masud et al. (2019)

$P_{\mu e}^{\rm SM+NSI}\simeq P_{\mu e}(\rm SM)+P_{\mu
e}(\epsilon_{e\mu})+P_{\mu e}(\epsilon_{e\tau}).$ (59)

Just like the LIV case, we can write the second and third terms in eq. (59)
as,

$\displaystyle P_{\mu e}(\epsilon_{e\beta})$ $\displaystyle=$ $\displaystyle
4|\epsilon_{e\beta}|\hat{A}\hat{\Delta}\sin\theta_{13}\sin
2\theta_{23}\sin\hat{\Delta}\left[Z_{e\beta}\sin(\delta_{\mathrm{CP}}+\phi_{e\beta})+W_{e\beta}\cos(\delta_{\mathrm{CP}}+\phi_{e\beta})\right]$
(60)

where $\beta=\mu,\,\tau$;

$\displaystyle Z_{e\beta}$ $\displaystyle=$
$\displaystyle-\cos\theta_{23}\sin\hat{\Delta},\,{\rm if}\,\beta=\mu$ (61)
$\displaystyle=$ $\displaystyle\sin\theta_{23}\sin\hat{\Delta},\,{\rm
if}\,\beta=\tau\,,$

and

$\displaystyle W_{e\beta}$ $\displaystyle=$
$\displaystyle\cos\theta_{23}\left(\frac{\sin^{2}{\theta_{23}}\sin\hat{\Delta}}{\cos^{2}{\theta_{12}}\hat{\Delta}}+\cos\hat{\Delta}\right),\,{\rm
if}\,\beta=\mu$ (62) $\displaystyle=$
$\displaystyle\sin\theta_{23}\left(\frac{\sin\hat{\Delta}}{\hat{\Delta}}-\cos\hat{\Delta}\right),\,{\rm
if}\,\beta=\tau.$

Because of the $\hat{A}$ term in eq. (60), the oscillation probability after
inclusion of NSI is dependent on matter effect unlike LIV. The first
oscillation maximum of NO$\nu$A (T2K) peaks at $1.4$ GeV ($0.6$ GeV).
Therefore, the matter effect is almost 3 times larger at NO$\nu$A
($\hat{A}\simeq 0.14$) than at T2K ($\hat{A}\simeq 0.05$). Hence, NO$\nu$A can
observe NSI effects which can be remain unseen by T2K. Therefore, NSI can be a
possible explanation behind the tension between the two experiments.

Figure 13: Allowed region in the $\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$
plane after analysing NO$\nu$A and T2K complete data set with NSI hypothesis.
Only the effect of $\epsilon_{e\mu}=|\epsilon_{e\mu}|e^{i\phi_{e\mu}^{\rm
NSI}}$ has been considered. The left (right) panel represents test hierarchy
to be 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 50 (88) bins is 49.08 (93.64) and it occurs
at NH. For the combined analysis, the minimum $\chi^{2}$ with 138 bins is
146.26.

In Ref. Chatterjee and Palazzo (2021); Denton et al. (2021) an effort to
resolve the tension with the NSI has been made. To do so, the authors first
analysed the combined data of NO$\nu$A and T2K. Finding the best-fit values of
the NSI parameters after the combined analysis, they analysed the individual
data from NO$\nu$A and T2K with the NSI parameter values fixed at the combined
best-fit values, and showed that the two experiments agree on their results on
$\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane. In this review article, we
have analysed the individual NO$\nu$A and T2K data as well as their combined
data independent of each other with NSI. To do so, at first we considered only
$\epsilon_{e\mu}=|\epsilon_{e\mu}|e^{i\phi_{e\mu}^{\rm NSI}}$, and fixed all
other NSI parameters to be $0$. We have presented the result on the
$\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane in Fig. 13. The minimum
$\chi^{2}$ for NO$\nu$A (T2K) with 50 (88) bins is 49.08 (93.64) and it occurs
at NH (IH). For the combined analysis, the minimum $\chi^{2}$ with 138 bins is
146.26 and it is at NH. Both the experiments lose hierarchy sensitivity after
analysing with NSI, and as a result NO$\nu$A (T2K) has a degenerate best-fit
point at IH (NH). The best-fit points of the two experiments are close to each
other for both the hierarchies. However, they exclude each other’s best-fit
point at the $1\,\sigma$ C.L. for both the hierarchies. There is a significant
overlap between the $1\,\sigma$ allowed regions of the two experiments. The
experiments lose their octant sensitivity as well after the inclusion of NSI.
At the best-fit points, the $\delta_{\mathrm{CP}}$ values are close to
$-90^{\circ}$ for both the experiments and both the hierarchies. The combined
analysis prefers NH, ${\theta_{23}}$ in HO and
$\delta_{\mathrm{CP}}\sim-90^{\circ}$ as the best-fit point. However, there is
a nearly degenerate best-fit point at IH, ${\theta_{23}}$ in HO and
$\delta_{\mathrm{CP}}\sim-90^{\circ}$.

Just like the non-unitary mixing and LIV, we have emphasized our argument
about resolution of the tension with the inclusion of NSI due to
$\epsilon_{e\mu}$ through bi-event plots presented in Fig. 14. It is clear
that for both the experiments, inclusion of NSI due to $\epsilon_{e\mu}$
brings the expected $\nu_{e}$ and $\bar{\nu}_{e}$ event numbers at the
combined best-fit points closer to the observed event numbers for NH, and thus
resolves the tension between the two experiments for NH. For IH, the change in
expected event numbers for T2K at the combined best-fit point after inclusion
of NSI due to $\epsilon_{e\mu}$ is negligible. For NO$\nu$A, expected
$\bar{\nu}_{e}$ event numbers at the combined IH best-fit point comes closer
to the observed event number after inclusion of NSI due to $\epsilon_{e\mu}$,
whereas change in the expected $\nu_{e}$ appearance event number due to the
same is negligible.

Figure 14: Bi-event plot for T2K (NO$\nu$A) in the left (right) panel. The
upper (lower) panel is for NH (IH). To generate the ellipses,
$\delta_{\mathrm{CP}}$ has been varied in the range
$[-180^{\circ}:180^{\circ}]$ while keeping all other parameters fixed. Among
the NSI parameters only the effect of $\epsilon_{e\mu}$ has been considered.
The black (red) ellipses represent the SM (NSI) case with the best-fit points
indicated by black square (red circle). The ellipses and the best-fit points
have been determined by fitting the combined data from NO$\nu$A and T2K. The
black circle with error bars represent the experimental data.

In the next step, we have considered the effects of
$\epsilon_{e\tau}=|\epsilon_{e\tau}|e^{i\phi_{e\tau}^{\rm NSI}}$. All other
NSI parameters have been kept fixed at $0$. The result has been displayed on
the $\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane in Fig. 15. The minimum
$\chi^{2}$ for NO$\nu$A (T2K) with 50 (88) bins is 48.59 (93.73) and it occurs
at IH. For the combined analysis, the minimum $\chi^{2}$ with 138 bins is
146.38 and it is at NH. Although NO$\nu$A and T2K both individually prefer IH
and ${\theta_{23}}$ in LO as their best-fit point, both of them lose hierarchy
and octant sensitivity after consideration of NSI parameter
$\epsilon_{e\tau}$, and therefore both of them have a degenerate best-fit
point at NH as well as ${\theta_{23}}$ in HO. The best-fit points of the two
experiments are close to each other. However, both of them exclude each
other’s best-fit point at $1\,\sigma$ C.L. for both the hierarchies.
Nonetheless, there is a significant overlap between the $1\,\sigma$ allowed
region of the two experiments after including the effect of the NSI parameter
$\epsilon_{e\tau}$. At the best-fit points the $\delta_{\mathrm{CP}}$ values
are close to $-90^{\circ}$ for both the experiments and both the hierarchies.
The combined analysis prefers NH, ${\theta_{23}}$ in HO and
$\delta_{\mathrm{CP}}\sim-90^{\circ}$ as the best-fit point. However, there is
a nearly degenerate best-fit point at IH, ${\theta_{23}}$ in HO and
$\delta_{\mathrm{CP}}\sim-90^{\circ}$.

Figure 15: Allowed region in the $\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$
plane after analysing NO$\nu$A and T2K complete data set with the NSI
hypothesis. Only the effect of
$\epsilon_{e\tau}=|\epsilon_{e\tau}|e^{i\phi_{e\tau}^{\rm NSI}}$ has been
considered. The left (right) panel represents test hierarchy to be 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 50 (88) bins is 48.59 (93.73) and it occurs at IH. For
the combined analysis, the minimum $\chi^{2}$ with 138 bins is 146.38 and it
is at NH.

As before, we have emphasized our argument about resolution of the tension
between NO$\nu$A and T2K with the inclusion of NSI due to $\epsilon_{e\tau}$
through bi-event plots in Fig. 16. At the best-fit point of the combined
analysis, the value of $|\epsilon_{e\tau}|$ is $0.73$. However, there is a
nearly degenerate best-fit point at $|\epsilon_{e\tau}|=0.19$ with
$\Delta\chi^{2}=0.13$. Because of the stronger constraint against
$\epsilon_{e\tau}$ to be large from IceCube data Ehrhardt (2019), we consider
the combined best-fit point at $|\epsilon_{e\tau}|=0.19$. It is clear that for
NH, inclusion of NSI due to $\epsilon_{e\tau}$ brings the expected $\nu_{e}$
and $\bar{\nu}_{e}$ appearance events for both the experiments at their
combined best-fit point closer to the observed event numbers. For IH, the
change in expected $\nu_{e}$ and $\bar{\nu}_{e}$ appearance events for T2K at
the combined best-fit point is negligible. For NO$\nu$A, after the inclusion
of NSI due to $\epsilon_{e\tau}$ the expected $\bar{\nu}_{e}$ event numbers at
the combined IH best-fit point comes closer to the observed event number,
whereas the change in expected $\nu_{e}$ appearance event number due to the
same is quite small.

Figure 16: Bi-event plot for T2K (NO$\nu$A) in left (right) panel. The upper
(lower) panel is for NH (IH). To generate the ellipses, $\delta_{\mathrm{CP}}$
has been varied in the range $[-180^{\circ}:180^{\circ}]$ while keeping all
other parameters fixed. Among the NSI parameters, only the effect of
$\epsilon_{e\tau}$ has been considered. The black (red) ellipses represent the
SM (NSI) case with the best-fit points indicated by black square (red circle).
The ellipses and the best-fit points have been determined by fitting the
combined data from NO$\nu$A and T2K. The black circle with error bars
represent the experimental data.

In Fig. 17, we have presented the precision plots for the NSI parameters. It
can be concluded that when we consider the effect of $\epsilon_{e\mu}$, the
present NO$\nu$A data cannot make any preference between SM and NSI. However,
both T2K and the combined data rule out SM at $1\,\sigma$ C.L. When the effect
of $\epsilon_{e\tau}$ is considered, all three cases – NO$\nu$A, T2K, and
their combined data – rule out SM at $1\,\sigma$ C.L.

Figure 17: $\Delta\chi^{2}$ as a function of individual NSI parameters.

The best-fit values for various non-standard parameters discussed in this
section have been listed in table 7. The $90\%$ C.L. limit for $1$ degree of
freedom (d.o.f.) of these parameters have been mentioned as well in the
parenthesis. When the $90\%$ limit falls beyond the studied range of a
parameter, NA has been mentioned instead of a number.

Parameters |  NO$\nu$A | T2K |  NO$\nu$A+T2K
---|---|---|---
| NH | IH | NH | IH | NH | IH
$\alpha_{00}$ |  $0.84$ (NA) |  $0.72$ (NA) |  $0.80$ (NA) |  $0.80$ (NA) |  $0.70$ (NA) |  $0.76$ (NA)
$|\alpha_{10}|$ |  $0.080$ (NA) |  $0.0.12^{+0.06}_{-0.12}$ (NA) |  $0.082^{+0.108}_{-0.082}$ ($<0.190$) |  $0.080^{+0.110}_{-0.060}$ (NA) |  $0.125^{+0.025}_{-0.085}$ ($<0.170$) |  $0.110^{+0.040}_{-0.070}$ ($<0.155$)
$\alpha_{11}$ |  $0.97^{+0.02}_{-0.03}$ ($>0.92$) |  $0.96^{+0.03}_{-0.03}$ ($>0.92$) |  $0.98$ ($>0.95$) |  $0.98^{+0.02}_{-0.03}$ ($>0.95$) |  $0.98^{+0.01}_{-0.02}$ ($>0.95$) |  $0.98^{+0.01}_{-0.02}$ ($>0.95$)
$\phi_{10}/^{\circ}$ | $-(125.68^{+54.32}_{-305.68})$ | $76.15^{+103.85}_{-86.4}$ | $54.77^{+97.10}_{-60.54}$ | $112.69^{+42.38}_{-79.93}$ | $120.41^{+59.57}_{-300.33}$ | $4.31^{+162.71}_{-181.51}$
$\frac{|a_{e\mu}|}{10^{-23}{\rm GeV}}$ |  $4.81$($<8.19$) |  $2.22$ ($<7.78$) |  $4.60$ ($<15.25$) |  $6.17_{-6.02}^{+4.86}$ ($<14.92$) |  $1.86^{+2.57}_{-1.86}$ ($<4.80$) |  $1.52$ ($<3.80$)
$\frac{|a_{e\tau}|}{10^{-23}{\rm GeV}}$ |  $2.52$ ($<3.18$) |  $5.33^{+9.20}_{-5.33}$ ($<15.71$) |  $11.14$ (NA) |  $8.06$ (NA) |  $0.57$ ($<6.70$) |  $4.16$ ($<9.50$)
$\phi_{e\mu}/^{\circ}$ | $-114.52$ | $141.18$ | $64.29$ | $-77.34$ | $-115.72$ | $110.08$
$\phi_{e\tau}/^{\circ}$ | $-145.02$ | $25.04$ | $-153.24$ | $158.07$ | $84.36$ | $-89.27$
$|\epsilon_{e\mu}|$ |  $0.23$ ($<0.40$) |  $0.06^{+0.10}_{-0.04}$ ($<0.42$) |  $0.24_{-0.21}^{+0.28}$ (NA) |  $0.46$ (NA) |  $0.19_{-0.10}^{+0.11}$ ($<0.36$) |  $0.09_{-0.05}^{+0.05}$ ($<0.24$)
$|\epsilon_{e\tau}|$ |  $0.16$ (NA) |  $0.95$ (NA) |  $0.15$ (NA) |  $0.33$ (NA) |  $0.19_{-0.09}^{+0.72}$ (NA) |  $0.12_{-0.08}^{+0.13}$ ($<0.91$)
$\phi_{e\mu}^{\rm NSI}/^{\circ}$ | $-160$ | $60$ | $-80$ | $50$ | $-150$ | $40$
$\phi_{e\tau}^{\rm NSI}/^{\circ}$ | $30$ | $90$ | $120$ | $-150$ | $-30$ | $120$
Table 7: Best-fit values of several BSM parameters discussed in section V
along with the $1\,\sigma$ error bar, where available. The $90\%$ C.L. limits
for 1 d.o.f. have been mentioned in the parenthesis. When the $90\%$ limit
falls beyond the studied range of a parameter, NA has been mentioned.

##  VI Summary and discussion

A tension between the best-fit points of T2K and NO$\nu$A existed from the
very beginning, which became only stronger with time. This tension arises
mostly from the $\nu_{e}$ appearance data of the two experiments. T2K observes
a large excess in the $\nu_{e}$ appearance events compared to the expected
event number at the reference point of vacuum oscillation,
${\theta_{23}}=\pi/4$, and $\delta_{\mathrm{CP}}=0$ (referred to as $000$).
This large excess dictates that $\delta_{\mathrm{CP}}$ be anchored around
$-90^{\circ}$ and that ${\theta_{23}}$ be in HO, for both the hierarchies
($+++$, and $-++$ with former being the best-fit point). The appearance events
observed by NO$\nu$A show a very different pattern. They are moderately larger
than the expectation from the reference point in the $\nu_{e}$ channel and are
consistent with it in the $\bar{\nu}_{e}$ channel. These two facts, when
combined together, lead to two possible degenerate solutions for NO$\nu$A : A.
NH - ${\theta_{23}}$ in HO - $\delta_{\mathrm{CP}}$ in UHP ($++-$), and B. IH
- ${\theta_{23}}$ in HO - $\delta_{\mathrm{CP}}$ in LHP ($-++$). A fit of the
combined T2K + NO$\nu$A data to standard three flavour oscillation framework,
has the best-fit point as IH - ${\theta_{23}}$ in HO - $\delta_{\mathrm{CP}}$
in LHP which is reasonably close to the IH best-fit points of T2K and
NO$\nu$A. If NH is assumed to be the true hierarchy, there is almost no
allowed region within $1\,\sigma$, even though the best-fit point of each
experiment picks NH. This is the essential tension between the two
experiments.

Several studies have been done to resolve this tension with BSM physics. Three
different BSM scenarios have been considered in the literarture: 1. non-
unitary mixing, 2. Lorentz invariance violation, and 3. non-standard
interaction during neutrino propagation. All these three scenarios bring the
expected event numbers at the combined best-fit point at NH closer to the
observed $\nu_{e}$, and $\bar{\nu}_{e}$ event numbers of both the experiments,
and thus reduce the tension between them.

T2K and NO$\nu$A individually prefer non-unitary mixing over unitary mixing at
$1\,\sigma$ C.L. Combined data from both of them prefer non-unitary mixing at
$2\,\sigma$ C.L. Both the experiments lose hierarchy and octant sensitivity
when analysed with non-unitary mixing. There is a large overlap between the
$1\,\sigma$ allowed regions on the
$\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane of the two experiments.

In the case of LIV, T2K data prefers LIV over standard 3-flavour oscillation
at $1\,\sigma$. NO$\nu$A data cannot make any preference between the two
hypothesis at $1\,\sigma$. The combined analysis rules out standard
oscillation at $1\,\sigma$ C.L. Just like, non-unitary mixing, both the
experiments lose hierarchy, and octant sensitivity in the case of LIV too. In
this case also, there is a large overlap between the $1\,\sigma$ allowed
region on the $\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane of the two
experiments.

In the case of NSI, we considered the effects of
$\epsilon_{e\mu}=|\epsilon_{e\mu}|e^{i\phi_{e\mu}^{\rm NSI}}$, and
$\epsilon_{e\tau}=|\epsilon_{e\tau}|e^{i\phi_{e\tau}^{\rm NSI}}$ one at a
time. In case of $\epsilon_{e\mu}$, T2K data rule out standard oscillation at
$1\,\sigma$ C.L., whereas NO$\nu$A data cannot make a preference between the
two hypothesis. The combined data from NO$\nu$A and T2K rule out standard
oscillation at more than $1.5\,\sigma$ C.L. In the case of $\epsilon_{e\tau}$,
data from both NO$\nu$A and T2K rule out standard oscillation at $1\,\sigma$,
whereas their combined data rule out standard oscillation at $1.5\,\sigma$
C.L. As before, in the case of NSI also, each of the two experiments loses
their hierarchy and octant sensitivity. A large overlap between the
$1\,\sigma$ allowed regions on the
$\sin^{2}{\theta_{23}}-\delta_{\mathrm{CP}}$ plane of the two experiments
exists.

T2K and NO$\nu$A continue to take data. The additional data may either sharpen
or reduce the tension. If the tension becomes sharper, then we need to explore
which new physics scenario can best relieve this tension. We also need to test
the predictions of the preferred new physics scenario at future neutrino
oscillation experiments, such as T2HK Abe et al. (2018c) and DUNE Abi et al.
(2018a, b, c).

###### Acknowledgements.

U.R. and S.R. were supported by a grant from the Research Council of the
University of Johannesburg.

## References

  * Cleveland et al. (1998) B. T. Cleveland, T. Daily, R. Davis, Jr., J. R. Distel, K. Lande, C. K. Lee, P. S. Wildenhain, and J. Ullman, Astrophys. J. 496, 505 (1998).
  * Hirata et al. (1991) K. S. Hirata et al. (Kamiokande-II), Phys. Rev. D 44, 2241 (1991), [Erratum: Phys.Rev.D 45, 2170 (1992)].
  * Hampel et al. (1999) W. Hampel et al. (GALLEX), Phys. Lett. B447, 127 (1999).
  * Abdurashitov et al. (2009) J. N. Abdurashitov et al. (SAGE), Phys. Rev. C80, 015807 (2009), eprint 0901.2200.
  * Altmann et al. (2005) M. Altmann et al. (GNO), Phys. Lett. B 616, 174 (2005), eprint hep-ex/0504037.
  * Abe et al. (2016a) K. Abe et al. (Super-Kamiokande), Phys. Rev. D 94, 052010 (2016a), eprint 1606.07538.
  * Aharmim et al. (2013) B. Aharmim et al. (SNO), Phys. Rev. C 87, 015502 (2013), eprint 1107.2901.
  * Agostini et al. (2019) M. Agostini et al. (Borexino), Phys. Rev. D 100, 082004 (2019), eprint 1707.09279.
  * Casper et al. (1991) D. Casper et al., Phys. Rev. Lett. 66, 2561 (1991).
  * Becker-Szendy et al. (1992) R. Becker-Szendy et al., Phys. Rev. D 46, 3720 (1992).
  * Hirata et al. (1992) K. Hirata et al. (Kamiokande-II), Phys. Lett. B 280, 146 (1992).
  * Fukuda et al. (1994) Y. Fukuda et al. (Kamiokande), Phys. Lett. B 335, 237 (1994).
  * Fukuda et al. (1998) Y. Fukuda et al. (Super-Kamiokande), Phys. Rev. Lett. 81, 1562 (1998), eprint hep-ex/9807003.
  * Steinberger (1991) J. Steinberger, Phys. Rept. 203, 345 (1991).
  * Maki et al. (1962) Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).
  * Bilenky and Pontecorvo (1978) S. M. Bilenky and B. Pontecorvo, Phys. Rept. 41, 225 (1978).
  * Wolfenstein (1978) L. Wolfenstein, Phys. Rev. D17, 2369 (1978).
  * Wolfenstein (1979) L. Wolfenstein, Phys. Rev. D 20, 2634 (1979).
  * Abe et al. (2018a) K. Abe et al. (Super-Kamiokande), Phys. Rev. D 97, 072001 (2018a), eprint 1710.09126.
  * Bahcall and Pinsonneault (2004) J. N. Bahcall and M. H. Pinsonneault, Phys. Rev. Lett. 92, 121301 (2004), eprint astro-ph/0402114.
  * Kobayashi and Maskawa (1973) M. Kobayashi and T. Maskawa, Prog. Theor. Phys. 49, 652 (1973).
  * Apollonio et al. (1998) M. Apollonio et al. (CHOOZ), Phys. Lett. B420, 397 (1998), eprint hep-ex/9711002.
  * Narayan et al. (1998) M. Narayan, G. Rajasekaran, and S. U. Sankar, Phys.Rev. D58, 031301 (1998), eprint hep-ph/9712409.
  * Araki et al. (2005) T. Araki et al. (KamLAND), Phys. Rev. Lett. 94, 081801 (2005), eprint hep-ex/0406035.
  * Adey et al. (2018) D. Adey et al. (Daya Bay), Phys. Rev. Lett. 121, 241805 (2018), eprint 1809.02261.
  * Bak et al. (2018) G. Bak et al. (RENO), Phys. Rev. Lett. 121, 201801 (2018), eprint 1806.00248.
  * Abe et al. (2016b) Y. Abe et al. (Double Chooz), JHEP 01, 163 (2016b), eprint 1510.08937.
  * Adamson et al. (2014) P. Adamson et al. (MINOS), Phys. Rev. Lett. 112, 191801 (2014), eprint 1403.0867.
  * Esteban et al. (2020) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09, 178 (2020), eprint 2007.14792.
  * de Salas et al. (2021) P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola, and J. W. F. Valle, JHEP 02, 071 (2021), eprint 2006.11237.
  * Jarlskog (1985) C. Jarlskog, Phys. Rev. Lett. 55, 1039 (1985).
  * Cervera et al. (2000) A. Cervera, A. Donini, M. B. Gavela, J. J. Gomez Cadenas, P. Hernandez, O. Mena, and S. Rigolin, Nucl. Phys. B579, 17 (2000), [Erratum: Nucl. Phys.B593,731(2001)], eprint hep-ph/0002108.
  * Lipari (2000) P. Lipari, Phys. Rev. D 61, 113004 (2000), eprint hep-ph/9903481.
  * Narayan and Sankar (2000) M. Narayan and S. U. Sankar, Phys. Rev. D 61, 013003 (2000), eprint hep-ph/9904302.
  * Itow et al. (2001) Y. Itow et al. (T2K), pp. 239–248 (2001), eprint hep-ex/0106019.
  * Ayres et al. (2007) D. Ayres et al. (NOvA) (2007).
  * McDonald, Kirk T (2001) McDonald, Kirk T (2001), eprint hep-ex/0111033.
  * Nunokawa et al. (2005) H. Nunokawa, S. J. Parke, and R. Zukanovich Funchal, Phys.Rev. D72, 013009 (2005), eprint hep-ph/0503283.
  * Barger et al. (2002) V. Barger, D. Marfatia, and K. Whisnant, Phys. Rev. D 65, 073023 (2002), eprint hep-ph/0112119.
  * Minakata and Sugiyama (2004) H. Minakata and H. Sugiyama, Phys. Lett. B580, 216 (2004), eprint hep-ph/0309323.
  * Mena and Parke (2004) O. Mena and S. J. Parke, Phys.Rev. D70, 093011 (2004), eprint hep-ph/0408070.
  * Prakash et al. (2012) S. Prakash, S. K. Raut, and S. U. Sankar, Phys.Rev. D86, 033012 (2012), eprint 1201.6485.
  * Agarwalla et al. (2012) S. K. Agarwalla, S. Prakash, S. K. Raut, and S. U. Sankar, JHEP 1212, 075 (2012), eprint 1208.3644.
  * Nichol (2012a) R. Nichol (MINOS) (2012a), talk given at the Neutrino 2012 Conference, June 3-9, 2012, Kyoto, Japan, http://neu2012.kek.jp/.
  * Forero et al. (2012) D. Forero, M. Tortola, and J. Valle, Phys.Rev. D86, 073012 (2012), eprint 1205.4018.
  * Fogli et al. (2012) G. Fogli, E. Lisi, A. Marrone, D. Montanino, A. Palazzo, et al., Phys.Rev. D86, 013012 (2012), eprint 1205.5254.
  * Gonzalez-Garcia et al. (2012) M. Gonzalez-Garcia, M. Maltoni, J. Salvado, and T. Schwetz, JHEP 1212, 123 (2012), eprint 1209.3023.
  * Prakash et al. (2014) S. Prakash, U. Rahaman, and S. U. Sankar, JHEP 07, 070 (2014), eprint 1306.4125.
  * Agarwalla et al. (2013) S. K. Agarwalla, S. Prakash, and S. U. Sankar (2013), eprint 1301.2574.
  * Huber et al. (2005) P. Huber, M. Lindner, and W. Winter, Comput.Phys.Commun. 167, 195 (2005), eprint hep-ph/0407333.
  * Huber et al. (2007) P. Huber, J. Kopp, M. Lindner, M. Rolinec, and W. Winter, Comput.Phys.Commun. 177, 432 (2007), eprint hep-ph/0701187.
  * Himmel (2020) A. Himmel (2020), talk given at the Neutrino 2020 meeting on July, 2nd, 2020,https://indico.fnal.gov/event/43209/contributions/187840/attachments/130740/159597/NOvA-Oscilations-NEUTRINO2020.pdf.
  * Acero et al. (2021) M. A. Acero et al. (NOvA) (2021), eprint 2108.08219.
  * Dunne (2020) P. Dunne (2020), talk given at the Neutrino 2020 meeting on July, 2nd, 2020,https://indico.fnal.gov/event/43209/contributions/187830/attachments/129636/159603/T2K_Neutrino2020.pdf.
  * Abe et al. (2021) K. Abe et al. (T2K), Phys. Rev. D 103, 112008 (2021), eprint 2101.03779.
  * Dohnal (2021) T. Dohnal (Daya Bay), PoS ICRC2021, 1175 (2021).
  * Esteban et al. (2019) I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni, and T. Schwetz, JHEP 01, 106 (2019), eprint 1811.05487.
  * Nichol (2012b) R. Nichol (MINOS) (2012b), talk given at the Neutrino 2012 Conference, June 3-9, 2012, Kyoto, Japan, http://neu2012.kek.jp/.
  * Acero et al. (2018) M. A. Acero et al. (NOvA), Phys. Rev. D 98, 032012 (2018), eprint 1806.00096.
  * Acero et al. (2019) M. A. Acero et al. (NOvA), Phys. Rev. Lett. 123, 151803 (2019), eprint 1906.04907.
  * Adamson et al. (2017) P. Adamson et al. (NOvA), Phys. Rev. Lett. 118, 231801 (2017), eprint 1703.03328.
  * Bharti et al. (2018) S. Bharti, S. Prakash, U. Rahaman, and S. Uma Sankar, JHEP 09, 036 (2018), eprint 1805.10182.
  * Lindner et al. (2018) M. Lindner, W. Rodejohann, and X.-J. Xu, Phys. Rev. D 97, 075024 (2018), eprint 1709.10252.
  * Sanchez (2018) M. Sanchez (2018), talk given at Neutrino 2018, https://doi.org/10.5281/zenodo.1286758.
  * Nizam et al. (2019) M. Nizam, S. Bharti, S. Prakash, U. Rahaman, and S. Uma Sankar, Mod. Phys. Lett. A 35, 06 (2019), eprint 1811.01210.
  * Rahaman and Raut (2021) U. Rahaman and S. Raut (2021), eprint 2112.13186.
  * Abe et al. (2014a) K. Abe et al. (T2K Collaboration), Phys.Rev.Lett. 112, 061802 (2014a), eprint 1311.4750.
  * Abe et al. (2014b) K. Abe et al. (T2K Collaboration), Phys.Rev.Lett. 112, 181801 (2014b), eprint 1403.1532.
  * Bharti et al. (2016) S. Bharti, S. Prakash, U. Rahaman, and S. U. Sankar (2016), eprint 1602.03513.
  * Abe et al. (2018b) K. Abe et al. (T2K), Phys. Rev. Lett. 121, 171802 (2018b), eprint 1807.07891.
  * Abe et al. (2020) K. Abe et al. (T2K), Nature 580, 339 (2020), [Erratum: Nature 583, E16 (2020)], eprint 1910.03887.
  * Berns (2021) L. Berns (T2K), in _55th Rencontres de Moriond on Electroweak Interactions and Unified Theories_ (2021), eprint 2105.06732.
  * Kelly et al. (2021) K. J. Kelly, P. A. N. Machado, S. J. Parke, Y. F. Perez-Gonzalez, and R. Z. Funchal, Phys. Rev. D 103, 013004 (2021), eprint 2007.08526.
  * Argüelles et al. (2020) C. A. Argüelles et al., Rept. Prog. Phys. 83, 124201 (2020), eprint 1907.08311.
  * Miranda et al. (2021) L. S. Miranda, P. Pasquini, U. Rahaman, and S. Razzaque, Eur. Phys. J. C 81, 444 (2021), eprint 1911.09398.
  * Rahaman (2021) U. Rahaman, Eur. Phys. J. C 81, 792 (2021), eprint 2103.04576.
  * Chatterjee and Palazzo (2021) S. S. Chatterjee and A. Palazzo, Phys. Rev. Lett. 126, 051802 (2021), eprint 2008.04161.
  * Denton et al. (2021) P. B. Denton, J. Gehrlein, and R. Pestes, Phys. Rev. Lett. 126, 051801 (2021), eprint 2008.01110.
  * Aguilar-Arevalo et al. (2001) A. Aguilar-Arevalo et al. (LSND), Phys. Rev. D64, 112007 (2001), eprint hep-ex/0104049.
  * Aguilar-Arevalo et al. (2007) A. A. Aguilar-Arevalo et al. (MiniBooNE), Phys. Rev. Lett. 98, 231801 (2007), eprint 0704.1500.
  * Abazajian et al. (2012) K. Abazajian et al. (2012), eprint 1204.5379.
  * Abratenko et al. (2021) P. Abratenko et al. (MicroBooNE) (2021), eprint 2110.00409.
  * Argüelles et al. (2021) C. A. Argüelles, I. Esteban, M. Hostert, K. J. Kelly, J. Kopp, P. A. N. Machado, I. Martinez-Soler, and Y. F. Perez-Gonzalez (2021), eprint 2111.10359.
  * Gupta et al. (2018) S. Gupta, Z. M. Matthews, P. Sharma, and A. G. Williams, Phys. Rev. D98, 035042 (2018), eprint 1804.03361.
  * Chatla et al. (2018) A. Chatla, S. Rudrabhatla, and B. A. Bambah, Adv. High Energy Phys. 2018, 2547358 (2018), eprint 1804.02818.
  * Choubey et al. (2018) S. Choubey, D. Dutta, and D. Pramanik, Eur. Phys. J. C78, 339 (2018), eprint 1711.07464.
  * Choubey et al. (2017) S. Choubey, D. Dutta, and D. Pramanik, Phys. Rev. D96, 056026 (2017), eprint 1704.07269.
  * Berryman et al. (2015) J. M. Berryman, A. de Gouvêa, K. J. Kelly, and A. Kobach, Phys. Rev. D92, 073012 (2015), eprint 1507.03986.
  * Escrihuela et al. (2015) F. J. Escrihuela, D. V. Forero, O. G. Miranda, M. Tortola, and J. W. F. Valle, Phys. Rev. D92, 053009 (2015), [Erratum: Phys. Rev.D93,no.11,119905(2016)], eprint 1503.08879.
  * Escrihuela et al. (2017) F. J. Escrihuela, D. V. Forero, O. G. Miranda, M. Tórtola, and J. W. F. Valle, New J. Phys. 19, 093005 (2017), eprint 1612.07377.
  * Ge et al. (2017) S.-F. Ge, P. Pasquini, M. Tortola, and J. W. F. Valle, Phys. Rev. D95, 033005 (2017), eprint 1605.01670.
  * Soumya and Rukmani (2018) C. Soumya and M. Rukmani, J. Phys. G45, 095003 (2018).
  * Fong et al. (2019) C. S. Fong, H. Minakata, and H. Nunokawa, JHEP 02, 015 (2019), eprint 1712.02798.
  * Verma and Bhardwaj (2018) S. Verma and S. Bhardwaj, Indian J. Phys. 92, 1161 (2018), eprint 1609.06412.
  * Antusch and Fischer (2014) S. Antusch and O. Fischer, Journal of High Energy Physics 2014 (2014), ISSN 1029-8479, URL http://dx.doi.org/10.1007/JHEP10(2014)094.
  * Kostelecky and Samuel (1989a) V. A. Kostelecky and S. Samuel, Phys. Rev. D 39, 683 (1989a).
  * Kostelecky and Samuel (1989b) V. A. Kostelecky and S. Samuel, Phys. Rev. Lett. 63, 224 (1989b).
  * Kostelecky and Potting (1991) V. A. Kostelecky and R. Potting, Nucl. Phys. B 359, 545 (1991).
  * Kostelecky and Potting (1995) V. A. Kostelecky and R. Potting, Phys. Rev. D 51, 3923 (1995), eprint hep-ph/9501341.
  * Kostelecky and Potting (1996) V. A. Kostelecky and R. Potting, Phys. Lett. B 381, 89 (1996), eprint hep-th/9605088.
  * Tanabashi et al. (2018) M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Ohlsson and Zhou (2015) T. Ohlsson and S. Zhou, Nucl. Phys. B 893, 482 (2015), eprint 1408.4722.
  * Kostelecky and Mewes (2004a) V. A. Kostelecky and M. Mewes, Phys. Rev. D 69, 016005 (2004a), eprint hep-ph/0309025.
  * Diaz and Kostelecky (2012) J. S. Diaz and A. Kostelecky, Phys. Rev. D 85, 016013 (2012), eprint 1108.1799.
  * Kostelecky and Mewes (2004b) V. A. Kostelecky and M. Mewes, Phys. Rev. D 70, 076002 (2004b), eprint hep-ph/0406255.
  * Katori et al. (2006) T. Katori, V. A. Kostelecky, and R. Tayloe, Phys. Rev. D 74, 105009 (2006), eprint hep-ph/0606154.
  * Dighe and Ray (2008) A. Dighe and S. Ray, Phys. Rev. D 78, 036002 (2008), eprint 0802.0121.
  * Barenboim and Lykken (2009) G. Barenboim and J. D. Lykken, Phys. Rev. D 80, 113008 (2009), eprint 0908.2993.
  * Rebel and Mufson (2013) B. Rebel and S. Mufson, Astropart. Phys. 48, 78 (2013), eprint 1301.4684.
  * de Gouvêa and Kelly (2017) A. de Gouvêa and K. J. Kelly, Phys. Rev. D 96, 095018 (2017), eprint 1709.06090.
  * Barenboim et al. (2018) G. Barenboim, C. A. Ternes, and M. Tórtola, Phys. Lett. B 780, 631 (2018), eprint 1712.01714.
  * Barenboim et al. (2019) G. Barenboim, M. Masud, C. A. Ternes, and M. Tórtola, Phys. Lett. B 788, 308 (2019), eprint 1805.11094.
  * Majhi et al. (2020) R. Majhi, S. Chembra, and R. Mohanta, Eur. Phys. J. C 80, 364 (2020), eprint 1907.09145.
  * Giunti and Laveder (2010) C. Giunti and M. Laveder, Phys. Rev. D 82, 113009 (2010), eprint 1008.4750.
  * Datta et al. (2004) A. Datta, R. Gandhi, P. Mehta, and S. U. Sankar, Phys. Lett. B 597, 356 (2004), eprint hep-ph/0312027.
  * Chatterjee et al. (2014) A. Chatterjee, R. Gandhi, and J. Singh, JHEP 06, 045 (2014), eprint 1402.6265.
  * Singh Koranga and Khurana (2014) B. Singh Koranga and P. Khurana, Int. J. Theor. Phys. 53, 3737 (2014).
  * Diaz and Schwetz (2016) J. S. Diaz and T. Schwetz, Phys. Rev. D 93, 093004 (2016), eprint 1603.04468.