# Probing mass orderings in presence of a very light sterile neutrino in a
liquid argon detector

Animesh Chatterjee1,2 0000-0002-2935-0958<EMAIL_ADDRESS>Srubabati
Goswami1 0000-0002-5614-4092<EMAIL_ADDRESS>Supriya Pan1,3
0000-0003-3556-8619<EMAIL_ADDRESS>1 Physical Research Laboratory,
Ahmedabad, Gujarat, 380009, India  2 European Organization for Nuclear
Research (CERN), 1211 Geneva 23, Switzerland  3 Indian Institute of
Technology, Gandhinagar, Gujarat, 382355, India

###### Abstract

Results from experiments like LSND and MiniBooNE hint towards the possible
presence of an extra eV scale sterile neutrino. The addition of such a
neutrino will significantly impact the standard three flavour neutrino
oscillations. In particular, it can give rise to additional degeneracies due
to additional sterile parameters. For an eV scale sterile neutrino, the
cosmological constraints dictate that the sterile state is heavier than the
three active states. However, for lower masses of sterile neutrinos, it can be
lighter than one and/or more of the three states. In such cases, the mass
ordering of the sterile neutrinos also becomes unknown along with the mass
ordering of the active states. In this paper, we explore the mass ordering
sensitivity in the presence of a sterile neutrino assuming the mass squared
difference $|\Delta_{41}|$ to be in the range $10^{-4}-1$ eV2. We study (i)
how the ordering of the active states, i.e. the determination of the sign of
$\Delta_{31}$ gets affected by the presence of a sterile neutrino in the above
mass range, (ii) the possible determination of the sign of $\Delta_{41}$ for
$\Delta_{41}$ in the range $10^{-4}-0.1$ eV2. This analysis is done in the
context of a liquid argon detector using both beam neutrinos traveling a
distance of 1300 km and atmospheric neutrinos which propagates through a
distance ranging from 10 - 10000 km allowing resonant matter effects. Apart
from presenting separate results from these sources, we also do a combined
study and probe the synergy between these two in giving an enhanced
sensitivity.

sterile neutrino and mass ordering and liquid argon

##  I Introduction

Neutrino oscillations in which one flavour of neutrinos gets converted to
others, have been discovered using diverse neutrino sources: from the sun to
the atmosphere, from reactors to accelerators, and a variety of detection
techniques in different terrestrial experiments. This has validated the three
neutrino oscillation picture and most of the parameters governing the
phenomenon have been determined with sizable precision. These parameters are
the mass squared differences $\Delta_{21}$, $|\Delta_{31}|$ ($\Delta_{ij}$ is
defined as $m_{i}^{2}-m_{j}^{2}$ for the mass eigenvalues $m_{i},m_{j}$), and
mixing angles $\theta_{12}$, $\theta_{23}$. The remaining unknowns of the
three flavour paradigm are the atmospheric mass ordering (the sign of
$\Delta_{31}$), the octant of the mixing angle $\theta_{23}$, and the CP phase
$\delta_{CP}$. The next generation beam based experiments like DUNE [1, 2],
T2HK [3], ESS$\nu$SB [4] are expected to resolve these issues. These
experiments plan to use a beam of higher intensity and larger volume detectors
which can increase the statistics as compared to the current generation
accelerator experiments T2K [5] and NO$\nu$A[6]. Apart from these there are
dedicated future atmospheric neutrino experiments like Hyper Kamiokande[7],
India-based Neutrino Observatory [8], IceCube and Pingu[9] which can also help
in throwing light on the final frontier of the three neutrino oscillation. The
liquid argon detector in DUNE will also be able to observe atmospheric
neutrinos [10, 11, 12]. Apart from the determination of parameters in the
standard three flavour oscillation, these next generation experiments also
open up the possibility of probing beyond standard model (BSM) physics, which
can occur at a sub-leading level. Several new physics scenarios like light
sterile neutrinos, non-standard interactions, long-range forces, neutrino
decay, and violation of fundamental symmetries like CPT, Lorentz-invariance,
etc have been explored in the context of neutrino oscillation experiments[13,
14].

Among these, the light sterile neutrino scenario is motivated by three long-
standing anomalies that have served as primary drivers in the development of a
vibrant short-baseline neutrino program over the last decade. Two of these
pieces of evidence come from the apparent oscillatory appearance of electron
(anti)neutrinos in muon-(anti)neutrino beams originating from charged-pion
decay-at-rest in the LSND experiment [15, 16], and charged-pion decay-in-
flight in the MiniBooNE experiment [17, 18]. There have also been an anomaly
associated with an overall normalization discrepancy of electron
(anti)neutrinos expected in the radioactive decay of Gallium-71[19, 20, 21].
One of the most theoretically motivated frameworks considered for the
interpretation of these anomalies is that of the presence of light ($\approx$
1 eV) sterile neutrino state[22].

A sterile neutrino is a neutral SU(2)$\times$U(1) singlet with no ordinary
weak interaction except those induced by the mixing. Very heavy sterile
neutrinos ($10^{14}-10^{16}$ GeV) are proposed as the mediators in the type I
seesaw model[23, 24, 25] which can give rise to small neutrino masses. Such
neutrinos also play a significant role in leptogenesis[26, 27]. Such neutrinos
are natural candidates in grand unified theories. Sterile neutrinos of TeV
energies have also been studied in the context of low-scale seesaw models[28,
29]. Sterile neutrinos of keV mass are especially interesting because the
sterile neutrinos would be a viable dark matter candidate[30].

The existence of an eV scale sterile neutrino motivated by the short-baseline
anomalies is in strong tension from cosmological bound on the effective number
of neutrinos $N_{eff}$ and the sum of total masses of light neutrinos. From
the recent measurement of Planck data[31] and combining together with the
Hubble parameter measurement [32] and Supernova Ia data from the Pantheon
sample [33], the extended fit to the parameters are

$\begin{split}N_{eff}&=3.11^{+0.37}_{-0.36}(95\%\text{ CL})\\\ &\sum
m_{\nu}<0.16\text{ eV}\end{split}$ (1)

In order to get around cosmological constraints, we need to introduce new
physics that directly affects the cosmological phenomenology of the light
sterile states. Since the main problem of the canonical light sterile neutrino
is that its thermalization in the early universe raises $N_{eff}$ to an
unacceptably large level for BBN and CMB/LSS constraints. All known new
physics solutions so far involve tampering with the thermalization process, in
order to maintain $N_{eff}$ as close to the SM value as possible. A number of
ideas have been proposed and explored throughout the years, e.g., large
chemical potentials or, equivalently, number density asymmetries for the
active neutrinos[34]; secret interactions of the sterile neutrinos[35, 36] and
low reheating temperature of the universe[37], etc.

Can there be sterile neutrinos lighter than the eV scale? In the presence of a
sterile neutrino, there is a new mass squared difference
$\Delta_{41}=m_{4}^{2}-m_{1}^{2}$. A very light sterile neutrino corresponding
to the mass-squared difference in ranges $10^{-4}-0.1$ eV2 is expected to be
consistent with cosmological mass bounds. It was suggested in ref. [38] that
the existence of a very light ($\approx 10^{-5}eV^{2}$) sterile neutrino can
provide the explanation for the lack of upturn in the solar neutrino
oscillation probability below $\approx$ 8 MeV. A recent study has probed the
possibility of alleviating the tension between the results of the ongoing beam
experiments, T2K and NO$\nu$A for the value $\delta_{cp}$ using very light
sterile neutrino with a wide mass difference range of $10^{-5}:0.1$ eV2[39].

We focus our study on only one sterile neutrino added to the three light
neutrinos, namely the 3+1 framework, and consider a wide mass range for
$|\Delta_{41}|$ varying in the range of $10^{-4}-1$ eV2. The cosmological
constraints on the sum of all the neutrino masses imply that the sign of
$\Delta_{41}$ can not be negative for $\Delta_{41}>0.1$ eV2. However, for
lower mass squared differences both signs of $\Delta_{41}$ are possible. In
this work, we investigate the possibility of determining (i) the sign of
$\Delta_{31}$ in the presence of a sterile neutrino corresponding to a)
$\Delta_{41}=1$ eV2, b)$\Delta_{41}$ in the range of $10^{-4}-0.1$ eV2; (ii)
the sign of $\Delta m^{2}_{41}$ for the mass range $10^{-4}-0.1$ eV2. To
answer these questions we use a liquid argon time projection chamber (LArTPC)
capable of detecting both beam and atmospheric neutrinos. The typical baseline
we have used for the beam neutrinos is $\sim$ 1300 km which is similar to the
DUNE experiment. We delineate the sensitivities to mass ordering by performing
a combined analysis of beam and atmospheric neutrinos, along with a separate
study for each. Additionally, we present the results including the charge
tagging capability of muon capture in liquid argon allowing one to
differentiate between $\mu^{+}$ and $\mu^{-}$ events in the context of
atmospheric neutrinos.

The mass ordering in the presence of a light sterile neutrino has been studied
in ref. [40] with the additional mass squared difference varying in a wide
range in the context of a magnetized iron calorimeter detector proposed by the
India-based Neutrino Observatory (INO) collaboration. Recently the sensitivity
of the sterile mass ordering in the same mass range has been studied in
reference [41] in the context of the DUNE experiment using beam neutrinos. We
perform our study in the context of a liquid argon time projection chamber
detector as in DUNE using both beam and atmospheric neutrino events separately
as well as in a combined analysis.

The plan of the paper is as follows. In section II, we present the 3+1
framework which is used for the analysis. In section III, we discuss the
probability level study of $P_{\mu e}$, $P_{\mu\mu}$ in the presence of a
sterile neutrino and explore the effect of sterile mixing and point out where
these effects will be significant. Simulation procedures used for the
neutrinos coming from the beam and atmosphere, detector specification, and
numerical analysis are given in section IV. Next, in section V, we present and
discuss the results. Finally, we conclude in section VI.

