aainstitutetext: Department of Physics and Astronomy, University of
Pittsburgh, Pittsburgh, PA 15260bbinstitutetext: Physical Research Laboratory,
Ahmedabad, Gujarat, 380009, Indiaccinstitutetext: Indian Institute of
Technology, Gandhinagar, Gujarat, 382355, India

# Matter effect in presence of a sterile neutrino and resolving the octant
degeneracy using liquid argon detector in DUNE

Animesh Chatterjee,111Corresponding author. b Srubabati Goswami b,c,1 Supriya
Pan<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Results from the 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 new sterile parameters. In our work, we investigate how the sensitivity to
determine the octant of $\theta_{23}$ is affected by introducing a sterile
neutrino to the standard neutrino oscillation framework. We compute the
oscillation probabilities in presence of a sterile neutrino, analytically,
using the approximation that $\Delta_{21}$, the smallest mass squared
difference, is zero. We use these probabilities to understand analytically the
degeneracies at different baselines. In our work, we investigate in the
context of the DUNE experiment, the sensitivity to determine the octant of
$\theta_{23}$ for 3+1 framework. We present our results of the sensitivity to
octant of $\theta_{23}$ for beam neutrinos using a liquid argon time
projection chamber (LArTPC) at the DUNE far detector. We also obtain octant
sensitivity using atmospheric neutrinos at DUNE using a LArTPC detector
without any charge identification capability. In addition, we obtain the
octant sensitivity including the charge tagging capability of muon capture in
argon which allows differentiating between muon neutrino and antineutrino
events. The combined sensitivity of beam and atmospheric neutrinos in DUNE is
also delineated. We observe that by combining simulated data from the beam and
atmospheric neutrinos (including charge-id for muons), the sensitivity to the
octant of $\theta_{23}$ for true values of
$\theta_{23}=41^{\circ}(49^{\circ})$ exceeds $4\sigma(3\sigma)$ for more than
$50\%$ values of true $\delta_{13}$.

##  1 Introduction

Coming a long way since the first observation of neutrino oscillations, the
standard three flavour neutrino oscillation paradigm is well established now
with most of the parameters being measured with considerable precisionEsteban
_et al._ (2017); de Salas  _et al._ (2021); Capozzi  _et al._ (2020). The
parameters describing the standard three flavour oscillations are the three
mixing angles $\theta_{12},\theta_{13},\theta_{23}$ corresponding to mixing
between the mass eigenstates with mass eigenvalues $m_{1},m_{2},m_{3}$, the
Dirac CP phase $\delta_{13}(\delta_{CP})$, the two mass squared differences
$\Delta_{21}=m_{2}^{2}-m_{1}^{2}$ driving the solar neutrino transitions and
$\Delta_{31}=m_{3}^{2}-m_{1}^{2}$ governing the atmospheric neutrino
oscillations. The oscillation probabilities are also dependent on the
experimental parameters like the baseline distance $L$ traversed by the
neutrinos from the source to the detector and the energy of the neutrinos
$E_{\nu}$. Currently, the unknowns in the standard oscillation sector are the
mass ordering among the three neutrino states, the octant of the atmospheric
mixing angle $\theta_{23}$, and the value of the CP violating phase
$\delta_{CP}$. The mass ordering refers to whether the sign of the atmospheric
mass squared difference $\Delta_{31}$ is positive (Normal Ordering/NH) or
negative (Inverted Ordering/IH). The octant of $\theta_{23}$ signifies if the
value of the angle lies above (Higher Octant/HO) or below (Lower Octant/LO)
$45^{\circ}$. One of the most impeding factors in the precise determination of
these three parameters is the occurrence of degeneracies, i.e., various sets
of different values of unknown parameters giving rise to equal probability
making an unambiguous determination of these parameters difficult. In view of
the current unknowns in the three flavour framework, the existing degeneracies
can be understood through a generalized hierarchy-octant-$\delta_{CP}$
degeneracyGhosh  _et al._ (2016a). Measuring these parameters with
considerable precision and alleviating the existing degeneracies are the focus
of the ongoing and upcoming experiments such as T2KAbe  _et al._ (2020),
NO$\nu$AAyres  _et al._ (2007), T2HK, DUNE Abi  _et al._ (2020a),
ESS$\nu$SBTolba (2022), etc. Planned atmospheric neutrino experiments like
HyperKamiokandeAbe  _et al._ (2015), KM3NeTBrunner (2018), PINGUAartsen  _et
al._ (2017), INOAhmed  _et al._ (2017), etc can also throw light on these
parameters. Synergy and complementarity between atmospheric and beam
experiments have been explored in the context of three generation framework in
Ghosh  _et al._ (2016b); Barger  _et al._ (2014).

Although neutrino oscillation is unequivocally the leading solution for the
flavour conversion to explain the solar, atmospheric, reactor, and accelerator
observations, the possibility of other beyond standard model (BSM) effects at
a sub-leading level cannot be precluded. Several new physics effects have been
discussed in the literature including sterile neutrinos, long rang forces,
neutrino decay, non-standard interactions, etc. Among these BSM scenarios, the
existence of a light sterile neutrino is one of the most promising new physics
hypotheses to explain the experimental anomalies observed by LSND and
MiniBooNE experiments.

The idea of the light sterile neutrino was initiated in view of the results
from the LSND experiment which reported a signature of
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ oscillation at $3.8\sigma$Hill
(1995). Later MiniBooNE experiment with the same $L/E$ ratio also confirmed
this result at $4.8\sigma$ significanceAguilar-Arevalo  _et al._ (2020). If we
want to interpret these results through effective two flavour oscillations
then corresponding to $L/E\sim 1$ GeV/km in these experiments there should be
a new mass-squared difference $\Delta_{s}\sim 1$ eV2. This new mass-squared
difference doesn’t fit in the standard three flavour oscillation scheme
requiring one to incorporate at least one extra neutrino with a mass of eV
scale. Result of the invisible decay width of the $Z$ boson at CERN suggests
that there can only be three different neutrinos below the mass range of half
of the Z bosonAbrams  _et al._ (1989). Hence, this extra neutrino has to be
inert with no standard model gauge interactions.

Additional support in favor of an additional light sterile neutrino came from
the observation of electron neutrino $\nu_{e}$ deficit in gallium-based radio-
chemical experiments SAGE and GALLEX (Gallium anomaly)Giunti and Laveder
(2011); Acero  _et al._ (2009) which has been reinforced with the recent
results from BEST experimentBarinov  _et al._ (2021) at $5\sigma$. There was
also the reactor antineutrino anomaly in which several reactor neutrino
experiments showed a deficit in the measured flux with an improved calculation
of the inverse beta decay cross-section Huber (2011); Mueller  _et al._
(2011). These could also be explained in terms of a sterile neutrino with a
mass of the order of eV. However, the results from reactor experiments such as
DANSSSvirida  _et al._ (2019); Danilov and Skrobova (2021), NEOSKo  _et al._
(2017), STEREOAlmazán  _et al._ (2020), and PROSPECTAndriamirado  _et al._
(2021) excluded most of the reactor antineutrino anomaly regionMinotti (2022)
at more than $90\%$ C.L. So far, the 3+1 framework including a light sterile
neutrino with a mass of 1 eV first introduced in Goswami (1997), offers the
most economical scenario to explain these anomalies. However, the 3+1
framework suffers from a tension between the $\nu_{\mu}$ disappearance and
appearance data. This tension originates from the non-observation of any
similar supportive signal in the accelerator based disappearance experiments
in $P_{\mu\mu}$ channel like CDHSW, MINOSAbazajian  _et al._ (2012); Adamson
_et al._ (2019), Super-KamiokandeAbe  _et al._ (2019), IceCube
DeepCoreVannerom  _et al._ (2022), MicroBooNEArgüelles  _et al._ (2022),
NO$\nu$AHausner (2022). Reactor based electron disappearance searches in the
experiments Bugey3Hu (2021) and Daya BayAdamson  _et al._ (2020) also didn’t
provide any evidence in support of sterile neutrino. Most recent results from
the MicroBooNE experiment did not report any evidence of electron neutrino
disappearance in their three years of dataMicroBooNE Collaboration  _et al._
(2021); Abratenko  _et al._ (2022). However, it was shown in Argüelles  _et
al._ (2021) that MicroBooNE data can not exclude the electron neutrino excess
observed in MiniBooNE in a model independent way. The joint analysis of
results from MiniBooNE, and MicroBooNE experiments preferred the 3+1 scenario
over no oscillationAguilar  _et al._ (2022). On the other hand, a recent study
by MicroBooNE performing a joint analysis of $\nu_{\mu}\rightarrow\nu_{e}$ and
$\nu_{e}\rightarrow\nu_{e}$ data finds that the three neutrino picture gives a
better fit as compared to the three active and one sterile scenario Abratenko
_et al._ (2022).

One of the sternest challenges for the existence of sterile neutrino comes
from cosmologyBöser  _et al._ (2020). The inclusion of an extra sterile
neutrino increases the effective no of neutrinos relevant for the Big Bang
Nucleosynthesis. It was proposed in Chu  _et al._ (2018) that secret
interaction between sterile neutrinos can remedy this situation but it was
later disfavoured by cosmic microwave background analysis Goswami  _et al._
(2022). Recently, a joint analysis of short baseline and cosmological data
showed that a sterile neutrino with a mass around 1 eV can be allowed for
interaction with a new light pseudo scalar. To summarize, the existence of
sterile neutrinos is still an open question, and more experimental efforts are
underway to resolve this. The upcoming TRISTAN detector at the KATRINAker  _et
al._ (2022) experiment, SBNAntonello  _et al._ (2015) at Fermilab, JSNS2
detectorAjimura  _et al._ (2021) at J-PARC in JAPAN are following up the
results of LSND, MiniBooNE and the results from these experiments are expected
to offer a clear indication on the existence of a eV scale sterile neutrino.

If we consider the sterile neutrino hypothesis to be true, then the standard
framework of neutrino oscillations is going to see some important
modifications. The addition of a light sterile neutrino comes with three extra
active-sterile sector mixing angles and two additional CP phases. These will
compound the effect of the parameter degeneracies which already exist in the
standard three flavour framework. In particular, it was shown in Agarwalla
_et al._ (2017) that for the 3+1 oscillation framework the octant degeneracy
is more pronounced due to the effect of an additional interference term in the
$\nu_{\mu}\rightarrow\nu_{e}$ vacuum oscillation probability relevant at long
baseline setups in the context of the DUNE detector. It is well known that the
addition of neutrino and anti-neutrino can evade the octant-$\delta_{13}$
degeneracy for three flavour caseKumar Agarwalla  _et al._ (2013); Ghosh  _et
al._ (2016a). However, in presence of a sterile neutrino, the
octant-$\delta_{14}$ degeneracy can’t be resolved even after the addition of
neutrino plus anti-neutrinoGhosh  _et al._ (2018). Implications of additional
octant degeneracies associated with the new phases in presence of a sterile
neutrino have also been studied in the context of the NO$\nu$AChoubey  _et
al._ (2019); Ghosh  _et al._ (2018) experiment. Other studies in the context
of long baseline experiments in presence of a sterile neutrino can be found
for instance in Dutta  _et al._ (2016).

The primary focus of our paper is to study the octant sensitivity if an
additional light sterile neutrino is present. We perform a comprehensive study
of the octant sensitivity in the DUNE experiment using a Liquid Argon Time
Projection Chamber(LArTPC) detector. Apart from the beam neutrinos at DUNE
which are detected at 1300 km baseline at the far detector, we also consider
the atmospheric neutrinos which provide large baselines as well as higher
energies. We aim to investigate if the effect of additional degeneracies
arising due to an extra light sterile neutrino can be reduced in presence of a
large matter effect encountered at higher baselines. The octant sensitivity
has also been studied for the combined analysis of beam neutrinos 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.

In order to have a proper understanding of the octant degeneracy seen from
numerical analysis, the study of the analytic expressions of neutrino
oscillation probabilities is important. We obtain analytic expressions of the
neutrino oscillation/survival probabilities assuming the solar mass squared
difference $\Delta_{21}$ to be negligible as compared to the mass squared
differences $\Delta_{31}$, and $\Delta_{41}=m_{4}^{2}-m_{1}^{2}$ driving the
atmospheric and sterile neutrino oscillations respectively. We use the
analytic expressions to understand the octant degeneracy both at 1300 km and
at a higher atmospheric baseline of 7000 km. There have been other analytical
calculations of neutrino oscillation probabilities in presence of sterile
neutrino in matter using the rotation methodsParke and Zhang (2020a), an exact
analytical methodLi  _et al._ (2018), but our calculations provide a compact
form to explain the octant sensitivity from an analytical point of view. We
discuss the region of validity and the error of the analytic expressions
compared to the numerical probabilities.

Studies related to sterile neutrinos in the context of atmospheric neutrino
observations at India-based Neutrino Observatory (INO) experiment have been
performed in Behera  _et al._ (2017); Gandhi and Ghoshal (2012). More recently
an analysis in Thakore  _et al._ (2018) consider sterile neutrinos in
atmospheric baselines for a wide $\Delta_{41}$ mass squared range
$10^{-5}:100$ eV2 in the context of the INO experiment. This paper obtained
bound on the active-sterile mixing angles as well as the sensitivity to the
neutrino mass ordering in the 3+1 oscillation framework. Our study in this
paper focuses on the impact of resonant matter effect on the probabilities at
very long baselines and its influence on the sensitivity to determine the
octant. To the best of our knowledge, this kind of study of the degeneracies
in presence of a light sterile neutrino under the influence of resonance
matter effect at very long baselines has not been looked into in past. We also
explore this aspect in the context of atmospheric and beam neutrinos both
separately as well as together using a LArTPC detector in the DUNE experiment
and examine the complementarities between these two. Such studies in the
context of a generic liquid argon detector have been performed for the three
generation case earlier in Chatterjee  _et al._ (2013).

The plan of the paper is as follows. To start with we establish the analytic
framework for neutrino oscillations in presence of sterile neutrino in section
2. The subsequent section 3 details the calculation of the probabilities.
Next, section 4 contains the discussion on octant degeneracy for different
baselines and energies as well as the dependence on the CP phases
$\delta_{13}$, and $\delta_{14}$. In section 5, we describe the experimental
details for the DUNE detector and outline the procedure of $\chi^{2}$ analysis
adopted. We discuss the results in section 6. Finally, we conclude in section
7.

##  2 3+1 Framework

The minimal scheme postulated to explain the results of LSND, and MiniBooNE is
the 3+1 framework in which one light sterile neutrino is added to the three
active neutrinos in the SM. In the 3+1 oscillation framework, the mixing
matrix $U$ depends on three additional mixing angles
$\theta_{14},\theta_{24},\theta_{34}$ corresponding to mixing between the
light sterile neutrino $\nu_{s}$ and the active sector neutrinos, two new CP
phases $\delta_{14},\delta_{34}$ along with the standard oscillation
parameters $\theta_{12},\theta_{13},\theta_{23},\delta_{13}$ and can be
expressed as,

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

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 the rotation 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. In the presence of an additional light sterile neutrino, there
is a new mass squared difference $\Delta_{41}$. The 3+1 picture considered
here is $m_{4}>>m_{3}>>m_{2}>>m_{1}$ which corresponds to $m_{4}$ being the
heaviest mass state. The case with $m_{4}$ as the lowest state is disfavoured
from cosmology. The mass ordering for three generation is considered to be NH.

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

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

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,
$G_{F}$ is the Fermi coupling constant, $N_{e}$, and $N_{n}$ correspond to
electron density and neutron density, respectively, of the medium in which
neutrinos travel. In presence of matter, the total Hamiltonian in flavour
basis is 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}$ (3)

where the propagation medium has been considered to be the earth 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}$ with neutrino energy
$E_{\nu}$ and the mass squared differences are given as
$\Delta_{ij}=m_{j}^{2}-m_{i}^{2}$ where $m_{i}$’s are mass eigenvalues.

##  3 Oscillation Probability

In order to obtain the probabilities in the matter, one has to solve the
neutrino propagation equation with the total Hamiltonian given by (3). This
would require diagonalization of the total Hamiltonian to go to the matter
mass basis. However, this possesses difficulty even in the three flavour case
and one has to resort to approximate methods. A comprehensive review of the
various approximations used in the three flavour case has been discussed
recently in Parke and Zhang (2020b). In the context of this work, we have
considered the two mass scale dominance (TMSD) approximation with
$\Delta_{21}$ set as zero similar to the OMSD approximation in three flavour
framework. This approximation allows us to obtain compact analytic expressions
for the probabilities which can facilitate the understanding of the underlying
physics.

###  3.1 TMSD approximation