##  II The 3+1 framework:

The minimal scheme postulated to explain the LSND and MiniBooNE results is the
addition of one light sterile neutrino of mass of the order of eV scale to the
three active neutrinos in the SM[42, 43, 44, 43, 15, 44, 45, 22, 46]. In this
case, there is one additional independent mass squared difference which we
take as $\Delta_{41}={m^{2}_{4}}-{m^{2}_{1}}$.

The possible mass orderings, in this case, are shown in fig. 1. We will refer
to the sign of $\Delta_{31}$ as the atmospheric mass ordering (AMO) and the
sign of $\Delta_{41}$ as sterile mass ordering (SMO) throughout the paper.
NO(IO) and SNO(SIO) will signify normal(inverted) ordering for AMO and SMO
respectively. There can be four possibilities as follows,
(i) SNO-NO: where $\Delta_{41}>0$ and $\Delta_{31}>0$. The positioning of the
4th state depends on the value of $\Delta_{41}$. For $\Delta_{41}>10^{-3}$ eV2
the 4th state lies above the 3rd state while if $\Delta_{41}<10^{-3}$ eV2 it
lies below the 3rd state.
(ii) SNO-IO: in this case $\Delta_{41}>0$ and $\Delta_{31}<0$ corresponding to
inverted ordering of the light active neutrino. The 4th state lies above the
three active states with positioning depending on the value of $\Delta_{41}$.
(iii) SIO-NO: this signifies to $\Delta_{41}<0$ and $\Delta_{31}>0$. The 4th
state will always lie below the lightest active states with the placement
depending on the value of $\Delta_{41}$.
(iv) SIO-IO: for this case both $\Delta_{41}$ and $\Delta_{31}$ are $<0$. For
$\Delta_{41}<10^{-3}$ eV2, the 4th state lies above $m_{3}$.

Note that the usual 3+1 picture corresponds to the cases (i) and (ii) with
$\Delta m^{2}_{41}\sim$ eV2. The Cases (iii) and (iv) with $\Delta
m^{2}_{41}\sim$ eV2 are disfavored from cosmology.

Figure 1: The 3+1 mass spectrum: ordering of mass states in the presence of an
extra sterile neutrino state $m_{4}$ (blue) corresponding to two different
sterile mass squared difference: A. $|\Delta_{41}|\sim 1$ eV2, B.
$|\Delta_{41}|\sim 10^{-4}$ eV2 when the standard mass ordering $\Delta_{31}$
lead by $m_{3}$ (red) can be both +ve and -ve.

For the $3+1$ oscillation framework, the new mixing matrix $U$ will have three
additional mixing angles $\theta_{14},\theta_{24},\theta_{34}$ corresponding
to mixing between the extra light sterile neutrino $\nu_{s}$ and the active
sector neutrinos along with two new CP phases $\delta_{14},\delta_{34}$ and
can be expressed as follows,

$\displaystyle U$
$\displaystyle=\tilde{R}_{34}(\theta_{34},\delta_{34})R_{24}(\theta_{24})\tilde{R}_{14}(\theta_{14},\delta_{14})\times$
$\displaystyle
R_{23}(\theta_{23})\tilde{R}_{13}(\theta_{13},\delta_{13})R_{12}(\theta_{12})$
(2)

where
$\tilde{R}_{ij}=U^{\delta}_{ij}(\delta_{ij})R_{ij}(\theta_{ij})U^{\dagger\delta}_{ij}(\delta_{ij})$,
$R_{ij}(\theta_{ij})$’s are rotational matrices in i-j plane and
$U^{\delta}_{ij}=diag(1,1,1,e^{\iota\delta_{ij}})$ with $\delta_{ij}$’s being
the CP phases.

The global analysis of data performed in ref. [47] gives the best-fit values
and sterile mixing angles for $\Delta_{41}=1.3$ eV2 as given in Table 1,

Parameters |  $3\sigma$ range | Best Fit
---|---|---
$\sin^{2}2\theta_{14}$ | 0.04 - 0.09 | 0.079
$\theta_{14}$ | $5.76^{\circ}-8.73^{\circ}$ | $8.15^{\circ}$
$\sin^{2}\theta_{24}$ | $6.7\times 10^{-3}-0.022$ | 0.015
$\theta_{24}$ | $4.68^{\circ}-8.6^{\circ}$ | $7.08^{\circ}$
Table 1: $3\sigma$ Levels and Best fit values extracted from [47]

Ref. [47] provides an upper bound of $\theta_{34}\leq 7.4^{\circ}$ at the
$90\%$ C.L. However, the analysis performed in [48] including the MINOS+ data
disfavoured the allowed regions in $\theta_{24}$ from above with a new bound
at 90% C.L. $\sin^{2}\theta_{24}\leq 0.006$, i.e., $\theta_{24}\leq
4.5^{\circ}$. Also, the analysis[48] of DayaBay and Bugey3 gives at 90% C.L.
$\sin^{2}2\theta_{14}\leq 0.046$. i.e., $\theta_{14}\leq 6.2^{\circ}$. The
bounds on these sterile mixing angles vary with $\Delta_{41}$. For values of
$\Delta_{41}$ in the range $10^{-4}:0.1$ eV2, the 90% C.L. bounds vary from
$\theta_{14}\leq 25^{\circ}:\theta_{14}\leq 2^{\circ}$, and $\theta_{24}\leq
56^{\circ}:\theta_{24}\leq 4^{\circ}$.

The effective interaction Hamiltonian in matter for the 3+1 framework in
flavour basis is given as follows,

$\displaystyle H_{int}$ $\displaystyle=diag(V_{CC},0,0,-V_{NC})$
$\displaystyle=diag(\sqrt{2}G_{F}N_{e},0,0,\sqrt{2}G_{F}N_{n}/2),$ (3)

where $V_{CC}=\sqrt{2}G_{F}N_{e}$ is the charge current interaction potential,
$V_{NC}=-\sqrt{2}G_{F}N_{n}/2$ is the neutral current interaction potential
and $G_{F}$ is the Fermi coupling constant with $N_{e},N_{n}$ corresponding to
the density of electron and neutron respectively of the medium in which
neutrinos travel. In the presence of matter, the total Hamiltonian in flavour
basis is expressed as follows,

$H=\frac{1}{2E_{\nu}}U\begin{bmatrix}0&0&0&0\\\ 0&\Delta_{21}&0&0\\\
0&0&\Delta_{31}&0\\\
0&0&0&\Delta_{41}\end{bmatrix}U^{\dagger}+\frac{1}{2E_{\nu}}\begin{bmatrix}A&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&\frac{A}{2}\end{bmatrix}$ (4)

where the propagation medium has been considered to be the Earth’s matter with
neutron density being equal to electron density, i.e, $N_{e}=N_{n}$ and the
matter potential term is $A=2\sqrt{2}G_{F}N_{e}E_{\nu}$ for neutrino with
energy $E_{\nu}$.

##  III Oscillation Probabilities in presence of a sterile neutrino

In this section, we present both the appearance probability $P_{e\mu}$ and
disappearance probability $P_{\mu\mu}$ at baselines of 1300 km and 7000 km for
different sterile parameters. As in the three flavour framework, there are two
useful approximate methods to calculate the oscillation probability in the
matter: (i) $\alpha-s_{13}$ approximation[41], (ii) $\Delta_{21}=0$
approximation[49]. In the first case
$\alpha=\frac{\Delta_{21}}{\Delta_{31}},\sin{\theta_{13}}$ are considered
small parameters whereas in the second case $\Delta_{21}$ is considered as
negligible. The probabilities evaluated in these approximations considering
$\theta_{34}=0^{\circ}$ are presented in the following section.

###  III.1 $\alpha-s_{13}$ Approximation

The appearance probability calculated using $\alpha-s_{13}$ approximation in
[41] is well suited for any $\Delta_{41}$ values and baseline of 1300 km and
can be expressed as,

$\begin{split}P_{\mu
e}^{m}&=4s_{13}^{2}s_{23}^{2}\frac{\sin^{2}[(A^{\prime}-1)\Delta]}{(A^{\prime}-1)^{2}}+8\alpha
s_{13}s_{12}c_{12}s_{23}c_{23}\frac{\sin[A^{\prime}\Delta]}{A^{\prime}}\frac{\sin[(A^{\prime}-1)\Delta]}{A^{\prime}-1}\cos(\Delta+\delta_{13})\\\
+&4s_{13}s_{14}s_{24}s_{23}\frac{\sin[(A^{\prime}-1)\Delta]}{A^{\prime}-1}[P_{14}^{s}\sin{\tilde{\delta}_{14}}+P_{14}^{c}\cos{\tilde{\delta}_{14}}]\end{split},$
(5)

where the terms corresponding to sterile neutrino are,

$\begin{split}P_{14}^{s}=R[\frac{1}{2}A^{\prime}c_{23}+(R-1)(1+s_{23}^{2})]\frac{\sin[(R-1+\frac{A^{\prime}}{2})\Delta]}{R-1+\frac{A^{\prime}}{2}}\frac{\sin[(R-\frac{A^{\prime}}{2})\Delta]}{R-\frac{A^{\prime}}{2}}+Rc_{23}^{2}\sin[(R-1-\frac{A^{\prime}}{2})\Delta]\frac{\sin[(R+\frac{A^{\prime}}{2})\Delta]}{R+\frac{A^{\prime}}{2}}\end{split}$
(6)
$\begin{split}P_{14}^{c}&=\frac{R}{R-\frac{1}{2}}\left([R-\frac{1}{2}s_{23}^{2}-\frac{1}{2}]\cos[(R-1-\frac{A^{\prime}}{2})\Delta]\frac{\sin[(R-\frac{A^{\prime}}{2})\Delta]}{R-\frac{A^{\prime}}{2}}+s_{23}^{2}(R-1)\cos[(R-\frac{A^{\prime}}{2})\Delta]\frac{\sin[(R-1+\frac{A^{\prime}}{2})\Delta]}{R-1+\frac{A^{\prime}}{2}}\right.\\\
&\left.+s_{23}^{2}\frac{\sin[(A^{\prime}-1)\Delta]}{A^{\prime}-1}\right)+Rc_{23}^{2}\cos[(R-1+\frac{A^{\prime}}{2})\Delta]\frac{\sin[(R+\frac{A^{\prime}}{2})\Delta]}{R+\frac{A^{\prime}}{2}}\end{split},$
(7)

where $A^{\prime}=\frac{A}{\Delta_{31}}$, $R=\frac{\Delta_{41}}{\Delta_{31}}$,
$\Delta=\frac{1.27\Delta_{31}L}{E}$,
$\tilde{\delta}_{14}=\delta_{13}+\delta_{14}$ and
$c_{ij}\sim\cos\theta_{ij},s_{ij}\sim\sin\theta_{ij}.$
at limit $R>>1$, for $R>>\frac{A^{\prime}}{2}$, approximately
$R-\frac{A^{\prime}}{2}\simeq R+\frac{A^{\prime}}{2}\simeq R$, also
$R-\frac{1}{2}\simeq R$, $R-\frac{1}{2}s_{23}^{2}-\frac{1}{2}\simeq R$

$\displaystyle P_{14}^{s}\simeq$
$\displaystyle\frac{1}{2}A^{\prime}c_{23}\frac{\sin[(R-1)\Delta]}{R-1}\sin[R\Delta]+2\sin[(R-1)\Delta]\sin[R\Delta]\simeq(\frac{1}{2R}A^{\prime}c_{23}+2)\sin[R\Delta]^{2}$
(8) $\displaystyle P_{14}^{c}\simeq$
$\displaystyle(1+c_{23}^{2})\cos[(R-1)\Delta]\sin[R\Delta]+s_{23}^{2}\cos[R\Delta]\sin[(R-1)\Delta]+s_{23}^{2}\frac{\sin[(A^{\prime}-1)\Delta]}{A^{\prime}-1}\simeq\sin[2R\Delta]+s_{23}^{2}\frac{\sin[(A^{\prime}-1)\Delta]}{A^{\prime}-1}$
(9)

###  III.2 $\Delta_{21}=0$ Approximation

The dominant term in the $\nu_{\mu}-\nu_{e}$ oscillation probability in OMSD
approximation, valid for $\Delta_{21}L/E<<1$, e.g., at 7000 km baseline around
the resonance energy (7 GeV), is given by [49],

$\begin{split}P_{\mu
e}^{1}&=4\cos^{2}\theta_{13m}\cos^{2}\theta_{14m}\sin^{2}\theta_{13m}(\cos^{2}\theta_{24m}\sin^{2}\theta_{23}-\sin^{2}\theta_{14m}\sin^{2}\theta_{24m})\sin^{2}\frac{1.27\Delta_{31}^{m}L}{E}\\\
&+2\cos^{3}\theta_{13m}\cos^{2}\theta_{14m}\sin\theta_{13m}\sin\theta_{14m}\sin
2\theta_{24m}\sin\theta_{23}\sin{\frac{1.27\Delta_{31}^{m}L}{E}}\sin(\frac{1.27\Delta_{31}^{m}L}{E}+\delta_{13}-\delta_{14})\\\
&-2\cos\theta_{13m}\cos^{2}\theta_{14m}\sin^{3}\theta_{13m}\sin\theta_{14m}\sin
2\theta_{24m}\sin\theta_{23}\sin{\frac{1.27\Delta_{31}^{m}L}{E}}\sin(\frac{1.27\Delta_{31}^{m}L}{E}-\delta_{13}+\delta_{14})\end{split}$
(10)

where $\theta_{13m},\theta_{14m},\theta_{24m}$ are modified mixing angles,
$\Delta_{ij}^{m}$’s are modified mass-squared differences as defined in ref.
[49].

###  III.3 Effect of sterile parameters on sign of $\Delta_{31}$

In this section, we discuss the dependence of the probabilities on the sterile
parameters. For our study, we chose two illustrative mixing angles
$\theta_{14}=\theta_{24}=4^{\circ}$ and/or $7^{\circ}$. The phases
$\delta_{13}$ and $\delta_{14}$ are varied in the range $-180^{\circ}$ to
$180^{\circ}$, unless otherwise mentioned. The probabilities for both
$\theta_{34}=0^{\circ}$ and non-zero $\theta_{34}=7^{\circ}$ are discussed. In
this section, GLoBES[50] package is used to generate the probabilities.

####  III.3.1 Effect of non-zero $\theta_{14},\theta_{24}$

Figure 2: Probabilities $P_{\mu e}$(left) and $P_{\mu\mu}$ (right) as a
function of energy $E_{\nu}$ due to variation of phases
$\delta_{13},\delta_{14}$ for NO and IO at 1300 km baseline for
$\Delta_{41}=1$ eV2. Blue and orange bands in top (bottom) panels refer to
varied phases for $\theta_{14}=4^{\circ}$ ($7^{\circ}$),
$\theta_{24}=4^{\circ}$ ($7^{\circ}$) corresponding to NO and IO respectively.
The regions between cyan(yellow) curves are due to the variation of
$\delta_{13}$ in $3\nu$ case for NO(IO).

We have plotted the appearance (left), and disappearance (right) probabilities
in fig. 2 as a function of the neutrino energy varying the phases
$\delta_{13},\delta_{14}$. The blue[NO] and orange[IO] bands at the top
(bottom) panels refer to mixing angles
$\theta_{14}=\theta_{24}=4^{\circ}(7^{\circ})$ in the 3+1 framework. Whereas,
the regions between cyan and yellow lines correspond to the variation of
$\delta_{13}$ in three generation framework in NO and IO respectively. In the
right panel, we also show $P_{\mu\mu}$ over $2-4$ GeV in a magnified inset.
The important observations are as follows,

  * •

In the 3+1 framework, the probability bands corresponding to NO and IO are
closer than those in the three generation framework. This suggests a reduced
hierarchy sensitivity in the 3+1 framework.

  * •

The gap between NO and IO bands increases with the decreasing values of the
sterile mixing angles.

  * •

In $P_{\mu e}$ channel, for the 3+1 framework the difference between the two
probability bands is seen in the energy range 1-3 GeV.

  * •

The disappearance channel probability $P_{\mu\mu}$ doesn’t depend on the
phases, as can be seen from the narrow bands in NO and IO cases, for both
three and 3+1 frameworks.

  * •

The $P_{\mu\mu}$ curves for opposite mass orderings are hard to separate from
each other at energies lower than 2 GeV. However, some demarcation is visible
at energies in the range of 2-7 GeV for both three generation and 3+1
generation, as shown in the right panels.

  * •

In the disappearance channel, we don’t see a significant effect of variation
of the sterile-active mixing angles $\theta_{14},\theta_{24}$ on the
probability bands.

Figure 3: Probabilities $P_{\mu e}$(left) and $P_{\mu\mu}$ (right) as a
function of energy $E_{\nu}$ at 7000 km baseline for $\Delta_{41}=1$ eV2. The
shaded bands of blue and orange refer to varied phases
$\delta_{13},\delta_{14}$ for $\theta_{14}=4^{\circ}$ ($7^{\circ}$) for
NO(IO). The regions between cyan(yellow) curves are due to the variation of
$\delta_{13}$ in $3\nu$ case for NO(IO).

In fig. 3, the appearance (left) and disappearance (right) probabilities have
been shown at 7000 km baseline. Such baselines will be relevant for
atmospheric neutrinos. The significant observations are as follows,

  * •

There is a prominent difference between the probability bands of NO and IO,
implying sensitivity to mass ordering at 7000 km baseline in the 3+1
framework. The difference decreases for sterile case w.r.t. the standard one.

  * •

Also, with lower values of $\theta_{14},\theta_{24}$, the gap between the
probability bands of opposite mass ordering increases.

  * •

For $P_{\mu\mu}$ channel, a significant gap exists between bands of opposite
mass orderings at energies higher than 4 GeV whereas, in $P_{\mu e}$ channel,
the gap is present even at much lower energies of 3 GeV.

In order to understand the features of the fig. 2, 3 from analytical
expression, we compute the minimum difference in the probability $P_{\mu e}$
for two different mass orderings by varying the phases,

$\Delta P\equiv[P_{\mu e}^{1NO}(\delta_{13}^{NO},\delta_{14}^{NO})-P_{\mu
e}^{1IO}(\delta_{13}^{IO},\delta_{14}^{IO})]_{min}$ (11)

where the other parameters have been kept fixed. Using only the dominant first
term of (10), the difference in probability can be expressed as,

$\begin{split}\Delta P&=\Delta P_{np}+\\\
&A_{1}[\sin^{2}(M-N)\cos\delta^{NO}-\sin^{2}(M+N)\cos\delta^{IO}]+\\\
&A_{2}[\sin 2(M-N)\sin\delta^{NO}+\sin 2(M+N)\sin\delta^{IO}]\end{split}$ (12)

where $\Delta P_{np}$ is the part with no phases involved, and is given as
follows,

$\displaystyle\Delta P_{np}$
$\displaystyle=\cos^{2}\theta_{14m}\sin^{2}2\theta_{13m}(\sin^{2}\theta_{24m}\sin^{2}\theta_{14m}-$
$\displaystyle\cos^{2}\theta_{24m}\sin^{2}\theta_{23})\sin 2M\sin 2N$ (13)

The other part containing the phases depends on the amplitude parameters
$A_{1},A_{2}$ and the frequency parameters $M,N$ that are defined as,