In the TMSD approximation we choose $\Delta_{21}=0$ since from the
experimental data $\Delta_{21}<<\Delta_{31}<<\Delta_{41}$. As a consequence,
the contribution of the solar angle $\theta_{12}$ drops out of mixing matrix
$U$ (1) as $R_{12}$ commutes with the mass matrix $M$ in the approximation
$\Delta_{21}=0$. This approximation holds well for
$\frac{\Delta_{21}L}{E_{\nu}}<1$. In our study, we further assume
$\theta_{34}=0$ which is allowed within current bounds Capozzi  _et al._
(2020); Esteban  _et al._ (2017). Thus we have only two additional non-zero
mixing angles $\theta_{14},\theta_{24}$ and a non-zero phase $\delta_{14}$.
This leads to the effective vacuum mixing matrix,

$\begin{split}\tilde{U}&=R_{24}(\theta_{24})\tilde{R}_{14}(\theta_{14},\delta_{14})R_{23}(\theta_{23})U_{\delta
13}R_{13}(\theta_{13})\\\
&=\begin{bmatrix}c_{13}c_{14}&0&c_{14}s_{13}&e^{-\iota\delta_{14}}s_{14}\\\
-e^{\iota\delta_{13}}c_{24}s_{13}s_{23}-e^{\iota\delta_{14}}c_{13}s_{14}s_{24}&c_{23}c_{24}&e^{\iota\delta_{13}}c_{13}c_{24}s_{23}-e^{\iota\delta_{14}}s_{13}s_{14}s_{24}&c_{14}s_{24}\\\
-e^{\iota\delta_{13}}c_{23}s_{13}&-s_{23}&e^{\iota\delta_{13}}c_{13}c_{23}&0\\\
-e^{\iota\delta_{14}}c_{13}c_{24}s_{14}+e^{\iota\delta_{13}}s_{13}s_{23}s_{24}&-c_{23}s_{24}&-e^{\iota\delta_{14}}c_{24}s_{13}s_{14}-e^{\iota\delta_{13}}c_{13}s_{23}s_{24}&c_{14}c_{24}\end{bmatrix}\end{split}$
(4)

where we have used notations $s_{ij}=\sin\theta_{ij},c_{ij}=\cos\theta_{ij}$.
Considering the current bounds and best-fit values the vacuum mixing angles
$\theta_{13},\theta_{14},\text{and }\theta_{24}$ are of a similar order. These
small parameters can be expressed in terms of $\mathcal{O}(\lambda^{n})$ with
$\lambda\sim 0.15$ as follows;

$\sin\theta_{13}\simeq\mathcal{O}(\lambda),\sin\theta_{14}\simeq\mathcal{O}(\lambda),\sin\theta_{24}\simeq\mathcal{O}(\lambda),\Delta_{21}\simeq\mathcal{O}(\lambda^{5}),\Delta_{31}\simeq\mathcal{O}(\lambda^{3}),A\simeq\mathcal{O}(\lambda^{3})$
(5)

We can split the total Hamiltonian $H$ into two parts as

$H=\frac{1}{2E_{\nu}}(H_{0}+H_{p})$ (6)

where $H_{p}$, the perturbed Hamiltonian, is proportional to the order of
$\Delta_{31},A[\mathcal{O}(\lambda^{3})]$ whereas the unperturbed Hamiltonian
$H_{0}$ is proportional to $\Delta_{41}$. These can be written as follows,

$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle\Delta_{41}\begin{bmatrix}s^{2}_{14}&e^{-\iota\delta_{14}}c_{14}s_{14}s_{24}&0&e^{-\iota\delta_{14}}c_{24}c_{14}s_{14}\\\
e^{\iota\delta_{14}}c_{14}s_{24}s_{14}&c^{2}_{14}s^{2}_{24}&0&c^{2}_{14}c_{24}s_{24}\\\
0&0&0&0\\\
e^{\iota\delta_{14}}c_{24}c_{14}s_{14}&c^{2}_{14}c_{24}s_{24}&0&c^{2}_{14}c^{2}_{24}\end{bmatrix}\mathrm{,}$
(7) $\displaystyle H_{p}$ $\displaystyle=$
$\displaystyle\tilde{U}\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\
0&0&\Delta_{31}&0\\\
0&0&0&0\end{bmatrix}\tilde{U}^{\dagger}+\begin{bmatrix}A&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&\frac{A}{2}\end{bmatrix}$ (8)

The unperturbed and perturbed Hamiltonian can be expressed in terms of the
small parameter $\lambda$ in the following manner,

$H_{0}\sim\begin{bmatrix}\lambda^{2}&\lambda^{2}&0&\lambda\\\
\lambda^{2}&\lambda^{2}&0&\lambda\\\ 0&0&0&0\\\
\lambda&\lambda&0&1\end{bmatrix}\mathrm{,}H_{p}\sim\begin{bmatrix}\lambda^{5}&\lambda^{4}&\lambda^{4}&-\lambda^{5}\\\
\lambda^{4}&\lambda^{3}&\lambda^{3}&-\lambda^{4}\\\
\lambda^{4}&\lambda^{3}&\lambda^{3}&-\lambda^{4}\\\
-\lambda^{5}&-\lambda^{4}&-\lambda^{4}&\lambda^{5}\end{bmatrix}$ (9)

The unperturbed Hamiltonian has the smallest terms proportional to
$\mathcal{O}(\lambda^{2})$ which is at least an order less than the largest
term in $H_{p}$, the perturbed Hamiltonian. The eigenvalues of $H_{0}$ are
$\lambda_{01}=0$, $\lambda_{02}=0$, $\lambda_{03}=0$,
$\lambda_{04}=\Delta_{41}$. This implies the need for degenerate perturbation
theory to determine the modified energy eigenvalues in the presence of the
matter potential. The modified energy eigenvalues evaluated using degenerated
perturbation theory in ascending order of energy are as follows,

$\begin{split}E_{1m}&=\frac{1}{2E_{\nu}}[\Delta_{31}\sin^{2}(\theta_{13}-\theta_{13m})+A^{{}^{\prime}}\cos^{2}\theta_{13m}(1+\cos^{2}\theta_{14}+\cos^{2}\theta_{14}\sin^{2}\theta_{24})\\\
&-A^{{}^{\prime}}\sin^{2}\theta_{24}\cos 2\theta_{13m}-A\sin
2\theta_{24}\sin\theta_{14}\sin\theta_{23}\sin
2\theta_{13m}\cos\delta/2]\mathrm{,}\\\ E_{2m}&=0\mathrm{,}\\\
E_{3m}&=\frac{1}{2E_{\nu}}[\Delta_{31}\cos^{2}(\theta_{13}-\theta_{13m})+A^{{}^{\prime}}\sin^{2}\theta_{13m}(1+\cos^{2}\theta_{14}+\cos^{2}\theta_{14}\sin^{2}\theta_{24})\\\
&+A^{{}^{\prime}}\sin^{2}\theta_{24}\cos 2\theta_{13m}+A\sin
2\theta_{24}\sin\theta_{14}\sin\theta_{23}\sin
2\theta_{13m}\cos\delta/2]\mathrm{,}\\\
E_{4m}&=\frac{1}{2E_{\nu}}[\Delta_{41}+A^{{}^{\prime}}(1+\sin^{2}\theta_{14}-\cos^{2}\theta_{14}\sin^{2}\theta_{24})]\end{split}$
(10)

where $A^{{}^{\prime}}=A/2=\sqrt{2}G_{F}N_{e}$, the modified angle
$\theta_{13m}$ in the matter is related to the original angles, and the new
phase $\delta=(\delta_{13}-\delta_{14})$ as,

$\displaystyle\sin 2\theta_{13m}$ $\displaystyle=[\Delta_{31}\sin
2\theta_{13}+A^{{}^{\prime}}\cos\delta\sin\theta_{14}\sin\theta_{23}\sin
2\theta_{24}]/f\mathrm{,}$ (11) $\displaystyle\cos 2\theta_{13m}$
$\displaystyle=[\Delta_{31}\cos
2\theta_{13}-A^{{}^{\prime}}(1+\cos^{2}\theta_{14}+\cos^{2}\theta_{14}\sin^{2}\theta_{24}-2\sin^{2}\theta_{24})]/f$
(12)

where $f$ is defined as,

$\begin{split}f&=\sqrt{[\Delta_{31}\sin
2\theta_{13}+A^{{}^{\prime}}s_{14}s_{23}\sin
2\theta_{24}\cos\delta]^{2}+[\Delta_{31}\cos
2\theta_{13}-A^{{}^{\prime}}(1+c^{2}_{14}+c^{2}_{14}s^{2}_{24}-2s^{2}_{24})]^{2}}\end{split}$
(13)

It is noteworthy that for the 3+1 framework, the modified angle $\theta_{13m}$
depends on cp phases unlike in three generation framework. Now if we put $\sin
2\theta_{13m}=1$, i.e., $\cos 2\theta_{13m}=0$, we will get maximum
$\theta_{13m}$, i.e., resonance in this sector for the matter. The
corresponding resonance energy is given by,

$E_{res}=\frac{\Delta_{31}\cos
2\theta_{13}}{\sqrt{2}G_{F}N_{e}(1+\cos^{2}\theta_{14}+\cos^{2}\theta_{14}\sin^{2}\theta_{24}-2\sin^{2}\theta_{24})}$
(14)

The resonance energy for 1300 km and 7000 km are $\sim 11$ GeV, and $8$ GeV
respectively corresponding to
$\theta_{14}=\theta_{24}=7^{\circ},\theta_{13}=8.57^{\circ},\Delta_{31}=2.515\times
10^{-3}\mathrm{eV}^{2}$. It only changes minimally from the three generation
case. The modified active-sterile mixing angles $\theta_{14m},\theta_{24m}$
are related to the vacuum angles as,

$\displaystyle\sin\theta_{14m}=\sin\theta_{14}[1+\frac{A^{{}^{\prime}}}{\Delta_{41}}\cos^{2}\theta_{14}(1+s^{2}_{24})]\mathrm{,}\cos\theta_{14m}=\cos\theta_{14}[1-\frac{A^{{}^{\prime}}}{\Delta_{41}}\sin^{2}\theta_{14}(1+s^{2}_{24})]\mathrm{,}$
(15)
$\displaystyle\sin\theta_{24m}=\sin\theta_{24}[1-\frac{A^{{}^{\prime}}}{\Delta_{41}}\cos^{2}\theta_{14}\cos^{2}\theta_{24}]\mathrm{,}\cos\theta_{24m}=\cos\theta_{24}[1+\frac{A^{{}^{\prime}}}{\Delta_{41}}\cos^{2}\theta_{14}\sin^{2}\theta_{24}]$
(16)

The mixing matrix in matter obtained from the modified eigenvectors using
degenerate perturbation theory is as follows,

$\begin{split}&\tilde{U}_{m}=R_{24}^{m}(\theta_{24m})\tilde{R}_{14}^{m}(\theta_{14m},\delta_{14})R_{23}(\theta_{23})U_{\delta
13}R_{13}^{m}(\theta_{13m})R_{12}^{m}(\theta_{12m})\\\
&=\begin{bmatrix}c_{13m}c_{14m}&(U_{m})_{12}&c_{14m}s_{13m}&e^{-\iota\delta_{14}}s_{14m}\\\
-e^{\iota\delta_{13}}c_{24m}s_{13m}s_{23}-e^{\iota\delta_{14}}c_{13m}s_{14m}s_{24m}&c_{23}c_{24m}&e^{\iota\delta_{13}}c_{13m}c_{24m}s_{23}-e^{\iota\delta_{14}}s_{13m}s_{14m}s_{24m}&c_{14m}s_{24m}\\\
-e^{\iota\delta_{13}}c_{23}s_{13m}&-s_{23}&e^{\iota\delta_{13}}c_{13m}c_{23}&0\\\
-e^{\iota\delta_{14}}c_{13m}c_{24m}s_{14m}+e^{\iota\delta_{13}}s_{13m}s_{23}s_{24m}&-c_{23}s_{24m}&-e^{\iota\delta_{14}}c_{24m}s_{13m}s_{14m}-e^{\iota\delta_{13}}c_{13m}s_{23}s_{24m}&c_{14m}c_{24m}\end{bmatrix}\end{split}$
(17)

where the original vacuum angles are replaced by modified angles as given by
(11), (12), (15), (16) and null value of the element $(\tilde{U})_{12}$ in
vacuum mixing matrix $\tilde{U}$(4) is modified as
$(U_{m})_{12}=\frac{A}{\Delta_{41}}e^{-\iota\delta_{14}}c_{14}c_{23}c_{24}s_{14}s_{24}\sim\mathcal{O}(\lambda^{5})$.
This is due to the fact that the matter effect introduces correction of mixing
angle $\theta_{12}$ which was absent before due to the approximation
$\Delta_{21}=0$. The other terms related to $\theta_{12}$ don’t show up as
they are $>\mathcal{O}(\lambda^{5})$. Now we can calculate the
oscillation(survival) probabilities using the elements of $\tilde{U}_{m}$ in
place of $U$ and $\Delta_{ij}^{m}=2E_{\nu}(E_{im}-E_{jm})$ replacing
$\Delta_{ij}$ in (18) assuming constant matter density,

$P_{\alpha\beta}=\delta_{\alpha\beta}-4\sum_{i>j}^{N}Re(U_{\alpha
i}^{\ast}U_{\beta i}U_{\alpha j}U_{\beta
j}^{\ast})\sin^{2}\frac{1.27\Delta_{ij}L}{E_{\nu}}+2\sum_{i>j}^{N}Im(U_{\alpha
i}^{\ast}U_{\beta i}U_{\alpha j}U_{\beta j}^{\ast})\sin
2\frac{1.27\Delta_{ij}L}{E_{\nu}}$ (18)

On the other hand, the exact numerical probability at constant matter density
with total Hamiltonian given by (3) can be evaluated as,

$P_{\alpha\beta}^{\mathrm{num}}=|[e^{-\iota HL}]_{\alpha\beta}|^{2}$ (19)

####  3.1.1 $P_{\mu e}$ Channel

The appearance channel, i.e., $\nu_{\mu}\rightarrow\nu_{e}$ oscillation
probability is given by,

$P_{\mu e}=P_{\mu e}^{1}+P_{\mu e}^{2}+P_{\mu e}^{3}+\mathcal{O}(\lambda^{6})$
(20)

where the different significant terms of the probability $P_{\mu e}$ are as
follows,

$\displaystyle\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)\\\
&-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)\mathrm{,}\end{split}$
(21) $\displaystyle\begin{split}P_{\mu e}^{2}&=\cos^{2}\theta_{14m}\sin
2\theta_{13m}\sin\theta_{14m}\sin\theta_{23}\sin
2\theta_{24m}\sin\frac{1.27\Delta_{41}^{m}L}{E}\sin(\frac{1.27\Delta_{41}^{m}L}{E}-\delta)\\\
&+\sin^{2}2\theta_{14m}\sin^{2}\theta_{24m}\cos^{2}\theta_{13m}\sin^{2}\frac{1.27\Delta_{41}^{m}L}{E}\mathrm{,}\end{split}$
(22) $\displaystyle\begin{split}P_{\mu e}^{3}&=-\cos^{2}\theta_{14m}\sin
2\theta_{13m}\sin\theta_{14m}\sin\theta_{23}\sin
2\theta_{24m}\sin\frac{1.27\Delta_{43}^{m}L}{E}\sin(\frac{1.27\Delta_{43}^{m}L}{E}-\delta)\\\
&+\sin^{2}2\theta_{14m}\sin^{2}\theta_{24m}\sin^{2}\theta_{13m}\sin^{2}\frac{1.27\Delta_{43}^{m}L}{E}\end{split}$
(23)

Figure 1: The total analytic probability $P_{\mu e}$ (orange) along with its
other dominant terms in the top panel and the absolute differences $|P_{\mu
e}-P_{\mu e}^{\mathrm{num}}|$ (red) and $|P_{\mu e}-P_{\mu e}^{\mathrm{GL}}|$
(cyan) in the bottom panel at 1300 km(left), and 7000 km(right) baselines.

The total analytic probability $P_{\mu e}$(orange) and the dominant terms
contributing to it are plotted at 1300 km and 7000 km baselines as a function
of neutrino energy $E_{\nu}$ in the top panel of Figure 1. For the plots, and
calculations of $P_{\mu e},P_{\mu\mu}$ in this section, we have considered
$\theta_{12}=33.44^{\circ},\theta_{13}=8.57^{\circ},\theta_{23}=49^{\circ},\theta_{14}=\theta_{24}=7^{\circ},\delta_{13}=195^{\circ},\delta_{14}=30^{\circ},\Delta_{31}=2.515\times
10^{-3}$eV2, and $\Delta_{41}=1$eV2. The analytic expression of $P_{\mu e}$
consists of three significant terms, although there are other higher order
terms [$\mathcal{O}(\lambda^{6})$] that are neglected. The first term in
(21)(blue curve) that is proportional to the modified mass squared difference
$\Delta_{31}^{m}$, is the most dominant one and provides the average curve of
the total probability as seen in Figure 1. The fast oscillations seen Figure 1
are a manifestation of the terms in (22) (green curve), (23) (violet curve)
which are proportional to the modified mass squared differences related to the
sterile neutrino mass states $\Delta_{41}^{m},\Delta_{43}^{m}$ respectively.
The fast oscillations are not reflected in experiments, as we can only get the
average probability. Also, these terms are relatively much smaller than the
$P_{\mu e}^{1}$ around probability maxima, so in the next section while
discussing the degeneracies we will only use the term $P_{\mu e}^{1}$. Putting
$\theta_{14},\theta_{24}$ angles to zero in equations (21), (22), (23) gives
the standard three flavour oscillation probability from the very first term of
the (21) as the other terms goes to zero due to presence of
$\sin\theta_{14m},\sin\theta_{24m}$.