$\displaystyle A_{1}=\cos^{2}\theta_{14m}\cos 2\theta_{13m}\sin
2\theta_{13m}\sin\theta_{14m}\sin 2\theta_{24m}\sin\theta_{23}$ (14)
$\displaystyle A_{2}=\cos^{2}\theta_{14m}\sin
2\theta_{13m}\sin\theta_{14m}\sin 2\theta_{24m}\sin\theta_{23}$ (15)
$\displaystyle M$ $\displaystyle=\Delta_{31}\cos
2(\theta_{13}-\theta_{13m})\times\frac{1.27L}{E}$ (16) $\displaystyle N$
$\displaystyle=A\cos
2\theta_{13m}(1+\cos^{2}\theta_{14}+\cos^{2}\theta_{14}\sin^{2}\theta_{24})\times\frac{1.27L}{E}$
(17)

We use $\Delta_{31}=2.5\times 10^{-3}$ eV2, $\theta_{24}=7^{\circ}$,
$\theta_{14}=7^{\circ}$, $\theta_{23}=45^{\circ}$ to calculate the $\Delta P$
at the first oscillation maxima. For $L=1300$ km at $E=2.5$ GeV, $\Delta P<0$,
which means the minima of the NO curve is below the maxima of the IO curve.
This implies the overlap of NO and IO bands suggesting the presence of
degeneracy. Whereas, for 7000 Km at $E=7$ GeV, $\Delta P>0$, i.e., mass
ordering can be determined even with varying the phases $\delta_{13}$ and
$\delta_{14}$.

####  III.3.2 Effect of non-zero $\theta_{34}$

Figure 4: $P_{\mu e}$ as a function of energy $E_{\nu}$ with the variation of
phases $\delta_{13},\delta_{14},\delta_{34}$ for NO (blue) and IO(orange),
$\Delta_{41}=1$ eV2 at 1300 km (left) and 7000 km (right). Shaded blue
(orange) bands refer to $\theta_{14},\theta_{24},\theta_{34}=7^{\circ}$ for
NO(IO). Regions between cyan(yellow) curves are due to variation of
$\delta_{13}$ in $3\nu$ case for NO(IO).

In fig. 4, the appearance probability is plotted as a function of neutrino
energy at 1300 km (left) and 7000 km (right) baseline for
$\theta_{14},\theta_{24},\theta_{34}=7^{\circ}$. The shaded blue (yellow)
bands are due to the variation of the phases
$\delta_{13},\delta_{14},\delta_{34}$ for NO(IO). The regions between the
solid cyan(yellow) curves correspond to the variation of $\delta_{13}$ in the
3$\nu$ framework for NO(IO). The most important observation is a notable
decrease in the gap between NO and IO bands at both baselines. In the case of
7000 km, there is still a gap between the opposite mass ordering bands, which
gets diminished for a non-zero $\theta_{34}$.

####  III.3.3 Effect of $|\Delta_{41}|$

In this section, the sensitivity to atmospheric mass ordering is studied at
the probability level, with the sterile mass squared difference $\Delta_{41}$
in the range $10^{-4}$ \- $0.1$ eV2.

To understand the effect of $\Delta_{41}$ on sensitivity to atmospheric mass
ordering, we probe the difference in the appearance probability $\Delta P_{\mu
e}$ as a function of $\Delta_{41}$.

$\displaystyle\Delta P_{\mu e}=|P_{\mu
e}^{true}(\Delta_{31},\delta_{13},\delta_{14})-P_{\mu
e}^{test}(-\Delta_{31}^{\prime},\delta_{13}^{\prime},\delta_{14}^{\prime})|_{min}$
(18)

We compute the minimum difference $\Delta P_{\mu e}$ (using GLoBES), by
considering a particular AMO with constant $\delta_{13},\delta_{14}$ in
$P_{\mu e}^{true}$ whereas $P_{\mu e}^{test}$ is calculated for the opposite
AMO by varying the phases and $|\Delta_{31}|$. The phases are varied over
their full range and for $\Delta_{31}$ variation over the current $3\sigma$
range in the opposite mass ordering is considered.

Figure 5: The difference in appearance channel probability $\Delta P_{\mu e}$
for opposite AMO as a function of $\Delta_{41}^{true}$ and $E_{\nu}$ at 1300
km baseline for SNO-NO (above) and SNO-IO (below). The plots in the left
(right) panel has $\delta_{13}^{true}=90^{\circ}(-90^{\circ})$ and fixed
$\delta_{14}=0^{\circ}$.

In fig. 5, we illustrate the probability difference $\Delta P_{\mu e}$ in the
$\Delta_{41}-E_{\nu}$ plane at 1300 km. The important observations are as
follows,

  * •

Around $\Delta_{41}=2.5\times 10^{-3}$ eV2, $\Delta P_{\mu e}$ is either very
high or low depending upon the values of $\delta_{13}$ and NO/IO.

  * •

We observe an oscillating pattern of $\Delta P_{\mu e}$ along $\Delta_{41}$
for a fixed energy. This oscillation becomes rapid at higher $\Delta_{41}$
values.

  * •

Significant contribution to $\Delta P_{\mu e}$ is seen for energies in the
range of $1.5-4$ GeV.

  * •

For the SNO-NO case (top panels), the occurrence for maxima and minima
reverses for $\delta_{13}=90^{\circ}$, and $-90^{\circ}$.

  * •

However, for the SNO-IO case, the maxima and minima occur at the same
$\Delta_{41}$ for $\delta_{13}=90^{\circ},-90^{\circ}$. Although, the
magnitude is higher for $\delta_{13}=90^{\circ}$.

###  III.4 Effect of the sign of $\Delta_{41}$ in $P_{\mu e}$ channel

As we consider $\Delta_{41}$ in the rage $5\times 10^{-4}:10^{-1}$ eV2, the
sterile mass ordering also becomes unknown giving us four possibilities
depending on the ordering of the three active states. The difference in the
probability for the opposite signs of $\Delta_{41}$ is defined as,

$\Delta P_{s}=|P_{\mu
e,\mu\mu}^{true}(+\Delta_{41},\delta_{13},\delta_{14})-P_{\mu
e,\mu\mu}^{test}(-\Delta_{41}^{\prime},\delta_{13}^{\prime},\delta_{14}^{\prime})|_{min}$
(19)

Figure 6: Difference in the appearance probability $\Delta P_{\mu e}$(left),
and disappearance probability $\Delta P_{\mu\mu}$(right) for different sterile
mass ordering in the $\Delta_{41}-E_{\nu}$ plane at 1300 km(top), 7000
km(bottom).

We plot $\Delta P_{s}$ (using GLoBES) by marginalizing over
$\Delta_{41}^{\prime}$ in opposite SMO and phases
$\delta_{13}^{\prime},\delta_{14}^{\prime}$ in fig. 6 for appearance (left)
and disappearance channel(right) over a wide range of the sterile mass squared
difference and neutrino energy at 1300 km baseline (top) and 7000 km (bottom).
We have taken constant values of $\theta_{14},\theta_{24}=7^{\circ}$,
$\delta_{13}=-90^{\circ},\delta_{14}=90^{\circ}$ in $P_{\mu e,\mu\mu}^{true}$.
It can be observed that the high values of $\Delta P_{s}$ are mostly
concentrated in the mass square range of $10^{-3}:10^{-2}$ eV2 range. In
$P_{\mu e}$ channel, the contribution is lower than $P_{\mu\mu}$. The
difference is observed to be larger at higher baselines.

In $P_{\mu e}$ channel, probability difference is higher around
$\Delta_{41}=1-2\times 10^{-3}$ eV2 while a dip is found immediately after
around $\Delta_{41}=2.5-3\times 10^{-3}$ eV2 . In $P_{\mu\mu}$ channel for
energy below 4 GeV, a similar pattern is observed. However, at higher
energies, a higher $\Delta P$ value is observed.

##  IV Simulation procedure and details of the experimental setup

The experimental setup under consideration consists of a megawatt-scale muon
neutrino beam source accompanied by a near detector (ND) and a far
detector(FD). The ND will be placed close to the source of the beam, while the
FD, comprising a 40 Kton LArTPC detector is placed at a distance of 1300 km
away from the neutrino source. The large LArTPC at an underground observatory
is also capable of observing atmospheric neutrinos. The proposed DUNE
experiment has a similar experimental configuration[1]. In this analysis, both
neutrino beam coming from the accelerator and the atmospheric neutrinos have
been considered.

###  IV.1 Events from accelerator beam

A beam-power of 1.2MW leading to a total exposure of $10\times 10^{21}$ pot
has been implemented for the numerical analysis. The neutrino beam simulation
has been carried out using the GLoBES[50] software. We assume the experiment
to be running for 3.5 years each in the neutrino mode and the antineutrino
mode.

###  IV.2 Atmospheric events

The atmospheric neutrinos produce muons and electrons. These event rates
corresponding to an energy bin of width $\mathrm{dE}$ and in a solid angle bin
of width ${\mathrm{d\Omega}}$ are expressed as,

$\rm{\frac{d^{2}N_{\mu}}{d\Omega\;dE}=\frac{\sigma_{c}D_{eff}}{2\pi}\left[\left(\frac{d^{2}\Phi_{\mu}}{d\cos\theta\;dE}\right)P_{\mu\mu}+\left(\frac{d^{2}\Phi_{e}}{d\cos\theta\;dE}\right)P_{e\mu}\right]},$
(20)
$\mathrm{\frac{d^{2}N_{e}}{d\Omega\;dE}=\frac{\sigma_{c}D_{eff}}{2\pi}\left[\left(\frac{d^{2}\Phi_{\mu}}{d\cos\theta\;dE}\right)P_{\mu
e}+\left(\frac{d^{2}\Phi_{e}}{d\cos\theta\;dE}\right)P_{ee}\right]}$ (21)

Here ${\mathrm{\Phi_{\mu}}}$ and ${\mathrm{\Phi_{e}}}$ stands for the
$\nu_{\mu}$ and ${\mathrm{\nu_{e}}}$ atmospheric fluxes [51] respectively, the
disappearance and appearance probabilities are given as $P_{\mu\mu}$ and
$P_{e\mu}$ respectively, $\rm{\sigma_{c}}$ is the total CC cross section and
$\rm{D_{eff}}$ is the detector efficiency. The energy and angular resolutions
for the LArTPC detector, implemented using Gaussian resolution function R, are
defined as follows,

$\mathrm{R_{E}(E_{t},E_{m})=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(E_{m}-E_{t})^{2}}{2\sigma^{2}}\right]}\,.$
(22)
$\rm{R_{\theta}(\Omega_{t},\Omega_{m})=N\exp\left[-\frac{(\theta_{t}-\theta_{m})^{2}+\sin^{2}\theta_{t}~{}(\phi_{t}-\phi_{m})^{2}}{2(\Delta\theta)^{2}}\right]}\,,$
(23)

where N is the normalization constant. Here, $\rm{E_{m}}$
(${\rm{\Omega_{m}}}$) and $\rm{E_{t}}$ (${\rm{\Omega_{t}}}$), denote the
measured and true values of energy (angle) respectively. The smearing width
$\sigma$ is a function of $\rm{E_{t}}$. The smearing function for the zenith
angle is a bit more complicated. The direction of the incident (measured)
neutrino is specified by two variables: the polar angle
${\rm{\theta_{t}}}({\rm{\theta_{m}}})$ and the azimuthal angle
${\rm{\phi_{t}}}({\rm{\phi_{m}}})$, denoted together by
${\rm{\Omega_{t}}}({\rm{\Omega_{m}}})$. The measured direction denoted by
${\rm{\Omega_{m}}}$ is expected to be within a cone of half angle
$\Delta\theta$ of the true direction. The far detector (LArTPC) parameters
assumed are mentioned in Table 2[52].

Parameter Uncertainty | Value
---|---
$\mu^{+/-}$ $\Delta\theta$ | $2.5^{\circ}$
$e^{+/-}$ $\Delta\theta$ | $3.0^{\circ}$
($\mu^{+/-}$, $e^{+/-}$) Energy | GLB files for each E bin [1]
Detection efficiency | GLB files for each E bin [1]
Flux normalization | 20$\%$
Zenith angle dependence | 5$\%$
Cross section | 10$\%$
Overall systematic | 5$\%$
Tilt | 5$\%$
Table 2: Assumptions of the LArTPC far detector parameters and uncertainties.

####  IV.2.1 Charge identification using muon capture in Argon

The charge id of the muon can be identified using the capture vs decay process
of the muon inside the argon[53]. The working principle of charge id of the
muon is as follows: a fraction of the $\mu^{-}$ like events that undergo the
capture process are identified using capture fraction efficiency and the rest
of the muons as well as all the $\mu^{+}$ undergo muon decay. We have
implemented these as mentioned in [49].

###  IV.3 $\chi^{2}$ analysis

We have evaluated the $\chi^{2}$ for a fixed set of parameters using the
method of pulls. We can tackle various statistical and systematic
uncertainties directly through this method. The flux, cross sections, and
other systematic uncertainties are taken into account by allowing these inputs
to vary from their standard values in the computation of the expected rate in
the $\text{i-j}^{\text{th}}$ bin, ${\mathrm{N^{th}_{ij}}}$. Let the
${\mathrm{k^{th}}}$ input deviate from its standard value by
${\mathrm{\sigma_{k}\;\xi_{k}}}$ with ${\mathrm{\sigma_{k}}}$ denoting the
uncertainty. The value of ${\mathrm{N^{th}_{ij}}}$ with the modified inputs is
given by,

$\mathrm{N^{th}_{ij}=N^{th}_{ij}(std)+\sum^{npull}_{k=1}c_{ij}^{k}\xi_{k}}\,,$
(24)

where ${\mathrm{N^{th}_{ij}(std)}}$ gives the expected rate obtained with the
standard values of the inputs in the $\text{i-j}^{\text{th}}$ bin and npull is
the total number of sources of uncertainty (5 for our case). The
${\mathrm{\xi_{k}}}$’s are the _pull_ variables and they determine the number
of ${\mathrm{\sigma}}$’s by which the ${\mathrm{k^{th}}}$ input deviates from
its standard value. In eq. (24), ${\mathrm{c_{ij}^{k}}}$ signifies the change
in ${\mathrm{N^{th}_{ij}}}$ when the ${\mathrm{k^{th}}}$ input is changed by
${\mathrm{\sigma_{k}}}$ (i.e. by 1 standard deviation). The uncertainties in
the inputs being not very large, we only consider changes in
${\mathrm{N^{th}_{ij}}}$ that are linear in ${\mathrm{\xi_{k}}}$. Hence we
evaluate the modified $\chi^{2}$ as,

$\displaystyle{\chi^{2}(\xi_{k})}=$
$\displaystyle\sum_{i,j}\;\frac{\left[~{}N_{ij}^{th}(std)+\sum^{npull}_{k=1}c_{ij}^{k}\xi_{k}-N_{ij}^{ex}~{}\right]^{2}}{N_{ij}^{ex}}$
$\displaystyle+\sum^{npull}_{k=1}\xi_{k}^{2}\,,$ (25)

where the additional ${\rm{\xi_{k}^{2}}}$-dependent term is added due to the
penalty imposed for moving the value of the ${\mathrm{k^{th}}}$ input away
from its standard value by ${\rm{\sigma_{k}\;\xi_{k}}}$. $\chi^{2}$ in the
case of the standard LArTPC detector and for a detector with change id for
muons are evaluated as,

$\displaystyle\chi^{2}_{standard}=\chi^{2}_{\mu^{-}+\mu^{+}}+\chi^{2}_{e^{-}+e^{+}}$
(26) $\displaystyle\chi^{2}_{charge-
id}=\chi^{2}_{\mu^{-}}+\chi^{2}_{\mu^{+}}+\chi^{2}_{e^{-}+e^{+}}$ (27)

The $\chi^{2}$ with pulls, which includes the effects of all theoretical and
systematic uncertainties, is obtained by minimizing ${\rm{\chi^{2}(\xi_{k})}}$
with respect to all the pulls ${\rm{\xi_{k}}}$.

${\mathrm{\chi^{2}_{pull}=Min_{\xi_{k}}\left[\chi^{2}(\xi_{k})\right]}}\,.$
(28)

Finally, We marginalize the $\chi^{2}$ over the allowed range of the
oscillation parameters as mentioned in Table 3. For the combined analysis we
add the chi-square for beam and atmospheric and then marginalize over the
oscillation parameters. The marginalization has been performed in
$\theta_{23},\theta_{14},\theta_{24},\delta_{13},\delta_{14}$ over the range
specified in Table 3 for all cases unless otherwise mentioned.

Parameter | True Values | Marginalization Range
---|---|---
$\theta_{12}$ | $33.47^{\circ}$ | N.A.
$\theta_{13}$ | $8.54^{\circ}$ | N.A.
$\theta_{23}$ | $45^{\circ}$ | $39^{\circ}:51^{\circ}$
$\theta_{14},\theta_{24}$ | $7^{\circ}$ | $0^{\circ}:10^{\circ}$
$\theta_{14},\theta_{24}$ | $4^{\circ}$ | $0^{\circ}:6^{\circ}$
$\theta_{34}$ | $0^{\circ},7^{\circ}$ | $0^{\circ}:17^{\circ}$
$\Delta_{21}$ |  $7.42\times 10^{-5}$ eV2 | N.A.
$\Delta_{31}$(NO) |  $2.5\times 10^{-3}$ eV2 |  $-(2.42:2.62)\times 10^{-3}$ eV2
$\Delta_{31}$(IO) |  $-2.5\times 10^{-3}$ eV2 |  $(2.42:2.62)\times 10^{-3}$ eV2
$\Delta_{41}$ (for AMO) |  $1$ eV2 | N.A.
$\Delta_{41}$ (for SMO) |  $0.0005:0.01$ eV2 |  $\pm 15\%$ of $-\Delta_{41}$
$\delta_{13}$ | many | $-180^{\circ}:180^{\circ}$
$\delta_{14}$ | $0^{\circ},90^{\circ},-90^{\circ}$ | $-180^{\circ}:180^{\circ}$
Table 3: The table depicts true values of all the parameters and their range
of marginalization as used in our analysis.

##  V Results and Discussion

In this section, we present the results for the analysis of beam, atmospheric,
and a combination of both data for the following cases,

  * •

determination of the sign of $\Delta_{31}$ (AMO) for $\Delta_{41}=1$ eV2

  * •

determination of the sign of $\Delta_{31}$ for $\Delta_{41}$ in the range of
$5\times 10^{-4}:0.1$ eV2

  * •

determination the sign of $\Delta_{41}$ (SMO) when it’s value lies in the
range of $5\times 10^{-4}:0.1$ eV2

###  V.1 Sensitivity to the AMO for $|\Delta_{41}|=1$ eV2

Figure 7: The sensitivity to the atmospheric mass ordering as a function of
true $\delta_{13}$ for various $\delta_{14}^{true}$ values at 1300 km baseline
considering normal (left), inverted (right) ordering.

In fig. 7, the sensitivity to the atmospheric mass ordering (AMO) is presented
as a function of $\delta_{13}^{true}$ in standard three flavor framework
(black) for normal (left) and inverted (right) ordering. We also present the
sensitivity in the presence of a sterile neutrino corresponding to SNO-NO
(left), and SNO-IO (right) for true values of $\delta_{14}=0^{\circ}$(blue),
$90^{\circ}$(green), $-90^{\circ}$(orange), $180^{\circ}$(red). We will call
$-180^{\circ}<\delta_{13}\leq 0^{\circ}$ as the lower half plane[LHP] and
$0^{\circ}<\delta_{13}\leq 180^{\circ}$ as the upper half plane[UHP]
throughout this section. The important points to be noted are,

  * •