We have shown the comparison of the absolute differences $|\Delta P|$ of the
analytic probability $P_{\mu e}$ (20) with the exact probability $P_{\mu
e}^{\mathrm{num}}$(19) (red) as well as with the probability $P_{\mu
e}^{\mathrm{GL}}$ (cyan) obtained using GLoBESHuber  _et al._ (2005) as a
function of neutrino energy at the bottom panel in Figure 1. We can see the
value of $|\Delta P|$ is around $10^{-3}$ for most of the energies. $|\Delta
P|$ values are smaller around the resonance energy of $11$ GeV for 1300 km and
$8$ GeV for 7000 km. Also, we can see the energies at which the value of
probability is smaller we get smaller values of $|\Delta P|$. Overall we can
conclude that the analytic probability $P_{\mu e}$ using TMSD approximation is
in good agreement with both numerical and exact ones, better with the exact
one $P_{\mu e}^{\mathrm{num}}$ for all energies (>0.5 GeV), especially around
the resonance. This is similar to OMSD probabilities matching well with
numerical ones in the standard three flavour case.

####  3.1.2 $P_{\mu\mu}$ Channel

The disappearance channel, i.e., $\nu_{\mu}\rightarrow\nu_{\mu}$ survival
probability is given by,

$P_{\mu\mu}=1-P_{\mu\mu}^{1}-P_{\mu\mu}^{2}-P_{\mu\mu}^{3}+\mathcal{O}(\lambda^{6})$
(24)

where the significant terms of the probability are as follows,

$\displaystyle\begin{split}P_{\mu\mu}^{1}&=\cos^{4}\theta_{24m}\sin^{2}2\theta_{13m}\sin^{4}\theta_{23}\sin^{2}\frac{1.27\Delta_{31}^{m}L}{E}+\sin^{4}\theta_{24m}\sin^{4}\theta_{14m}\sin^{2}2\theta_{13m}\sin^{2}\frac{1.27\Delta_{31}^{m}L}{E}\\\
&+\sin 2\theta_{24m}\sin\theta_{14m}\sin
4\theta_{13m}\sin\theta_{23}\cos\delta(\cos^{2}\theta_{24m}\sin^{2}\theta_{23}-\sin^{2}\theta_{24m}\sin^{2}\theta_{13m})\sin^{2}\frac{1.27\Delta_{31}^{m}L}{E}\\\
&+4\cos^{2}\theta_{24m}\sin^{2}\theta_{24m}\sin^{2}\theta_{14m}\sin^{2}\theta_{23}(1-\frac{\sin^{2}2\theta_{13m}}{2}-\sin^{2}2\theta_{13m}\cos^{2}\delta)\sin^{2}\frac{1.27\Delta_{31}^{m}L}{E}\mathrm{,}\end{split}$
(25)
$\displaystyle\begin{split}P_{\mu\mu}^{2}&=\cos^{4}\theta_{24m}\cos^{2}\theta_{13m}\sin^{2}2\theta_{23}\sin^{2}\frac{1.27\Delta_{32}^{m}L}{E}\\\
&+4\cos^{2}\theta_{24m}\sin^{2}\theta_{24m}\sin^{2}\theta_{14m}\sin^{2}\theta_{13m}\cos^{2}\theta_{23}\sin^{2}\frac{1.27\Delta_{32}^{m}L}{E}\\\
&-4\cos^{3}\theta_{24m}\sin\theta_{24m}\sin\theta_{14m}\sin
2\theta_{13m}\cos^{2}\theta_{23}\sin\theta_{23}\cos\delta\sin^{2}\frac{1.27\Delta_{32}^{m}L}{E}\mathrm{,}\end{split}$
(26)
$\displaystyle\begin{split}P_{\mu\mu}^{3}&=\cos^{4}\theta_{24m}\sin^{2}\theta_{13m}\sin^{2}2\theta_{23}\sin^{2}\frac{1.27\Delta_{21}^{m}L}{E}\\\
&+4\cos^{2}\theta_{24m}\sin^{2}\theta_{24m}\sin^{2}\theta_{14m}\cos^{2}\theta_{13m}\cos^{2}\theta_{23}\sin^{2}\frac{1.27\Delta_{21}^{m}L}{E}\\\
&+4\cos^{3}\theta_{24m}\sin\theta_{24m}\sin\theta_{14m}\sin
2\theta_{13m}\cos^{2}\theta_{23}\sin\theta_{23}\cos\delta\sin^{2}\frac{1.27\Delta_{21}^{m}L}{E}\end{split}$
(27)

Figure 2: The total analytic probability $P_{\mu\mu}$ (orange) along with its
other dominant terms in the top panel and the absolute differences
$|P_{\mu\mu}-P_{\mu\mu}^{\mathrm{num}}|$ (red) and
$|P_{\mu\mu}-P_{\mu\mu}^{\mathrm{GL}}|$ (cyan) in the bottom panel at 1300
km(left), and 7000 km(right) baselines.

We show the total analytic probability $P_{\mu\mu}$ (orange) and the different
terms contributing significantly to it at 1300 km and 7000 km baselines in the
top panel of Figure 2 as a function of neutrino energy. The analytic
expression of $P_{\mu\mu}$ consists of three significant terms (25), (26), and
(27), although there are three other fast oscillating terms that are
neglected. Here, the fast oscillating terms are proportional to the sterile
mass squared differences $\Delta_{41}^{m},\Delta_{42}^{m},\Delta_{43}^{m}$ and
are of higher orders [$\mathcal{O}(\lambda^{6})$]. The first term in (25)
(blue curve) which is proportional to the modified mass squared difference
$\Delta_{31}^{m}$, has a dependence on octant of $\theta_{23}$ in the leading
order due to the presence of $\sin^{4}\theta_{23}$. $P_{\mu\mu}^{1}$ grows
with energy initially and decreases after resonance energy. The second, and
third terms in (26) (green curve), (27) (violet curve) which are proportional
to the modified mass squared differences $\Delta_{32}^{m},\Delta_{21}^{m}$
respectively, show no octant dependence in the leading order due to the
presence of $\sin^{2}2\theta_{23}$. The second term is the most dominant one
before resonance energy but almost becomes zero after resonance energy,
whereas the third term only grows after the resonance energy. In the case of
7000 km at oscillation maxima of 7.5 GeV,
$P_{\mu\mu}^{1},P_{\mu\mu}^{2},P_{\mu\mu}^{3}$ all have significant
contributions. Putting the $\theta_{14},\theta_{24}$ angles to zero, we will
get back the three flavour oscillation probability from the first term of the
equations (25), (26), and (27).

It has been also shown in the bottom panel of Figure 2, the absolute
differences $|\Delta P|$ of the analytical probability $P_{\mu\mu}$(24) with
the exact probability $P_{\mu\mu}^{\mathrm{num}}$(19) (red) and the
probability $P_{\mu\mu}^{\mathrm{GL}}$ (cyan) obtained using GLoBES at 1300 km
and 7000 km baselines. We observe that value of $|\Delta P|$ is mostly around
$10^{-3}$. The $|\Delta P|$ values are seen to be lower around resonance
energies. We can also see the $|\Delta P|$ value going down at the minima or
at the regions where the value of probability is less. The $|\Delta
P_{\mu\mu}|$ for 7000 km is increasing after resonance energy as the dominant
term in those energies is $P_{\mu\mu}^{3}$ that is proportional to
$\Delta_{21m}$ and hence is affected by the $\Delta_{21}=0$ approximation222In
the appendix we have shown that with non-zero $\Delta_{21}$ in the Cayley
Hamilton method we get better fit at these regions as well as at very low
energies.. Hence, we can conclude that the analytical probability $P_{\mu\mu}$
using TMSD approximation is in agreement with exact and numerical
probabilities to a good extent, matching better with the exact one
$P_{\mu\mu}^{\mathrm{num}}$.

##  4 Octant Degeneracy

The degeneracy in the determination of the octant of $\theta_{23}$ can arise
from both the survival/oscillation probabilities as follows:

  * •

When the probability is a function of $sin^{2}{2\theta_{23}}$, it is not
possible to differentiate between the probabilities arising due to
$\theta_{23}$ and $\frac{\pi}{2}-\theta_{23}$. This is called intrinsic octant
degeneracyFogli  _et al._ (1994).

  * •

When the probability is a function of $sin^{2}{\theta_{23}}$ or
$cos^{2}{\theta_{23}}$, the degeneracy of the octant arises due to the
uncertainties in the Dirac CP phase $\delta_{CP}$.

$P(\theta_{23}^{\rm{right}},\delta_{13})=P(\theta_{23}^{\rm{wrong}},\delta_{13}^{{}^{\prime}})$
(28)

  * •

Addition of a light sterile neutrino brings an extra phase $\delta_{14}$ which
will also affect the determination of octant just like in the above case
through additional degeneracies.

$P(\theta_{23}^{\rm{right}},\delta_{14})=P(\theta_{23}^{\rm{wrong}},\delta_{14}^{{}^{\prime}})$
(29)

  * •

Considering known hierarchy and two unknown phases there will be a new 8 fold
octant-$\delta_{13}$-$\delta_{14}$ degeneracy.

$P(\theta_{23}^{\rm{right}},\delta_{13},\delta_{14})=P(\theta_{23}^{\rm{wrong}},\delta_{13}^{{}^{\prime}},\delta_{14}^{{}^{\prime}})$
(30)

We consider the normal hierarchy ($\Delta_{31}=2.515\times 10^{-3}$ eV2) for
our octant degeneracy study. Therefore, we have a 8-fold
octant-$\delta_{13}$-$\delta_{14}$ degeneracy in presence of a sterile
neutrino as depicted in Table 1. For unknown hierarchy, this will become a
16-fold degeneracy.

Solution with right octant | Solution with wrong octant
---|---
RO-R$\delta_{13}$-R$\delta_{14}$ | WO-R$\delta_{13}$-R$\delta_{14}$
RO-R$\delta_{13}$-W$\delta_{14}$ | WO-R$\delta_{13}$-W$\delta_{14}$
RO-W$\delta_{13}$-R$\delta_{14}$ | WO-W$\delta_{13}$-R$\delta_{14}$
RO-W$\delta_{13}$-W$\delta_{14}$ | WO-W$\delta_{13}$-W$\delta_{14}$
Table 1: New degeneracies in presence of unknown octant and phases with fixed
hierarchy.

In order to understand the degeneracy analytically we follow the method
outlined in Agarwalla  _et al._ (2017) and use the TMSD probabilities derived
in the earlier section. The current $3\sigma$ range of $\theta_{23}$ is
$[39.7^{\circ},50.9^{\circ}]$ Esteban  _et al._ (2017) for normal hierarchy.
We can express $\theta_{23}$ w.r.t. $\pi/4$ as,

$\theta_{23}=\frac{\pi}{4}\pm\eta$ (31)

where the deviation in value of $\theta_{23}$ from current global analysis fit
is given by $\eta\sim 0.1$ with the plus and minus sign in (31) indicating
higher octant(HO), and lower octant (LO) of $\theta_{23}$ respectively. The
octant sensitivity will be there if there is a difference between
probabilities of the two opposite octants even when the phases
$\delta_{13},\delta_{14}$ vary in the range $[-\pi,\pi]$. The octant
sensitivity from the appearance channel probability $P_{\mu e}$ is defined as,

$\Delta P_{oct,1}\equiv P_{\mu
e}^{1HO}(\delta_{13}^{HO},\delta_{14}^{HO})-P_{\mu
e}^{1LO}(\delta_{13}^{LO},\delta_{14}^{LO})>0$ (32)

As $\eta$ is small we can have the following expansion

$\sin^{2}\theta_{23}\simeq\frac{1}{2}\pm\eta,\sin\theta_{23}\simeq\frac{1}{\sqrt{2}}(1\pm\eta),\cos\theta_{23}\simeq\frac{1}{\sqrt{2}}(1\mp\eta)$
(33)

Putting $P_{\mu e}^{1}$ from (21) in (32) and using the above expressions of
(33), we get three contributions to $\Delta P_{oct,1}$ corresponding to the
three terms in $P_{\mu e}^{1}$,

$\begin{split}\Delta
P_{0}&=8\eta\cos^{2}\theta_{13m}\cos^{2}\theta_{14m}\cos^{2}\theta_{24m}\sin^{2}\theta_{13m}\sin^{2}D_{31}^{m}\mathrm{,}\\\
\Delta
P_{1}&=X_{1}[\sin(D_{31}^{m}+\delta^{HO})-\sin(D_{31}^{m}+\delta^{LO})]+\eta
X_{1}[\sin(D_{31}^{m}+\delta^{HO})+\sin(D_{31}^{m}+\delta^{LO})]\mathrm{,}\\\
\Delta
P_{2}&=-Y_{1}[\sin(D_{31}^{m}-\delta^{HO})-\sin(D_{31}^{m}-\delta^{LO})]-\eta
Y_{1}[\sin(D_{31}^{m}-\delta^{HO})+\sin(D_{31}^{m}-\delta^{LO})]\end{split}$
(34)

The contribution of the fast oscillation terms $P_{\mu e}^{2},P_{\mu e}^{3}$
to the octant sensitivity is,

$\Delta
P_{fast}=\sum_{k=1,3}Z_{k}[\sin(D_{4k}^{m}-\delta^{HO})-\sin(D_{4k}^{m}-\delta^{LO})]+\eta
Z_{k}[\sin(D_{4k}^{m}-\delta^{HO})+\sin(D_{4k}^{m}-\delta^{LO})]$ (35)

Where $D_{ij}^{m}=\frac{1.27\Delta_{ij}^{m}}{E}$. Now we can rewrite (32) for
octant sensitivity as,

$\Delta P_{oct,1}=\Delta P_{0}+\Delta P_{1}+\Delta P_{2}+\Delta P_{fast}$ (36)

Among the terms of $\Delta P_{oct,1}$ (34), $\Delta P_{0}$ has no dependence
on phase and is positive whereas the values of $\Delta P_{1},\Delta
P_{2},\Delta P_{fast}$ can be both positive and negative as they contain
phases. Thus degeneracy can occur when $\Delta P_{1}+\Delta P_{2}+\Delta
P_{fast}$ is negative and is of the same order as $\Delta P_{0}$ making
$\Delta P$ zero. $X_{1},Y_{1}$ the positive definite amplitudes of $\Delta
P_{1}$, $\Delta P_{2}$ respectively as well as the amplitudes $Z_{1},Z_{3}$ of
$\Delta P_{fast}$ are as follows,

$\begin{split}X_{1}&=\sqrt{2}\cos^{3}\theta_{13m}\cos^{2}\theta_{14m}\sin\theta_{13m}\sin\theta_{14m}\sin
2\theta_{24m}\sin D_{31}^{m}\mathrm{,}\\\
Y_{1}&=\sqrt{2}\cos\theta_{13m}\cos^{2}\theta_{14m}\sin^{3}\theta_{13m}\sin\theta_{14m}\sin
2\theta_{24m}\sin D_{31}^{m}\mathrm{,}\\\ Z_{1}&=\cos^{2}\theta_{14m}\sin
2\theta_{13m}\sin\theta_{14m}\sin 2\theta_{24m}\sin
D_{41}^{m}/\sqrt{2}\mathrm{,}\\\ Z_{3}&=-\cos^{2}\theta_{14m}\sin
2\theta_{13m}\sin\theta_{14m}\sin 2\theta_{24m}\sin
D_{43}^{m}/\sqrt{2}\end{split}$ (37)

Now, if we inspect the possibility of the octant degeneracy for the DUNE
experiment with a baseline of 1300 km. We use the following values of the
oscillation parameters:
$\theta_{12}=33.47^{\circ},\theta_{13}=8.54^{\circ},\theta_{14}=7^{\circ},\theta_{24}=7^{\circ},\Delta_{31}=2.515\times
10^{-3}\rm{eV^{2}},\Delta_{41}=1\rm{eV^{2}}$. For 1300 km at oscillation
maxima of $2.5$ GeV, the values of various terms of $\Delta P_{oct,1}$ are,

$\Delta P_{0}=0.0279,X_{1}=0.0073,Y_{1}=0.0003,Z_{1}=-0.0056,Z_{3}=0.0064$
(38)