The sensitivity decreases in the presence of a sterile neutrino compared to
the three flavor case.

  * •

The sensitivity for the sterile cases depends on the true values of
$\delta_{14},\delta_{13}$.

  * •

For $\theta_{14},\theta_{24}=4^{\circ}$, the sensitivity is higher than that
of $\theta_{14},\theta_{24}=7^{\circ}$ and also closer to the standard 3$\nu$
case. This is due to the fact that the smaller the sterile mixing angles are
the more the sensitivity of the sterile case gets closer to the results of the
standard $3\nu$ case.

  * •

For NO, in the LHP of true $\delta_{13}$ the highest sensitivity is observed
for $\delta_{14}^{true}=-90^{\circ}$, and the lowest sensitivity for
$\delta_{14}^{true}=90^{\circ}$ whereas in the UHP, this order gets reversed.

  * •

For IO, $\delta_{14}^{true}=0^{\circ}$ (blue) shows the lowest sensitivity in
the LHP of true $\delta_{13}$ and also the highest sensitivity in UHP.

Next, in fig. 8, the AMO sensitivity is shown as a function of true
$\delta_{13}$ corresponding to the analysis of only beam(red), only
atmospheric (blue) and a combination of both beam and atmospheric(green)
simulated data. The cases with charge identification in atmospheric only
(violet) and combined analysis (orange) are also depicted. The representative
sensitivity curves are obtained for $\theta_{14},\theta_{24}=7^{\circ}$,
$\delta_{14}=0^{\circ}$ corresponding to true hierarchy considered as normal
(left) and inverted (right). The observations from fig. 8 are following,

  * •

Sensitivity for atmospheric neutrinos doesn’t have significant dependence on
$\delta_{13}$.

  * •

Combined sensitivity of beam and atmospheric is greater than the sum of
individual sensitivities of these two, which demonstrates the synergy between
them.

  * •

We observe slightly higher sensitivity when we use partial charge
identification for atmospheric neutrinos. This also leads to higher
sensitivity for combined analysis.

Figure 8: Atmospheric mass ordering sensitivity as a function of
$\delta_{13}^{true}$ corresponding to the analysis of only beam (red), only
atmospheric (violet), combined atmospheric+beam (green) neutrinos for normal
(left) and inverted (right) hierarchy with 400 kt-yr exposure of LArTPC

In fig. 9, we present the effect of $\theta_{34}$ on the sensitivity to the
atmospheric mass ordering. In this plot, the sensitivity is shown as a
function of true $\delta_{13}$ for various combinations of true values of
$\theta_{34},\delta_{34}$. We consider for beam only analysis,
$\theta_{34}=0^{\circ}$, $\delta_{34}=0^{\circ}$(green dotted), and
$\theta_{34}=7^{\circ}$,$\delta_{34}=90^{\circ}$(red dotted) along with the
sensitivity curve for standard three flavors (black). We also plot the
sensitivity with combined beam and atmospheric analysis
$\theta_{34}=7^{\circ}$,$\delta_{34}=90^{\circ}$ (red-dashed). Other sterile
parameters are fixed as $\delta_{14}=90^{\circ}$,
$\theta_{14}=\theta_{24}=7^{\circ}$. The observations are as follows,

  * •

The sensitivity decreases as $\theta_{34}$ becomes higher.

  * •

In the combined analysis of beam and atmospheric data, the sensitivity is
higher than the beam analysis compensating for the decrease due to non-zero
$\theta_{34}$.

Figure 9: The sensitivity to mass ordering as a function of
$\delta_{13}^{true}$ for various $\theta_{34}^{true}$ values for beam
neutrinos with 1300 km baseline for normal (left), inverted (right) ordering

###  V.2 Sensitivity of the sign of $\Delta_{31}$ for
$\Delta_{41}=10^{-4}:10^{-1}$ eV2

In this section, we study how the sensitivity to the sign of $\Delta_{31}$
behaves with $\Delta_{41}$ where the latter varies in the range of
$10^{-4}:10^{-1}$ eV2. Note that for $\Delta_{41}\sim 1$ eV2 only SNO-NO, SNO-
IO cases are cosmologically allowed. However, for
$\Delta_{41}=10^{-4}:10^{-1}$ eV2 all the four possibilities depicted in fig.
1 are admissible.

Figure 10: Sensitivity to atmospheric mass ordering as a function of
$\delta_{13}^{true}$ in SNO(top), SIO (bottom) scenarios with true
$\Delta_{31}$ for different values of $\Delta_{41}^{true}$ at 1300 km baseline

In fig. 10, the AMO sensitivity is shown as a function of $\delta_{13}^{true}$
at various true values of $\Delta_{41}$. The upper(lower) panels correspond to
the true value in SNO (SIO) cases while the left (right) panels are for
NO(IO). During the computation of $\chi^{2}$, the $|\Delta_{41}|$ is fixed in
true and test cases for this plot. The observations of significance in fig. 10
are as follows,

  * •

The nature of variation of sensitivity with $\delta_{13}^{true}$ doesn’t
change significantly for different true values of $\Delta_{41}$.

  * •

Sensitivity gets notably reduced for $\Delta_{41}=0.001$ eV2(blue) for the
most of values of $\delta_{13}^{true}$ in the UHP in SNO-NO and SIO-IO case.
In SNO-IO and SIO-NO cases, blue curves give minimum sensitivity over the full
range of $\delta_{13}$

  * •

However, sensitivity for $\Delta_{41}=0.001$ eV2 is very high in the LHP of
$\delta_{13}^{true}$ in SNO-NO, SIO-IO.

  * •

The maximum sensitivity is observed for $\Delta_{41}=0.01$ eV2(violet) in SIO-
NO and SIO-IO case over almost the full range of $\delta_{13}^{true}$.

Figure 11: Sensitivity to atmospheric mass ordering as a function of
$\Delta_{41}^{true}$ with marginalisation in $\Delta_{31}$ for
$\theta_{23}^{true}=45^{\circ}$, $\delta_{14}^{true}=0^{\circ}$,
$\delta_{13}^{true}=-90^{\circ}$(red), $90^{\circ}$(blue) at 1300 km baseline.
The violet curve shows $\Delta P_{14}^{s}$ at 2.5 GeV.

We note that for the difference in the probability between the opposite AMO,
the dependence on $\Delta_{41}$ will come from the last term in eq. (5). The
probability difference for fixed values of
$\theta_{13},\theta_{23},\theta_{14},\theta_{24},\delta_{13},\delta_{14}$ and
a given SMO can be expressed as,

$\displaystyle\Delta P_{\mu e}^{st}$ $\displaystyle=P_{\mu
e}^{m,true}(i_{1}\Delta_{41},j_{1}\Delta_{31})-P_{\mu
e}^{m,test}(i_{2}\Delta_{41},j_{2}\Delta_{31})$ (29) $\displaystyle\propto
4s_{13}s_{14}s_{24}s_{23}[\Delta P_{14}^{s}\sin{\tilde{\delta}_{14}}+\Delta
P_{14}^{c}\cos{\tilde{\delta}_{14}]}\text{,}$

where we define;

$\displaystyle\Delta P_{14}^{s,c}$
$\displaystyle=\frac{\sin[(A^{\prime}-1)\Delta]}{(A^{\prime}-1)}P_{14}^{s,c}(i_{1}\Delta_{41},j_{1}\Delta_{31})$
(30)
$\displaystyle-\frac{\sin[(A^{\prime}+1)\Delta]}{-(A^{\prime}+1)}P_{14}^{s,c}(i_{2}\Delta_{41},j_{2}\Delta_{31})$

Depending upon the true value considered in a certain mass ordering, there can
be four scenarios as follows

  * •

SNO-NO: $i_{1}=i_{2}=$+ve, $j_{1}=$+ve and $j_{2}=$-ve

  * •

SNO-IO: $i_{1}=i_{2}=$+ve, $j_{1}=$-ve and $j_{2}=$+ve

  * •

SIO-NO: $i_{1}=i_{2}=$-ve, $j_{1}=$+ve and $j_{2}=$-ve

  * •

SIO-IO: $i_{1}=i_{2}=$+ve, $j_{1}=$+ve and $j_{2}=$-ve

Note that $\Delta P_{14}^{s,c}$ for different cases are connected as,

$\displaystyle\Delta P_{14}^{s,c}|_{\rm SNO-NO}$ $\displaystyle=-\Delta
P_{14}^{s,c}|_{\rm SNO-IO},$ (31) $\displaystyle\Delta P_{14}^{s,c}|_{\rm SIO-
NO}$ $\displaystyle=-\Delta P_{14}^{s,c}|_{\rm SIO-IO}$ (32)

In fig. 11, we have depicted the sensitivity to AMO as a function of
$|\Delta_{41}|$ with marginalization performed only over $\Delta_{31}$ with
all other parameters being fixed. The red (blue) curve refers to
$\delta_{13}=-90^{\circ}(90^{\circ})$. We also show the difference in
probability term $\Delta P_{14}^{s}$(30) evaluated at 2.5 GeV by the violet
curve. The understandings from fig. 11 are as follows,

  * •

Since we have chosen $\delta_{14}=0^{\circ}$ for $\delta_{13}=90^{\circ}$ and
$-90^{\circ}$, the value of
$\sin\tilde{\delta}_{14}=\sin[\delta_{13}+\delta_{14}]$ is +1, -1 respectively
and $\cos\tilde{\delta}_{14}=0$. For the fixed phases and mixing angles in
both true and test cases, the difference in probability(29) between NO and IO
will only depend on $\Delta P_{14}^{s}$ as,

$\Delta P_{\mu e}^{st}=4s_{13}s_{14}s_{24}s_{23}\Delta
P_{14}^{s}\sin\tilde{\delta}_{14}$ (33)

This means $\Delta P_{\mu e}^{st}$ will be opposite for
$\delta_{13}=90^{\circ}$ and $-90^{\circ}$ leading to the opposite nature of
chi-square.

  * •

In the case of SNO-NO, for $\delta_{13}=90^{\circ}$, we have $\Delta P_{\mu
e}^{st}\propto\Delta P_{14}^{s}|_{\rm SNO-NO}$. This can be seen from the top-
left panel which shows that the nature of the blue and violet curves are
similar. For $\delta_{13}=-90^{\circ}$, $\Delta P_{\mu e}^{st}\propto-\Delta
P_{14}^{s}|_{\rm SNO-NO}$ which is reflected in the red curve being opposite
to the violet curve.

  * •

The sensitivity for SNO-IO is just opposite in nature to SNO-NO as shown in
eq. (31). Therefore, in the top-right panel, the red curve for
$\delta_{13}=-90^{\circ}$ is similar to the violet curve here, and the blue
curve for $\delta_{13}=+90^{\circ}$ is the opposite of the violet curve.

  * •

In SIO-NO, and SIO-IO cases, we also observe similar patterns as SNO-NO and
SNO-IO, respectively.

  * •

The sensitivity is seemed to be almost constant for $R>>1$, i.e.,
$\Delta_{41}>>\Delta_{31}$. This can be understood analytically as follows.
The difference in $P_{\mu e}$ for fixed energy and phases will only depend on
$\Delta P_{14}^{s}$ and can be evaluated using eq. (30) for $\Delta_{41}>>1$
limit as,

$\displaystyle\Delta P_{14}^{s}$
$\displaystyle=\left(\frac{\sin[(A^{\prime}-1)\Delta]}{(A^{\prime}-1)}+\frac{\sin[(A^{\prime}+1)\Delta]}{(A^{\prime}+1)}\right)$
(34)
$\displaystyle\times\Big{\\{}(\frac{A^{\prime}}{2R}c_{23}+2)\sin[R\Delta]^{2}\Big{\\}},$

where the term in the braces is dependent on $\Delta_{41}$. For $R>>1$ the
term $\sin[R\Delta]$ shows fast oscillation. Summing over energies,
$\sin[R\Delta]$, will give a constant value that is reflected in the
sensitivity curves.

In fig. 12, we depict the sensitivity to the sign of $\Delta_{31}$ as a
function of true $\Delta_{41}$ for
$\delta_{13}=-90^{\circ}(\rm{red}),90^{\circ}(\rm{blue})$, and
$\delta_{14}=0^{\circ}$ at 1300 km. For this figure, we have marginalized over
the parameters as mentioned in table‘ 3.

Figure 12: Sensitivity to atmospheric mass ordering as a function of
$\Delta_{41}^{true}$ for $\theta_{23}^{true}=45^{\circ}$,
$\delta_{14}^{true}=0^{\circ}$, $\delta_{13}^{true}=-90^{\circ}$ (red),
$90^{\circ}$(blue) at 1300 km baseline

Some interesting features of sensitivity to the AMO as seen from fig. 12 are
as follows,

  * •

For SNO-NO (top-left panel) and SIO-IO (bottom-right panel) cases we observe a
contrasting nature of the sensitivity between $\delta_{13}=90^{\circ}$ and
$-90^{\circ}$. Note that a similar contrasting nature has been observed in the
top two panels of fig. 5 showing the oscillogram of $\Delta P_{\mu e}$. For
instance, in SNO-NO, at $\Delta_{41}=3\times 10^{-3}$ eV2 a maxima of
sensitivity occurs for $\delta_{13}=90^{\circ}$ whereas minima occurs for
$\delta_{13}=-90^{\circ}$.

  * •

However, the nature of the sensitivity curves for $\delta_{13}=90^{\circ}$ and
$-90^{\circ}$ is similar in SNO-IO (top-right panel) and SIO-NO (bottom-left
panel) cases. This behavior was also observed in $\Delta P_{\mu e}$
oscillogram in the bottom panels of fig. 5 for SNO-IO case.

  * •

In all the cases the minima and maxima are observed in the range of
$0.001-0.01$ eV2. Beyond that, the sensitivity is relatively flat with
$\Delta_{41}$.

###  V.3 Sensitivity to the sign of $\Delta_{41}$ (SMO)

In this section, we present the sensitivity of the sign of $\Delta_{41}$
considering it’s values in the range of $[5\times 10^{-4}:10^{-1}]$ eV2 for
which all four possibilities depicted in fig 1 will be viable. In fig. 13, the
sensitivity of the sign of sterile mass squared differences are depicted as a
function of the true value of $\Delta_{41}$ for various true values of
$\delta_{13},\delta_{14}$111This plot has been presented in ref [41] for fixed
values of $\delta_{13}=0^{\circ},\delta_{14}=0^{\circ}$.

Figure 13: Sensitivity to sterile mass ordering as a function of
$\Delta_{41}^{true}$ for $\Delta_{31}=+ve$(left), $-ve$(right) and
$\Delta_{41}=+ve$(top), $-ve$(bottom) using different values of
$\delta_{14}^{true}$, $\delta_{13}^{true}$ at 1300 km.

The salient features of the fig. 13 are as follows,

  * •

The sensitivity curves for SNO-NO and SIO-IO show two prominent maxima around
true values of $\Delta_{41}=1\times 10^{-3}$ eV2, $5\times 10^{-3}$ eV2 cases.
There is a dip in sensitivity when true $\Delta_{41}$ is around $2.5\times
10^{-3}$ eV2 due to its proximity to atmospheric mass squared difference as
discussed in ref. [41].

  * •

In the case of SNO-IO and SIO-NO, the maxima occurs around
$\Delta_{41}=2.5\times 10^{-3}$ eV2, i.e, when the sterile mass squared
difference is close to the atmospheric mass squared difference $\Delta_{31}$.

  * •

The SNO-IO and SIO-IO cases provide relatively higher sensitivity than the
SNO-NO and SIO-IO cases.

  * •

The features of sensitivity in SNO-NO can be qualitatively understood from the
plot of probability difference in top panels of fig. 6.

Figure 14: Sensitivity to sterile mass hierarchy as a function of
$\Delta_{41}^{true}$ using combined beam and atmospheric neutrinos at 1300 km
baseline.

In fig. 14, we have used the simulated data from atmospheric neutrino analysis
to perform a combined analysis of beam and atmospheric neutrinos and get the
sensitivity of the sterile mass ordering as a function of the true value of
$\Delta_{41}$ for the true value of phases $\delta_{13}=-90^{\circ}$(left),
$\delta_{14}=90^{\circ}$(right). In these four panels, it is observed that the
sensitivity to SMO is bettered in combined analysis than the beam neutrinos.
The nature of the sensitivity is almost similar for beam and atmospheric
analysis. This can be understood from the similar profile of difference in
probabilities $\Delta P_{\nu e},\Delta P_{\mu\mu}$ at 1300 km and 7000 km as
shown in fig. 5. The combined sensitivity is above $3\sigma$ for most of the
parameter space up to $\Delta_{41}=10^{-2}$ eV2. For the values of
$\Delta_{41}$ greater than that, even with the addition of atmospheric
neutrinos, we get fixed sensitivity of $1.5\sigma-2\sigma$.

##  VI Conclusions

Our work focuses on the effect of an additional light sterile neutrino with a
mass squared difference in the range of $10^{-4}:1$ eV2 on the determination
of atmospheric mass ordering and sterile mass ordering. For our study, we use
a liquid argon detector for a) beam neutrinos (baseline 1300 km), and b)
atmospheric neutrinos.

The main new aspect of our study is to do a combined analysis of beam and
atmospheric neutrinos to determine the sign of $\Delta_{31}$ and $\Delta_{41}$
when the latter is varied in the range $10^{-4}-1$ eV2. We also study in
detail the impact of various sterile parameters.

The sensitivity to the sign of $\Delta_{31}$(AMO) for $\Delta_{41}=1$ eV2 gets
diminished w.r.t. to the 3$\nu$ case in the presence of a sterile neutrino.
The decrement in the sensitivity is higher for larger values of
$\theta_{14},\theta_{24}$. The values of $\chi^{2}$ also depend on
$\delta_{13},\delta_{14}$. The sensitivity to AMO decreases further in the
presence of a non-zero $\theta_{34}$. However, with the combined analysis of
beam and atmospheric neutrino, we are able to recover the sensitivity over
10$\sigma$, irrespective of the choice of true values of
$\delta_{13},\delta_{14}$. The presence of a light sterile neutrino gives the
possibility of both positive and negative values of $\Delta_{41}$. Our study
also demonstrates for the first time the dependence of the atmospheric mass
ordering sensitivity on the absolute value of $\Delta_{41}$ as well as the on
the nature of the 3+1 mass spectrum (SNO-NO, SNO-IO, SIO-NO, SIO-IO).

We also study the sensitivity to the sign of $\Delta_{41}$ (SMO) for the
$\Delta_{41}=5\times 10^{-4}:0.1$ eV2 in different scenarios of the 3+1 mass
spectrum. The sensitivity gets reduced when $\Delta_{41}$ is in the proximity
of $\Delta_{31}$ for the SNO-NO and SIO-IO cases whereas in the SNO-IO and
SIO-NO cases, the sensitivity gets enhanced. The addition of atmospheric
neutrinos boosts the sensitivity over 3$\sigma$ for $\Delta_{41}<10^{-2}$ eV2.
However, for higher values of $\Delta_{41}$, the sensitivity falls off to
$\sim 1.5\sigma-2\sigma$ for the combined analysis.