Therefore, $\Delta P_{2}$ is negligible compared to $\Delta P_{0}$, $\Delta
P_{fast}$, and $\Delta P_{1}$ in DUNE due to presence of extra
$\sin^{2}\theta_{13m}$ in $Y_{1}$(37). It can be seen from (34) the square
bracketed terms multiplying $X_{1},Z_{1},Z_{3}$ can vary from $-2:+2$ and
therefore for certain phase choices cancellation can occur resulting in loss
of octant sensitivity in 1300 km baseline when fast oscillations considered.
However, in the absence of fast oscillations, there is octant sensitivity.

Next, we use the analytic expressions in (36) to understand the octant
sensitivity at 7000 km. In case of 7000 km at oscillation maxima of $E=6.5$
GeV, the values of the different terms contributing to $\Delta P$ are,

$\Delta P_{0}=0.1453,X_{1}=0.0133,Y_{1}=0.0040,Z_{1}=0.0001,Z_{3}=-0.0164$
(39)

It shows that $\Delta P_{0},X_{1},Z_{3}$ are the dominant contributions and
any combination of phases can not make $\Delta P_{oct,1}=0$ as the value of
$P_{0}$ is one order greater than $X_{1}$. It shows probabilities $(P_{\mu
e})$ corresponding to two different octants will be always well separated from
each other, i.e., the octant-$\delta_{13}$-$\delta_{14}$ degeneracy will be
removed. This suggests unlike in 1300 km here even with the variation of
phases in both octants, we can have significant octant sensitivity at higher
baselines. This is mainly because at the higher baselines $\Delta P_{0}$ has
much higher values than others due to higher matter effects.

The octant sensitivity from the disappearance channel probability $P_{\mu\mu}$
is defined as,

$\Delta P_{oct,2}\equiv
P_{\mu\mu}^{1HO}(\delta_{13}^{HO},\delta_{14}^{HO})-P_{\mu\mu}^{1LO}(\delta_{13}^{LO},\delta_{14}^{LO})>0$
(40)

As we have seen earlier the largest octant sensitive term in $P_{\mu\mu}$
comes from (25), we put that in the above (40) to get the difference in
opposite octant probabilities as,

$\Delta
P_{oct,2}=\cos^{2}\theta_{24m}[2\eta\cos^{2}\theta_{24m}\sin^{2}2\theta_{13m}+\sin
2\theta_{24m}\sin\theta_{14m}\sin
4\theta_{13m}(\cos\delta^{HO}-\cos\delta^{LO})\frac{1+3\eta}{2\sqrt{2}}]\sin^{2}D^{m}_{31}$
(41)

###  4.1 Degeneracy in $\cos\theta_{\nu}-E_{\nu}$ Plane

Figure 3:  $\Delta P_{\mu e}$(left), $\Delta P_{\mu\mu}$(right), i.e, the
absolute differences in probabilities for $\theta_{23}$ values from opposite
octant with fixed value of $\delta_{13},\delta_{14}$ in
$\cos\theta_{\nu}-E_{\nu}$ plane.

To see the octant sensitivity spanning over all the baselines and energies, we
present the oscillogram plots of the differences in probabilities
corresponding to the value of $\theta_{23}=40^{\circ}$(LO) and
$\theta_{23}=50^{\circ}$(HO) in $\cos\theta_{\nu}-E_{\nu}$ plane for normal
hierarchy in Figure 3. The phases are kept fixed at same
$\delta_{13}=195^{\circ},\delta_{14}=30^{\circ}$ for both the octants. From
the figure, it can be seen that the maximum difference is obtained at the
energy range of $5:10$ GeV for $\cos\theta_{\nu}$ in the range of $-0.5:-0.8$
which roughly translates to baselines around 5000-10000 km. This figure serves
as a reference to show at which baselines and energies the octant sensitivity
can be maximum and motivates us to add the contribution from atmospheric
neutrinos to obtain better octant sensitivity in our analysis.

###  4.2 Degeneracy with variation of $\delta_{13},\delta_{14}$ at fixed
baseline

In this section, we study the probabilities as a function of the phases to
understand the dependency of the degeneracy on these parameters. The variation
of the probability $P_{\mu e}$(left), $P_{\bar{\mu}\bar{e}}$(right) as a
function of phases $\delta_{13}$ (top), and $\delta_{14}$ (bottom) at 2.5 GeV
is displayed in Figure 4 for values of
$\theta_{23}=39^{\circ}(grey),42^{\circ}(orange),48^{\circ}(violet),51^{\circ}(blue)$
spanning over both octants. The curves for other values of $\theta_{23}$ will
lie in between these ranges. The bands correspond to variation over the non-
displayed phase $\delta_{14}(top)/\delta_{13}(bottom)$ over the range
$-180^{\circ}:180^{\circ}$ respectively. Three horizontal iso-probability
lines are drawn in Figure 4 to indicate the values of
$\delta_{13}/\delta_{14}$ for which there are degeneracies (dot-dashed line)
and there are no degeneracies (dotted, dashed lines) between the two octants.
Note that in the probability vs $\delta_{13}$ plots for the three-generation
case there is a single curve for each $\theta_{23}$ whereas in presence of
sterile neutrino there are bands due to $\delta_{14}$ variation for a fixed
$\theta_{23}$.

We can infer the following points from Figure 4,

  * •

The regions above the dotted line in the top panels indicate the values of
$\delta_{13}$ for which there is no degeneracy in HO. This is around
$\delta_{13}=-90^{\circ}(90^{\circ})$ in $P_{\mu e}(P_{\bar{\mu}\bar{e}})$
channel. However, some portions of the blue and violet bands extend below the
dotted lines in both figures and sometimes also overlap with the orange and
violet bands indicating that for these values of $\delta_{13}$ there are still
degeneracies for certain values of $\delta_{14}$.

  * •

Similarly, the regions below the dashed lines in the top panels signify the
$\delta_{13}$ values devoid of degeneracy for $\theta_{23}$ in LO. This region
for $P_{\mu e}(P_{\bar{\mu}\bar{e}})$ channel is around
$\delta_{13}=90^{\circ}(-90^{\circ})$. Here also, the portions of grey and
orange bands above the dashed lines as well as the portions coinciding with
the blue and violet bands indicate the existence of degeneracies at these
values of $\delta_{13}$.

  * •

From the top panels, we can clearly see a synergy between neutrino and anti-
neutrino channels for octant degeneracy in both HO and LO. For instance, for
HO (LO), the degeneracy is present around
$\delta_{13}=90^{\circ}(-90^{\circ})$ at $P_{\mu e}$ channel but absent for
$P_{\bar{\mu}\bar{e}}$.

  * •

In the bottom panels, the regions above the dotted line indicate that the no
degeneracy region in HO lies around $\delta_{14}=-60^{\circ}(60^{\circ})$ for
$P_{\mu e}(P_{\bar{\mu}\bar{e}})$ channel. Note that the region has a larger
spread in $\delta_{14}$ over
$-180^{\circ}:95^{\circ}(-70^{\circ}:140^{\circ})$ for
$\theta_{23}=51^{\circ}$, and over
$-180^{\circ}:65^{\circ}(-50^{\circ}:120^{\circ})$ for
$\theta_{23}=48^{\circ}$ in $P_{\mu e}(P_{\bar{\mu}\bar{e}})$ channel.
Corresponding nondegenerate regions have a smaller spread in $\delta_{13}$ as
seen from the top panel plots.

  * •

There are regions below the dashed line signifying no degeneracy in LO for the
plots in the bottom panels. These regions occur around
$\delta_{14}=-60^{\circ}(60^{\circ})$ for $P_{\mu e}(P_{\bar{\mu}\bar{e}})$
channel. However, it is to be noted that unlike in the top panel, the non-
degenerate region in LO is over the similar range of $\delta_{14}$ w.r.t HO as
mentioned in the previous point. Therefore we see in the neutrino (anti-
neutrino) channel maximum sensitivity for both HO and LO is around
$\delta_{14}=60^{\circ}(-60^{\circ})$.

  * •

In the bottom panels, the probability bands are wider and the extent of
overlap is higher around $-60^{\circ}(60^{\circ})$ in $P_{\mu
e}(P_{\bar{\mu}\bar{e}})$ channel. These give rise to WO-R$\delta_{14}$
degeneracies which are hard to resolve using neutrino plus anti-neutrino. The
synergy between neutrino and anti-neutrino channels for octant degeneracy is
less pronounced here.

  * •

In the bottom panels for $P_{\mu e}$ channel, around
$\delta_{14}=130^{\circ}$, there is a small region where there is no
WO-R$\delta_{14}$ degeneracy between HO and LO for all values of
$\delta_{13}$. For $P_{\bar{\mu}\bar{e}}$ channel there is a similar region
with minimum degeneracy around $\delta_{14}=-130^{\circ}$.

  * •

When the probability bands from HO (blue and violet) coincide with bands from
LO (orange and grey) at the same $\delta_{13}/\delta_{14}$ values, those are
examples of WO-R$\delta_{13}$/R$\delta_{14}$ degeneracies. While the regions
of bands from opposite octants connected through iso-probability lines show
WO-W$\delta_{13}$/W$\delta_{14}$ degeneracies.

Figure 4: $P_{\mu e}$ (left), and $P_{\bar{\mu}\bar{e}}$ (right) as a function
of $\delta_{13}$ (top), $\delta_{14}$ (bottom) for variation of the respective
another phase at neutrino energy 2.5 GeV at 1300 km baseline for NH.

While performing $\chi^{2}$ analysis, we take fixed true values of parameters
in one octant and marginalize over the relevant parameters in the opposite
octant. Therefore, a better understanding of the octant degeneracy can be
achieved if we keep $\theta_{23}$, and the phases constant in one octant and
vary them in the opposite one. We replicate this in Figure 5 where the
probabilities in neutrino (left) and anti-neutrino (right) channels are drawn
as a function of phase $\delta_{13}$. In the top [bottom] panel the green
[red] solid(dashed) line corresponds to $\theta_{23}=49^{\circ}[41^{\circ}]$
and $\delta_{14}=0^{\circ}(90^{\circ})$. The grey and orange [violet and blue]
bands correspond to $\theta_{23}=39^{\circ},42^{\circ}[48^{\circ},51^{\circ}]$
in LO[HO] for $\delta_{14}$ varying over $-180^{\circ}:180^{\circ}$. The
horizontal iso-probability lines in the plots demarcate different degenerate
and non-degenerate regions. The important points from Figure 5 are as follows,

  * •

In the top panel, the region above the dotted line corresponds to no
degeneracy. This region is around $\delta_{13}=-90^{\circ}(90^{\circ})$ at
$P_{\mu e}(P_{\bar{\mu}\bar{e}})$ channel for green solid
($\delta_{14}=0^{\circ}$) curve. However, the green dashed
($\delta_{14}=90^{\circ}$) curve has a non-degenerate region only in
$P_{\bar{\mu}\bar{e}}$ channel around $\delta_{13}=90^{\circ}$. This suggests
that for $\delta_{14}=0^{\circ}$, the octant sensitivity comes from both
$P_{\mu e}$, and $P_{\bar{\mu}\bar{e}}$ channel around $\delta_{13}\sim
0^{\circ}$ whereas for $\delta_{14}=90^{\circ}$ sensitivity comes only from
$P_{\bar{\mu}\bar{e}}$ channel around $\delta_{13}\sim 90^{\circ}$.

  * •

For the bottom panel, the non-degenerate regions are below the dashed
horizontal line. In $P_{\mu e}$ channel this region is around
$\delta_{13}=120^{\circ}$ for $\delta_{14}=0^{\circ}$. A very small region for
$\delta_{14}=90^{\circ}$ also extends below the dashed line. In
$P_{\bar{\mu}\bar{e}}$ channel the region of no degeneracy lies around
$\delta_{13}=-120^{\circ}(-60^{\circ})$ for
$\delta_{14}=0^{\circ}(90^{\circ})$.

Figure 5: $P_{\mu e}$ (left), and $P_{\bar{\mu}\bar{e}}$ (right) as a function
of $\delta_{13}$ for variation of the phase $\delta_{14}$ at neutrino energy
2.5 GeV at 1300 km baseline for NH.

Figure 6: $P_{\mu\mu}$ (left), and $P_{\bar{\mu}\bar{\mu}}$ (right) vs
$\delta_{13}$ (top), $\delta_{14}$ (bottom) for variation of the respective
another phase at neutrino energy 2.5 GeV at 1300 km baseline for NH.

Now we focus on the disappearance channel probabilities $P_{\mu\mu}$ (left),
and $P_{\bar{\mu}\bar{\mu}}$ (right) as a function of phases $\delta_{13}$(top
panel), $\delta_{14}$(bottom panel) at 2.5 GeV in Figure 6. The following
points may be noted,

  * •

The bands due to variation of $\delta_{13}/\delta_{14}$ are narrower than the
ones for appearance channel. Hence these bands are well separated from each
other.

  * •

The bands corresponding to $\theta_{23}=51^{\circ}$(blue) in HO comes in
between the bands corresponding to $\theta_{23}=39^{\circ}$(grey) and
$\theta_{23}=42^{\circ}$(yellow) in LO. On the other hand, the violet band
corresponding to $\theta_{23}=48^{\circ}$ is outside the whole region of LO
between the grey and yellow band. This implies the presence (absence) of the
octant degeneracy for $\theta_{23}=51^{\circ}(48^{\circ})$ in $P_{\mu\mu}$
channel.

  * •

Similarly, $\theta_{23}=39^{\circ}$ (grey) in LO demonstrates octant
sensitivity since it lies outside the HO region between the blue and violet
bands but $\theta_{23}=42^{\circ}$ (yellow) lies within the HO region and
therefore is not sensitive to the octant. A similar feature can also be seen
from probability vs $\delta_{14}$ plots in the bottom panel.

  * •

We can conclude that even for the $P_{\mu\mu}$ channel, there are certain
trues values of $\theta_{23}$ in both the octants for which there can be
contributions to the octant sensitivity at 1300 km.

Next, we study the behavior of the probabilities at a higher baseline of 7000
km where the resonant matter effect comes into play. We observe the appearance
probability $P_{\mu e}$ as a function of the phase $\delta_{13}$ (left), and
$\delta_{14}$ while the respective other phase variation creates band at
different values of $\theta_{23}=39^{\circ},42^{\circ},48^{\circ},51^{\circ}$
spanning over both octants at energy maxima of 6.5 GeV in Figure 7. We see
similar variations of the disappearance channel probability $P_{\mu\mu}$ at
maxima energy of 7 GeV in Figure 8. The following facts can be noted,

  * •

Unlike at 1300 km, the $P_{\mu e}$ probability bands of opposite octant at
7000 km are clearly separated. It suggests that even with the variation of
phases and $\theta_{23}$ in both octants, the octant degeneracy can be clearly
removed at higher baselines.

  * •

In $P_{\mu\mu}$ channel, the LO and HO bands are mostly separated apart from
occurrence of WO-W$\delta_{13}$ (left panel), WO-R$\delta_{14}$/W$\delta_{14}$
(right panel) degeneracies respectively around $\delta_{13}$, $\delta_{14}$
values of $\pm 150^{\circ}$ in a tiny region. This suggests contributions to
the octant sensitivity also come from the $P_{\mu\mu}$ channel. The
sensitivity of the octant in $P_{\mu\mu}$ comes from the first term in (25)
which has a more significant contribution at 7000 km than 1300 km as noted in
Figure 2 due to larger matter effect.

Figure 7: $P_{\mu e}$ vs $\delta_{13}$(left), and $\delta_{14}$(right) for
variation of the respective another phase at neutrino energy 6.5 GeV at 7000
km baseline for NH. Figure 8: $P_{\mu\mu}$ vs $\delta_{13}$(left), and
$\delta_{14}$(right) for variation of the respective another phase at energy 7
GeV at 7000 km baseline for NH.

##  5 Experimental and Simulation Details of the DUNE experiment

The Deep Underground Neutrino Experiment (DUNE)Abi  _et al._ (2020b) will
consist of a near detector (ND) and far detector (FD) exposed to a megawatt-
scale muon neutrino beam produced by Long Baseline Neutrino Facility
(LBNF)Acciarri  _et al._ (2016) at Fermilab. The ND will be placed close to
the source of the beam, while the FD, comprising four LArTPC detectors of 10
kt each will be installed 1300 km away and 1.5 km underground from the
neutrino source in the Sanford Underground Research Facility in South Dakota.
The large LArTPC detector at this depth will also collect atmospheric
neutrinos. In this analysis, we have used neutrino beams coming from the
accelerator as well as atmospheric neutrinos.

###  5.1 Events from beam neutrinos

In the DUNE experiment, we use a beam-power of 1.2MW leading to a total
exposure of $10\times 10^{-21}$ POT. The neutrino beam simulation for the DUNE
experiment has been carried out using the GLoBESHuber  _et al._ (2005)
software and with the most recent DUNE configuration file provided by the
CollaborationAlion  _et al._ (2016). We assume DUNE to be running for 3.5
years each in the neutrino, and the antineutrino mode with a total exposure of
280 kt-yr.

We have plotted the electron and muon events spectrum for 1300 km baseline at
DUNE considering normal hierarchy with sterile mixing angle of
$\theta_{14},\theta_{24}=7^{\circ}$ at fixed phases
$\delta_{13}=-90^{\circ},\delta_{14}=90^{\circ}$ in Figure 9. There are
differences between the spectra of the events for the true value of
$\theta_{23}=41^{\circ}$(green) in LO with the values of $\theta_{23}$ in HO
for $46^{\circ}$(orange), $50^{\circ}$(blue). This is indicative of the octant
sensitivity. It should be noted that although the green spectrum is closer to
the orange one($46^{\circ}$) for electron events (left panels), for muon
events (right panels) the green one is closer to the blue one($50^{\circ}$).
This indicates that the maximum sensitivity occurs at different $\theta_{23}$
values in the opposite octant for electron and muon events. This will lead to
the synergy between electron and muon events when we compute the combined
octant sensitivity at $\chi^{2}$ level. The maximum difference in events is
observed in the energy region of 2-4 GeV where the DUNE events spectra have
maxima in the case of both electrons and muons.

Figure 9: Electron (left) and muon neutrino (right) event spectrum for
neutrinos (top) and anti-neutrinos (bottom) as a function of energy for true
$\theta_{23}=41^{\circ}$(green) with true phases
$\delta_{13}=-90^{\circ},\delta_{14}=90^{\circ}$ at 1300 km for test values of
$\theta_{23}=46^{\circ}$(blue) and $\theta_{23}=50^{\circ}$(orange) for NH.

We present a bi-events plot considering electron neutrino and anti-neutrino
events for varying phases for normal hierarchy in Figure 10. This figure shows
that in the case of three flavour oscillation framework the ellipses for
$\theta_{23}$ being in two different octants are well separated showing no
octant degeneracy with combined $\nu_{e}+\bar{\nu}_{e}$ events of 3.5+3.5
years with 40 kt LArTPC detector at DUNE. Now if we add a sterile neutrino,
these ellipses turn into blobs which is a combination of many
ellipsesAgarwalla  _et al._ (2017). Now the interesting point to be noted is
these colorful green (LO) and purple (HO) blobs coincide when the sterile
angles have higher values of $10^{\circ}$, but blobs are slightly separated
when the angles are of $7^{\circ}$. Higher sterile mixing angles change
neutrino and anti-neutrino events spread more to make the two blobs coincide.
So we can conclude that the extent of octant degeneracy depends on values of
the sterile-active mixing angles, i.e., higher values of angles will give a
larger extent of degeneracy. The current upper bounds on
$\theta_{14},\theta_{24}$ being around $7^{\circ}$, we can observe that it is
possible to distinguish between events due to $\theta_{23}$ value $40^{\circ}$
and $50^{\circ}$ from bi-events plot.

Figure 10: Bi-events plot in $\nu_{e}-\bar{\nu}_{e}$ plane for
$\theta_{23}=41^{\circ}$(red, green), $49^{\circ}$(blue, violet) at 1300 km
with variation of phases $\delta_{13},\delta_{14}$ corresponding to
$\theta_{14},\theta_{24}=7^{\circ}$ (left), $10^{\circ}$ (right) for NH.

###  5.2 Events from atmospheric neutrinos

The atmospheric neutrino and anti-neutrino events are obtained by folding the
relevant incident fluxes with the appropriate disappearance and appearance
probabilities, charge current (CC) cross sections, detector efficiency,
resolution, detector mass, and exposure time. The $\mu^{-}$, and $e^{-}$ event
rates in an energy bin of width $\mathrm{dE_{\nu}}$ and in a solid angle bin
of width ${\mathrm{d\Omega_{\nu}}}$ are as follows,

$\rm{\frac{d^{2}N_{\mu}}{d\Omega_{\nu}\;dE_{\nu}}=\frac{1}{2\pi}\left[\left(\frac{d^{2}\Phi_{\mu}}{d\cos\theta_{\nu}\;dE_{\nu}}\right)P_{\mu\mu}+\left(\frac{d^{2}\Phi_{e}}{d\cos\theta_{\nu}\;dE_{\nu}}\right)P_{e\mu}\right]\sigma_{CC}D_{eff}}\,.$
(42)
$\mathrm{\frac{d^{2}N_{e}}{d\Omega_{\nu}\;dE_{\nu}}=\frac{1}{2\pi}\left[\left(\frac{d^{2}\Phi_{\mu}}{d\cos\theta_{\nu}\;dE_{\nu}}\right)P_{\mu
e}+\left(\frac{d^{2}\Phi_{e}}{d\cos\theta_{\nu}\;dE_{\nu}}\right)P_{ee}\right]\sigma_{CC}D_{eff}}\,,$
(43)

Here ${\mathrm{\Phi_{\mu}}}$ and ${\mathrm{\Phi_{e}}}$ are the
$\mathrm{\nu_{\mu}}$ and ${\mathrm{\nu_{e}}}$ atmospheric fluxes respectively
obtained from Honda et.al. Gaisser and Honda (2002); $P_{\mu\mu}(P_{ee})$ and
$P_{\mu e}$ are disappearance and appearance probabilities; $\rm{\sigma_{CC}}$
is the total charge current (CC) cross-section and $\rm{D_{eff}}$ is the
detector efficiency. The $\mu^{+}$, and $e^{+}$ event rates are similar to the
above expression with the fluxes, probabilities, and cross sections replaced
by those for $\mathrm{\bar{\nu}_{\mu}}$ and ${\mathrm{\bar{\nu}_{e}}}$
respectively. For the LArTPC detector at DUNE, the energy and angular
resolution are implemented using the Gaussian resolution function as follows,

$\mathrm{R_{E_{\nu}}(E_{t},E_{m})=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(E_{m}-E_{t})^{2}}{2\sigma^{2}}\right]}\,.$
(44)
$\rm{R_{\theta_{\nu}}(\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]}\,,$
(45)

where N is a normalization constant. Here, $\rm{E_{m}}$ (${\rm{\Omega_{m}}}$),
and $\rm{E_{t}}$ (${\rm{\Omega_{t}}}$) denote the measured and true values of
energy (zenith angle) respectively. The smearing width $\sigma$ is a function
of the energy $\rm{E_{t}}$. The smearing function for the zenith angle is a
bit more complicated because the direction of the incident neutrino is
specified by two variables: the polar angle ${\rm{\theta_{t}}}$ and the
azimuthal angle ${\rm{\phi_{t}}}$. We denote both these angles together by
${\rm{\Omega_{t}}}$. The measured direction of the neutrino, with polar angle
${\rm{\theta_{m}}}$ and azimuthal angle ${\rm{\phi_{m}}}$, which together we
denote by ${\rm{\Omega_{m}}}$, is expected to be within a cone of half-angle
$\Delta\theta$ of the true direction. Assumptions of the DUNE far detector
(LArTPC) parameters are mentioned in Table 2Barger  _et al._ (2016).

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

####  5.2.1 Charge identification using muon capture in liquid argon

Magnetizing the large 40 kt (10 kt each module) LArTPC detector is difficult
and expensive, but the charge id of the muon can be identified using the
capture vs decay process of the muon inside the argon as studied previously
for the DUNE detectorTernes  _et al._ (2019). We have implemented the charge
id of the muon as follows: some 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. The
lifetime of the muon resulting from the capture and decay processes can be
written as,

$\tau=\bigl{(}\frac{1}{\tau_{cap}}+\frac{Q}{\tau_{free}})^{-1}$ (46)

where $\tau_{cap}$ is the lifetime in the capture process, $\tau_{free}$ is
the decay lifetime, and Q is the Huff correction factorSuzuki  _et al._
(1987). We can define $\mu^{-}$ capture fraction as,

$\epsilon^{cap}=\frac{\tau}{\tau_{cap}}=1-\frac{\tau}{\tau_{free}}$ (47)

We use the most precise value of $\mu^{-}$ lifetime in argonKlinskikh  _et
al._ (2008), $\mu^{-}$ capture fraction becomes $\epsilon^{cap}$= 71.9$\%$.
Electron charge identification is impossible at GeV energies and electron
events are summed for each energy and angular bin. For the sensitivity
calculation, the $\mu^{-}$ and $\mu^{+}$ are separated as follows: the
$\mu^{-}$ events selected that undergo muon capture are given by,

${N_{i,j,\mu^{-}}}^{cap}={\epsilon}^{cap}\times N_{\mu^{-}}$ (48)

and the remaining $\mu^{-}$ events are included within the $\mu^{+}$ event bin
as follows,

$N_{i,j,\mu^{+}}^{rest}=\bigl{(}1-{\epsilon}^{cap})N_{i,j,\mu^{-}}+N_{i,j,\mu^{+}}$
(49)

In Figure 11, we show the absolute differences of atmospheric events between
HO $\&$ LO in $E_{\nu}$-$\cos\theta_{\nu}$ plane for $\mu^{+}+\mu^{-}$ (left),
and $e^{+}+e^{-}$ (right). This clearly shows that the difference is larger at
the matter-resonance region as observed from the probability oscillogram plot
in Figure 3. The electron event spectrum shows a significant difference in the
energy range $2-8$ GeV for $\cos\theta_{\nu}$ range of $-0.5:-0.9$. The muon
events also contribute especially in a few parts of the energy range $3-8$ GeV
for $\cos\theta_{\nu}$ range of $-0.5:-0.9$.

Figure 11: The difference of atmospheric events between HO and LO has been
plotted in $E_{\nu}-\cos\theta_{\nu}$ plane for $e^{+}+e^{-}$(left), and
$\mu^{+}+\mu^{-}$(right) events.

###  5.3 $\chi^{2}$ analysis

The computation of $\chi^{2}$ is performed using the method of pulls. This
method allows us to take into account the various statistical and systematic
uncertainties in a straightforward way. The flux, cross sections and other
systematic uncertainties are included by allowing these inputs to deviate 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}}}$, where ${\mathrm{\sigma_{k}}}$ is its
uncertainty. Then 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}}\,,$
(50)

where ${\mathrm{N^{th}_{ij}(std)}}$ is the expected rate in the
$\text{i-j}^{\text{th}}$ bin calculated with the standard values of the inputs
and npull=5 is the number of sources of uncertainty. The
${\mathrm{\xi_{k}}}$’s are called 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. (50), ${\mathrm{c_{ij}^{k}}}$ is 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). Since the
uncertainties in the inputs are not very large, we only consider changes in
${\mathrm{N^{th}_{ij}}}$ that is linear in ${\mathrm{\xi_{k}}}$. Thus we have
the modified $\chi^{2}$ as,

$\mathrm{{\chi^{2}(\xi_{k})}=\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}}+\sum^{npull}_{k=1}\;\xi_{k}^{2}}\,,$
(51)

where the additional ${\rm{\xi_{k}^{2}}}$-dependent term is the penalty
imposed for moving the value of the ${\mathrm{k^{th}}}$ input away from its
standard value by ${\rm{\sigma_{k}\;\xi_{k}}}$. The $\chi^{2}$ with pulls,
which includes the effects of all theoretical and systematic uncertainties (as
mentioned in Table 2), is obtained by minimizing ${\rm{\chi^{2}(\xi_{k})}}$
with respect to all the pulls ${\rm{\xi_{k}}}$ as follows,

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

In the case of a LArTPC detector without charge-id and with change-id,
$\chi^{2}$ is defined as,

$\chi^{2}_{w/o\text{ }charge-
id}=\chi^{2}_{\mu^{-}+\mu^{+}}+\chi^{2}_{e^{-}+e^{+}}$ (53) $\chi^{2}_{charge-
id}=\chi^{2}_{\mu^{-}}+\chi^{2}_{\mu^{+}}+\chi^{2}_{e^{-}+e^{+}}$ (54)

Finally, $\Delta{\chi^{2}}$ is marginalized over the oscillation parameters as
mentioned in Table 3.

Parameter | True Value | Marginalization Range
---|---|---
$\theta_{12}$ | $33.47^{\circ}$ | N.A.
$\theta_{13}$ | $8.54^{\circ}$ | N.A.
$\theta_{23}$ | $49^{\circ}(41^{\circ})$ | $39^{\circ}:44^{\circ}(46^{\circ}:51^{\circ})$
$\theta_{14},\theta_{24}$ | $7^{\circ}$ | $3^{\circ}:9^{\circ}$
$\Delta_{21}$ |  $7.42\times 10^{-5}$ eV2 | N.A.
$\Delta_{31}$ |  $2.515\times 10^{-3}$ eV2 | N.A.
$\Delta_{41}$ |  $1$ eV2 | N.A.
$\delta_{13},\delta_{14}$ | many | $-180^{\circ}:180^{\circ}$
Table 3: True values of all the oscillation parameters and their range of
marginalization.

##  6 Results and Discussion

Here, we present sensitivities to the octant of $\theta_{23}$ at the DUNE
experiment, assuming the 3+1 scenario. The results are demonstrated for beam
only, atmospheric only, and a combination of both of these simulated data for
fixed NH. We also explain the underlying degeneracies through the contour
plots of octant sensitivity in $\delta_{13}-\delta_{14}$ test plane.

Figure 12: Sensitivity to the octant of $\theta_{23}$ with beam only analysis
as a function of $\delta_{13}^{true}$ for true values of
$\delta_{14}=0^{\circ}$ (left), $90^{\circ}$ (right). The representative plots
are shown for HO (red), and LO (blue) for 7 years (solid) and 10 years of
simulated data (dashed). Figure 13: Sensitivity to the octant of $\theta_{23}$
as a function of $\delta_{13}^{\rm{true}}$ at
$\delta_{14}^{\rm{true}}=0^{\circ}$ for $\theta_{23}^{\rm{true}}=41^{\circ}$
(left) and $49^{\circ}$ (right). The representative plots are shown for
simulated data from beam only(red), atmospheric only w/o charge-id (blue),
atmospheric only with charge-id (violet), beam+atmospheric w/o charge-
id(green), and beam+atmospheric with charge-id(yellow) analysis with 280 kt-yr
exposure.

In the Figure 12, the sensitivity to the octant of $\theta_{23}$ degeneracy
$(\Delta\chi^{2})$ has been plotted as a function of true $\delta_{13}$ for
NH. The marginalised $\Delta\chi^{2}$ value for true
$\theta_{23}=41^{\circ}$(blue), $49^{\circ}$(red) has been shown for true
$\delta_{14}=0^{\circ}$ (left panel), $90^{\circ}$ (right panel). We present
the results for two sets of run-time; 3.5+3.5 years (solid) and 5+5 years
(dashed) in neutrino and anti-neutrino mode. We can observe the following
points from Figure 12,

  * •

The sensitivity for $\theta_{23}$ is prominently higher for LO as compared to
HO for most of the $\delta_{13}^{\rm{true}}$ values.

  * •

The $\Delta\chi^{2}$ vs $\delta_{13}$ curve has strikingly different features
for different $\delta_{14}^{true}$ values as can be seen from the two panels
in Figure 12.

  * •

For $\delta_{14}^{true}=0^{\circ}$ and LO the highest sensitivity comes around
$\delta_{13}=\pm 120^{\circ}$. This feature can be understood from Figure 5
which shows that there is no degeneracy in $P_{\mu e}(P_{\bar{\mu}\bar{e}})$
channel at $\delta_{13}=120^{\circ}(-120^{\circ})$.

  * •

On the other hand for $\delta_{14}^{true}=90^{\circ}$ the maximum sensitivity
occurs for $\delta_{13}=-90^{\circ}$. From the red dashed curves depicted in
the bottom panels of Figure 5, we can see that this sensitivity comes from
$P_{\bar{\mu}\bar{e}}$ channel.

  * •

For HO and $\delta_{14}^{true}=0^{\circ}$ the octant sensitivity is higher
around the range $\delta_{13}=-60^{\circ}:60^{\circ}$. From the solid green
curve drawn in the top panels of Figure 5, we can see that there is no
degeneracy in the range $-120^{\circ}:0^{\circ}(0^{\circ}:120^{\circ})$ comes
from $P_{\mu e}(P_{\bar{\mu}\bar{e}})$ channel with a maximum difference
between the HO curve and the LO band occurring at
$\delta_{13}=-60^{\circ}(60^{\circ})$.

  * •

In case of $\delta_{14}=90^{\circ}$ in HO, the highest sensitivity is at
$\delta_{13}=90^{\circ}$. From the top panel in Figure 5, it can be seen that
is no degeneracy in $P_{\bar{\mu}\bar{e}}$ around $\delta_{13}=90^{\circ}$.

In the above discussion, we try to explain the salient features of Figure 12
in terms of the probabilities plotted in Figure 5 for an energy of 2.5 GeV.
However, it should be borne in mind that DUNE has a broadband beam, and
contributions from other energy bins also influence the $\Delta\chi^{2}$.