###### Acknowledgements.

AC acknowledges the University of Pittsburgh, where the initial phase of the
work has been done. AC also acknowledges the Ramanujan Fellowship
(RJF/2021/000157), of the Science and Engineering Research Board of the
Department of Science and Technology, Government of India. AC also
acknowledges CERN for providing facilities. SG acknowledges the J.C. Bose
Fellowship (JCB/2020/000011) of the Science and Engineering Research Board of
the Department of Science and Technology, Government of India. It is to be
noted that this work has been done solely by the authors and is not
representative of the DUNE collaboration.

## References

  * Alion  _et al._ [2016] T. Alion _et al._ (DUNE), (2016), arXiv:1606.09550 [physics.ins-det] .
  * Abi  _et al._ [2020] B. Abi _et al._ (DUNE), (2020), arXiv:2002.03005 [hep-ex] .
  * Hyp [2016] (2016).
  * Alekou  _et al._ [2022] A. Alekou _et al._ ,  (2022), 10.1140/epjs/s11734-022-00664-w, arXiv:2206.01208 [hep-ex] .
  * Abe  _et al._ [2020] K. Abe _et al._ (T2K), Nature 580, 339 (2020), [Erratum: Nature 583, E16 (2020)], arXiv:1910.03887 [hep-ex] .
  * Ayres  _et al._ [2007] D. S. Ayres _et al._ (NOvA),  (2007), 10.2172/935497.
  * Abe  _et al._ [2015] K. Abe _et al._ (Hyper-Kamiokande Proto-), PTEP 2015, 053C02 (2015), arXiv:1502.05199 [hep-ex] .
  * Ahmed  _et al._ [2017] S. Ahmed _et al._ (ICAL), Pramana 88, 79 (2017), arXiv:1505.07380 [physics.ins-det] .
  * Aartsen  _et al._ [2016] M. G. Aartsen _et al._ (IceCube Collaboration), Phys. Rev. Lett. 117, 071801 (2016).
  * Barger  _et al._ [2014] V. Barger, A. Bhattacharya, A. Chatterjee, R. Gandhi, D. Marfatia, and M. Masud, Phys. Rev. D 89, 011302 (2014), arXiv:1307.2519 [hep-ph] .
  * Kelly  _et al._ [2022] K. J. Kelly, P. A. N. Machado, I. Martinez-Soler, and Y. F. Perez-Gonzalez, JHEP 05, 187 (2022), arXiv:2110.00003 [hep-ph] .
  * Argüelles  _et al._ [2022] C. A. Argüelles, P. Fernández, I. Martínez-Soler, and M. Jin, (2022), arXiv:2211.02666 [hep-ph] .
  * Abi  _et al._ [2021] B. Abi _et al._ (DUNE), Eur. Phys. J. C 81, 322 (2021), arXiv:2008.12769 [hep-ex] .
  * Argüelles  _et al._ [2020] C. A. Argüelles _et al._ , Rept. Prog. Phys. 83, 124201 (2020), arXiv:1907.08311 [hep-ph] .
  * Aguilar  _et al._ [2001] A. Aguilar _et al._ (LSND), Phys. Rev. D 64, 112007 (2001), arXiv:hep-ex/0104049 .
  * Armbruster  _et al._ [2002] B. Armbruster _et al._ (KARMEN), Phys. Rev. D 65, 112001 (2002), arXiv:hep-ex/0203021 .
  * Aguilar-Arevalo  _et al._ [2009] A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Nucl. Instrum. Meth. A 599, 28 (2009), arXiv:0806.4201 [hep-ex] .
  * Aguilar-Arevalo  _et al._ [2007] A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. Lett. 98, 231801 (2007), arXiv:0704.1500 [hep-ex] .
  * Hampel  _et al._ [1998] W. Hampel _et al._ (GALLEX), Phys. Lett. B 420, 114 (1998).
  * Abdurashitov  _et al._ [1999] J. N. Abdurashitov _et al._ (SAGE), Phys. Rev. C 59, 2246 (1999), arXiv:hep-ph/9803418 .
  * Barinov  _et al._ [2022] V. V. Barinov _et al._ , Phys. Rev. C 105, 065502 (2022), arXiv:2201.07364 [nucl-ex] .
  * Abazajian  _et al._ [2012a] K. N. Abazajian _et al._ , (2012a), arXiv:1204.5379 [hep-ph] .
  * Minkowski [1977] P. Minkowski, Phys. Lett. B 67, 421 (1977).
  * Mohapatra and Senjanovic [1980] R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).
  * Schechter and Valle [1980] J. Schechter and J. W. F. Valle, Phys. Rev. D 22, 2227 (1980).
  * Fukugita and Yanagida [1986] M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45 (1986).
  * Davidson  _et al._ [2008] S. Davidson, E. Nardi, and Y. Nir, Phys. Rept. 466, 105 (2008), arXiv:0802.2962 [hep-ph] .
  * Appelquist and Shrock [2002] T. Appelquist and R. Shrock, Phys. Lett. B 548, 204 (2002), arXiv:hep-ph/0204141 .
  * Appelquist and Shrock [2003] T. Appelquist and R. Shrock, Phys. Rev. Lett. 90, 201801 (2003), arXiv:hep-ph/0301108 .
  * Dodelson and Widrow [1994] S. Dodelson and L. M. Widrow, Phys. Rev. Lett. 72, 17 (1994), arXiv:hep-ph/9303287 .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Riess  _et al._ [2018] A. G. Riess, S. Casertano, W. Yuan, L. Macri, J. Anderson, J. W. MacKenty, J. B. Bowers, K. I. Clubb, A. V. Filippenko, D. O. Jones, and B. E. Tucker, The Astrophysical Journal 855, 136 (2018).
  * Scolnic  _et al._ [2018] D. M. Scolnic _et al._ (Pan-STARRS1), Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Foot and Volkas [1995] R. Foot and R. R. Volkas, Phys. Rev. Lett. 75, 4350 (1995), arXiv:hep-ph/9508275 .
  * Dasgupta and Kopp [2014] B. Dasgupta and J. Kopp, Phys. Rev. Lett. 112, 031803 (2014), arXiv:1310.6337 [hep-ph] .
  * Chu  _et al._ [2018] X. Chu, B. Dasgupta, M. Dentler, J. Kopp, and N. Saviano, JCAP 11, 049 (2018), arXiv:1806.10629 [hep-ph] .
  * Yaguna [2007] C. E. Yaguna, JHEP 06, 002 (2007), arXiv:0706.0178 [hep-ph] .
  * de Holanda and Smirnov [2011] P. C. de Holanda and A. Y. Smirnov, Phys. Rev. D 83, 113011 (2011), arXiv:1012.5627 [hep-ph] .
  * de Gouvêa  _et al._ [2022] A. de Gouvêa, G. Jusino Sánchez, and K. J. Kelly, Phys. Rev. D 106, 055025 (2022), arXiv:2204.09130 [hep-ph] .
  * Thakore  _et al._ [2018] T. Thakore, M. M. Devi, S. Kumar Agarwalla, and A. Dighe, JHEP 08, 022 (2018), arXiv:1804.09613 [hep-ph] .
  * Chattopadhyay  _et al._ [2023] D. S. Chattopadhyay, M. M. Devi, A. Dighe, D. Dutta, D. Pramanik, and S. K. Raut, JHEP 02, 044 (2023), arXiv:2211.03473 [hep-ph] .
  * Goswami [1997] S. Goswami, Phys. Rev. D 55, 2931 (1997), arXiv:hep-ph/9507212 .
  * Abazajian  _et al._ [2012b] K. N. Abazajian _et al._ , (2012b), arXiv:1204.5379 [hep-ph] .
  * Donini  _et al._ [2002] A. Donini, M. Lusignoli, and D. Meloni, Nucl. Phys. B 624, 405 (2002), arXiv:hep-ph/0107231 .
  * Meloni  _et al._ [2010] D. Meloni, J. Tang, and W. Winter, Phys. Rev. D 82, 093008 (2010), arXiv:1007.2419 [hep-ph] .
  * Kopp  _et al._ [2013] J. Kopp, P. A. N. Machado, M. Maltoni, and T. Schwetz, JHEP 05, 050 (2013), arXiv:1303.3011 [hep-ph] .
  * Gariazzo  _et al._ [2017] S. Gariazzo, C. Giunti, M. Laveder, and Y. F. Li, JHEP 06, 135 (2017), arXiv:1703.00860 [hep-ph] .
  * Adamson  _et al._ [2020] P. Adamson _et al._ (MINOS+, Daya Bay), Phys. Rev. Lett. 125, 071801 (2020), arXiv:2002.00301 [hep-ex] .
  * Chatterjee  _et al._ [2022] A. Chatterjee, S. Goswami, and S. Pan, (2022), arXiv:2212.02949 [hep-ph] .
  * Huber  _et al._ [2005] P. Huber, M. Lindner, and W. Winter, Comput. Phys. Commun. 167, 195 (2005), arXiv:hep-ph/0407333 [hep-ph] .
  * Honda  _et al._ [2015] M. Honda, M. S. Athar, T. Kajita, K. Kasahara, and S. Midorikawa, Phys. Rev. D 92, 023004 (2015).
  * Barger  _et al._ [2016] V. Barger, A. Bhattacharya, A. Chatterjee, R. Gandhi, D. Marfatia, and M. Masud, Int. J. Mod. Phys. A 31, 1650020 (2016), arXiv:1405.1054 [hep-ph] .
  * Ternes  _et al._ [2019] C. A. Ternes, S. Gariazzo, R. Hajjar, O. Mena, M. Sorel, and M. Tórtola, Phys. Rev. D 100, 093004 (2019).
  * Note [1] This plot has been presented in ref [41] for fixed values of $\delta_{13}=0^{\circ},\delta_{14}=0^{\circ}$.