In Figure 13, we have shown the sensitivity to the octant of $\theta_{23}$ for
atmospheric neutrinos without and with partial charge id of muon events(blue
and violet curves respectively) as well as combining both beam and atmospheric
data (green and orange curves) using the 40 kt DUNE far detector. In the
figure, we also present the $\Delta\chi^{2}$ for beam only data (red curve).
These plots are obtained for true values of $\theta_{23}=41^{\circ}$(left),
$49^{\circ}$(right) respectively. Here are the observations from Figure 13,

  * •

The sensitivity for atmospheric data is less than $2\sigma$ for HO and
slightly higher than $2\sigma$ for LO for whole $\delta_{13}^{\rm{true}}$
parameter space.

  * •

For the case including charge id, the sensitivity increases slightly. In
matter $P_{\mu\mu}$, and $P_{\bar{\mu}\bar{\mu}}$ probabilities are very
different due to the presence of resonant matter effect in $P_{\mu\mu}$ since
we are considering normal hierarchy. This leads to a synergy when neutrino and
anti-neutrino $\chi^{2}$ are added separately and enhances the sensitivity.

  * •

Combining atmospheric and beam data, the sensitivity increases up to more than
$4\sigma(3\sigma)$ for LO(HO) depending on the values of $\delta_{13}^{true}$.

  * •

The $\Delta\chi^{2}$ for atmospheric data has very less dependence on
$\delta_{13}^{true}$. Therefore in the combined case, the nature of
$\Delta\chi^{2}$ is mostly dictated by the beam data.

$\theta_{23}$ | $\delta_{14}$ | Above $2\sigma$ | Above $3\sigma$ | Above $2\sigma$ | Above $3\sigma$ | Above $2\sigma$ | Above $3\sigma$
---|---|---|---|---|---|---|---
Beam+Atmospheric w/o(with) charge-id | Beam
True Value | 3.5+3.5 Years | 5+5 Years
$41^{\circ}$ | $0^{\circ}$ |  $100\%$(100%) |  $100\%$(100%) | $100\%$ | $46\%$ | $100\%$ | $62\%$
$49^{\circ}$ | $0^{\circ}$ |  $100\%$(100%) |  $38\%$(53%) | $42\%$ | $0\%$ | $65\%$ | $8\%$
$41^{\circ}$ | $90^{\circ}$ |  $100\%$(100%) |  $100\%$(100%) | $100\%$ | $32\%$ | $100\%$ | $75\%$
$49^{\circ}$ | $90^{\circ}$ |  $100\%$(100%) |  $30\%$(48%) | $100\%$ | $0\%$ | $100\%$ | $23\%$
Table 4: The percentages of $\delta_{13}^{\rm{true}}$ parameter space that has
$\chi^{2}$ value above $2\sigma,3\sigma$ for various combination of true
values of $\theta_{23},\delta_{14}$ as seen in Figure 12, Figure 13.

The percentage of values of $\delta_{13}^{\rm{true}}$ for which
$\Delta\chi^{2}$ for octant sensitivity is above $2\sigma$, and $3\sigma$ are
shown in the above Table 4.

  * •

The percentage of values of the $\delta_{13}^{\rm{true}}$ for which $3\sigma$
sensitivity is achieved, is higher for $\theta_{23}^{\rm{true}}$ in lower
octant than in higher octant.

  * •

The sensitivity for $\theta_{23}^{\rm{true}}=41^{\circ}$ (LO) is more than
$3\sigma$ for $46\%(32\%)$ values of the $\delta_{13}^{\rm{true}}$ for
$\delta_{14}^{\rm{true}}=0^{\circ}(90^{\circ})$ with 7 years of beam only
data. However, in case of $\theta_{23}^{\rm{true}}=49^{\circ}$ (HO) $3\sigma$
sensitivity isn’t observed for any values of $\delta_{13}^{\rm{true}}$ as
$2\sigma$ sensitivity is achieved for $42\%(100\%)$ values of the
$\delta_{13}^{true}$ for $\delta_{14}^{\rm{true}}=0^{\circ}(90^{\circ})$.

  * •

The octant sensitivity shows a significant increment for 10 years of beam only
data while the nature of the sensitivity w.r.t. $\delta_{13}^{\rm{true}}$
remains the same. Therefore, $3\sigma$ sensitivity is observed for a
$8\%(23\%)$ parameter space of $\delta_{13}^{\rm{true}}$ even for
$\theta_{23}^{\rm{true}}=49^{\circ}$ as well as $62\%(75\%)$ of parameter
space for $\theta_{23}^{\rm{true}}=41^{\circ}$ corresponding to
$\delta_{14}^{\rm{true}}=0^{\circ}(90^{\circ})$.

  * •

For the combination of 7 years of both the beam and the atmospheric data (w/o
charge-id), the sensitivity for $\theta_{23}=49^{\circ}$ increases to more
than 3$\sigma$ for 38$\%(30\%)$ values of the $\delta_{13}^{\rm{true}}$ while
for $41^{\circ}$ the whole $\delta_{13}^{\rm{true}}$ parameter space is
allowed.

  * •

When we use the combined data for beam, and atmospheric neutrinos with charge-
id, the sensitivity improves further to provide more than 3$\sigma$ for all
$\delta_{13}^{\rm{true}}$ values $\theta_{23}=41^{\circ}$ and for $53\%(48\%)$
of $\delta_{13}^{\rm{true}}$ values corresponding to
$\delta_{14}^{\rm{true}}=0^{\circ}(90^{\circ})$ for $\theta_{23}=49^{\circ}$.

Figure 14: Octant sensitivity as a function of $\theta_{23}^{\rm{test}}$ from
beam (left), and atmospheric (right) neutrinos using 280 kt-yr exposure of
LArTPC detector with
$\theta_{23}^{tr}=41^{\circ},\delta_{13}^{tr}=-90^{\circ},\delta_{14}^{tr}=90^{\circ}$.

One of the noteworthy features of a liquid argon detector is its sensitivity
to both electron and muon events. In order to explore if there is any synergy
between these, we show in Figure 14 how $\chi^{2}$ for octant sensitivity from
muon (red) and electron events (blue) varies with $\theta_{23}^{\rm{test}}$
for true values of
$\theta_{23}=41^{\circ},\delta_{13}=-90^{\circ},\delta_{14}=90^{\circ}$ for
beam (left) and atmospheric (right) neutrinos. The observations from Figure 14
are as follows,

  * •

In the case of beam neutrinos, the octant sensitivity for appearance channel
increases with $\theta_{23}^{\rm{test}}$ whereas the sensitivity for
disappearance channel mimics the nature of $\sin^{2}2\theta_{23}$ with minima
at $41^{\circ}$, and $50^{\circ}$. This different feature of octant
sensitivity for $P_{\mu e},P_{\mu\mu}$ channels can be seen in Figure 9. When
we combine these two channels, the position of minimum sensitivity at
$\theta_{23}^{\rm{test}}=50^{\circ}$ is still guided by muon events but due
rising nature of electron $\chi^{2}$ a large octant sensitive contribution
gets added and increases the overall value of the $\chi^{2}$.

  * •

For atmospheric neutrinos, both muon, and electron $\chi^{2}$ are similar. The
muon $\chi^{2}$ is dictated by probabilities $P_{\mu\mu},P_{e\mu}$, and the
octant sensitivity coming from these channels is opposite which dilutes the
sensitivity for muons. On the other hand, for electron events, the octant
sensitivity comes from only $P_{\mu e}$ since $P_{ee}$ doesn’t depend on
$\theta_{23}$. Therefore, even though atmospheric $\nu_{\mu}$ flux is almost
twice as $\nu_{e}$ flux, both muon and electron events can give similar values
of $\chi^{2}$. These features were also noted in three flavour case in
Chatterjee  _et al._ (2013).

In order to understand the $\theta_{23}$-$\delta_{13}$-$\delta_{14}$
degeneracies listed in Table 1, we have provided the contour plots in
$\delta_{13}$-$\delta_{14}$ plane showing the regions with octant sensitivity
more than $3\sigma$. In Figure 15, the 3$\sigma$ contours are shown for the
true value of sterile CP phase $\delta_{14}=0^{\circ}$ with four different
true values of $\delta_{13}=-90^{\circ},0^{\circ},90^{\circ},150^{\circ}$. In
each panel, solid (dashed) lines represent the RO (WO) solutions. The blue,
yellow (violet, red) correspond to contours from beam only (beam and
atmospheric combined) analysis for $\theta_{23}^{true}=41^{\circ},49^{\circ}$
respectively.

Figure 15: $3\sigma$ contour plot of sensitivity to the octant of $\theta_{23}$ in test $\delta_{13}-\delta_{14}$ plane with 7 years of data for $\delta_{14}^{true}=0^{\circ}$ and $\delta_{13}^{true}=-90^{\circ},0^{\circ},90^{\circ},150^{\circ}$ in panels a,b,c,d respectively. The representative plots are shown for the true value of $\theta_{23}=41^{\circ}$ in LO (blue and violet) and $49^{\circ}$ in HO (yellow and red) for right octant solutions(solid) and wrong octant solutions(dashed) for simulated beam only (B) and beam+atmospheric (B+A) data. True $\delta_{13}$ | True $\delta_{14}$ | Present Degeneracies
---|---|---
$-90^{\circ}$ | $0^{\circ}$ | WO-R$\delta_{13}$-W$\delta_{14}$
$0^{\circ}$ | $0^{\circ}$ | WO-R$\delta_{13}$-R$\delta_{14}(49^{\circ})$, WO-R$\delta_{13}$-W$\delta_{14}(41^{\circ})$
$90^{\circ}$ | $0^{\circ}$ | WO-R$\delta_{13}$-W$\delta_{14}(49^{\circ})$
$150^{\circ}$ | $0^{\circ}$ | WO-R$\delta_{13}$-R$\delta_{14}(49^{\circ})$, WO-R$\delta_{13}$-W$\delta_{14}(49^{\circ})$
Table 5: The degeneracies for different true value of $\delta_{13}$ with true $\delta_{14}=0^{\circ}$ as seen in Figure 15. Figure 16: $3\sigma$ contour plot of sensitivity to the octant of $\theta_{23}$ in test $\delta_{13}-\delta_{14}$ plane with 7 years of beam only simulated data for $\delta_{14}^{true}=90^{\circ}$ and $\delta_{13}^{true}=-90^{\circ},0^{\circ},90^{\circ},150^{\circ}$ in panels a,b,c,d respectively. The representative plots shown for true value of $\theta_{23}=41^{\circ}$ in LO (blue) and $49^{\circ}$(yellow) in HO for right octant solutions(solid) and wrong octant solutions(dashed). True $\delta_{13}$ | True $\delta_{14}$ | Present Degeneracies
---|---|---
$-90^{\circ}$ | $90^{\circ}$ | WO-R$\delta_{13}$-R$\delta_{14}(49^{\circ})$, WO-R$\delta_{13}$-W$\delta_{14}(49^{\circ})$
$0^{\circ}$ | $90^{\circ}$ | WO-R$\delta_{13}$-W$\delta_{14}$
$90^{\circ}$ | $90^{\circ}$ | WO-W$\delta_{13}$-W$\delta_{14}(49^{\circ})$, WO-R$\delta_{13}$-W$\delta_{14}(41^{\circ})$
$150^{\circ}$ | $90^{\circ}$ | WO-R$\delta_{13}$-W$\delta_{14}$
Table 6: The degeneracies for different true value of $\delta_{13}$ with true
$\delta_{14}=90^{\circ}$ as seen in Figure 16.

The noteworthy observations from Figure 15 are as follows,

  * •

In panel "a", the solid contours spanning the full range of $\delta_{14}$
indicate true solutions with poor precision in $\delta_{14}$ for both
$\theta_{23}^{true}=41^{\circ},49^{\circ}$. We also observe dashed contours
indicating WO-R$\delta_{13}$-W$\delta_{14}$ solutions for both
$\theta_{23}^{true}$.

  * •

In panel "b", the precision of the true solutions improves significantly. A
small region of WO solutions for $\theta_{23}^{true}=49^{\circ}$ occurs
adjacent to the true value. We also find WO-R$\delta_{13}$-W$\delta_{14}$
solutions for $\theta_{23}^{true}=41^{\circ}$.

  * •

Comparing the true solutions in panels "c", and "d" but the precision of
$\delta_{14}$ is notably better in "d". In these panels, WO solutions are
present for only $\theta_{23}^{true}=49^{\circ}$. For
$\theta_{23}^{true}=41^{\circ}$, the octant can be determined at more than
$3\sigma$ sensitivity as seen from the solid blue curve in the left panel of
Figure 12 and hence WO solutions are not observed. In panel "c" we find
WO-R$\delta_{13}$-W$\delta_{14}$ solution wheres the WO-R$\delta_{13}$
solutions are observed in panel "d".

  * •

Inclusion of atmospheric analysis shrinks all the contours improving octant
sensitivity. The choice of $\delta_{13}^{true}$ affects the precision of RO
solutions as well as the occurrence of degeneracies.

Similarly, we have plotted the $3\sigma$ contours in Figure 16 showing WO
(dashed), and RO (solid) solutions w.r.t. true values of
$\theta_{23}=41^{\circ}$ (blue), and $49^{\circ}$ (yellow) for the true value
of $\delta_{14}=90^{\circ}$ with
$\delta_{13}=-90^{\circ},0^{\circ},90^{\circ},150^{\circ}$ using beam-only
analysis. The observations from Figure 16 are as follows,

  * •

In panel "a", we see the WO-R$\delta_{13}$ solutions spanning the full range
of $\delta_{14}$ for only $\theta_{23}^{true}=49^{\circ}$. We also find true
solutions with notable precision in $\delta_{14}$ for both
$\theta_{23}^{true}=41^{\circ},49^{\circ}$ as compared to panel "a" in Figure
15.

  * •

In panel "b", the precision of $\delta_{14}$ in true solutions deteriorates
w.r.t panel "a" covering the full $\delta_{14}$ range. We observe a small
region of WO-R$\delta_{13}$-W$\delta_{14}$ solution for
$\theta_{23}^{true}=49^{\circ}$, along with a bigger region of
WO-R$\delta_{13}$-W$\delta_{14}$ solution for $\theta_{23}^{true}=41^{\circ}$.

  * •

In panels "c" and "d", the true solutions show better precision in
$\delta_{14}$ as compared to the same panels in Figure 15. We can also observe
for $\theta_{23}^{true}=49^{\circ}$ a tine region of
WO-W$\delta_{13}$-W$\delta_{14}$ in panel "c" while in panel "d"
WO-R$\delta_{13}$-W$\delta_{14}$ solutions occur. There are
WO-R$\delta_{13}$-W$\delta_{14}$ solutions for $\theta_{23}^{true}=41^{\circ}$
in both panel "c", and "d" but the region is smaller in "c".

  * •

Overall, we see the precision of the RO true solutions along with the size and
type of WO contours depend on $\delta_{13}^{true}$ for fixed
$\delta_{14}^{true}$.

The most common degeneracies seen in Figure 15, Figure 16 are
WO-R$\delta_{13}$-R$\delta_{14}$, WO-R$\delta_{13}$-W$\delta_{14}$. It
indicates that the presence of $\delta_{14}$ creates more problems in precise
measurement of the octant of $\theta_{23}$. We also observe true solutions
with poor precision in $\delta_{14}$.

The regions under $3\sigma$ sensitivity in the contour plots of Figure 15,
Figure 16 can be understood using the difference in the probability plots in
$\delta_{13}-\delta_{14}$ plane. We will mainly focus on the dominant $P_{\mu
e}$ channel to understand the effect.

Figure 17: Contour plot in test $\delta_{13}-\delta_{14}$ plane showing the
difference in probability $\Delta P_{\mu e}$ with $\theta_{23}$ being fixed at
one octant while $\theta_{23}$ varies in the opposite octant for WO solutions
(right) and in the same octant for RO solutions (left) at true values of
$\delta_{13}=-90^{\circ},\delta_{14}=0^{\circ},\theta_{23}=49^{\circ}$(top),
$41^{\circ}$(bottom). Black and dark red show the least differences while blue
and white show the highest.

In Figure 17, the contour plot in test $\delta_{13}$-$\delta_{14}$ plane
represents the difference between the probabilities $P_{\mu e}$ of opposite
octants while varying the $\theta_{23}$ value only in same (left) /opposite
(right) octant for the true $\theta_{23}=49^{\circ}$(top),
$41^{\circ}$(bottom) with
$\delta_{13}^{tr}=-90^{\circ},\delta_{14}^{tr}=0^{\circ}$ corresponding to
panel "a" of Figure 15. The understandings are as follows,

  * •

First, we consider the right octant solutions in the panels at the left side
column. It can be clearly seen that the black and darker red regions around
the true value on the left side of $\delta_{13}-\delta_{14}$ plane where the
difference in the probability is minimum in Figure 17 is similar to the
$3\sigma$ regions under the solid curves in panel "a" of Figure 15. These
darker regions also indicate poor precision of $\delta_{14}$.

  * •

For $49^{\circ}$-WO solution, minima arise in the darker red region including
the true value in the top-right panel of Figure 17 similar to the yellow-
dashed contour in panel "a" of Figure 15. Similarly, for $41^{\circ}$-WO
solutions in the bottom-right panel, the minimum difference is observed in the
darker red region just above the true value similar to the blue-dashed contour
in panel "a" of Figure 15. These darker red regions clearly show precise
WO-R$\delta_{13}$-W$\delta_{14}$ degenerate solutions.

##  7 Conclusions

We study the octant sensitivity in presence of an additional sterile neutrino
using beam and atmospheric neutrinos separately as well as in a combined study
of both in the context of the LArTPC detector at DUNE. We provide the analytic
expressions for oscillation probabilities in presence of an extra sterile
neutrino using the approximation that the mass squared difference
$\Delta_{21}$ is zero. It has been shown that these expressions are a good
match with the numerical probability, especially in the resonance region. We
perform a $\chi^{2}$ analysis and show that for a set of true values of
sterile parameters one can achieve more than $3\sigma$ octant sensitivity
depending upon the true value of $\delta_{13}$ using beam neutrinos. For
instance using $3.5+3.5$ years of combined run time, for true value of
$\theta_{23}=41^{\circ}(49^{\circ}),\delta_{14}=90^{\circ}$ one gets more than
$3\sigma$ sensitivity for $51\%(18\%)$ of the $\delta_{13}^{true}$ space. If
the statistics are increased to $5+5$ years, one gets more than $3\sigma$
sensitivity for $69\%(28\%)$ of the $\delta_{13}^{true}$ space.

For only atmospheric neutrinos with $280$ kt-yr exposure, for true values of
$\theta_{23}=41^{\circ},\delta_{14}=0^{\circ}$, more than $3\sigma$
sensitivity can be achieved for all values of the true $\delta_{13}$. The
inclusion of the charge tagging capability of muon capture in argon helps us
in getting an increased sensitivity. But for true values of
$\theta_{23}=49^{\circ},\delta_{14}=0^{\circ}$ we can not achieve $3\sigma$
sensitivity even with charge id using the atmospheric neutrinos only.

Combining the beam and the atmospheric neutrinos, we can obtain $3\sigma$
sensitivity in the entire true $\delta_{13}$ space for true values of
$\theta_{23}=41^{\circ},\delta_{14}=0^{\circ}$. In case of true values of
$\theta_{23}=49^{\circ},\delta_{14}=0^{\circ}$, $3\sigma$ sensitivity is
achieved for $55\%$ of the true $\delta_{13}$ values.

At fixed hierarchy, there can be a total of 8-fold degeneracies (Table 1) with
at least one of the octant, $\delta_{13}$, $\delta_{14}$ assuming a wrong
value. We also identify the extra degeneracies due to the presence of
$\delta_{14}$ assuming the normal hierarchy and summarise these in Table 5,
and Table 6. We can also conclude that presence of the phase $\delta_{14}$
leads to the occurrence of new degeneracies that hinders the discovery of the
octant of $\theta_{23}$ precisely.

In summary, the combination of beam and atmospheric neutrinos provides
promising results in the LArTPC detector at DUNE in the presence of sterile
neutrinos. We can achieve $3\sigma$ sensitivity for octant of $\theta_{23}$
corresponding to true values of $41^{\circ},49^{\circ}$ in both octants.

###### Acknowledgements.

We would like to express special acknowledgment to Jaydeep Dutta, Sushant
Raut, and Samiron Roy for their involvement in the earlier stages. We
acknowledge discussions with Monojit Ghosh. SP acknowledges Kaustav
Chakraborty for help in learning numerical analysis using GLoBES. AC
acknowledges the University of Pittsburgh where initial work has been done. AC
also acknowledges Ramanujan Fellowship (grant no: RJF2021/000157), of the
Science and Engineering Research Board of the Department of Science and
Technology, Government of India. 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.

##  Appendix A Probability calculation using Cayley Hamilton formalism

We will now find out the analytic probability using the Cayley-Hamilton
formalismOhlsson and Snellman (2000); Kamo  _et al._ (2003); Akhmedov  _et
al._ (2004). We calculate the time evolution operator and do not introduce
auxiliary matter mixing angles.

The flavour eigenstates $\psi_{\alpha}$ and mass eigenstates $\psi_{i}$ are
related as

$\psi_{i}=\sum_{j=e,\mu,\tau,s}U^{\star}_{\alpha j}\psi_{j}$ (55)

where $U_{\alpha j}$ is component of unitary mixing matrix corresponding to
mixing between $\psi_{\alpha},\psi_{j}$,

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

The Schrodinger equation in mass basis is given as,

$\iota\dfrac{d}{dt}\psi_{m}(t)=\mathcal{H}_{m}\psi_{m}(t)$ (57)

where total Hamiltonian $\mathcal{H}_{m}$ in mass basis, and interaction
Hamiltonian $V_{f}$ in flavour basis are given as follows,

$\displaystyle\mathcal{H}_{m}=H_{m}+U^{-1}V_{f}U$ (58) $\displaystyle
V_{f}=H_{int}=diag(2A^{\prime},0,0,A^{\prime})$ (59)

Equation (57) gives the solution with time evolution operator
$e^{-\iota\mathcal{H}_{m}t}$ as,

$\psi_{m}(t)=e^{-\iota\mathcal{H}_{m}t}\psi_{m}(0)$ (60)

We get the solution in terms of distance $L$ traveled by neutrinos in time $t$
as,

$\psi_{m}(L)=\psi_{m}(t=L)=e^{-\iota\mathcal{H}_{m}t}\psi_{m}(0)\equiv
U_{m}(L)\psi_{m}(0)$ (61)

Solution in flavour state $\psi_{f}$ is expressed at $t=L$ as,

$\psi_{f}(L)=U\phi_{m}(L)=Ue^{-\iota\mathcal{H}_{m}t}U^{-1}U\psi_{m}(0)=Ue^{-\iota\mathcal{H}_{m}t}U^{-1}\psi_{f}(0)\equiv
U_{f}(L)\psi_{f}(0)$ (62)

We will calculate the time evolution operator, i.e., the exponential of the
matrix $\mathcal{H}_{m}$ using Cayley-Hamilton theorem. We construct traceless
matrix out of $\mathcal{H}_{m}$ as ,

$\mathcal{H}_{m}=T+\dfrac{1}{4}(tr\mathcal{H}_{m})I$ (63)

The time evolution operator is then redefined as,

$U_{m}(L)=e^{-\iota\mathcal{H}_{m}L}=\phi e^{-\iota TL}$ (64)

The elements of traceless matrix $T$ in mass basis are as follows,

$\displaystyle T_{11}$ $\displaystyle=$ $\displaystyle
A\left[-\cos^{2}\theta_{12}\left(2\sin\theta_{13}\cos\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\cos(\delta_{13}-\delta_{14})+\cos
2\theta_{23}\sin^{2}\theta_{24}\right)\right.$ (65)
$\displaystyle\left.+2\sin\theta_{12}\cos\theta_{12}\cos\theta_{23}\sin\theta_{24}(\cos\delta_{13}\sin\theta_{13}\sin\theta_{23}\sin\theta_{24}-\cos\delta_{14}\cos\theta_{13}\sin\theta_{14}\cos\theta_{24})\right.$
$\displaystyle\left.+\cos^{2}\theta_{12}\cos^{2}\theta_{13}\left(2-\sin^{2}\theta_{24}\left(\sin^{2}\theta_{14}+\sin^{2}\theta_{23}\right)-\sin^{2}\theta_{14}\right)\right.$
$\displaystyle\left.+\cos^{2}\theta_{23}\sin^{2}\theta_{24}\right]-\frac{3A}{4}+\frac{1}{4}(-\Delta_{21}-\Delta_{31}-\Delta_{41})$
$\displaystyle T_{12}$ $\displaystyle=$ $\displaystyle
A\left[-\sin\theta_{12}\cos\theta_{12}\left(2\sin\theta_{13}\cos\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\cos(\delta_{13}-\delta_{14})+\cos
2\theta_{23}\sin^{2}\theta_{24}\right)\right.$ (66)
$\displaystyle\left.-\sin\theta_{13}\sin\theta_{23}\cos\theta_{23}\sin^{2}\theta_{24}\left(e^{-i\delta_{13}}\cos^{2}\theta_{12}-e^{i\delta_{13}}\sin^{2}\theta_{12}\right)\right.$
$\displaystyle\left.+\cos\theta_{13}\sin\theta_{14}\cos\theta_{23}\sin\theta_{24}\cos\theta_{24}\left(e^{-i\delta_{14}}\cos^{2}\theta_{12}-e^{i\delta_{14}}\sin^{2}\theta_{12}\right)\right.$
$\displaystyle\left.+\sin\theta_{12}\cos\theta_{12}\cos^{2}\theta_{13}\left(2-\sin^{2}\theta_{24}\left(\sin^{2}\theta_{14}+\sin^{2}\theta_{23}\right)-\sin^{2}\theta_{14}\right)\right]$
$\displaystyle T_{13}$ $\displaystyle=$ $\displaystyle
A\left[-e^{i\delta_{14}-2i\delta_{13}}\cos\theta_{12}\sin^{2}\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\right.$
(67)
$\displaystyle\left.-e^{i\delta_{14}-i\delta_{13}}\sin\theta_{12}\sin\theta_{13}\sin\theta_{14}\cos\theta_{23}\sin\theta_{24}\cos\theta_{24}\right.$
$\displaystyle\left.+e^{-i\delta_{14}}\cos\theta_{12}\cos^{2}\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\right.$
$\displaystyle\left.+e^{-i\delta_{13}}\cos\theta_{12}\sin\theta_{13}\cos\theta_{13}\left(2-\sin^{2}\theta_{24}\left(\sin^{2}\theta_{14}+\sin^{2}\theta_{23}\right)-\sin^{2}\theta_{14}\right)\right]$
$\displaystyle-A\sin\theta_{12}\cos\theta_{13}\sin\theta_{23}\cos\theta_{23}\sin^{2}\theta_{24}$
$\displaystyle T_{14}$ $\displaystyle=$ $\displaystyle
A\left[e^{-i\delta_{13}}\cos\theta_{12}\sin\theta_{13}\cos\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\right.$
(68)
$\displaystyle\left.+e^{-i\delta_{14}}\cos\theta_{12}\cos\theta_{13}\sin\theta_{14}\cos\theta_{14}\left(2-\cos^{2}\theta_{24}\right)+\sin\theta_{12}\cos\theta_{14}\cos\theta_{23}\sin\theta_{24}\cos\theta_{24}\right]$
$\displaystyle T_{22}$ $\displaystyle=$ $\displaystyle
A\left[-\sin^{2}\theta_{12}\left(2\sin\theta_{13}\cos\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\cos(\delta_{13}-\delta_{14})+\cos
2\theta_{23}\sin^{2}\theta_{24}\right)\right.$ (69)
$\displaystyle\left.+2\sin\theta_{12}\cos\theta_{12}\cos\theta_{23}\sin\theta_{24}(\cos\delta_{14}\cos\theta_{13}\sin\theta_{14}\cos\theta_{24}-\cos\delta_{13}\sin\theta_{13}\sin\theta_{23}\sin\theta_{24})\right.$
$\displaystyle\left.+\sin^{2}\theta_{12}\cos^{2}\theta_{13}\left(2-\sin^{2}\theta_{24}\left(\sin^{2}\theta_{14}+\sin^{2}\theta_{23}\right)-\sin^{2}\theta_{14}\right)\right.$
$\displaystyle\left.+\cos^{2}\theta_{23}\sin^{2}\theta_{24}\right]-\frac{3A}{4}+\frac{1}{4}(3\Delta_{21}-\Delta_{31}-\Delta_{41})$
$\displaystyle T_{23}$ $\displaystyle=$ $\displaystyle
A\left[-e^{i\delta_{14}-2i\delta_{13}}\sin\theta_{12}\sin^{2}\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\right.$
(70)
$\displaystyle\left.+e^{i\delta_{14}-i\delta_{13}}\cos\theta_{12}\sin\theta_{13}\sin\theta_{14}\cos\theta_{23}\sin\theta_{24}\cos\theta_{24}\right.$
$\displaystyle\left.+e^{-i\delta_{14}}\sin\theta_{12}\cos^{2}\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\right.$
$\displaystyle\left.+e^{-i\delta_{13}}\sin\theta_{12}\sin\theta_{13}\cos\theta_{13}\left(2-\sin^{2}\theta_{24}\left(\sin^{2}\theta_{14}+\sin^{2}\theta_{23}\right)-\sin^{2}\theta_{14}\right)\right]$
$\displaystyle+A\cos\theta_{12}\cos\theta_{13}\sin\theta_{23}\cos\theta_{23}\sin^{2}\theta_{24}$
$\displaystyle T_{24}$ $\displaystyle=$ $\displaystyle
A\left[e^{-i\delta_{13}}\sin\theta_{12}\sin\theta_{13}\cos\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\right.$
(71)
$\displaystyle\left.+e^{-i\delta_{14}}\sin\theta_{12}\cos\theta_{13}\sin\theta_{14}\cos\theta_{14}\left(2-\cos^{2}\theta_{24}\right)\right.$
$\displaystyle\left.-\cos\theta_{12}\cos\theta_{14}\cos\theta_{23}\sin\theta_{24}\cos\theta_{24}\right]$
$\displaystyle T_{33}$ $\displaystyle=$
$\displaystyle\,A\left[2\sin\theta_{13}\cos\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\cos(\delta_{13}-\delta_{14})\right.$
(72)
$\displaystyle\left.+\sin^{2}\theta_{13}\left(2-\sin^{2}\theta_{24}\left(\sin^{2}\theta_{14}+\sin^{2}\theta_{23}\right)-\sin^{2}\theta_{14}\right)+\sin^{2}\theta_{23}\sin^{2}\theta_{24}\right]$
$\displaystyle-\frac{3A}{4}+\frac{1}{4}(-\Delta_{21}+3\Delta_{31}-\Delta_{41})$
$\displaystyle T_{34}$ $\displaystyle=$ $\displaystyle
A\left[e^{i\delta_{13}-i\delta_{14}}\sin\theta_{13}\sin\theta_{14}\cos\theta_{14}\left(2-\cos^{2}\theta_{24}\right)\right.$
(73)
$\displaystyle\left.-\cos\theta_{13}\cos\theta_{14}\sin\theta_{23}\sin\theta_{24}\cos\theta_{24}\right]$
$\displaystyle T_{44}$ $\displaystyle=$ $\displaystyle
A\left[\cos^{2}\theta_{14}\cos^{2}\theta_{24}+2A\sin^{2}\theta_{14}\right]-\frac{3A}{4}+\frac{1}{4}(-\Delta_{21}-\Delta_{31}+3\Delta_{41})$
(74)

Cayley-Hamilton theorem is used to get the form of the time evolution operator
$e^{-\iota TL}$. We need to solve the characteristic equation of matrix $T$
given by,

$\lambda^{4}+c_{3}\lambda^{3}+c_{2}\lambda^{2}+c_{1}\lambda+c_{0}=0$ (75)

to obtain the energy eigenvalues $\lambda$ where the constants are defined as
follows,

$\displaystyle c_{0}$ $\displaystyle=$
$\displaystyle\frac{A^{2}}{128}\Delta_{41}^{2}\left(8\sin^{2}\theta_{14}+29\right)+$
(76)
$\displaystyle\frac{A}{64}\left[\left(-\Delta_{31}^{3}+2\Delta_{31}^{2}\Delta_{41}+3\Delta_{31}\Delta_{41}^{2}\right)\sin
2\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin
2\theta_{24}\cos(\delta_{13}-\delta_{14})\right.$
$\displaystyle\left.+\Delta_{31}^{3}\left(3-4\cos^{2}\theta_{13}\sin^{2}\theta_{23}\sin^{2}\theta_{24}-4Q\sin^{2}\theta_{13}\right)\right.$
$\displaystyle\left.+\Delta_{31}^{2}\Delta_{41}\left(8\cos^{2}\theta_{13}\sin^{2}\theta_{23}\sin^{2}\theta_{24}+12\cos^{2}\theta_{24}-4Q\left(2\cos^{2}\theta_{13}+1\right)+9\right)\right.$
$\displaystyle\left.+\Delta_{31}\Delta_{41}^{2}\left(12\cos^{2}\theta_{13}\sin^{2}\theta_{23}\sin^{2}\theta_{24}+8\cos^{2}\theta_{24}+4Q\left(1-3\cos^{2}\theta_{13}\right)+1\right)\right.$
$\displaystyle\left.-\Delta_{41}^{3}\left(4\cos^{2}\theta_{24}-4Q+5\right)\right]+\frac{\Delta_{21}}{64}\left(\Delta_{31}^{3}-5\Delta_{31}^{2}\Delta_{41}-5\Delta_{31}\Delta_{41}^{2}+\Delta_{41}^{3}\right)$
$\displaystyle+\left(-\frac{3\Delta_{31}^{4}}{256}+\frac{\Delta_{31}^{3}\Delta_{41}}{64}+\frac{7\Delta_{31}^{2}\Delta_{41}^{2}}{128}+\frac{\Delta_{31}\Delta_{41}^{3}}{64}-\frac{3\Delta_{41}^{4}}{256}\right)$
$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{8}A^{2}\Delta_{41}\left(5-7\sin^{2}\theta_{14}\right)+\frac{A}{8}\Delta_{31}^{2}\left(3-4\sin^{2}\theta_{23}\sin^{2}\theta_{24}\cos^{2}\theta_{13}-4Q\sin^{2}\theta_{13}\right)$
(77)
$\displaystyle-\frac{A}{8}\Delta_{41}^{2}\left(5+4\cos^{2}\theta_{24}-4Q\right)+\frac{A}{16}\Delta_{31}\Delta_{41}\left(4+8\cos^{2}\theta_{24}-5P\cos^{2}\theta_{13}\right)$
$\displaystyle+\frac{A}{4}\left(-\text{$\Delta_{31}$}^{2}+\text{$\Delta_{31}$}\text{$\Delta_{41}$}\right)\sin
2\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin
2\theta_{24}\cos(\delta_{13}-\delta_{14})$
$\displaystyle+\frac{\Delta_{21}}{8}\left(\Delta_{41}-\Delta_{31}\right)^{2}+\frac{1}{8}\left(-\Delta_{31}^{3}+\Delta_{31}^{2}\Delta_{41}+\Delta_{31}\Delta_{41}^{2}-\Delta_{41}^{3}\right)$
$\displaystyle c_{2}$ $\displaystyle=$
$\displaystyle\frac{A}{4}\Delta_{31}\left(3-4\sin^{2}\theta_{23}\sin^{2}\theta_{24}\cos^{2}\theta_{13}-4Q\sin^{2}\theta_{13}\right)-\frac{A}{4}\Delta_{41}\left(5+4\cos^{2}\theta_{24}-4Q\right)$
(78) $\displaystyle-\frac{11}{8}A^{2}-\frac{A}{2}\Delta_{31}\sin
2\theta_{13}\sin\theta_{14}\sin\theta_{23}\sin
2\theta_{24}\cos(\delta_{13}-\delta_{14})$
$\displaystyle+\frac{\Delta_{21}}{4}\left(\Delta_{41}+\Delta_{31}\right)+\frac{1}{8}\left(-3\Delta_{31}^{2}+2\Delta_{31}\Delta_{41}-3\Delta_{41}^{2}\right)$
$\displaystyle c_{3}=\text{Trace(T)}=0$ (79) $\displaystyle P$
$\displaystyle=$ $\displaystyle
2-\sin^{2}\theta_{14}-\sin^{2}\theta_{24}\left(\sin^{2}\theta_{14}+\sin^{2}\theta_{23}\right)$
(80) $\displaystyle Q$ $\displaystyle=$ $\displaystyle
2-\sin^{2}\theta_{14}-\sin^{2}\theta_{24}\sin^{2}\theta_{14}$ (81)

The energy eigenvalues are as follows,

$\displaystyle\lambda_{1,2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left[\sqrt{-c_{2}+t_{0}}\pm\sqrt{-c_{2}-t_{0}-2\sqrt{-4c_{0}+t_{0}^{2}}}\right]$
(82) $\displaystyle\lambda_{3,4}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left[-\sqrt{-c_{2}+t_{0}}\pm\sqrt{-c_{2}-t_{0}+2\sqrt{-4c_{0}+t_{0}^{2}}}\right]$
(83)

where $t_{0}$ is a real root of the following equation,

$t^{3}-c_{2}t^{2}-4c_{0}t+4c_{0}c_{2}-c_{1}^{2}=0$ (84)

The general form of probability is given by

$\displaystyle
P_{\alpha\beta}=\sum_{a=1}^{4}\sum_{b=1}^{4}(\tilde{B}_{a})_{\alpha\beta}(\tilde{B}_{b})_{\alpha\beta}^{\star}e^{-\iota
L(\lambda_{a}-\lambda_{b})}$ (85)

where,

$(\tilde{B}_{a})_{\alpha\beta}=\frac{(c_{1}+c_{2}\lambda_{a}+\lambda_{a}^{3})\delta_{\alpha\beta}+(c_{2}+\lambda_{a}^{2})\tilde{T}_{\alpha\beta}+\lambda_{a}\tilde{T}^{2}_{\alpha\beta}+\tilde{T}_{\alpha\beta}^{3}}{4\lambda_{a}^{3}+c_{1}+2c_{2}\lambda_{a}}$
(86)

where components of $T,T^{2},T^{3}$ in flavour basis are given as following,

$\tilde{T}_{\alpha\beta}=<\alpha|UTU^{-1}|\beta>,\tilde{T^{2}}_{\alpha\beta}=<\alpha|UT^{2}U^{-1}|\beta>,\tilde{T^{3}}_{\alpha\beta}=<\alpha|UT^{3}U^{-1}|\beta>$
(87)

Figure 18: Comparison of the probability using GLoBES $P_{\mu e}^{GL}$ (red),
the Cayley Hamilton probability $P_{\mu e}^{CH}$ (green), and TMSD probability
$P_{\mu e}$ (blue) at 1300 km(left), 7000 km(right) baseline.

In Figure 18, we see that the Cayley Hamilton probabilities at lower energies
show a better match with numerical probabilities evaluated using GLoBES,
whereas at higher energies especially at resonance region the TMSD
probabilities match better as was also seen in Figure 1, and Figure 2.

## References

  * Esteban  _et al._ (2017) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, I. Martinez-Soler, and T. Schwetz, JHEP 01, 087 (2017), arXiv:1611.01514 [hep-ph] .
  * 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), arXiv:2006.11237 [hep-ph] .
  * Capozzi  _et al._ (2020) F. Capozzi, E. D. Valentino, E. Lisi, A. Marrone, A. Melchiorri, and A. Palazzo, Physical Review D 101 (2020), 10.1103/physrevd.101.116013.
  * Ghosh  _et al._ (2016a) M. Ghosh, P. Ghoshal, S. Goswami, N. Nath, and S. K. Raut, Phys. Rev. D 93, 013013 (2016a), arXiv:1504.06283 [hep-ph] .
  * 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.
  * Abi  _et al._ (2020a) B. Abi _et al._ (DUNE), (2020a), arXiv:2002.03005 [hep-ex] .
  * Tolba (2022) T. Tolba (ESSnuSB), PoS PANIC2021, 284 (2022).
  * Abe  _et al._ (2015) K. Abe _et al._ (Hyper-Kamiokande Proto-), PTEP 2015, 053C02 (2015), arXiv:1502.05199 [hep-ex] .
  * Brunner (2018) J. Brunner (KM3Net), PoS NEUTEL2017, 057 (2018).
  * Aartsen  _et al._ (2017) M. G. Aartsen _et al._ (IceCube), J. Phys. G 44, 054006 (2017), arXiv:1607.02671 [hep-ex] .
  * Ahmed  _et al._ (2017) S. Ahmed _et al._ (ICAL), Pramana 88, 79 (2017), arXiv:1505.07380 [physics.ins-det] .
  * Ghosh  _et al._ (2016b) M. Ghosh, S. Goswami, and S. K. Raut, Eur. Phys. J. C 76, 114 (2016b), arXiv:1412.1744 [hep-ph] .
  * 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] .
  * Hill (1995) J. E. Hill, Phys. Rev. Lett. 75, 2654 (1995), arXiv:hep-ex/9504009 .
  * Aguilar-Arevalo  _et al._ (2020) A. Aguilar-Arevalo _et al._ (MiniBooNE), (2020), arXiv:2006.16883 [hep-ex] .
  * Abrams  _et al._ (1989) G. S. Abrams _et al._ , Phys. Rev. Lett. 63, 2173 (1989).
  * Giunti and Laveder (2011) C. Giunti and M. Laveder, Phys. Rev. C 83, 065504 (2011), arXiv:1006.3244 [hep-ph] .
  * Acero  _et al._ (2009) M. A. Acero, C. Giunti, and M. Laveder, Nucl. Phys. B Proc. Suppl. 188, 211 (2009).
  * Barinov  _et al._ (2021) V. V. Barinov _et al._ , “Results from the baksan experiment on sterile transitions (best),”  (2021).
  * Huber (2011) P. Huber, Physical Review C 84 (2011), 10.1103/physrevc.84.024617.
  * Mueller  _et al._ (2011) T. A. Mueller, D. Lhuillier, M. Fallot, A. Letourneau, S. Cormon, M. Fechner, L. Giot, T. Lasserre, J. Martino, G. Mention, A. Porta, and F. Yermia, Physical Review C 83 (2011), 10.1103/physrevc.83.054615.
  * Svirida  _et al._ (2019) D. Svirida _et al._ , PoS NOW2018, 066 (2019).
  * Danilov and Skrobova (2021) M. Danilov and N. Skrobova, “New results from the danss experiment,”  (2021).
  * Ko  _et al._ (2017) Y. J. Ko _et al._ (NEOS Collaboration), Phys. Rev. Lett. 118, 121802 (2017).
  * Almazán  _et al._ (2020) H. Almazán _et al._ (STEREO Collaboration), Phys. Rev. D 102, 052002 (2020).
  * Andriamirado  _et al._ (2021) M. Andriamirado _et al._ (PROSPECT Collaboration), Phys. Rev. D 103, 032001 (2021).
  * Minotti (2022) A. Minotti, PoS NuFact2021, 246 (2022).
  * Goswami (1997) S. Goswami, Phys. Rev. D 55, 2931 (1997), arXiv:hep-ph/9507212 .
  * Abazajian  _et al._ (2012) K. N. Abazajian _et al._ , (2012), arXiv:1204.5379 [hep-ph] .
  * Adamson  _et al._ (2019) P. Adamson _et al._ (MINOS+), Phys. Rev. Lett. 122, 091803 (2019), arXiv:1710.06488 [hep-ex] .
  * Abe  _et al._ (2019) K. Abe _et al._ (T2K), Phys. Rev. D 99, 071103 (2019), arXiv:1902.06529 [hep-ex] .
  * Vannerom  _et al._ (2022) D. Vannerom, L. Fischer, J. Conrad, S. Blot, and C. Arguelles, PoS PANIC2021, 299 (2022).
  * Argüelles  _et al._ (2022) 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, Phys. Rev. Lett. 128, 241802 (2022), arXiv:2111.10359 [hep-ph] .
  * Hausner (2022) H. Hausner,  _Sterile Neutrino Search with the NOvA Detectors_ , Ph.D. thesis, Wisconsin U., Madison, SAL (2022).
  * Hu (2021) Z. Hu (MINOS, MINOS+, Daya Bay, Bugey-3), PoS ICHEP2020, 201 (2021).
  * Adamson  _et al._ (2020) P. Adamson _et al._ (MINOS+, Daya Bay), Phys. Rev. Lett. 125, 071801 (2020), arXiv:2002.00301 [hep-ex] .
  * MicroBooNE Collaboration  _et al._ (2021) MicroBooNE Collaboration, P. Abratenko, _et al._ , “Search for an excess of electron neutrino interactions in microboone using multiple final state topologies,”  (2021).
  * Abratenko  _et al._ (2022) P. Abratenko _et al._ (MicroBooNE), (2022), arXiv:2210.10216 [hep-ex] .
  * 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, “Microboone and the $\nu_{e}$ interpretation of the miniboone low-energy excess,”  (2021).
  * Aguilar  _et al._ (2022) A. Aguilar _et al._ , “Miniboone and microboone joint fit to a 3+1 sterile neutrino scenario,”  (2022).
  * Böser  _et al._ (2020) S. Böser, C. Buck, C. Giunti, J. Lesgourgues, L. Ludhova, S. Mertens, A. Schukraft, and M. Wurm, Progress in Particle and Nuclear Physics 111, 103736 (2020).
  * Chu  _et al._ (2018) X. Chu, B. Dasgupta, M. Dentler, J. Kopp, and N. Saviano, JCAP 11, 049 (2018), arXiv:1806.10629 [hep-ph] .
  * Goswami  _et al._ (2022) S. Goswami, V. K. N., A. Mukherjee, and N. Narendra, Phys. Rev. D 105, 095040 (2022), arXiv:2111.14719 [hep-ph] .
  * Aker  _et al._ (2022) M. Aker _et al._ (KATRIN), (2022), arXiv:2207.06337 [nucl-ex] .
  * Antonello  _et al._ (2015) M. Antonello _et al._ (MicroBooNE, LAr1-ND, ICARUS-WA104), (2015), arXiv:1503.01520 [physics.ins-det] .
  * Ajimura  _et al._ (2021) S. Ajimura _et al._ (JSNS2), Nucl. Instrum. Meth. A 1014, 165742 (2021), arXiv:2104.13169 [physics.ins-det] .
  * Agarwalla  _et al._ (2017) S. K. Agarwalla, S. S. Chatterjee, and A. Palazzo, Phys. Rev. Lett. 118, 031804 (2017), arXiv:1605.04299 [hep-ph] .
  * Kumar Agarwalla  _et al._ (2013) S. Kumar Agarwalla, S. Prakash, and S. Uma Sankar, PoS EPS-HEP2013, 534 (2013).
  * Ghosh  _et al._ (2018) M. Ghosh, S. Gupta, Z. M. Matthews, P. Sharma, and A. Williams, PoS NuFact2017, 133 (2018), arXiv:1712.06714 [hep-ph] .
  * Choubey  _et al._ (2019) S. Choubey, D. Dutta, and D. Pramanik, The European Physical Journal C 79 (2019), 10.1140/epjc/s10052-019-7479-8.
  * Dutta  _et al._ (2016) D. Dutta, R. Gandhi, B. Kayser, M. Masud, and S. Prakash, JHEP 11, 122 (2016), arXiv:1607.02152 [hep-ph] .
  * Parke and Zhang (2020a) S. J. Parke and X. Zhang, Phys. Rev. D 101, 056005 (2020a).
  * Li  _et al._ (2018) W. Li, J. Ling, F. Xu, and B. Yue, Journal of High Energy Physics 2018 (2018), 10.1007/jhep10(2018)021.
  * Behera  _et al._ (2017) S. P. Behera, A. Ghosh, S. Choubey, V. M. Datar, D. K. Mishra, and A. K. Mohanty, Eur. Phys. J. C 77, 307 (2017), arXiv:1605.08607 [hep-ph] .
  * Gandhi and Ghoshal (2012) R. Gandhi and P. Ghoshal, Phys. Rev. D 86, 037301 (2012), arXiv:1108.4360 [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] .
  * Chatterjee  _et al._ (2013) A. Chatterjee, P. Ghoshal, S. Goswami, and S. K. Raut, JHEP 06, 010 (2013), arXiv:1302.1370 [hep-ph] .
  * Parke and Zhang (2020b) S. J. Parke and X. Zhang, Phys. Rev. D 101, 056005 (2020b), arXiv:1905.01356 [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] .
  * Fogli  _et al._ (1994) G. L. Fogli, E. Lisi, and D. Montanino, Phys. Rev. D 49, 3626 (1994).
  * Abi  _et al._ (2020b) B. Abi _et al._ (DUNE), JINST 15, T08008 (2020b), arXiv:2002.02967 [physics.ins-det] .
  * Acciarri  _et al._ (2016) R. Acciarri _et al._ (DUNE), (2016), arXiv:1601.05471 [physics.ins-det] .
  * Alion  _et al._ (2016) T. Alion _et al._ (DUNE), (2016), arXiv:1606.09550 [physics.ins-det] .
  * Gaisser and Honda (2002) T. K. Gaisser and M. Honda, Ann. Rev. Nucl. Part. Sci. 52, 153 (2002), arXiv:hep-ph/0203272 .
  * 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).
  * Suzuki  _et al._ (1987) T. Suzuki, D. F. Measday, and J. P. Roalsvig, Phys. Rev. C 35, 2212 (1987).
  * Klinskikh  _et al._ (2008) A. Klinskikh, S. Brianson, V. Brudanin, _et al._ , Bull. Russ. Acad. Sci. Phys. 72 (2008), https://doi.org/10.3103.
  * Ohlsson and Snellman (2000) T. Ohlsson and H. Snellman, J. Math. Phys. 41, 2768 (2000), [Erratum: J.Math.Phys. 42, 2345 (2001)], arXiv:hep-ph/9910546 .
  * Kamo  _et al._ (2003) Y. Kamo, S. Yajima, Y. Higasida, S.-I. Kubota, S. Tokuo, and J.-I. Ichihara, Eur. Phys. J. C 28, 211 (2003), arXiv:hep-ph/0209097 .
  * Akhmedov  _et al._ (2004) E. K. Akhmedov, R. Johansson, M. Lindner, T. Ohlsson, and T. Schwetz, JHEP 04, 078 (2004), arXiv:hep-ph/0402175 .