# Search for sterile neutrinos by shower events at a future neutrino telescope

Yabin Wang and Osamu Yasuda  Department of Physics, Tokyo Metropolitan
University,
Hachioji, Tokyo 192-0397, Japan

###### Abstract

It is pointed out that searching for sterile neutrinos at high energy regions
at a future IceCube-like facility has advantages compared with that of reactor
or short baseline accelerator neutrino experiments, in which the size of the
detector and energy resolution make it difficult to gain good sensitivity for
a large mass-squared difference. In this study we show that it is possible to
improve constraints on sterile neutrino mixing $\theta_{14}$ for 1
eV${}^{2}~{}\lesssim\Delta m_{41}^{2}\lesssim$ 100 eV2 by looking for a dip in
the shower events at an IceCube-like neutrino telescope whose volume is at
least 10 times as large as that of IceCube and duration is 10 years. We also
give an analytic expression for the oscillation probabilities in two cases
where the condition of one mass scale dominance is satisfied.

##  1 Introduction

Since Super-Kamiokande first gave firm evidence of neutrino oscillations in
1998, studying neutrino oscillations due to non-zero neutrino mass has become
a promising direction of exploring new physics beyond the Standard Model of
particle physics. According to LEP ParticleDataGroup:2020ssz , there are only
three light active neutrinos $(\nu_{e},\nu_{\mu},\nu_{\tau})$. Based on
atmospheric and solar neutrino experiments as well as reactor and accelerator
experiments, there are three mass eigenstates $(\nu_{1},\nu_{2},\nu_{3})$ with
corresponding masses $(m_{1},m_{2},m_{3})$. The three-flavor eigenstates of
neutrinos are mixed with the three mass eigenstates via a $3\times 3$ unitary
matrix, which is called the PMNS matrix. In the standard three-flavor
framework of neutrino oscillations, there are six parameters: three mixing
angles, $\theta_{12}$, $\theta_{13}$, and $\theta_{23}$, two independent mass-
squared differences $\Delta m_{21}^{2}$ and $\Delta m_{32}^{2}$, where $\Delta
m_{ij}^{2}=m_{i}^{2}-m_{j}^{2}\ (i>j;i,j=1,2,3)$, and one CP phase, $\delta$.
As for the problem of the neutrino mass hierarchy, if $m_{3}>m_{2}>m_{1}$ we
call it a normal ordering (NO), whereas, if $m_{3}<m_{1}<m_{2}$, it is called
an inverted ordering (IO). Furthermore, based on atmospheric and solar
neutrino experiments, there are two mass-squared differences:
$\Delta{m}^{2}_{\mbox{\rm\scriptsize
atm}}\equiv|\Delta{m}^{2}_{32}|\simeq|\Delta{m}^{2}_{31}|$ and
$\Delta{m}^{2}_{\mbox{\rm\scriptsize sol}}\equiv\Delta{m}^{2}_{21}$.
Therefore, we have a hierarchy in the mass-squared differences:

$\displaystyle|\Delta{m}^{2}_{31}|\simeq|\Delta{m}^{2}_{32}|\gg\Delta{m}^{2}_{21}\,.$
(1)

The updated values of the above mentioned parameters are as follows
ParticleDataGroup:2020ssz :

$\displaystyle\sin^{2}{\theta_{12}}$ $\displaystyle=0.307\pm 0.013\,,$ (2)
$\displaystyle\sin^{2}{\theta_{13}}$ $\displaystyle=(2.18\pm 0.07)\times
10^{-2}\,,$ (3) $\displaystyle\sin^{2}{\theta_{23}}$ $\displaystyle=0.545\pm
0.021\ (\mathrm{NO})\,,$ (4) $\displaystyle\sin^{2}{\theta_{23}}$
$\displaystyle=0.547\pm 0.021\ (\mathrm{IO})\,,$ (5)
$\displaystyle\Delta{m}^{2}_{21}$ $\displaystyle=(7.53\pm 0.18)\times 10^{-5}\
\mathrm{eV}^{2}\,,$ (6) $\displaystyle\Delta{m}^{2}_{32}$
$\displaystyle=(2.453\pm 0.034)\times 10^{-3}\ \mathrm{eV}^{2}\
(\mathrm{NO})\,,$ (7) $\displaystyle\Delta{m}^{2}_{32}$
$\displaystyle=(-2.546^{+0.034}_{-0.040})\times 10^{-3}\ \mathrm{eV}^{2}\
(\mathrm{IO})\,,$ (8) $\displaystyle\delta$ $\displaystyle=1.36\pm 0.17\ \pi\
\mathrm{rad}\,.$ (9)

The first experimental anomaly that cannot be explained by the standard three-
flavor scheme was given by LSND Aguilar:2001ty . LSND was looking for an
oscillation signal in the channel of $\bar{\nu}_{\mu}\to\bar{\nu}_{e}$ and
observed an excess of $\bar{\nu}_{e}$ events above the background. This result
was supported by the MiniBooNE experiment MiniBooNE:2018esg . This anomaly can
be solved if there is at least one extra neutrino, called a sterile neutrino
(see, e.g., Abazajian:2012ys ), with a mass-squared difference of $O(1)\
\mathrm{eV}^{2}$, which is much larger than
$\Delta{m}^{2}_{\mbox{\rm\scriptsize atm}}$ and
$\Delta{m}^{2}_{\mbox{\rm\scriptsize sol}}$.

Moreover, solar and reactor neutrino experiments showed that there are the so-
called “gallium neutrino anomaly” Giunti:2010zu  and “reactor antineutrino
anomaly” Mention:2011rk ; Huber:2011wv , in which the observed $\nu_{e}$ and
reactor $\bar{\nu}_{e}$ event numbers are smaller than the theoretical
prediction, and these anomalies may be solved by introducing the existence of
light sterile neutrinos whose mass is in the order of 1 eV, which also
motivates us to search for sterile neutrinos. 111The result in Ref.
Adey:2019ywk  disfavors the new flux of reactor antineutrinos Mention:2011rk ;
Huber:2011wv , which provided the original motivation for reactor antineutrino
anomaly from the shape analysis of the energy spectrum. However, as far as the
overall normalization of reactor antineutrinos is concerned, there is
uncertainty which is as large as that of the difference between the old and
new flux Hayes:2013wra  and it is not clear whether the reactor antineutrino
anomaly is disfavored from the result of Ref. Adey:2019ywk .

While the combined results on the parameter space of $(\sin^{2}{2\theta_{\mu
e}},\Delta{m}^{2}_{41})$, where $\sin^{2}{2\theta_{\mu e}}=4|U_{\mu
4}|^{2}|U_{e4}|^{2}$, $U_{e4}$ and $U_{\mu 4}$ are the elements of the PMNS
matrix, given by DayaBay+Bugey-3 ($\bar{\nu}_{e}\to\bar{\nu}_{e}$) and MINOS
($\nu_{\mu}\to\nu_{\mu}$) MINOS:2020iqj  conflict with that of LSND, because
of the uncertainty in the theoretical calculation of the reactor antineutrino
flux Hayes:2013wra , the constraint in the $(\sin^{2}{2\theta_{\mu
e}},\Delta{m}^{2}_{41})$-plane in Ref. MINOS:2020iqj  may not be as strong as
they claimed. This encourages us to improve constraints on sterile neutrino
mixing.

Observing high energy neutrinos at IceCube is advantageous to search for
sterile neutrinos. Because of the matter effect, there can be a “dip” in the
neutrino oscillation probability at certain neutrino energy for the
disappearance channel of muon neutrinos. 222See Refs. Nicolaidis:1999be ;
Yasuda:2000xs ; Nunokawa:2003ep  for earlier references. By looking for the
dip in the energy spectrum, we can constrain the upper and lower bounds of the
mass-squared difference, $\Delta{m}^{2}_{41}$. Searching for sterile neutrinos
at high energy regions at IceCube has advantages compared with that of lower
energy experiments, such as reactor or short baseline accelerator neutrino
experiments, in which the size of the detectors and energy resolution make it
difficult to gain good sensitivity for a large mass-squared difference. In
this article, assuming the absence of the dip in the energy spectrum, we study
the sensitivity to sterile neutrino oscillations using shower events at
IceCube rather than track events studied in Ref. IceCube:2020tka .

The rest of this article is organized as follows. In Section 2, we briefly
introduce the (3+1)-scheme and IceCube experiment. In Section 3, we present
some analytical formulae for oscillation probabilities. In Section 4, we
perform a numerical analysis of shower events expected to be observed at a
future IceCube-like facility. Finally, we draw our conclusions in Section 5.
In Appendices A \- E, we give some details of deriving analytical expressions
of oscillation probabilities.

##  2 Preliminaries

###  2.1 (3+1)-scheme

The simplest extension of the standard three-flavor framework assumes that
there is only one sterile neutrino. The mixing between neutrino flavor
eigenstates $(\nu_{e},\nu_{\mu},\nu_{\tau},\nu_{s})$ and mass eigenstates
$(\nu_{1},\nu_{2},\nu_{3},\nu_{4})$ is formulated as

$\displaystyle\left(\begin{array}[]{rl}\nu_{e}\\\ \nu_{\mu}\\\ \nu_{\tau}\\\
\nu_{s}\end{array}\right)=\left(\begin{array}[]{cccc}U_{e1}&U_{e2}&U_{e3}&U_{e4}\\\
U_{\mu 1}&U_{\mu 2}&U_{\mu 3}&U_{\mu 4}\\\ U_{\tau 1}&U_{\tau 2}&U_{\tau
3}&U_{\tau 4}\\\
U_{s1}&U_{s2}&U_{s3}&U_{s4}\end{array}\right)\left(\begin{array}[]{rl}\nu_{1}\\\
\nu_{2}\\\ \nu_{3}\\\ \nu_{4}\end{array}\right)\,,$ (22)

where $\nu_{s}$ represents the flavor eigenstate of sterile neutrinos,
$\nu_{4}$ denotes the fourth mass eigenstate, and $U_{\alpha i}\
(\alpha=e,\mu,\tau,s;\ i=1,2,3,4)$ denotes the element of the PMNS matrix $U$.
As the (2+2)-scheme has been excluded by the atmospheric and solar neutrino
data Maltoni:2004ei , we only consider the (3+1)-scheme in this study. In this
scheme, the mass $m_{4}$ of the fourth mass eigenstate is much larger than
those of other mass eigenstates. Therefore, we have

$\displaystyle\Delta{m}^{2}_{41}\simeq\Delta{m}^{2}_{42}\simeq\Delta{m}^{2}_{43}\gg|\Delta{m}^{2}_{31}|\simeq|\Delta{m}^{2}_{32}|\gg\Delta{m}^{2}_{21}\,.$
(23)

Further, because the fourth mass eigenstate $\nu_{4}$ is weekly mixed with the
active three-flavor eigenstates $(\nu_{e},\nu_{\mu},\nu_{\tau})$, we have the
following constraint condition

$\displaystyle|U_{\alpha 4}|^{2}\ll 1\qquad(\alpha=e,\mu,\tau)\,.$ (24)

In this study, we adopt the following parametrization of U:

$\displaystyle U$ $\displaystyle={\cal
R}_{34}(\theta_{34},\,\delta_{3})\;{\cal R}_{24}(\theta_{24},\,0)\;{\cal
R}_{14}(\theta_{14},\,\delta_{2})\;{\cal R}_{23}(\theta_{23},\,0)\;{\cal
R}_{13}(\theta_{13},\,\delta_{1})\;{\cal R}_{12}(\theta_{12},\,0)$ (25)
$\displaystyle[{\cal R}_{ij}(\theta_{ij},\
\delta_{l})]_{pq}=\left\\{\begin{array}[]{ll}\cos\theta_{ij}&p=q=i,j\\\
1&p=q\not=i,j\\\ \sin\theta_{ij}\ e^{-i\delta_{l}}&p=i;q=j\\\
-\sin\theta_{ij}\ e^{i\delta_{l}}&p=j;q=i\\\ 0&\mbox{\rm
otherwise.}\end{array}\right.$ (31)

Propagation of the neutrino and ani-neutrino flavor eigenstates $\Psi$ and
$\bar{\Psi}$ in matter can be described by

$\displaystyle i\frac{d\Psi}{dt}$ $\displaystyle=\left(U{\cal E}U^{-1}+{\cal
A}\,\right)\,\Psi$ (32) $\displaystyle i\frac{d\bar{\Psi}}{dt}$
$\displaystyle=\left(U^{\ast}{\cal
E}U^{\ast-1}+{\cal\bar{A}}\,\right)\,\bar{\Psi}$ (33) $\displaystyle\Psi$
$\displaystyle\equiv\left(\begin{array}[]{rl}\nu_{e}\\\ \nu_{\mu}\\\
\nu_{\tau}\\\
\nu_{s}\end{array}\right)\,,\quad\bar{\Psi}\equiv\left(\begin{array}[]{rl}\bar{\nu}_{e}\\\
\bar{\nu}_{\mu}\\\ \bar{\nu}_{\tau}\\\ \bar{\nu}_{s}\end{array}\right)$ (42)
$\displaystyle{\cal E}$ $\displaystyle\equiv\mbox{\rm diag}\left(0,\Delta
E_{21},\Delta E_{31},\Delta E_{41}\right)$ (43) $\displaystyle{\cal A}$
$\displaystyle\equiv\mbox{\rm diag}\left(A_{e}-iB,-iB,-iB,-A_{n}\right)$
$\displaystyle\equiv\mbox{\rm diag}\left(A_{e}-iB,-iB,-iB,A_{s}\right)$ (44)
$\displaystyle{\cal\bar{A}}$ $\displaystyle\equiv\mbox{\rm
diag}\left(-A_{e}-iB,-iB,-iB,A_{n}\right)$ $\displaystyle\equiv\mbox{\rm
diag}\left(-A_{e}-iB,-iB,-iB,-A_{s}\right)$ (45) $\displaystyle A_{e}$
$\displaystyle=\sqrt{2}\,G_{F}N_{e}=\left[\frac{\rho}{2.6~{}\mathrm{(g\cdot
cm^{-3})}}\right]\cdot\left(\frac{Y_{e}}{0.5}\right)\cdot 1.0\times 10^{-13}\
\mathrm{eV}$ (46) $\displaystyle A_{s}$ $\displaystyle\equiv-
A_{n}=\frac{1}{\sqrt{2}}\,G_{F}N_{n}{\color[rgb]{0,0,0}=\left[\frac{\rho}{2.6~{}\mathrm{(g\cdot
cm^{-3})}}\right]\cdot\left(\frac{1-Y_{e}}{0.5}\right)\cdot 5.0\times
10^{-14}\ \mathrm{eV}}$ (47) $\displaystyle B$
$\displaystyle=\frac{1}{2}\,\sigma_{tot}\,(N_{p}+N_{n})\,,$ (48)

where $\Delta E_{jk}{\color[rgb]{0,0,0}\equiv
E_{j}-E_{k}\equiv\sqrt{p^{2}+m^{2}_{j}}-\sqrt{p^{2}+m^{2}_{k}}}\simeq\Delta{m}^{2}_{jk}/(2E)$,
$E$ ($p$) represents the energy (momentum) of neutrinos, $A_{e}$ represents
the matter effect induced by the charged current (CC) interaction between
$\nu_{e}$ and electrons, and $A_{n}\equiv-A_{s}$ represents the matter effect
induced by the neutral current (NC) interaction of $\nu_{\alpha}\
(\alpha=e,\mu,\tau)$ with electrons, $u$ and $d$ quarks. Note that a $4\times
4$ matrix $(E_{1}+A_{n})\,{\bf 1}_{4}$ ($(E_{1}-A_{n})\,{\bf 1}_{4}$), which
is proportional to the $4\times 4$ identity matrix ${\bf 1}_{4}$, was
subtracted from the right hand side of Eq. (32) (Eq. (33)), because it affects
only the phase of the probability amplitude $A(\nu_{\alpha}\to\nu_{\beta})$
($A(\bar{\nu}_{\alpha}\to\bar{\nu}_{\beta})$) and does not affect the
oscillation probability $P(\nu_{\alpha}\to\nu_{\beta})$
($P(\bar{\nu}_{\alpha}\to\bar{\nu}_{\beta})$). $G_{F}=1.17\times 10^{-5}\
\mathrm{GeV}^{-2}$ is the Fermi constant, $N_{e}$ is the number of electrons
in the Earth per unit volume, which is approximately equal to the number
$N_{n}$ of neutrons in matter per unit volume, $\rho$ is the density of
matter, and $Y_{e}=N_{p}/(N_{p}+N_{n})$ is the relative number density of
electrons in matter, where $N_{p}$, which is equal to $N_{e}$ because an atom
is electrically neutral, represents the number of protons per unit volume. $B$
represents the absorption effect Naumov:2001ci  and is related by the total
neutrino–nucleon cross section $\sigma_{tot}$ by Eq. (48), which leads to the
expected result
$P(\nu_{\alpha}\to\nu_{\alpha})=\exp[-\sigma_{tot}\,(N_{p}+N_{n})\,L]$ in the
absence of neutrino mixing. The numerical value of the total neutrino–nucleon
cross section is $\sigma_{tot}\sim 10^{-35}$cm2 at $E=1$ TeV, and the
contribution of the CC and NC interactions to
$\sigma_{tot}=\sigma_{CC}+\sigma_{NC}$ is shown in Figure 1.

In vacuum, if we adopt the condition of “one mass scale dominance”, which
represents the case where only the terms including $\Delta{m}^{2}_{41}$
dominate in the oscillation probability, we have the following formula of
oscillation probability:

$\displaystyle P(\nu_{\alpha}\to\nu_{\beta})\simeq
P(\bar{\nu}_{\alpha}\to\bar{\nu}_{\beta})\simeq\left|\delta_{\alpha\beta}-\sin^{2}{2\theta_{\alpha\beta}}\,\sin^{2}\left(\frac{\Delta{m}^{2}_{41}L}{4E}\right)\right|\,,$
(49)

where

$\displaystyle\sin^{2}{2\theta_{\alpha\beta}}=4|U_{\alpha
4}|^{2}\left|\delta_{\alpha\beta}-|U_{\beta
4}|^{2}\right|\qquad(\alpha,\beta=e,\mu,\tau,s)\,,$ (50)

and $L$ represents the baseline length.

In the discussions on neutrino oscillations at high energy ($E_{\nu}\gtrsim 1$
TeV) with the baseline length of $O(10^{4})$ km,
$|\Delta{m}^{2}_{\mbox{\rm\scriptsize atm}}L/E|$ and
$|\Delta{m}^{2}_{\mbox{\rm\scriptsize sol}}L/E|$ are negligible, and the
condition of the one mass scale dominance is satisfied. As we will see below,
the effect of CP violation is negligible in this case. In the analytical
treatment in Section 3, we discuss oscillation probability in the limit of one
mass scale dominance.

###  2.2 IceCube experiment

IceCube is a neutrino observatory located at the South Pole. IceCube observes
atmospheric and astrophysical neutrinos, whose energy ranges approximately
from 100 GeV to above 109 GeV. Among the neutrino events, there are track and
shower events. Track events are recorded by the Cherenkov light of muons
produced by the CC interaction between muon neutrinos and nuclei in ice.
Shower events include electromagnetic and hadronic shower events, which are
produced by the CC interactions of $\nu_{e}$ and $\nu_{\tau}$ as well as the
NC interactions of any flavor of neutrinos. Track events have good angular
resolution but poor energy resolution, whereas shower events have poor angular
resolution but good energy resolution.

When up-going neutrinos arrive at IceCube after passing through the Earth, the
effective mixing angle in matter can be enhanced due to the matter effect and
a dip can appear in the neutrino oscillation probability for the disappearance
channel. By searching for the dip in the oscillation probability, we can give
the lower and upper bounds on the third mass-squared difference
$\Delta{m}^{2}_{41}$. This was done for the disappearance channel of muon
neutrinos IceCube:2020tka .

At high energy ($\gtrsim 10$ GeV) muon neutrinos constitute most atmospheric
neutrinos because most high energy atmospheric muons do not decay into
electrons and neutrinos before they reach the Earth. Therefore, at high energy
the flux of atmospheric $\nu_{e}+\bar{\nu}_{e}$, which mainly come from kaon
decays, is much smaller than that of $\nu_{\mu}+\bar{\nu}_{\mu}$. Despite this
fact, motivated by the reactor antineutrino and gallium neutrino anomalies, we
study the sensitivity of shower events to neutrino oscillations, anticipating
the large detector volume may improve this disadvantage in the future.

In the energy range of 1-100 TeV, prompt neutrinos, which come from the decay
of charmed mesons, may have fluxes of $\nu_{\mu}$ and $\nu_{e}$, comparable to
that of the conventional atmospheric neutrinos. Although prompt neutrinos can
be essential components of atmospheric neutrinos, the result of Ref.
IceCube:2015mgt  agrees with a hypothesis of zero prompt neutrino event.
Therefore, in the analysis in Section 4, we do not consider the contribution
of prompt neutrinos in the search for sterile neutrinos.

##  3 Some analytical treatment

In the numerical analysis in Section 4, it is useful to have analytic
expressions for oscillation probabilities to understand the qualitative
behavior of the energy spectrum etc. So in this section we discuss some
analytic expressions for oscillation probabilities in matter with constant
density. For simplicity we discuss a case without absorption, i.e., $B=0$ in
Eqs. (32) and (33), since the analytic expressions for oscillation
probabilities with $B\neq 0$ are complicated.

###  3.1 Atmospheric neutrinos

In this subsection, we discuss the oscillation probabilities of atmospheric
neutrinos. Upward-going atmospheric neutrinos pass through the Earth and are
observed at the IceCube detector. The baseline length $L$ is given by

$\displaystyle L=-R\cos\Theta+\sqrt{(R+h)^{2}-(R\sin\Theta)^{2}}\,,$ (51)

where $R$ is the Earth’s radius, $\Theta$ is the zenith angle, and $h\sim$ 20
km is the thickness of the atmosphere.

In this study, we discuss neutrinos at high energy ($E_{\nu}\gtrsim 1$ TeV)
with a baseline length $L\lesssim 13000)$ km. Since
$|\Delta{m}^{2}_{\mbox{\rm\scriptsize atm}}L/E|\ll 1$ and
$|\Delta{m}^{2}_{\mbox{\rm\scriptsize sol}}L/E|\ll 1$, the condition of the
one mass scale dominance is satisfied. In this case, the energy difference
matrix in Eq. (43) can be approximated as

$\displaystyle{\cal E}\simeq\mbox{\rm diag}\left(0,0,0,\Delta E\right)\,,$
(52)

where we have defined

$\displaystyle\Delta E\equiv\Delta E_{41}$ (53)

for simplicity, and the matrices ${\cal R}_{23}(\theta_{23},\,0)$, ${\cal
R}_{13}(\theta_{13},\,\delta_{1})$, ${\cal R}_{12}(\theta_{12},\,0)$ disappear
from $U$ in Eqs. (32) and (33).

It turns out that the two cases $\theta_{14}\neq 0$,
$\theta_{24}=\theta_{34}=0$ and $\theta_{14}=0$, $\theta_{24}\neq 0$,
$\theta_{34}\neq 0$ in the limit of one mass scale dominance can be reduced to
a two-flavor problem. The former (latter) has a resonance in the neutrino
(antineutrino) mode, so we discuss these two cases separately. In both cases,
CP phases do not appear in the oscillation probability, as in the two-flavor
case.

####  3.1.1 The case with $\theta_{14}\neq 0$, and $\theta_{24}=\theta_{34}=0$

We first consider a case where there is only one non-zero mixing angle
$\theta_{14}$, whereas $\theta_{24}$ and $\theta_{34}$ vanish. In this case,
there is mixing only between electron and sterile neutrinos, and Eq. (32)
reduces to the one in the two-flavor case. From Eq. (129) in Appendix A, we
can calculate the oscillation probability of $\nu_{e}\to\nu_{e}$:

$\displaystyle
P(\nu_{e}\to\nu_{e})=1-\sin^{2}2\tilde{\theta}_{14}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)\,,$
(54)

where $\tilde{\theta}_{14}$ and $\Delta\tilde{E}_{1}$ are, respectively,
defined by

$\displaystyle\tan{2\tilde{\theta}_{14}}$ $\displaystyle=\frac{\Delta
E\,\sin{2\theta_{14}}}{\Delta E\cos{2\theta_{14}}-A_{e}+A_{s}}\,.$ (55)
$\displaystyle\Delta\tilde{E}_{1}$ $\displaystyle=\sqrt{(\Delta
E\cos{2\theta_{14}}-A_{e}+A_{s})^{2}+(\Delta E\sin{2\theta_{14}})^{2}}\,.$
(56)

Since $A_{e}-A_{s}\simeq A_{e}/2>0$ and $\Delta E=\Delta{m}^{2}_{41}/(2E)>0$,
Eq. (55) indicates that a resonance occurs only in the neutrino mode at the
neutrino energy
$E={\color[rgb]{0,0,0}\Delta{m}^{2}_{41}\cos{2\theta_{14}}/\\{2(A_{e}-A_{s})\\}\simeq}$
$\Delta{m}^{2}_{41}/\\{2(A_{e}-A_{s})\\}$, where we have assumed
$|\theta_{14}|\ll 1$.

Moreover, it is useful to study the behavior of oscillation probabilities when
considering the following:

$\displaystyle\Delta{m}^{2}_{41}\to+\infty\,.$ (57)

If $\Delta{m}^{2}_{41}\to+\infty$, then the vacuum oscillation contribution
$\Delta E$ dominates over the matter term $A_{e}-A_{s}$ in Eq. (55), so we
have $\tilde{\theta}_{14}\to\theta_{14}$, whereas $\Delta\tilde{E}\to+\infty$
leads to rapid oscillation in Eq. (54). Thus we obtain

$\displaystyle P(\nu_{e}\to\nu_{e})\to 1-\frac{1}{2}\,\sin^{2}2\theta_{14}\,.$
(58)

From Eq. (58), we observe that if $\Delta{m}^{2}_{41}\to+\infty$, the
oscillation probability $P(\nu_{e}\to\nu_{e})$ becomes independent of the
zenith angle $\Theta$ and neutrino energy $E$. The fact that
$P(\nu_{e}\to\nu_{e})$ becomes a constant suppression as
$\Delta{m}^{2}_{41}\to+\infty$ implies that the significance of the
(3+1)-scheme with $\theta_{14}\neq 0$ at very large $\Delta{m}^{2}_{41}$
depends entirely on the systematic errors.

####  3.1.2 The case with $\theta_{14}=0$, $\theta_{24}\neq 0$ and
$\theta_{34}\neq 0$

Next, we consider a case where there are two non-zero mixing angles
$\theta_{24}$ and $\theta_{34}$, whereas $\theta_{14}$ vanishes. Although the
(3+1)-scheme for which we perform the numerical analysis is not this case, it
is instructive to see how a resonance occurs in the antineutrino mode in this
case.

As explained in Appendix B, there can be a resonance in the antineutrino mode
(see Eq. (206).). Therefore, we discuss the disappearance channel
$\bar{\nu}_{\mu}\to\bar{\nu}_{\mu}$, introducing a new set of angles

$\displaystyle\tan{\varphi_{12}}\equiv\displaystyle\frac{t_{34}}{s_{24}}\,,$
(59) $\displaystyle\sin{\varphi_{13}}\equiv c_{24}c_{34}\,,$ (60)
$\displaystyle\tan{\varphi_{23}}\equiv\displaystyle\frac{s_{34}}{t_{24}}\,,$
(61)

where $c_{jk}\equiv\cos\theta_{jk}$, $s_{jk}\equiv\sin\theta_{jk}$,
$t_{jk}\equiv\tan\theta_{jk}$, and $\theta_{jk}~{}(j,k=1,\cdots,4)$ is the
original mixing angle in the parametrization (25). The oscillation
probabilities $P(\bar{\nu}_{\mu}\to\bar{\nu}_{\mu})$ and
$P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})$ can be expressed as:

$\displaystyle P(\bar{\nu}_{\mu}\to\bar{\nu}_{\mu})={}$ $\displaystyle
1-C_{23}^{4}\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)$
$\displaystyle~{}~{}-\sin^{2}2\varphi_{23}\,\left\\{\tilde{S}_{13}^{2}\sin^{2}\left(\frac{\varepsilon_{2-}L}{2}\right)+\tilde{C}_{13}^{2}\sin^{2}\left(\frac{\varepsilon_{2+}L}{2}\right)\right\\}\,,$
(62) $\displaystyle P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})={}$
$\displaystyle\sin^{2}2\varphi_{23}\left\\{\tilde{S}_{13}^{2}\sin^{2}\left(\frac{\varepsilon_{2-}L}{2}\right)+\tilde{C}_{13}^{2}\sin^{2}\left(\frac{\varepsilon_{2+}L}{2}\right)\right.$
$\displaystyle\hskip
56.9055pt-\left.\frac{1}{4}\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\right\\}\,,$
(63)

where $\tilde{\varphi}_{13}$, $\Delta\tilde{E}_{2}$ and $\varepsilon_{2\mp}$
are, respectively, given by

$\displaystyle\tan{2\tilde{\varphi}_{13}}$ $\displaystyle=\frac{\Delta
E\sin{2\varphi_{13}}}{\Delta E\cos{2\varphi_{13}}+A_{s}}\,,$ (64)
$\displaystyle\Delta\tilde{E}_{2}$ $\displaystyle=\sqrt{(\Delta
E\cos{2\varphi_{13}}+A_{s})^{2}+(\Delta E\sin{2\varphi_{13}})^{2}}\,,$ (65)
$\displaystyle\varepsilon_{2\mp}$ $\displaystyle=\frac{\Delta
E-A_{s}}{2}\mp\frac{1}{2}\Delta\tilde{E}_{2}\,,$ (66)

and we introduced the notations $\tilde{S}_{13}=\sin{\tilde{\varphi}_{13}}$,
$\tilde{C}_{13}=\cos{\tilde{\varphi}_{13}}$, $S_{23}=\sin{\varphi_{23}}$ and
$C_{23}=\cos{\varphi_{23}}$. If $\theta_{24}$ and $\theta_{34}$ are both small
so that $0<s_{24}\ll 1$ and $0<s_{34}\ll 1$ are satisfied, from Eq. (60), we
have $0<\pi/2-\varphi_{13}\ll 1$, which implies that $2\varphi_{13}\simeq\pi$
and $\cos 2\varphi_{13}\simeq-1$, so that the resonance condition, in which
the denominator of Eq. (64) vanishes, is possible only for antineutrinos.
Notice that the denominator for neutrinos would be $\Delta E\cos
2\varphi_{13}-A_{s}$, which is negative for any value of the neutrino energy
$E>0$.

Now, let us consider the behavior of the oscillations probabilities
$\bar{\nu}_{\mu}\to\bar{\nu}_{\mu}$ and $\bar{\nu}_{\mu}\to\bar{\nu}_{\tau}$
as $\Delta{m}^{2}_{41}\to+\infty$, where $\Delta E$ dominates over the matter
term $A_{s}$ in Eq. (64), so we have $\tilde{\varphi}_{13}\to\varphi_{13}$.
Moreover, for the three energy differences, two of them become large whereas
one remains finite:

$\displaystyle\Delta\tilde{E}_{2}$ $\displaystyle\sim\Delta E\to+\infty\,.$
(67) $\displaystyle\varepsilon_{2+}$ $\displaystyle=\frac{\Delta
E-A_{s}}{2}+\frac{1}{2}\Delta\tilde{E}_{2}\sim\Delta E\to+\infty\,.$ (68)
$\displaystyle\varepsilon_{2-}$
$\displaystyle=\frac{\varepsilon_{2+}\varepsilon_{2-}}{\varepsilon_{2+}}$
$\displaystyle=\frac{A_{s}\Delta E\cos^{2}\varphi_{13}}{\varepsilon_{2+}}\to
A_{s}\cos^{2}{\varphi_{13}}\,,$ (69)

From this, $\sin^{2}(\Delta\tilde{E}_{2}L/2)$ and
$\sin^{2}(\varepsilon_{2+}L/2)$ are averaged to yield $1/2$, whereas
$\sin^{2}(\varepsilon_{2+}L/2)$ becomes
$\sin^{2}(A_{s}\cos^{2}{\varphi_{13}}L/2)$, which does not depend on the
neutrino energy $E$ but depends on the baseline length $L$. Thus, we obtain

$\displaystyle P(\bar{\nu}_{\mu}\to\bar{\nu}_{\mu})$ $\displaystyle\to
1-\frac{C_{23}^{4}}{2}\,\sin^{2}2\varphi_{13}$
$\displaystyle~{}~{}-\sin^{2}2\varphi_{23}\,\left\\{{S}_{13}^{2}\sin^{2}\left(\frac{A_{s}\cos^{2}{\varphi_{13}}L}{2}\right)+\frac{{C}_{13}^{2}}{2}\right\\}\,,$
(70) $\displaystyle P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})$
$\displaystyle\to\sin^{2}2\varphi_{23}\left\\{{S}_{13}^{2}\sin^{2}\left(\frac{A_{s}\cos^{2}{\varphi_{13}}L}{2}\right)+\frac{{C}_{13}^{2}}{2}-\frac{1}{8}\sin^{2}2\varphi_{13}\right\\}\,,$
(71)

Eq. (70) shows that, even as $\Delta{m}^{2}_{41}\to+\infty$, the disappearance
probability has some dependence on the zenith angle $\Theta$ (not only
$P(\bar{\nu}_{\mu}\to\bar{\nu}_{\mu})$ but also $P(\nu_{\mu}\to\nu_{\mu})$
because Eq. (70) is an even function of $A_{s}$.), which is why the IceCube
has some sensitivity to $\theta_{24}$ even at $\Delta m^{2}_{41}=100$ eV2 (cf.
the upper panel of FIG.19 in Ref. IceCube:2020tka ). This phenomenon has also
been known in the standard three-flavor scenario Fogli:1997ff ; Foot:1998pv ,
where the atmospheric neutrino data was accounted for to some extent, even for
a low value of the mass-squared difference
$|\Delta{m}^{2}_{\mbox{\rm\scriptsize atm}}|$.

###  3.2 Astrophysical neutrinos

Next, we discuss neutrino oscillations of astrophysical neutrinos. We consider
the following propagation process: astrophysical neutrinos start from the
source and travel a long distance in outer space (in vacuum), then reach and
pass through the Earth (in matter), and are finally observed at IceCube. In
the following, we develop formulae for calculating the oscillation
probabilities in this process.

Calculating these oscillation probabilities is similar to that of the day-
night effect of solar neutrinos. We denote $L_{0}$ as the astronomical travel
distance of astrophysical neutrinos in outer space, and $L$ as the propagation
distance of astrophysical neutrinos inside the Earth. Because $|\Delta
m^{2}_{ij}L_{0}/E|\gg 1~{}(i,j=1,\cdots,4)$, to a good approximation, we can
assume $L_{0}\to+\infty$. Further, we denote $T$ as the transition matrix of
the neutrino flavor state inside the Earth. Then we have the formulae of
oscillation probabilities in the entire propagation process:

$\displaystyle P(\nu_{\alpha}\to\nu_{\beta})$
$\displaystyle=\lim_{L_{0}\to+\infty}\sum_{\rho
j}T_{\beta\rho}\,U_{\rho\,j}\,\exp(-iL_{0}E_{j})\,U^{-1}_{j\alpha}$
$\displaystyle\qquad\qquad\times\left(\sum_{\sigma
n}T_{\beta\sigma}\,U_{\sigma
n}\,\exp(-iL_{0}E_{n})\,U^{-1}_{n\alpha}\right)^{\ast}$
$\displaystyle=\lim_{L_{0}\to+\infty}\sum_{\rho\,\sigma
jn}T_{\beta\rho}\,T_{\beta\sigma}^{\ast}U_{\rho j}\,U_{\sigma
n}^{\ast}\exp(iL_{0}\Delta{E_{nj}})U_{\alpha j}^{\ast}U_{\alpha n}$
$\displaystyle=\sum_{\rho\,\sigma
n}T_{\beta\rho}\,T_{\beta\sigma}^{\ast}U_{\rho n}\,U_{\sigma
n}^{\ast}\,|U_{\alpha n}|^{2}\qquad(\alpha,\beta=e,\mu,\tau,s)\,,$ (72)

where we have used the unitary condition of the PMNS matrix $U$.

####  3.2.1 The case with $\theta_{14}\neq 0$, and $\theta_{24}=\theta_{34}=0$

The flavor ratio of neutrinos in the case of $\theta_{14}\neq 0$,
$\theta_{24}=\theta_{34}=0$ is given by the following (for a detailed
description, see Appendix C):

$\displaystyle F(\nu_{e}):F(\nu_{\mu}):F(\nu_{\tau}):F(\nu_{s})$
$\displaystyle\simeq 1-|T_{14}|^{2}:1:1:|T_{14}|^{2}\,,$ (73)

where $|T_{14}|^{2}$ is defined as

$\displaystyle|T_{14}|^{2}=\sin^{2}2\tilde{\theta}_{14}\sin^{2}\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)\,,$
(74)

$\tilde{\theta}_{14}$ and $\Delta\tilde{E}_{1}$ are defined by Eqs. (55) and
(56), respectively, and $L$ is the path length of astrophysical neutrinos
inside the Earth. Note that the factor $|T_{14}|^{2}$ and hence the flavor
ratio (73) have a nontrivial dependence on the neutrino energy $E$ especially
near the resonance region $E=\Delta m^{2}_{41}\cos
2\theta_{14}/\\{2(A_{e}-A_{s})\\}\sim\Delta m^{2}_{41}/\\{2(A_{e}-A_{s})\\}$.

####  3.2.2 The case with $\theta_{14}=0$, $\theta_{24}\neq 0$ and
$\theta_{34}\neq 0$

In the case with $\theta_{14}=0$, $\theta_{24}\neq 0$, and $\theta_{34}\neq
0$, as described in detail in Appendix D, the flavor ratio of antineutrinos is
expressed as follows:

$\displaystyle
F(\bar{\nu}_{e}):F(\bar{\nu}_{\mu}):F(\bar{\nu}_{\tau}):F(\bar{\nu}_{s})$
$\displaystyle\simeq
1:1-C_{23}^{2}\,|T_{44}|^{2}:1-S_{23}^{2}\,|T_{44}|^{2}:1-|T_{44}|^{2}\,,$
(75)

where $|T_{44}|^{2}$ is defined as

$\displaystyle|T_{44}|^{2}$
$\displaystyle=1-\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\,,$
(76)

$C_{23}\equiv\cos\varphi_{23}$, $S_{23}\equiv\sin\varphi_{23}$, the new mixing
angles $\varphi_{ij}~{}(i,j=1,2,3)$ are defined in Eqs. (59), (60), (61),
$\tilde{\varphi}_{13}$ and $\Delta\tilde{E}_{2}$ are defined by Eqs. (64) and
(65), respectively, and $L$ is the path length of astrophysical neutrinos
inside the Earth. Notice that the factor $|T_{44}|^{2}$ and therefore the
flavor ratio (75) also have a nontrivial dependence on the neutrino energy $E$
especially near the resonance region $E{\color[rgb]{0,0,0}=-\Delta
m^{2}_{41}\cos 2\varphi_{13}/(2A_{s})}\sim\Delta m^{2}_{41}/(2A_{s})$, where
we have assumed $|\pi/2-\varphi_{13}|\ll 1$.

##  4 Numerical analysis of shower events

In this section, for simplicity, we perform numerical analysis to obtain the
sensitivity to $\theta_{14}$ as a function of $\Delta{m}^{2}_{41}$ in the
(3+1)-scheme, assuming $\theta_{24}=\theta_{34}=0$. As shown in Section 3, the
case with $\theta_{14}\neq 0,\theta_{24}=\theta_{34}=0$ has a resonance in the
neutrino mode, whereas the case with $\theta_{14}=0,\theta_{24}\neq
0,\theta_{34}\neq 0$ has one in the antineutrino mode. If we turn on small
$\theta_{24}$ and $\theta_{34}$, the resonance in the neutrino mode has small
perturbation with respect to $\theta_{24}$ and $\theta_{34}$, and we expect
that our conclusion does not change significantly.

We assume that the numbers of events at an IceCube-like neutrino telescope are
given by extrapolating those in FIG.6 in Ref. IceCube:2015mgt , i.e., they are
given by a factor (volume of the IceCube-like detector) $\times$ (number of
days of measurement) / [(volume of the IceCube) $\times$ (332 days)]. The
figure shows the energy spectrum of various numbers of shower events for 332
days at IceCube. They comprise conventional atmospheric
$\nu_{e}+\bar{\nu}_{e}$, conventional atmospheric $\nu_{\mu}+\bar{\nu}_{\mu}$,
(which look like shower events even after the experimental cuts),
astrophysical $\nu_{\beta}+\bar{\nu}_{\beta}~{}(\beta=e,\mu,\tau)$, and cosmic
ray muons. Note that the IceCube collaboration concluded that zero prompt
neutrino event is consistent with their data, so we also assume that there is
no prompt neutrino event in our analysis.

###  4.1 Numbers of events

The neutrinos observed at an IceCube-like detector are either atmospheric or
astrophysical, and are expected to be produced mainly from pion and kaon
decays. So, it is expected that there are only $\nu_{e}+\bar{\nu}_{e}$ and
$\nu_{\mu}+\bar{\nu}_{\mu}$ at the production point. However, during neutrino
propagation, the processes $\nu_{\alpha}\to\nu_{\tau}~{}(\alpha=e,\mu)$ and
$\bar{\nu}_{\alpha}\to\bar{\nu}_{\tau}~{}(\alpha=e,\mu)$ exist in general
neutrino mixing. For atmospheric neutrinos, in the standard three-flavor
mixing scenario or (3+1)-scheme where the only non-zero sterile neutrino
mixing angle is $\theta_{14}$, there are very few
$\nu_{\tau}+\bar{\nu}_{\tau}$ from the conventional atmospheric neutrinos
through the oscillations $\nu_{\mu}\to\nu_{\tau}$ and
$\bar{\nu}_{\mu}\to\bar{\nu}_{\tau}$ because the oscillation probability is
very small: $P(\nu_{\mu}\to\nu_{\tau})\simeq
P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})\simeq\sin^{2}2\theta_{23}\sin^{2}(\Delta
m^{2}_{31}L/4E)\ll 1$ for $E>1$ TeV, $L\lesssim 2R$, where $R$ is the Earth’s
radius. Meanwhile, for astrophysical neutrinos, we have
$\nu_{\tau}+\bar{\nu}_{\tau}$ because of oscillations with an astronomical
baseline length Learned:1994wg .

Thus, the number of events ($N_{jk}(\nu_{\beta})$ for $\nu_{\beta}$ and
$N_{jk}(\bar{\nu}_{\beta})$ for $\bar{\nu}_{\beta}$) for the energy bin
($E_{j}<E<E_{j+1}$) and zenith angle bin
($\cos\Theta_{k}<\cos\Theta<\cos\Theta_{k+1}$) can be expressed in terms of
the flux $\Phi(\nu_{\alpha},E^{\prime},\Theta^{\prime})$, the cross section
Gandhi:1998ri  ($\sigma_{X}(E^{\prime})$ for neutrinos and
$\bar{\sigma}_{X}(E^{\prime})$ for antineutrinos, for CC ($X=CC$) and NC
($X=NC$) interactions, respectively) and the detection efficiency
$\epsilon_{X}(E,\nu_{\beta})~{}(X=CC,NC)$ as follows:

$\displaystyle{\ }\left\\{\begin{array}[]{c}N_{jk}(\nu_{\beta})\\\
N_{jk}(\bar{\nu}_{\beta})\end{array}\right\\}$ (79) $\displaystyle=$
$\displaystyle\int_{\log_{10}E_{j}}^{\log_{10}E_{j+1}}d\log_{10}E\,\int_{\cos\Theta_{j}}^{\cos\Theta_{j+1}}d\cos\Theta\,\int_{-\infty}^{\infty}d\log_{10}E^{\prime}\,\int_{-1}^{1}d\cos\Theta^{\prime}\,$
$\displaystyle\times\,R_{e}(E,E^{\prime},\Delta_{e}(E))\,R_{z}(\Theta,\Theta^{\prime},\Delta_{z}(E))$
$\displaystyle\times\sum_{\alpha}\sum_{X=CC,NC}\,\left\\{\begin{array}[]{c}\Phi(\nu_{\alpha},E^{\prime},\Theta^{\prime})\,P(\nu_{\alpha}\to\nu_{\beta})\,\sigma_{X}(E^{\prime})\,\epsilon_{X}(E^{\prime},\nu_{\beta})\\\
\Phi(\bar{\nu}_{\alpha},E^{\prime},\Theta^{\prime})\,P(\bar{\nu}_{\alpha}\to\bar{\nu}_{\beta})\,\bar{\sigma}_{X}(E^{\prime})\,\epsilon_{X}(E^{\prime},\bar{\nu}_{\beta})\end{array}\right\\}$
(82)

$R_{e}$ and $R_{z}$ represent the resolution functions in the neutrino energy
and zenith angle, respectively:

$\displaystyle R_{e}(E,E^{\prime},\Delta_{e}(E))$
$\displaystyle\equiv\frac{1}{\sqrt{\pi}\Delta_{e}(E)}\,\exp\left\\{-\frac{(\log_{10}E-\log_{10}E^{\prime})^{2}}{\Delta_{e}(E)^{2}}\right\\}$
$\displaystyle R_{z}(\Theta,\Theta^{\prime},\Delta_{z}(E))$
$\displaystyle\equiv
N(\Theta)\,\exp\left\\{-\frac{(\cos\Theta-\cos\Theta^{\prime})^{2}}{\Delta_{z}(E)^{2}}\right\\}$
$\displaystyle N(\Theta)$
$\displaystyle\equiv\left[\int_{-1}^{1}d\cos\Theta^{\prime}\,\exp\left\\{-\frac{(\cos\Theta-\cos\Theta^{\prime})^{2}}{\Delta_{z}(E)^{2}}\right\\}\right]^{-1}$
(83)

where $E^{\prime}$ and $\Theta^{\prime}$ are the true neutrino energy and
zenith angle, respectively, whereas $E$ and $\Theta$ are the reconstructed
ones, respectively. We use the energy resolution $\Delta_{e}(E)$ and the
zenith angle resolution $\Delta_{z}(E)$ which are given by the blue curve
“iterative CREDO” in FIG.5 in Ref. IceCube:2015mgt .

####  4.1.1 Atmospheric neutrinos

We use the atmospheric neutrino flux in Ref. Honda:2006qj  in our analysis.
However, the flux in Ref. Honda:2006qj  is given only up to 10 TeV, and we
need the flux up to ${\cal O}$(100 TeV) for our analysis. Therefore, we adopt
the following empirical formula in Ref. Chirkin:2004ic  for the atmospheric
neutrino flux, and extrapolate it up to ${\cal O}$(100 TeV):

$\displaystyle\left\\{\begin{array}[]{c}\Phi^{\mbox{\scriptsize
atm}}(\nu_{\mu},E,\Theta)\\\ \Phi^{\mbox{\scriptsize
atm}}(\bar{\nu}_{\mu},E,\Theta)\end{array}\right\\}=$
$\displaystyle\frac{2.85\cdot 10^{-2}}{\mbox{ cm}^{2}\mbox{ sr}\mbox{ s}\mbox{
GeV}}\cdot\mbox{A}\cdot\left(\frac{E}{\mbox{GeV}}\right)^{-\gamma}$ (86)
$\displaystyle\times\left(\frac{1}{1+6\,E\cos\Theta^{\ast}/\mathrm{121\
GeV}}+\frac{0.213}{1+1.44\,E\cos\Theta^{\ast}/\mathrm{897\ GeV}}\right)$ (87)
$\displaystyle\left\\{\begin{array}[]{c}\Phi^{\mbox{\scriptsize
atm}}(\nu_{e},E,\Theta)\\\ \Phi^{\mbox{\scriptsize
atm}}(\bar{\nu}_{e},E,\Theta)\end{array}\right\\}=$
$\displaystyle\frac{2.4\cdot 10^{-3}}{\mbox{cm}^{2}\mbox{ sr}\mbox{ s}\mbox{
GeV}}\cdot\mbox{A}\cdot\left(\frac{E}{\mbox{GeV}}\right)^{-\gamma}$ (90)
$\displaystyle\times\left(\frac{0.05}{1+1.5\,E\cos\Theta^{\ast}/\mathrm{897\
GeV}}+\frac{0.185}{1+1.5\,E\cos\Theta^{\ast}/\mathrm{194\ GeV}}\right.$
$\displaystyle\left.+\frac{11.4\,E^{\,\zeta(\Theta)}}{1+1.21\,E\cos\Theta^{\ast}/\mathrm{121\
GeV}}\right)$ (91)

where the cosine of the zenith angle $\Theta$ is modified as that of the
effective one $\Theta^{\ast}$ by considering the curvature of the atmosphere:

$\displaystyle\cos\Theta^{\ast}$
$\displaystyle=\sqrt{\frac{\cos^{2}\Theta+p_{1}^{2}+p_{2}\cos^{\,p_{3}}\Theta+p_{4}\cos^{\,p_{5}}\Theta}{1+p_{1}^{2}+p_{2}+p_{4}}}$
$\displaystyle(p_{1},p_{2},p_{3},p_{4},p_{5})$
$\displaystyle=(0.102573,-0.068287,0.958633,0.0407253,0.817285).$

$A$ is the normalization, given by $A=0.646$ for $\nu_{\mu}$,
$\bar{\nu}_{\mu}$ and $A=0.828$ for $\nu_{e}$, $\bar{\nu}_{e}$. $\gamma$ is
the spectral index for atmospheric neutrinos, and the flux for each flavor
($\nu_{\mu}$, $\bar{\nu}_{\mu}$, $\nu_{e}$ and $\bar{\nu}_{e}$) in Ref.
Honda:2006qj  is best reproduced by $\gamma=2.680$ (for $\nu_{\mu}$), 2.696
(for $\bar{\nu}_{\mu}$), 2.719 (for $\nu_{e}$), and 2.741 (for
$\bar{\nu}_{e}$). Hence, we use Eqs. (87) and (91) with these values for
$\gamma$ up to the neutrino energy $E\lesssim$ 1 PeV. $\zeta$ in Eq. (91) is
given by

$\displaystyle\zeta$ $\displaystyle=a+b\log_{10}(E/\mathrm{GeV})$
$\displaystyle a$
$\displaystyle=\max(0.11-2.4\cos\Theta,-0.46-0.54\cos\Theta)$ $\displaystyle
b$ $\displaystyle=\min(-0.22+0.69\cos\Theta,-0.01+0.01\cos\Theta)$

####  4.1.2 Detection efficiency for shower events

The detection efficiency
$\epsilon_{X}(E_{\color[rgb]{0,0,0}\nu},\nu_{\beta})~{}(\beta=e,\mu;X=CC,NC)$
in Eq. (82) for shower events can be deduced from the number of shower events
of atmospheric neutrinos given in Ref. IceCube:2015mgt , if we make the
following assumptions:

  1. (i)

The detection efficiencies for neutrinos and antineutrinos are the same.

$\displaystyle\epsilon_{X}(E_{\color[rgb]{0,0,0}\nu},\nu_{\beta})=\epsilon_{X}(E_{\color[rgb]{0,0,0}\nu},\bar{\nu}_{\beta})\quad(\beta=e,\mu,\tau;X=CC,NC)$
(92)

  2. (ii)

The energy dependencies of $\epsilon_{NC}(E_{\nu}^{NC},\nu_{e})$ and
$\epsilon_{CC}(E_{\nu}^{CC},\nu_{e})$ are related by the following formula:

$\displaystyle\color[rgb]{0,0,0}\epsilon_{NC}(E_{\nu}^{NC},\nu_{e})=\mathrm{const.}\,\epsilon_{CC}(E_{\nu}^{CC},\nu_{e})$
(93)

where the neutrino energy $E_{\nu}^{CC}$ ($E_{\nu}^{NC}$) in the charged
(neutral) current interaction are related through the observed energy
$E_{obs}$ by the following: 333 Since the contribution of the neutral current
interaction to the numbers of events is not dominant, the precise value
$E_{\nu}-E_{obs}$ in the relation between the neutrino energy $E_{\nu}^{CC}$
($E_{\nu}^{NC}$) and the observed energy $E_{obs}$ does not affect the
sensitivity to $\theta_{14}$ very much.

$\displaystyle\color[rgb]{0,0,0}E_{obs}\simeq E_{\nu}^{CC}-1~{}\mbox{\rm
GeV}\simeq E_{\nu}^{NC}-10~{}\mbox{\rm GeV}\,.$ (94)

  3. (iii)

The detection efficiencies $\epsilon_{NC}(E,\nu_{\beta})$ for all flavors
$\beta=e,\mu$, and $\tau$ are the same.

$\displaystyle\epsilon_{NC}(E_{\nu},\nu_{e})=\epsilon_{NC}(E_{\nu},\nu_{\mu})=\epsilon_{NC}(E_{\nu},\nu_{\tau})$
(95)

  4. (iv)

The detection efficiency for hadronic shower events from $\nu_{\tau}$ and
$\bar{\nu}_{\tau}$ is the same as that for electromagnetic shower events from
$\nu_{e}$ and $\bar{\nu}_{e}$. Since 65% (18%) of tau decay is hadronic
(electromagnetic) ParticleDataGroup:2020ssz , this assumption implies

$\displaystyle\epsilon_{CC}(E_{\nu},\nu_{\tau})=0.83\,\epsilon_{CC}(E_{\nu},\nu_{e})\,.$
(96)

From TABLE I in Ref. IceCube:2015mgt , the numbers of events for 332 days at
IceCube are 198 (CC) + 17 (NC) for $\nu_{e}+\bar{\nu}_{e}$, 387 (CC) + 258
(NC) for $\nu_{\mu}+\bar{\nu}_{\mu}$, 115 for cosmic ray muons, and 103 for
astrophysical $\nu_{\beta}+\bar{\nu}_{\beta}~{}(\beta=e+\mu+\tau)$. Since we
know the flux of conventional atmospheric neutrinos and the cross section at
each energy, the CC and NC numbers of events for $\nu_{e}+\bar{\nu}_{e}$ along
with assumptions (i) and (ii) allow us to deduce the detection efficiency for
$\epsilon_{CC}(E,\nu_{e})$ and $\epsilon_{NC}(E,\nu_{e})$, where they are
extrapolated up to the neutrino energy $E\lesssim$ 10 PeV. Then, from the CC
and NC numbers of events for $\nu_{\mu}+\bar{\nu}_{\mu}$ with assumptions
(iii), we can deduce $\epsilon_{CC}(E,\nu_{\mu})$. Thus, we obtain
$\epsilon_{CC}(E,\nu_{e})$, $\epsilon_{NC}(E,\nu_{e})$ and
$\epsilon_{CC}(E,\nu_{\mu})$, and they are depicted in an arbitrary unit in
Figure 1. Assumption (iv) can be used to deduce $\epsilon_{CC}(E,\nu_{\tau})$
if there are $\nu_{\tau}+\bar{\nu}_{\tau}$.

Figure 1:  $\epsilon_{CC}(E,\nu_{e})$, $\epsilon_{NC}(E,\nu_{e})$ and
$\epsilon_{CC}(E,\nu_{\mu})$, which are deduced from FIG.6 in Ref.
IceCube:2015mgt , the flux in Refs. Honda:2006qj  and Chirkin:2004ic , the
cross sections in Ref. Gandhi:1998ri  with the CTEQ4–DIS parton distributions,
and assumptions (i) - (iii) in this article. The absolute magnitude of the
detection efficiency is arbitrary, and we fix the normalization so that the
total numbers of events coincide with those in Ref. IceCube:2015mgt . The
cross sections for the charged (neutral) current neutrino-nucleon
(antineutrino-nucleon) interactions in Ref. Gandhi:1998ri  are also shown in
the figure.

####  4.1.3 Astrophysical neutrinos

The flux of astrophysical neutrinos is assumed to be isotropical, and the
ratios $\Phi^{\mbox{\scriptsize ast}}(\nu_{\mu})/\Phi^{\mbox{\scriptsize
ast}}(\nu_{e})$ and $\Phi^{\mbox{\scriptsize
ast}}(\bar{\nu}_{\mu})/\Phi^{\mbox{\scriptsize ast}}(\bar{\nu}_{e})$ at the
production point are assumed to be 2:

$\displaystyle{\ }\Phi^{\mbox{\scriptsize
ast}}(\nu_{\mu},E,\Theta)=\Phi^{\mbox{\scriptsize
ast}}(\bar{\nu}_{\mu},E,\Theta)=2\,\Phi^{\mbox{\scriptsize
ast}}(\nu_{e},E,\Theta)=2\,\Phi^{\mbox{\scriptsize
ast}}(\bar{\nu}_{e},E,\Theta)$ $\displaystyle=$
$\displaystyle\,\frac{1}{3}\,\Phi^{\mbox{\scriptsize
ast}}_{0}\left(\frac{E}{100\mathrm{TeV}}\right)^{-\gamma}$ (97)

where $\gamma$ is the astrophysical spectral index, and
$\Phi^{\mbox{\scriptsize ast}}_{0}$ is the total flux of astrophysical
neutrinos at $E=$100 TeV.444 If there is no enhancement of the sterile mixing
angle $\theta_{14}$ due to the matter effect in the Earth, then each flux at
$E=$100 TeV is $\Phi^{\mbox{\scriptsize
ast}}(\nu_{e})\simeq\Phi^{\mbox{\scriptsize
ast}}(\bar{\nu}_{e})\simeq\Phi^{\mbox{\scriptsize
ast}}(\nu_{\mu})\simeq\Phi^{\mbox{\scriptsize
ast}}(\bar{\nu}_{\mu})\simeq\Phi^{\mbox{\scriptsize
ast}}(\nu_{\tau})\simeq\Phi^{\mbox{\scriptsize
ast}}(\bar{\nu}_{\tau})\simeq\Phi^{\mbox{\scriptsize ast}}_{0}/6$, so that the
total flux is approximately $\Phi^{\mbox{\scriptsize ast}}_{0}$.  We adopt the
reference value $\gamma=2.4$ used in Ref. IceCube:2015mgt .

If there is no enhancement of the sterile mixing angle $\theta_{14}$ due to
the matter effect in the Earth, the oscillation probabilities are
approximately given by those in vacuum in the limit of an infinite baseline
length, and

$\displaystyle\sum_{\alpha=e,\mu}\,\Phi^{\mbox{\scriptsize
ast}}(\nu_{\alpha},E^{\prime},\Theta^{\prime})\,P(\nu_{\alpha}\to\nu_{e})$
$\displaystyle\simeq$
$\displaystyle\sum_{\alpha=e,\mu}\,\Phi^{\mbox{\scriptsize
ast}}(\nu_{\alpha},E^{\prime},\Theta^{\prime})\,P(\nu_{\alpha}\to\nu_{\mu})$
$\displaystyle\simeq$
$\displaystyle\sum_{\alpha=e,\mu}\,\Phi^{\mbox{\scriptsize
ast}}(\nu_{\alpha},E^{\prime},\Theta^{\prime})\,P(\nu_{\alpha}\to\nu_{\tau})\,,$
(98)

as was pointed out in Ref. Learned:1994wg . Meanwhile, if there is an
enhancement of the sterile mixing angle $\theta_{14}$, Eq. (98) is modified
for the energy range in which the enhancement of the effective mixing angle
$\theta_{14}$ occurs. Particularly, if the sterile mixing angle $\theta_{14}$
and mass-squared difference lie in a suitable range, then the flavor ratio
$(\nu_{e}+\bar{\nu}_{e})$ : $(\nu_{\mu}+\bar{\nu}_{\mu})$ :
$(\nu_{\tau}+\bar{\nu}_{\tau})$ $=1-|T_{14}|^{2}/2:1:1$, (where the flavor
ratio of neutrinos is given by Eq. (73), and the flavor ratio of antineutrinos
is 1:1:1 in the absence of resonance), may be observed, instead of the
standard flavor ratio Learned:1994wg  $(\nu_{e}+\bar{\nu}_{e})$ :
$(\nu_{\mu}+\bar{\nu}_{\mu})$ : $(\nu_{\tau}+\bar{\nu}_{\tau})$ =1:1:1. This
deviation depends on the neutrino energy, and it can be a signal to search for
sterile neutrino oscillations. Although we discuss only the case with
$\theta_{14}\neq 0$, $\theta_{24}=\theta_{34}=0$ in this study, if the
astrophysical neutrino events can be separated from those of atmospheric
neutrinos, the energy dependence of the flavor ratio is in general one of the
signatures of sterile neutrino oscillations for the range
$1~{}\mathrm{eV}^{2}\lesssim\Delta m^{2}_{41}\lesssim 100~{}\mathrm{eV}^{2}$.

###  4.2 $\chi^{2}$

In this study, we examine whether the numbers of shower events predicted by
the (3+1)-scheme with $\theta_{14}\neq 0$ can be distinguished from those
without sterile neutrino oscillation, i.e., those with the standard three-
flavor framework with the central value of the oscillation parameters given in
Eqs. (2)-(8) with NO. To evaluate the significance between the number of
events with $\theta_{14}\neq 0$ and the one without non-zero $\theta_{14}$, we
use the following Poisson likelihood chi-square:

$\displaystyle\chi^{2}_{0}=2\sum\limits_{j,k}\left[Y_{jk}-N_{jk}-N_{jk}\log\left(\frac{Y_{jk}}{N_{jk}}\right)\right]\,,$
(99)

where $j$ represents the $j$-th zenith angle bin
($(j-11)/10\leq\cos\Theta\leq(j-10)/10$, $1\leq j\leq 20$), and $k$ represents
the $k$-th energy bin
($({\color[rgb]{0,0,0}k}+7)/4\leq\log_{10}(E/\mathrm{GeV})\leq({\color[rgb]{0,0,0}k}+8)/4$,
$1\leq k\leq 16$). Moreover,

$\displaystyle Y_{jk}=$
$\displaystyle{}(1+\alpha_{1})\left[N^{\mbox{\scriptsize
atm}}_{jk}(\nu_{e}+\bar{\nu}_{e};\theta_{14})+N^{\mbox{\scriptsize
atm}}_{jk}(\nu_{\mu}+\bar{\nu}_{\mu};0)\right]$
$\displaystyle+(1+\alpha_{2})N^{\mbox{\scriptsize
cr}-\mu}_{jk}+(1+\alpha_{3})N^{\mbox{\scriptsize ast}}_{jk}(\theta_{14})$
(100)

and

$\displaystyle N_{jk}=$ $\displaystyle N^{\mbox{\scriptsize
atm}}_{jk}(\nu_{e}+\bar{\nu}_{e};0)+N^{\mbox{\scriptsize
atm}}_{jk}(\nu_{\mu}+\bar{\nu}_{\mu};0)$ $\displaystyle+N^{\mbox{\scriptsize
cr}-\mu}_{jk}+N^{\mbox{\scriptsize ast}}_{jk}(0)\,,$ (101)

where $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ are three pull parameters;
$N^{\mbox{\scriptsize atm}}_{jk}(\nu_{e}+\bar{\nu}_{e};\theta_{14})$ and
$N^{\mbox{\scriptsize atm}}_{jk}(\nu_{e}+\bar{\nu}_{e};0)$ represent the
number of events of atmospheric electron neutrinos with and without sterile
neutrino oscillation (i.e., $\theta_{14}\neq 0$ or $\theta_{14}=0$),
respectively, while assuming all three-flavor oscillation parameters as the
central value given in Eqs. (2)-(8) with NO in both cases;
$N^{\mbox{\scriptsize atm}}_{jk}(\nu_{\mu}+\bar{\nu}_{\mu};0)$ represents the
number of events of atmospheric muon neutrinos without sterile neutrino
oscillation, $N^{\mbox{\scriptsize ast}}_{jk}(\theta_{14})$ and
$N^{\mbox{\scriptsize ast}}_{jk}(0)$ represent the number of events of
astrophysical neutrinos with and without sterile neutrino oscillation,
respectively, $N^{\mbox{\scriptsize cr}-\mu}_{jk}$ stands for the number of
events of cosmic ray muons. Notice that in all numbers of events
$N^{\ast}_{jk}(\ast,\theta_{14})$ and $N^{\ast}_{jk}(\ast,0)$, the effects of
the standard three-flavor oscillation are considered, since we are interested
in the difference between the numbers of shower events with the standard
three-flavor oscillation scenario and those with the (3+1)-scheme.

The total $\chi^{2}$ is defined as follows:

$\displaystyle\chi^{2}=\min_{\alpha_{1},\alpha_{2},\alpha_{3}}\left[\chi_{0}^{2}+\left(\frac{\alpha_{1}}{\sigma_{1}}\right)^{2}+\left(\frac{\alpha_{2}}{\sigma_{2}}\right)^{2}+\left(\frac{\alpha_{3}}{\sigma_{3}}\right)^{2}\right]\,,$
(102)

where $\sigma_{1}$, $\sigma_{2}$ and $\sigma_{3}$ are the three systematic
errors for atmospheric neutrinos, cosmic ray muons, and astrophysical
neutrinos, respectively. In this study, we set
$\sigma_{1}=\sigma_{2}=\sigma_{3}=0.4$ for simplicity. 555 Although various
systematic errors are given in Ref. IceCube:2020tka , in our analysis, for
simplicity, we consider only the overall normalization for the atmospheric
neutrinos ($\sigma_{1}$), astrophysical neutrinos ($\sigma_{3}$), and cosmic
ray muons ($\sigma_{2}$). The reference values quoted in Ref. IceCube:2020tka
are $\sigma_{1}=0.4$ and $\sigma_{3}=0.36$, but we could not find the one for
$\sigma_{2}$. For simplicity, we assumed that all of them are 0.4.
Implications of this choice will be discussed in Section 5.  In (102) we fix
the three-flavor oscillation parameters and assume the central value given in
Eqs. (2)-(8) with Normal Ordering. The reason that we do not marginalize
$\chi^{2}$ with respect to the three-flavor oscillation parameters is that we
are interested in fitting the energy spectrum with a dip (i.e., the one with
non-zero $\theta_{14}$) to that without a dip, and variations with respect to
the three-flavor oscillation parameters are not expected to affect this
fitting. Thus, $\chi^{2}$ becomes a function of $\theta_{14}$ and $\Delta
m^{2}_{41}$.

###  4.3 Sensitivity

We use the number of events mentioned in (4.1) and evaluate $\chi^{2}$ in Eq.
(102). The region with $\chi^{2}>4.61$ represents the region in which the
hypothesis with $\theta_{14}\neq 0$ can be excluded at 90%CL, and it is
depicted in the ($\theta_{14}$, $\Delta{m}^{2}_{41}$)-plane in Figure 2. From
Figure 2, if we use the number of events observed for more than 100 years with
the present IceCube detector, the sensitivity can be better than the one given
by DayaBay+Bugey at 90% CL MINOS:2020iqj  for
$1~{}\mathrm{eV}^{2}\lesssim\Delta{m}^{2}_{41}\lesssim 100~{}\mathrm{eV}^{2}$.
To achieve this sensitivity in the near future, we must increase the volume of
the detector to get more numbers of events. We also observe from Figure 2 that
the sensitivity to $\Delta{m}^{2}_{41}$ is saturated for
$\Delta{m}^{2}_{41}\gtrsim$ 100 eV2. Figure 3 shows the dependence of the
oscillation probability $P(\nu_{e}\to\nu_{e})$, which is calculated
numerically by considering the absorption effect $B$ in Eq. (32), as a
function of the neutrino energy $E$. To see the dip, realized by the neutrino
energy $E_{\mbox{\rm\scriptsize
dip}}=\Delta{m}^{2}_{41}/\\{2(A_{e}-A_{s})\\}$, $E_{\mbox{\rm\scriptsize
dip}}\lesssim$ 500 TeV must be satisfied, and this condition yields
$\Delta{m}^{2}_{41}\lesssim$ 100 eV2.

Figure 2:  The expected excluded region from the shower events of the IceCube
experiment for data size 10, 20, 100, 500, and 1000 years. The bound at 90% CL
from the combined analysis of the Bugey3 and Daya Bay experiments
MINOS:2020iqj  is also shown. Figure 3:  Oscillation probability
$P(\nu_{e}\to\nu_{e})$ as a function of the neutrino energy $E$ for various
values of $\Delta m^{2}_{41}$ in the case of $\sin^{2}2\theta_{14}=0.05$ and
the baseline length $L$ = 12756 km (= the diameter of the Earth). For
$E\gtrsim$ 100 TeV, the probability decreases because of the absorption effect
of neutrino.

Figure 4 shows the contribution of atmospheric and astrophysical neutrinos to
the sensitivity expected by observation for 500 years at the present IceCube.
We observe that atmospheric neutrinos contribute more to the sensitivity than
astrophysical neutrinos because the number of events of the former is larger.
Moreover, the peak of the sensitivity is obtained at $\Delta{m}^{2}_{41}\sim$
0.1 eV2 ($\Delta{m}^{2}_{41}\sim$ 1 eV2) for atmospheric (astrophysical)
neutrinos, which is because the number of events of atmospheric neutrinos
(astrophysical neutrinos) has a peak at $E\sim$ 1 TeV ($E\sim$ 10 TeV) in the
energy spectrum (cf. FIG.6 in Ref. IceCube:2015mgt ), and the condition for
$\Delta{m}^{2}_{41}$ to have the dip (the denominator of Eq. (123) vanishes)
leads to $\Delta{m}^{2}_{41}\sim$ 0.1 eV2 ($\Delta{m}^{2}_{41}\sim$ 1 eV2).

Figure 4:  Hypothetical excluded region from the shower events of the IceCube
experiment for 500 years, assuming only atmospheric neutrinos (red curve) or
only astrophysical neutrinos (green curve), where the thin blue curve stands
for the region with both atmospheric and astrophysical neutrinos for 500
years. The bound from the combined analysis of the Bugey3 and the Daya Bay
experiments MINOS:2020iqj  is also shown (the thin black curve).

We note in passing that the curve from astrophysical neutrinos in Figure 4 has
two intersections with a straight line $\Delta{m}^{2}_{41}$=constant for
$0.06~{}\mbox{\rm eV}^{2}\lesssim\Delta{m}^{2}_{41}\lesssim O(1)$ eV2, i.e.,
there are two solutions of $\chi^{2}=4.61$ for a fixed value of
$\Delta{m}^{2}_{41}$ in this region of $\Delta{m}^{2}_{41}$, which can be
qualitatively understood as follows. As described in Appendix E, as
$\theta_{14}\to\pi/4$, if $\Delta
E\equiv\Delta{m}^{2}_{41}/(2E)\ll|A_{e}-A_{s}|$, the $\nu_{e}$ flux becomes
approximately independent of the neutrino energy $E$ and baseline length $L$
in the Earth. In this case, if the systematic error is not very small, the
suppression due to neutrino oscillation is constant and the data may be fit by
considering the overall normalization uncertainty. Meanwhile, if
$\Delta{m}^{2}_{41}$ becomes large, the condition
$\Delta{m}^{2}_{41}/(2E)\ll|A_{e}-A_{s}|$ is no longer satisfied, and the
above argument does not hold. This argument qualitatively explains the
behavior of the curve from astrophysical neutrinos in Figure 4.

##  5 Conclusions

In this study we showed that a future neutrino telescope facility, whose
volume is at least 10 times as large as that of IceCube, can improve
sensitivity to $\theta_{14}$ by running for 10 years for 1
eV${}^{2}\lesssim\Delta m^{2}_{41}\lesssim$100 eV2 in the (3+1)-scheme, if no
dip is observed in the energy spectrum of the shower events.

We have assumed that $\theta_{24}=\theta_{34}=0$ in our numerical analysis. If
we turn on small $\theta_{24}$ and $\theta_{34}$, the resonance in the
neutrino mode has small perturbation with respect to $\theta_{24}$ and
$\theta_{34}$, and we expect that our conclusion do not change very much.

Although the variety of systematic errors in this study is smaller than that
in Ref. IceCube:2020tka , it is not expected to significantly affect the
sensitivity to $\sin^{2}2\theta_{14}$ at the peak, since the sensitivity is
mainly given by the energy spectrum’s shape. What is significantly affected by
the details of the systematic errors is the behavior for $\Delta
m^{2}_{41}\gtrsim$ 100 eV2, because the oscillation probability in the region
is expected to be averaged out to constant, and the details of the overall
normalization significantly differ. So, the sensitivity to
$\sin^{2}2\theta_{14}$ for $\Delta m^{2}_{41}\gtrsim$100 eV2 is expected to
depend on the details of the systematic errors, and more careful analysis
would be required.

In our analysis, we assumed that there is no prompt neutrino event. If the
volume of the future detector is increased, there can be a small contribution
of prompt neutrino events for the neutrino energy $E\gtrsim$ 104.5 GeV.
However, from the result of Ref. IceCube:2015mgt  in which the data agree with
a hypothesis of zero prompt neutrino event, the contribution from prompt
neutrinos is expected to be much smaller than that of astrophysical neutrinos.
Therefore, we expect that our results are not significantly altered even if
the contribution from prompt neutrinos is considered.

We studied only the shower events for simplicity. If we combine the track and
shower events, we can even give stronger constraints on the sterile mixing
angles. In the presence of all non-zero sterile mixing angles, it is expected
that we can also see a dip in the antineutrino disappearance channel
$\bar{\nu}_{\mu}\to\bar{\nu}_{\mu}$ at the same energy region as that in
neutrino disappearance channel $\nu_{e}\to\nu_{e}$. This dip can be observed
in the track events, and it gives us a constraint on
$\varphi_{13}=\sin^{-1}(c_{24}c_{34})$, so we could obtain more information on
the sterile mixing angles than in this study.

In this article we obtain the sensitivity to $\theta_{14}$ on the assumption
that no dip is observed. If the dip is discovered, however, the situation
becomes different from that of the present work. By looking at the shower
events only, it is impossible to conclude whether the dip comes from the
neutrino mode which is caused mainly by the enhancement due to $\theta_{14}$,
from the antineutrino mode which is caused mainly by the enhancement due to
$\theta_{24}$ and $\theta_{34}$, or from both in the presence of all
$\theta_{j4}~{}(j=1,2,3)$. To solve such degeneracy, we need information from
the track events. Discussions on such degeneracy is complicated, and it will
be a future work.

In summary, we emphasize that measurements of high energy neutrinos at a
neutrino telescope have a potential advantage over low energy short baseline
experiments in investigating sterile neutrino oscillations for 1
eV${}^{2}~{}\lesssim\Delta m_{41}^{2}\lesssim$ 100 eV2, if the volume of the
detector is at least 10 times as large as that of IceCube.

## Acknowledgements

This research was partly supported by a Grant-in-Aid for Scientific Research
of the Ministry of Education, Science and Culture, under Grants No. 18K03653,
No. 18H05543 and No. 21K03578.

Appendix

##  Appendix A Derivation of Eq. (54) in the case of $\theta_{14}\neq 0$, and
$\theta_{24}=\theta_{34}=0$

In this Appendix, we discuss oscillations of atmospheric neutrinos for
$\theta_{14}\neq 0$, and $\theta_{24}=\theta_{34}=0$. Since $\nu_{\mu}$ and
$\nu_{\tau}$ are decoupled from $\nu_{e}$ and $\nu_{s}$, the system can be
effectively described by the two-flavor framework. The only remaining CP phase
$\delta_{2}$ in Eq. (25) disappears by absorbing in the flavor eigenstate
$\nu_{\alpha}~{}(\alpha=e,s)$. Thus, we express the evolution equation as
follows:

$\displaystyle i\frac{d}{dt}\left(\begin{array}[]{rl}\nu_{e}\\\
\nu_{s}\end{array}\right)=\left[U_{1}\left(\begin{array}[]{cc}0&0\\\ 0&\Delta
E\\\ \end{array}\right)U_{1}^{-1}+\left(\begin{array}[]{cc}A_{e}&0\\\
0&A_{s}\\\ \end{array}\right)\right]\left(\begin{array}[]{rl}\nu_{e}\\\
\nu_{s}\end{array}\right)\,,$ (111)

where $\Delta
E\equiv\sqrt{p^{2}+m^{2}_{2}}-\sqrt{p^{2}+m^{2}_{1}}\simeq(m_{2}^{2}-m_{1}^{2})/(2p)\equiv\Delta
m^{2}/(2p)\simeq\Delta m^{2}/(2E)$,

$\displaystyle U_{1}=\exp(i\theta_{14}\sigma_{2})\,,$ (112)

and $\sigma_{2}$ is one of the three Pauli matrices

$\displaystyle\sigma_{1}=\left(\begin{array}[]{cc}0&1\\\ 1&0\\\
\end{array}\right)\,,~{}~{}\sigma_{2}=\left(\begin{array}[]{cc}0&-i\\\ i&0\\\
\end{array}\right)\,,~{}~{}\sigma_{3}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\\\
\end{array}\right)\,.$ (119)

We denote the Hamiltonian of Eq. (111) as $H_{1}$ and rewrite it as follows:

$\displaystyle H_{1}$ $\displaystyle=\exp(i\theta_{14}\sigma_{2})\frac{\Delta
E}{2}({\bf
1}_{2}-\sigma_{3})\exp(-i\theta_{14}\sigma_{2})+\frac{A_{e}}{2}({\bf
1}_{2}+\sigma_{3})+\frac{A_{s}}{2}({\bf 1}_{2}-\sigma_{3})$
$\displaystyle=\frac{\Delta E+A_{e}+A_{s}}{2}{\bf
1}_{2}-\frac{1}{2}(\alpha\sigma_{3}-\beta\sigma_{1})\,,$ (120)

where

$\displaystyle\alpha=\Delta E\cos{2\theta_{14}}-A_{e}+A_{s},\ \beta=\Delta
E\sin{2\theta_{14}}\,,$ (121)

and ${\bf 1}_{2}$ is the $2\times 2$ identity matrix. Denoting the mixing
angle in matter as $\tilde{\theta}_{14}$, we have

$\displaystyle\exp(-i\tilde{\theta}_{14}\sigma_{2})(\alpha\sigma_{3}-\beta\sigma_{1})\exp(i\tilde{\theta}_{14}\sigma_{2})$
$\displaystyle=$
$\displaystyle(\alpha\cos{2\tilde{\theta}_{14}}+\beta\sin{2\tilde{\theta}_{14}})\sigma_{3}+(\alpha\sin{2\tilde{\theta}_{14}}-\beta\cos{2\tilde{\theta}_{14}})\sigma_{1}\,.$
(122)

If we define

$\displaystyle\tan{2\tilde{\theta}_{14}}=\frac{\beta}{\alpha}=\frac{\Delta
E\sin{2\theta_{14}}}{\Delta E\cos{2\theta_{14}}-A_{e}+A_{s}}\,,$ (123)

which implies

$\displaystyle\sin{2\tilde{\theta}_{14}}=\frac{\beta}{\sqrt{\alpha^{2}+\beta^{2}}}\,,\
\cos{2\tilde{\theta}_{14}}=\frac{\alpha}{\sqrt{\alpha^{2}+\beta^{2}}}\,,$
(124)

we obtain

$\displaystyle\alpha\sigma_{3}-\beta\sigma_{1}=\sqrt{\alpha^{2}+\beta^{2}}\exp(i\tilde{\theta}_{14}\sigma_{2})\sigma_{3}\exp(-i\tilde{\theta}_{14}\sigma_{2})\,.$
(125)

Further, if we define

$\displaystyle\Delta\tilde{E}_{1}=\sqrt{\alpha^{2}+\beta^{2}}=\sqrt{(\Delta
E\cos{2\theta_{14}}-A_{e}+A_{s})^{2}+(\Delta E\sin{2\theta_{14}})^{2}}\,,$
(126)

Eq. (120) becomes

$\displaystyle H_{1}$
$\displaystyle=\exp(i\tilde{\theta}_{14}\sigma_{2})\left(\frac{\Delta
E+A_{e}+A_{s}}{2}{\bf
1}_{2}-\frac{1}{2}\Delta\tilde{E}_{1}\sigma_{3}\right)\exp(-i\tilde{\theta}_{14}\sigma_{2})$
$\displaystyle=\exp(i\tilde{\theta}_{14}\sigma_{2})\left(\begin{array}[]{cc}\varepsilon_{1-}&0\\\
0&\varepsilon_{1+}\end{array}\right)\exp(-i\tilde{\theta}_{14}\sigma_{2})\,,$
(129)

where

$\displaystyle\varepsilon_{1\mp}=\frac{\Delta
E+A_{e}+A_{s}}{2}\mp\frac{1}{2}\Delta\tilde{E}_{1}\,.$ (130)

From these discussions, the oscillation probability in this case is obtained
by replacing $\theta_{14}$ by $\tilde{\theta}_{14}$ and $\Delta E$ by
$\Delta\tilde{E}_{1}$ in the vacuum oscillation probability, and we obtain Eq.
(54).

##  Appendix B Derivation of Eqs. (62) and (63) in the case of
$\theta_{14}=0$, $\theta_{24}\neq 0$, and $\theta_{34}\neq 0$

In this Appendix, we discuss oscillations of atmospheric neutrinos for
$\theta_{14}=0$, $\theta_{24}\neq 0$, and $\theta_{34}\neq 0$. In this case,
it turns out that resonance occurs in the antineutrino mode. So, we start with
Eq. (33) for antineutrinos. There is mixing among $\bar{\nu}_{\mu}$,
$\bar{\nu}_{\tau}$ and $\bar{\nu}_{s}$, and Eq. (33) reduces to that in the
three-flavor case

$\displaystyle i\frac{d}{dt}\left(\begin{array}[]{rl}\bar{\nu}_{\mu}\\\
\bar{\nu}_{\tau}\\\
\bar{\nu}_{s}\end{array}\right)=\left[U_{2}^{\ast}\left(\begin{array}[]{ccc}0&0&0\\\
0&0&0\\\ 0&0&\Delta E\\\
\end{array}\right)U_{2}^{\ast-1}+\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&0\\\
0&0&-A_{s}\\\
\end{array}\right)\right]\left(\begin{array}[]{rl}\bar{\nu}_{\mu}\\\
\bar{\nu}_{\tau}\\\ \bar{\nu}_{s}\end{array}\right)\,.$ (143)

In the discussions of the three-flavor mixing case, it is useful to introduce
the following three Gell-Mann matrices:

$\displaystyle\lambda_{2}=\left(\begin{array}[]{ccc}0&-i&0\\\ i&0&0\\\
0&0&0\end{array}\right)\,,~{}~{}\lambda_{5}=\left(\begin{array}[]{ccc}0&0&-i\\\
0&0&0\\\
i&0&0\end{array}\right)\,,~{}~{}\lambda_{7}=\left(\begin{array}[]{ccc}0&0&0\\\
0&0&-i\\\ 0&i&0\end{array}\right)\,.$ (153)

Now using the property of $3\times 3$ matrices

$\displaystyle\mbox{\rm
diag}\left(1,e^{-i\delta_{3}/2},e^{i\delta_{3}/2}\right)=\mbox{\rm
diag}\left(1,e^{i\delta_{3}/2},e^{i\delta_{3}/2}\right)\cdot\mbox{\rm
diag}\left(1,e^{-i\delta_{3}},1\right),$
$\displaystyle\exp(i\theta_{34}\lambda_{7})\cdot\mbox{\rm
diag}\left(1,e^{i\delta_{3}/2},e^{i\delta_{3}/2}\right)=\mbox{\rm
diag}\left(1,e^{i\delta_{3}/2},e^{i\delta_{3}/2}\right)\cdot\exp(i\theta_{34}\lambda_{7})\,,$
$\displaystyle\mbox{\rm
diag}\left(1,e^{-i\delta_{3}},1\right)\cdot\exp(i\theta_{24}\lambda_{5})=\exp(i\theta_{24}\lambda_{5})\cdot\mbox{\rm
diag}\left(1,e^{-i\delta_{3}},1\right)\,,$

in the limit of one mass scale dominance, we can show that the diagonal matrix
$\mbox{\rm diag}(1,e^{i\delta_{3}/2},e^{i\delta_{3}/2})$ ($\mbox{\rm
diag}(1,e^{-i\delta_{3}},1)$) with a CP phase $\delta_{3}$ can be moved to the
left of $\exp(i\theta_{34}\lambda_{7})$ (the right of
$\exp(i\theta_{24}\lambda_{5})$), respectively:

$\displaystyle U_{2}^{\ast}$ $\displaystyle=\mbox{\rm
diag}\left(1,e^{i\delta_{3}/2},e^{-i\delta_{3}/2}\right)\cdot\exp\left(i\theta_{34}\lambda_{7}\right)\cdot\mbox{\rm
diag}\left(1,e^{-i\delta_{3}/2},e^{i\delta_{3}/2}\right)\cdot\exp\left(i\theta_{24}\lambda_{5}\right)$
$\displaystyle=\mbox{\rm
diag}\left(1,e^{i\delta_{3}/2},e^{-i\delta_{3}/2}\right)\cdot\mbox{\rm
diag}\left(1,e^{i\delta_{3}/2},e^{i\delta_{3}/2}\right)\cdot\exp\left(i\theta_{34}\lambda_{7}\right)$
$\displaystyle\quad\times\exp\left(i\theta_{24}\lambda_{5}\right)\cdot\mbox{\rm
diag}\left(1,e^{-i\delta_{3}},1\right)$ $\displaystyle=\mbox{\rm
diag}\left(1,e^{i\delta_{3}},1\right)\cdot\exp\left(i\theta_{34}\lambda_{7}\right)\cdot\exp\left(i\theta_{24}\lambda_{5}\right)\cdot\mbox{\rm
diag}\left(1,e^{-i\delta_{3}},1\right)\,.$ (154)

The phase factor $\mbox{\rm diag}(1,e^{-i\delta_{3}},1)$ cancels in
$U_{2}^{\ast}\,\mbox{\rm diag}(0,0,\Delta E)\,U_{2}^{\ast-1}$, and by
redefining $(\bar{\nu}_{\mu},\bar{\nu}_{\tau},\bar{\nu}_{s})^{T}$ by
$\mbox{\rm
diag}(1,e^{-i\delta_{3}},1)(\bar{\nu}_{\mu},\bar{\nu}_{\tau},\bar{\nu}_{s})^{T}$,
the CP phase $\delta_{3}$ disappears in the one mass scale dominance limit.
Therefore, in the following discussions, we set $\delta_{3}=0$ for simplicity.
In this case, there is no difference between $U_{2}$ and $U_{2}^{\ast}$, so we
denote $U_{2}^{\ast}$ as $U_{2}$ for simplicity in the following discussions.
To transform Eq. (143) into a more familiar form, we introduce the following
matrix

$\displaystyle R=\left(\begin{array}[]{ccc}0&0&1\\\ 0&1&0\\\
1&0&0\end{array}\right)\quad(R^{2}={\bf 1}_{3})\,,$ (158)

where ${\bf 1}_{3}$ is the $3\times 3$ identity matrix. Then the Hamiltonian
$H_{2}$ on the right hand side of Eq. (143) can be transformed into

$\displaystyle H_{2}=R\left[RU_{2}\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&0\\\
0&0&\Delta
E\end{array}\right)U_{2}^{-1}R+\left(\begin{array}[]{ccc}-A_{s}&0&0\\\
0&0&0\\\ 0&0&0\end{array}\right)\right]R\,.$ (165)

Next we define

$\displaystyle H_{2}^{{}^{\prime}}=RU_{2}\left(\begin{array}[]{ccc}0&0&0\\\
0&0&0\\\ 0&0&\Delta
E\end{array}\right)U_{2}^{-1}R+\left(\begin{array}[]{ccc}-A_{s}&0&0\\\
0&0&0\\\ 0&0&0\end{array}\right)\,,$ (172)

and find that $H_{2}^{{}^{\prime}}$ is similar to the Hamiltonian of the
evolution equation of flavor eigenstates of neutrinos in the standard three-
flavor framework including the matter effect, when $\Delta
E_{21}=\Delta{m}^{2}_{21}/(2E)$ vanishes. Therefore, in comparison with the
standard PMNS matrix (the CP phase is defined as $\eta$), which is marked as
$U_{PMNS}$ in the following, we have

$\displaystyle
RU_{2}=U_{PMNS}=\exp(i\varphi_{23}\lambda_{7})\cdot\Gamma(\eta)\cdot\exp(i\varphi_{13}\lambda_{5})\cdot\Gamma(\eta)^{-1}\cdot\exp(i\varphi_{12}\lambda_{2})\,,$
(173)

where

$\displaystyle\Gamma(\eta)=\mbox{\rm
diag}\left(e^{-i\eta/2},0,e^{i\eta/2}\right)\,.$ (174)

Note that we have used $\varphi_{23}$, $\varphi_{13}$ and $\varphi_{12}$ to
avoid confusion with the original mixing angles $\theta_{23}$, $\theta_{13}$,
and $\theta_{12}$ in the $4\times 4$ matrix $U$ in Eq. (25). To make the sign
of the first row consistent with the standard representation $U_{PMNS}$, we
have to introduce a phase matrix, and we have

$\displaystyle RU_{2}$
$\displaystyle=\left(\begin{array}[]{ccc}-s_{24}c_{34}&-s_{34}&c_{24}c_{34}\\\
-s_{24}s_{34}&c_{34}&c_{24}s_{34}\\\ c_{24}&0&s_{24}\end{array}\right)$ (178)
$\displaystyle=\mbox{\rm diag}\left(-1,1,1\right)\cdot U_{PMNS}$
$\displaystyle=\left(\begin{array}[]{ccc}-C_{12}C_{13}&-S_{12}C_{13}&-e^{-i\eta}\,S_{13}\\\
-S_{12}C_{23}-e^{i\eta}\,C_{12}S_{13}S_{23}&C_{12}C_{23}-e^{i\eta}\,S_{12}S_{13}S_{23}&C_{13}S_{23}\\\
S_{12}S_{23}-e^{i\eta}\,C_{12}S_{13}C_{23}&-C_{12}S_{23}-e^{i\eta}\,S_{12}S_{13}C_{23}&C_{13}C_{23}\end{array}\right)\,,$
(182)

where $\mbox{\rm diag}\,(-1,1,1)$ is a diagonal matrix with a phase that does
not affect the oscillation probability, and we have defined
$s_{ij}=\sin{\theta_{ij}}$, $c_{ij}=\cos{\theta_{ij}}$,
$S_{ij}=\sin{\varphi_{ij}}$ and $C_{ij}=\cos{\varphi_{ij}}$. From Eq. (182),
we observe that the CP phase must satisfy $\eta=\pi$ for variables $S_{ij}$
and $C_{ij}$, which are semipositive definite, to satisfy
$[RU_{2}]_{32}=-C_{12}S_{23}-e^{i\eta}\,S_{12}S_{13}C_{23}=0$. In this case,
we have

$\displaystyle\tan\varphi_{12}=\frac{\tan\varphi_{23}}{S_{13}}\ ;$ (183)

from this we can derive the following three components of $U_{PMNS}$:

$\displaystyle[RU_{2}]_{21}=$ $\displaystyle-
S_{12}C_{23}-e^{i\eta}\,C_{12}S_{13}S_{23}$ $\displaystyle=$ $\displaystyle-
S_{12}C_{23}+C_{12}S_{13}S_{23}$ $\displaystyle=$ $\displaystyle
C_{12}\left(-C_{23}\tan\varphi_{12}+S_{13}S_{23}\right)$ $\displaystyle=$
$\displaystyle
C_{12}\left(-C_{23}\frac{\tan\varphi_{23}}{S_{13}}+S_{13}S_{23}\right)$
$\displaystyle=$ $\displaystyle-\frac{C_{12}S_{23}C_{13}^{2}}{S_{13}}<0\,,$
(184) $\displaystyle[RU_{2}]_{22}=$ $\displaystyle
C_{12}C_{23}-e^{i\eta}\,S_{12}S_{13}S_{23}$ $\displaystyle=$ $\displaystyle
C_{12}C_{23}+S_{12}S_{13}S_{23}$ $\displaystyle=$ $\displaystyle
C_{12}\left(C_{23}+S_{13}S_{23}\tan\varphi_{12}\right)$ $\displaystyle=$
$\displaystyle
C_{12}\left(C_{23}+S_{13}S_{23}\frac{\tan\varphi_{23}}{S_{13}}\right)$
$\displaystyle=$
$\displaystyle\frac{C_{12}}{C_{23}}\left(C_{23}^{2}+S_{23}^{2}\right)=\frac{C_{12}}{C_{23}}>0\,,$
(185) $\displaystyle[RU_{2}]_{31}=$ $\displaystyle
S_{12}S_{23}-e^{i\eta}\,C_{12}S_{13}C_{23}$ $\displaystyle=$ $\displaystyle
S_{12}S_{23}+C_{12}S_{13}C_{23}$ $\displaystyle=$ $\displaystyle
C_{12}\left(S_{23}\tan\varphi_{12}+S_{13}C_{23}\right)$ $\displaystyle=$
$\displaystyle
C_{12}\left(S_{23}\frac{\tan\varphi_{23}}{S_{13}}+S_{13}C_{23}\right)$
$\displaystyle=$
$\displaystyle\frac{C_{12}}{C_{23}S_{13}}\left(S_{23}^{2}+S_{13}^{2}C_{23}^{2}\right)>0\,,$
(186)

and obtain

$\displaystyle\mbox{\rm diag}\left(-1,1,1\right)\cdot U_{PMNS}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccc}-C_{12}C_{13}&-S_{12}C_{13}&S_{13}\\\
\displaystyle-\frac{C_{12}S_{23}C_{13}^{2}}{S_{13}}&\displaystyle\frac{C_{12}}{C_{23}}&C_{13}S_{23}\\\
\displaystyle\frac{C_{12}}{C_{23}S_{13}}\left(S_{23}^{2}+S_{13}^{2}C_{23}^{2}\right)&0&C_{13}C_{23}\end{array}\right)\,.$
(190)

From Eqs. (182) and (190) we obtain

$\displaystyle\tan{\varphi_{12}}$ $\displaystyle=\frac{t_{34}}{s_{24}}$ (59)
$\displaystyle\sin{\varphi_{13}}$ $\displaystyle=c_{24}c_{34}$ (60)
$\displaystyle\tan{\varphi_{23}}$ $\displaystyle=\frac{s_{34}}{t_{24}}\,,$
(61)

where $t_{ij}=\tan{\theta_{ij}}~{}(i,j=1,\cdots,4)$. Eqs. (59), (60), (61) are
consistent with all elements of $RU_{2}$ in Eq. (182) for the region
$0\leq\theta_{i4}\leq\pi/2~{}(i=2,3)$,
$0\leq\varphi_{ij}\leq\pi/2~{}(i,j=1,2,3)$. Based on Eqs. (165), (172) and
(173), we can rewrite the Hamiltonian $H_{2}$ as follows:

$\displaystyle H_{2}={}$ $\displaystyle R\exp(i\varphi_{23}\lambda_{7})$
$\displaystyle\times\left[\exp(i\varphi_{13}\lambda_{5})\left(\begin{array}[]{ccc}0&0&0\\\
0&0&0\\\ 0&0&\Delta
E\end{array}\right)\exp(-i\varphi_{13}\lambda_{5})+\left(\begin{array}[]{ccc}-A_{s}&0&0\\\
0&0&0\\\ 0&0&0\end{array}\right)\right]$ (197)
$\displaystyle\times\exp(-i\varphi_{23}\lambda_{7})R\,.$ (198)

Similarly to the case in Subsubsection 3.1.1, we have

$\displaystyle
H_{2}={}R\exp(i\varphi_{23}\lambda_{7})\exp(i\tilde{\varphi}_{13}\lambda_{5})$
(202) $\displaystyle\times\left(\begin{array}[]{ccc}\varepsilon_{2-}&0&0\\\
0&0&0\\\ 0&0&\varepsilon_{2+}\end{array}\right)$
$\displaystyle\times\exp(-i\tilde{\varphi}_{13}\lambda_{5})\exp(-i\varphi_{23}\lambda_{7})R\,.$
(203)

where

$\displaystyle\tan{2\tilde{\varphi}_{13}}=\frac{\Delta
E\sin{2\varphi_{13}}}{\Delta E\cos{2\varphi_{13}}+A_{s}}\,,$ (64)
$\displaystyle\varepsilon_{2\mp}=\frac{\Delta
E-A_{s}}{2}\mp\frac{1}{2}\Delta\tilde{E}_{2}\,,$ (65)

and

$\displaystyle\Delta\tilde{E}_{2}=\sqrt{(\Delta
E\cos{2\varphi_{13}}+A_{s})^{2}+(\Delta E\sin{2\varphi_{13}})^{2}}\,.$ (66)

Since the mixing between active and sterile neutrinos are expected to be very
weak in the (3+1)-scheme, we can consider

$\displaystyle|\theta_{24}|\ll 1,~{}~{}\ |\theta_{34}|\ll 1\,.$ (204)

Then, from (60) we obtain

$\displaystyle\varphi_{13}\simeq\frac{\pi}{2}\,.\hskip 85.35826pt$ (205)

Thus, Eq. (64) becomes

$\displaystyle\tan{2\tilde{\varphi}_{13}}\simeq\frac{\Delta
E\sin{2\varphi_{13}}}{-\Delta E+A_{s}}\,.$ (206)

Eq. (206) indicates that there is resonance in the disappearance channel
$\bar{\nu}_{\mu}\to\bar{\nu}_{\mu}$ and the appearance channel
$\bar{\nu}_{\mu}\to\bar{\nu}_{\tau}$. From these, the oscillation
probabilities $P(\bar{\nu}_{\mu}\to\bar{\nu}_{\mu})$ and
$P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})$ can be, respectively, expressed as
follows:

$\displaystyle P(\bar{\nu}_{\mu}\to\bar{\nu}_{\mu})={}$ $\displaystyle
1-C_{23}^{4}\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)$
$\displaystyle~{}~{}-\sin^{2}2\varphi_{23}\,\left\\{\tilde{S}_{13}^{2}\sin^{2}\left(\frac{\varepsilon_{2-}L}{2}\right)+\tilde{C}_{13}^{2}\sin^{2}\left(\frac{\varepsilon_{2+}L}{2}\right)\right\\}\,,$
(62) $\displaystyle P(\bar{\nu}_{\mu}\to\bar{\nu}_{\tau})={}$
$\displaystyle\sin^{2}2\varphi_{23}\left\\{\tilde{S}_{13}^{2}\sin^{2}\left(\frac{\varepsilon_{2-}L}{2}\right)+\tilde{C}_{13}^{2}\sin^{2}\left(\frac{\varepsilon_{2+}L}{2}\right)\right.$
$\displaystyle\hskip
56.9055pt-\left.\frac{1}{4}\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\right\\}\,,$
(63)

##  Appendix C Derivation of Eq. (73) in the case of $\theta_{14}\neq 0$, and
$\theta_{24}=\theta_{34}=0$

In this Appendix, we discuss oscillations of astrophysical neutrinos for
$\theta_{14}\neq 0$, and $\theta_{24}=\theta_{34}=0$. As in the discussion in
Subsubsection 3.1.1, we consider the case of the neutrino mode. Based on Eq.
(129) in Subsubsection 3.1.1, it is not difficult to explicitly obtain the
transition matrix $T$:

$\displaystyle T$ $\displaystyle=\exp(i\tilde{\theta}_{14}{\cal
T}_{14})\cdot\left(\begin{array}[]{cccc}e^{-iL\varepsilon_{1-}}&0&0&0\\\
0&1&0&0\\\ 0&0&1&0\\\
0&0&0&e^{-iL\varepsilon_{1+}}\end{array}\right)\cdot\exp(-i\tilde{\theta}_{14}{\cal
T}_{14})$ (211)
$\displaystyle=\left(\begin{array}[]{cccc}\tilde{c}_{14}^{2}e^{-iL\varepsilon_{1-}}+\tilde{s}_{14}^{2}e^{-iL\varepsilon_{1+}}&0&0&\tilde{s}_{14}\tilde{c}_{14}\left(e^{-iL\varepsilon_{1+}}-e^{-iL\varepsilon_{1-}}\right)\\\
0&1&0&0\\\ 0&0&1&0\\\
\tilde{s}_{14}\tilde{c}_{14}\left(e^{-iL\varepsilon_{1+}}-e^{-iL\varepsilon_{1-}}\right)&0&0&\tilde{s}_{14}^{2}e^{-iL\varepsilon_{1-}}+\tilde{c}_{14}^{2}e^{-iL\varepsilon_{1+}}\end{array}\right)\,,$
(216)

where

$\displaystyle{\cal T}_{14}=\left(\begin{array}[]{cccc}0&0&0&-i\\\ 0&0&0&0\\\
0&0&0&0\\\ i&0&0&0\end{array}\right).$ (221)

Here, we can write down the non-zero matrix elements of $T$: 666The two
indices of $T$ and one of the indices of $U$ must be that of flavor
$\alpha=e,\mu,\tau,s$. For simplicity, we denote $T_{ij}~{}(i,j=1,\cdots,4)$
and $U_{ij}~{}(i,j=1,\cdots,4)$ instead of
$T_{\alpha\beta}~{}(\alpha,\beta=e,\mu,\tau,s)$ and $U_{\alpha
i}~{}(\alpha=e,\mu,\tau,s;i=1,\cdots,4)$ in this appendix.

$\displaystyle\begin{aligned}
T_{11}&=\tilde{c}_{14}^{2}e^{-iL\varepsilon_{1-}}+\tilde{s}_{14}^{2}e^{-iL\varepsilon_{1+}}\,,\\\
T_{14}&=T_{41}=\tilde{s}_{14}\tilde{c}_{14}\left(e^{-iL\varepsilon_{1+}}-e^{-iL\varepsilon_{1-}}\right)\,,\\\
T_{22}&=T_{33}=1\,,\\\
T_{44}&=\tilde{s}_{14}^{2}e^{-iL\varepsilon_{1-}}+\tilde{c}_{14}^{2}e^{-iL\varepsilon_{1+}}\,,\end{aligned}$
(222)

where $\tilde{s}_{14}=\sin{\tilde{\theta}_{14}}$,
$\tilde{c}_{14}=\cos{\tilde{\theta}_{14}}$ and $\varepsilon_{1\mp}$ is given
by Eq. (130). From Eq. (222) we have the following:

$\displaystyle|T_{11}|^{2}$
$\displaystyle=1-\sin^{2}2\tilde{\theta}_{14}\sin^{2}\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)$
(223) $\displaystyle|T_{14}|^{2}$
$\displaystyle=\sin^{2}2\tilde{\theta}_{14}\sin^{2}\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)$
(74) $\displaystyle|T_{41}|^{2}$
$\displaystyle=\sin^{2}2\tilde{\theta}_{14}\sin^{2}\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)$
(224) $\displaystyle|T_{44}|^{2}$
$\displaystyle=1-\sin^{2}2\tilde{\theta}_{14}\sin^{2}\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)$
(225)

For later purposes, we keep the dependence of $U$ on $\theta_{14}$ but we use
the following approximation for the three-flavor mixing angles in $U$:

$\displaystyle\theta_{13}\simeq 0\,.$ (226)

We introduce the product of matrices $T$ and $U$:

$\displaystyle V$ $\displaystyle\equiv TU$
$\displaystyle\simeq\left(\begin{array}[]{cccc}T_{11}&T_{12}&T_{13}&T_{14}\\\
T_{21}&T_{22}&T_{23}&T_{24}\\\ T_{31}&T_{32}&T_{33}&T_{34}\\\
T_{41}&T_{42}&T_{43}&T_{44}\end{array}\right)\left(\begin{array}[]{cccc}c_{14}&0&0&s_{14}\\\
0&1&0&0\\\ 0&0&1&0\\\ -s_{14}&0&0&c_{14}\end{array}\right)$ (235)
$\displaystyle\quad\times\left(\begin{array}[]{cccc}c_{12}&s_{12}&0&0\\\
-s_{12}c_{23}&c_{12}c_{23}&s_{23}&0\\\ s_{12}s_{23}&-c_{12}s_{23}&c_{23}&0\\\
0&0&0&1\end{array}\right)$ (240)
$\displaystyle=\left(\begin{array}[]{cccc}c_{12}\,T_{1-}&s_{12}\,T_{1-}&0&T_{1+}\\\
-s_{12}c_{23}&c_{12}c_{23}&s_{23}&0\\\ s_{12}s_{23}&-c_{12}s_{23}&c_{23}&0\\\
c_{12}\,T_{4-}&s_{12}\,T_{4-}&0&T_{4+}\end{array}\right)$ (245)

where

$\displaystyle T_{1-}$ $\displaystyle\equiv c_{14}\,T_{11}-s_{14}\,T_{14}$
(246) $\displaystyle T_{1+}$ $\displaystyle\equiv
s_{14}\,T_{11}+c_{14}\,T_{14}$ (247) $\displaystyle T_{4-}$
$\displaystyle\equiv c_{14}\,T_{41}-s_{14}\,T_{44}$ (248) $\displaystyle
T_{4+}$ $\displaystyle\equiv s_{14}\,T_{41}+c_{14}\,T_{44}$ (249)

Then, from Eq. (72), we obtain the formulae of the oscillation probabilities:

$\displaystyle\quad P(\nu_{\alpha}\to\nu_{\beta})$
$\displaystyle=\sum_{i}|V_{\beta i}|^{2}|U_{\alpha i}|^{2}$
$\displaystyle\simeq\left[\left(\begin{array}[]{cccc}c_{12}^{2}|T_{1-}|^{2}&s_{12}^{2}|T_{1-}|^{2}&0&|T_{1+}|^{2}\\\
s_{12}^{2}c_{23}^{2}&c_{12}^{2}c_{23}^{2}&s_{23}^{2}&0\\\
s_{12}^{2}s_{23}^{2}&c_{12}^{2}s_{23}^{2}&c_{23}^{2}&0\\\
c_{12}^{2}|T_{4-}|^{2}&s_{12}^{2}|T_{4-}|^{2}&0&|T_{4+}|^{2}\end{array}\right)\left(\begin{array}[]{cccc}c_{12}^{2}c_{14}^{2}&s_{12}^{2}c_{23}^{2}&s_{12}^{2}s_{23}^{2}&c_{12}^{2}s_{14}^{2}\\\
s_{12}^{2}c_{14}^{2}&c_{12}^{2}c_{23}^{2}&c_{12}^{2}s_{23}^{2}&s_{12}^{2}s_{14}^{2}\\\
0&s_{23}^{2}&c_{23}^{2}&0\\\
s_{14}^{2}&0&0&c_{14}^{2}\end{array}\right)\right]_{\beta\alpha}\,.$ (258)

Each probability is as follows:

$\displaystyle P(\nu_{e}\to\nu_{e})$
$\displaystyle=c_{14}^{2}\left(s_{12}^{4}+c_{12}^{4}\right)\,|T_{1-}|^{2}+s_{14}^{2}\,|T_{1+}|^{2}\,,$
(259) $\displaystyle P(\nu_{e}\to\nu_{\mu})$
$\displaystyle=2c_{14}^{2}s_{12}^{2}c_{12}^{2}c_{23}^{2}\,,$ (260)
$\displaystyle P(\nu_{e}\to\nu_{\tau})$
$\displaystyle=2c_{14}^{2}s_{12}^{2}c_{12}^{2}s_{23}^{2}\,,$ (261)
$\displaystyle P(\nu_{e}\to\nu_{s})$
$\displaystyle=c_{14}^{2}\left(s_{12}^{4}+c_{12}^{4}\right)\,|T_{4-}|^{2}+s_{14}^{2}\,|T_{4+}|^{2}\
;$ (262)

and

$\displaystyle
P(\nu_{\mu}\to\nu_{e})=2s_{12}^{2}c_{12}^{2}c_{23}^{2}\,|T_{1-}|^{2}\,,$ (263)
$\displaystyle
P(\nu_{\mu}\to\nu_{\mu})=\left(s_{12}^{4}+c_{12}^{4}\right)c_{23}^{4}+s_{23}^{4}\,,$
(264) $\displaystyle
P(\nu_{\mu}\to\nu_{\tau})=\left(s_{12}^{4}+c_{12}^{4}+1\right)s_{23}^{2}c_{23}^{2}\,,$
(265) $\displaystyle
P(\nu_{\mu}\to\nu_{s})=2s_{12}^{2}c_{12}^{2}c_{23}^{2}\,|T_{4-}|^{2}\,.$ (266)

It is known that at the production source of astrophysical neutrinos, such as
active galactic nuclei, the ratios of the flux of $\nu_{e}$, $\nu_{\mu}$ and
$\nu_{\tau}$ is 1:2:0. Moreover, it is expected that there is no sterile
neutrino at the source. Denoting the flux of astrophysical neutrinos at the
source as $F^{0}(\nu_{\alpha})\ (\alpha=e,\mu,\tau,s)$, we have
$F^{0}(\nu_{e}):F^{0}(\nu_{\mu}):F^{0}(\nu_{\tau}):F^{0}(\nu_{s})=1:2:0:0$.
From Eq. (4), to a good approximation, we can assume that
$\theta_{23}\simeq\pi/4$, so that $s_{23}\simeq c_{23}\simeq 1/\sqrt{2}$. Thus
we can estimate the flux of $\nu_{e}$, $\nu_{\mu}$, $\nu_{\tau}$ and
$\nu_{s}$, which is observed at IceCube is given in the unit of
$F^{0}(\nu_{e})$ by:

$\displaystyle F(\nu_{e})$
$\displaystyle=P(\nu_{e}\to\nu_{e})+2P(\nu_{\mu}\to\nu_{e})$
$\displaystyle\simeq\left\\{c_{14}^{2}\left(s_{12}^{4}+c_{12}^{4}\right)+2s_{12}^{2}c_{12}^{2}\right\\}\,|T_{1-}|^{2}+s_{14}^{2}|T_{1+}|^{2}\,,$
(267) $\displaystyle F(\nu_{\mu})$
$\displaystyle=P(\nu_{e}\to\nu_{\mu})+2P(\nu_{\mu}\to\nu_{\mu})$
$\displaystyle\simeq
c_{14}^{2}s_{12}^{2}c_{12}^{2}+\frac{1}{2}\left(s_{12}^{4}+c_{12}^{4}+1\right)\,,$
(268) $\displaystyle F(\nu_{\tau})$
$\displaystyle=P(\nu_{e}\to\nu_{\tau})+2P(\nu_{\mu}\to\nu_{\tau})$
$\displaystyle\simeq
c_{14}^{2}s_{12}^{2}c_{12}^{2}+\frac{1}{2}\left(s_{12}^{4}+c_{12}^{4}+1\right)\,,$
(269) $\displaystyle F(\nu_{s})$
$\displaystyle=P(\nu_{e}\to\nu_{s})+2P(\nu_{\mu}\to\nu_{s})$
$\displaystyle\simeq\left\\{c_{14}^{2}\left(s_{12}^{4}+c_{12}^{4}\right)+2s_{12}^{2}c_{12}^{2}\right\\}\,|T_{4-}|^{2}+s_{14}^{2}|T_{4+}|^{2}\,.$
(270)

If we use approximation $\theta_{14}\simeq 0$ except in $T_{i\pm}~{}(i=1,4)$
while keeping $\tilde{\theta}_{14}\neq 0$ in each flux, from Eqs. (246),
(248), (223) and (74), we have

$\displaystyle|T_{1-}|^{2}$
$\displaystyle\simeq|T_{11}|^{2}=1-\sin^{2}2\tilde{\theta}_{14}\sin^{2}\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)$
(271) $\displaystyle|T_{4-}|^{2}$
$\displaystyle\simeq|T_{41}|^{2}=\sin^{2}2\tilde{\theta}_{14}\sin^{2}\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)$
(272)

Hence, the ratio of each flavor of neutrino can be approximated as

$\displaystyle F(\nu_{e}):F(\nu_{\mu}):F(\nu_{\tau}):F(\nu_{s})$
$\displaystyle\simeq 1-|T_{14}|^{2}:1:1:|T_{14}|^{2}\,,$ (73)

where $|T_{14}|^{2}$ is defined in Eq. (74) and $|T_{14}|^{2}$ has nontrivial
dependence on the neutrino energy $E$, especially near the resonance region
$E{\color[rgb]{0,0,0}=\Delta{m}^{2}_{41}\cos{2\theta_{14}}/\\{2(A_{e}-A_{s})\\}}\sim\Delta
m^{2}_{41}/\\{2(A_{e}-A_{s})\\}$, where we have assumed $|\theta_{14}|\ll 1$.

##  Appendix D Derivation of Eq. (75) in the case of $\theta_{14}=0$,
$\theta_{24}\neq 0$, and $\theta_{34}\neq 0$

In this Appendix, we discuss oscillations of astrophysical neutrinos for
$\theta_{14}=0$, $\theta_{24}\neq 0$ and $\theta_{34}\neq 0$. As in
Subsubsection 3.1.2, it is useful to discuss in terms of the new mixing angles
in Eqs. (59), (60), (61). According to Eq. (203) in Subsubsection 3.1.2, we
can explicitly express the transition matrix $T$:

$\displaystyle T=\left(\begin{array}[]{cccc}e^{-iLA_{e}}&0&0&0\\\
0&T_{22}&T_{23}&T_{24}\\\ 0&T_{32}&T_{33}&T_{34}\\\
0&T_{42}&T_{43}&T_{44}\end{array}\right)\,.$ (277)

Here, we denote the submatrix of $T$ by $T_{0}$:

$\displaystyle T_{0}$
$\displaystyle=\left(\begin{array}[]{ccc}T_{22}&T_{23}&T_{24}\\\
T_{32}&T_{33}&T_{34}\\\ T_{42}&T_{43}&T_{44}\end{array}\right)$ (281)
$\displaystyle=R\,e^{i\varphi_{23}\lambda_{7}}e^{i\tilde{\varphi}_{13}\lambda_{5}}\left(\begin{array}[]{ccc}e^{-iL\varepsilon_{2-}}&0&0\\\
0&1&0\\\
0&0&e^{-iL\varepsilon_{2+}}\end{array}\right)e^{-i\tilde{\varphi}_{13}\lambda_{5}}e^{-i\varphi_{23}\lambda_{7}}R\,,$
(285)

where

$\displaystyle\begin{aligned}
T_{22}&=\left(\tilde{S}_{13}^{2}e^{-iL\varepsilon_{2-}}+\tilde{C}_{13}^{2}e^{-iL\varepsilon_{2+}}\right)C_{23}^{2}+S_{23}^{2}\,,\\\
T_{23}&=T_{32}=\left(\tilde{S}_{13}^{2}e^{-iL\varepsilon_{2-}}+\tilde{C}_{13}^{2}e^{-iL\varepsilon_{2+}}-1\right)S_{23}C_{23}\,,\\\
T_{24}&=T_{42}=\left(e^{-iL\varepsilon_{2+}}-e^{-iL\varepsilon_{2-}}\right)\tilde{S}_{13}\tilde{C}_{13}C_{23}\,,\\\
T_{33}&=\left(\tilde{S}_{13}^{2}e^{-iL\varepsilon_{2-}}+\tilde{C}_{13}^{2}e^{-iL\varepsilon_{2+}}\right)S_{23}^{2}+C_{23}^{2}\,,\\\
T_{34}&=T_{43}=\left(e^{-iL\varepsilon_{2+}}-e^{-iL\varepsilon_{2-}}\right)\tilde{S}_{13}\tilde{C}_{13}S_{23}\,,\\\
T_{44}&=\tilde{C}_{13}^{2}e^{-iL\varepsilon_{2-}}+\tilde{S}_{13}^{2}e^{-iL\varepsilon_{2+}}\,.\end{aligned}$
(286)

From Eqs. (286), we have the following:

$\displaystyle|T_{22}|^{2}$
$\displaystyle=S_{23}^{4}+C_{23}^{4}\left\\{1-\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\right\\}$
$\displaystyle\quad+2C_{23}^{2}S_{23}^{2}\,\mbox{\rm Re}\left[e^{-i(\Delta
E-A_{s})L/2}\,\left\\{\cos\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)-i\cos
2\tilde{\varphi}_{13}\sin\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\right\\}\right]\,,$
(287) $\displaystyle|T_{23}|^{2}$ $\displaystyle=|T_{32}|^{2}$
$\displaystyle=C_{23}^{2}S_{23}^{2}\,\left\\{2-\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\right\\}$
$\displaystyle\quad-2C_{23}^{2}S_{23}^{2}\,\mbox{\rm Re}\left[e^{-i(\Delta
E-A_{s})L/2}\,\left\\{\cos\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)-i\cos
2\tilde{\varphi}_{13}\sin\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\right\\}\right]\,,$
(288) $\displaystyle|T_{24}|^{2}$ $\displaystyle=|T_{42}|^{2}$
$\displaystyle=C_{23}^{2}\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\,,$
(289) $\displaystyle|T_{33}|^{2}$
$\displaystyle=C_{23}^{4}+S_{23}^{4}\left\\{1-\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\right\\}$
$\displaystyle\quad+2C_{23}^{2}S_{23}^{2}\,\mbox{\rm Re}\left[e^{-i(\Delta
E-A_{s})L/2}\,\left\\{\cos\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)-i\cos
2\tilde{\varphi}_{13}\sin\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\right\\}\right]\,,$
(290) $\displaystyle|T_{34}|^{2}$ $\displaystyle=|T_{43}|^{2}$
$\displaystyle=S_{23}^{2}\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\,,$
(291) $\displaystyle|T_{44}|^{2}$
$\displaystyle=1-\sin^{2}2\tilde{\varphi}_{13}\,\sin^{2}\left(\frac{\Delta\tilde{E}_{2}L}{2}\right)\,.$
(76)

Here, we have defined $\tilde{S}_{13}=\sin{\tilde{\varphi}_{13}}$,
$\tilde{C}_{13}=\cos{\tilde{\varphi}_{13}}$, $S_{23}=\sin{\varphi_{23}}$,
$C_{23}=\cos{\varphi_{23}}$ and $\varepsilon_{2\mp}$ is given by Eq. (65).

Then, according to Eq. (72), using approximation $\theta_{13}\simeq 0$, we can
obtain the formulae of oscillation probabilities:

$\displaystyle P(\nu_{e}\to\nu_{e})$
$\displaystyle=1-2s_{12}^{2}c_{12}^{2}\,,$ (292) $\displaystyle
P(\nu_{e}\to\nu_{\mu})$ $\displaystyle=2|T_{22}|^{2}\cdot
s_{12}^{2}c_{12}^{2}c_{23}^{2}-4\mbox{\rm Re}\,(T_{22}T_{23}^{\ast})\cdot
s_{12}^{2}c_{12}^{2}s_{23}c_{23}+2|T_{23}|^{2}\cdot
s_{12}^{2}c_{12}^{2}s_{23}^{2}\,,$ (293) $\displaystyle
P(\nu_{e}\to\nu_{\tau})$ $\displaystyle=2|T_{32}|^{2}\cdot
s_{12}^{2}c_{12}^{2}c_{23}^{2}-4\mbox{\rm Re}\,(T_{32}T_{33}^{\ast})\cdot
s_{12}^{2}c_{12}^{2}s_{23}c_{23}+2|T_{33}|^{2}\cdot
s_{12}^{2}c_{12}^{2}s_{23}^{2}\,,$ (294) $\displaystyle P(\nu_{e}\to\nu_{s})$
$\displaystyle=2|T_{42}|^{2}\cdot s_{12}^{2}c_{12}^{2}c_{23}^{2}-4\mbox{\rm
Re}\,(T_{42}T_{43}^{\ast})\cdot
s_{12}^{2}c_{12}^{2}s_{23}c_{23}+2|T_{43}|^{2}\cdot
s_{12}^{2}c_{12}^{2}s_{23}^{2}\ ;$ (295)

and

$\displaystyle P(\nu_{\mu}\to\nu_{e})={}$ $\displaystyle
2s_{12}^{2}c_{12}^{2}c_{23}^{2}\,,$ (296) $\displaystyle
P(\nu_{\mu}\to\nu_{\mu})={}$
$\displaystyle|T_{22}|^{2}\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)c_{23}^{4}+s_{23}^{4}\right]$
$\displaystyle-2\mbox{\rm
Re}\,(T_{22}T_{23}^{\ast})\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}c_{23}^{3}-s_{23}^{3}c_{23}\right]$
$\displaystyle+|T_{23}|^{2}\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}^{2}c_{23}^{2}+s_{23}^{2}c_{23}^{2}\right]\,,$
(297) $\displaystyle P(\nu_{\mu}\to\nu_{\tau})={}$
$\displaystyle|T_{32}|^{2}\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)c_{23}^{4}+s_{23}^{4}\right]$
$\displaystyle-2\mbox{\rm
Re}\,(T_{32}T_{33}^{\ast})\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}c_{23}^{3}-s_{23}^{3}c_{23}\right]$
$\displaystyle+|T_{33}|^{2}\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}^{2}c_{23}^{2}+s_{23}^{2}c_{23}^{2}\right]\,,$
(298) $\displaystyle P(\nu_{\mu}\to\nu_{s})={}$
$\displaystyle|T_{42}|^{2}\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)c_{23}^{4}+s_{23}^{4}\right]$
$\displaystyle-2\mbox{\rm
Re}\,(T_{42}T_{43}^{\ast})\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}c_{23}^{3}-s_{23}^{3}c_{23}\right]$
$\displaystyle+|T_{43}|^{2}\cdot\left[\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}^{2}c_{23}^{2}+s_{23}^{2}c_{23}^{2}\right]\,.$
(299)

As in Subsubsection 3.2.1, using approximation $\theta_{23}\simeq\pi/4$, we
obtain the flux of $\nu_{e}$, $\nu_{\mu}$ and $\nu_{s}$ observed at IceCube:

$\displaystyle F(\nu_{e})={}$ $\displaystyle
P(\nu_{e}\to\nu_{e})+2P(\nu_{\mu}\to\nu_{e})$ $\displaystyle\simeq{}$
$\displaystyle s_{12}^{4}+c_{12}^{4}+4s_{12}^{2}c_{12}^{2}c_{23}^{2}$
$\displaystyle\simeq{}$ $\displaystyle 1\,,$ (300) $\displaystyle
F(\nu_{\mu})={}$ $\displaystyle
P(\nu_{e}\to\nu_{\mu})+2P(\nu_{\mu}\to\nu_{\mu})$ $\displaystyle\simeq{}$
$\displaystyle|T_{22}|^{2}\cdot\left[2s_{12}^{2}c_{12}^{2}c_{23}^{2}+2\left(s_{12}^{4}+c_{12}^{4}\right)c_{23}^{4}+2s_{23}^{4}\right]$
$\displaystyle-4\mbox{\rm
Re}\,(T_{22}T_{23}^{\ast})\cdot\left[s_{12}^{2}c_{12}^{2}s_{23}c_{23}+\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}c_{23}^{3}-s_{23}^{3}c_{23}\right]$
$\displaystyle+|T_{23}|^{2}\cdot\left[2s_{12}^{2}c_{12}^{2}s_{23}^{2}+2\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}^{2}c_{23}^{2}+2s_{23}^{2}c_{23}^{2}\right]$
$\displaystyle\simeq{}$ $\displaystyle|T_{22}|^{2}+|T_{23}|^{2}\,,$ (301)
$\displaystyle F(\nu_{\tau})={}$ $\displaystyle
P(\nu_{e}\to\nu_{\tau})+2P(\nu_{\mu}\to\nu_{\tau})$ $\displaystyle\simeq{}$
$\displaystyle|T_{32}|^{2}\cdot\left[2s_{12}^{2}c_{12}^{2}c_{23}^{2}+2\left(s_{12}^{4}+c_{12}^{4}\right)c_{23}^{4}+2s_{23}^{4}\right]$
$\displaystyle-4\mbox{\rm
Re}\,(T_{32}T_{33}^{\ast})\cdot\left[s_{12}^{2}c_{12}^{2}s_{23}c_{23}+\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}c_{23}^{3}-s_{23}^{3}c_{23}\right]$
$\displaystyle+|T_{33}|^{2}\cdot\left[2s_{12}^{2}c_{12}^{2}s_{23}^{2}+2\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}^{2}c_{23}^{2}+2s_{23}^{2}c_{23}^{2}\right]$
$\displaystyle\simeq{}$ $\displaystyle|T_{32}|^{2}+|T_{33}|^{2}\,,$ (302)
$\displaystyle F(\nu_{s})={}$ $\displaystyle
P(\nu_{e}\to\nu_{s})+2P(\nu_{\mu}\to\nu_{s})$ $\displaystyle\simeq{}$
$\displaystyle|T_{42}|^{2}\cdot\left[2s_{12}^{2}c_{12}^{2}c_{23}^{2}+2\left(s_{12}^{4}+c_{12}^{4}\right)c_{23}^{4}+2s_{23}^{4}\right]$
$\displaystyle-4\mbox{\rm
Re}\,(T_{42}T_{43}^{\ast})\cdot\left[s_{12}^{2}c_{12}^{2}s_{23}c_{23}+\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}c_{23}^{3}-s_{23}^{3}c_{23}\right]$
$\displaystyle+|T_{43}|^{2}\cdot\left[2s_{12}^{2}c_{12}^{2}s_{23}^{2}+2\left(s_{12}^{4}+c_{12}^{4}\right)s_{23}^{2}c_{23}^{2}+2s_{23}^{2}c_{23}^{2}\right]$
$\displaystyle\simeq{}$ $\displaystyle|T_{42}|^{2}+|T_{43}|^{2}\,.$ (303)

Hence the flavor ratio is given by

$\displaystyle F(\nu_{e}):F(\nu_{\mu}):F(\nu_{\tau}):F(\nu_{s})$
$\displaystyle\simeq{}$ $\displaystyle
1:1-|T_{24}|^{2}:1-|T_{34}|^{2}:1-|T_{44}|^{2}$ $\displaystyle={}$
$\displaystyle
1:1-C_{23}^{2}\,(1-|T_{44}|^{2}):1-S_{23}^{2}\,(1-|T_{44}|^{2}):1-|T_{44}|^{2}\,,$
(75)

where $|T_{44}|^{2}$ is defined in Eq. (76) and $|T_{44}|^{2}$ has nontrivial
dependence on the neutrino energy $E$, especially near the resonance region
$E{\color[rgb]{0,0,0}=-\Delta m^{2}_{41}\cos 2\varphi_{13}/(2A_{s})}\sim\Delta
m^{2}_{41}/(2A_{s})$, where we have assumed $|\pi/2-\varphi_{13}|\ll 1$.

##  Appendix E The behavior of astrophysical $\nu_{e}$ as
$\theta_{14}\to\pi/4$

If $\theta_{14}$ is very close to $\pi/4$ and $|\Delta E|\ll|A_{e}-A_{s}|$,
from Eqs. (55) and (56), we obtain

$\displaystyle\cos 2\tilde{\theta}_{14}=\frac{\Delta
E\cos{2\theta_{14}}-A_{e}+A_{s}}{\Delta\tilde{E}_{1}}\simeq-1\,,$ (304)

which indicates

$\displaystyle\tilde{\theta}_{14}\simeq\frac{\pi}{2}\,.$ (305)

From Eq. (246), we obtain

$\displaystyle T_{1-}$
$\displaystyle\simeq\frac{1}{\sqrt{2}}\,\left(T_{11}-T_{14}\right)$ (306)

We can now calculate the factor $|T_{1-}|^{2}$ in Eq. (267) as
$\theta_{14}\to\pi/4$ as follows.

$\displaystyle|T_{1-}|^{2}$
$\displaystyle\simeq\frac{1}{2}\,\left|T_{11}-T_{14}\right|^{2}$
$\displaystyle=\frac{1}{2}\,\left\\{|T_{11}|^{2}+|T_{14}|^{2}-2\,\mbox{\rm
Re}\,\left(T_{11}T_{14}^{\ast}\right)\right\\}$
$\displaystyle\simeq\frac{1}{2}$ (307)

where we used the following properties and Eqs. (223) and (74):

$\displaystyle T_{11}$ $\displaystyle=e^{-i(\Delta
E+A_{e}+A_{s})L/2}\left\\{\cos\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)+i\cos
2\tilde{\theta}_{14}\sin\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)\right\\}$
(308) $\displaystyle T_{14}$ $\displaystyle=-ie^{-i(\Delta
E+A_{e}+A_{s})L/2}\sin
2\tilde{\theta}_{14}\sin\left(\frac{\Delta\tilde{E}_{1}L}{2}\right)\simeq
0\,.$ (309)

Eq. (307) implies that the oscillation probability of $\nu_{e}$ flux is
independent of the neutrino energy $E$ and baseline length $L$ in the Earth.

## References

  * (1)
  * (2) P. A. Zyla et al. [Particle Data Group], PTEP 2020 (2020) no.8, 083C01 doi:10.1093/ptep/ptaa104
  * (3) A. Aguilar-Arevalo et al. [LSND], Phys. Rev. D 64 (2001), 112007 doi:10.1103/PhysRevD.64.112007 [arXiv:hep-ex/0104049 [hep-ex]].
  * (4) A. A. Aguilar-Arevalo et al. [MiniBooNE], Phys. Rev. Lett. 121 (2018) no.22, 221801 doi:10.1103/PhysRevLett.121.221801 [arXiv:1805.12028 [hep-ex]].
  * (5) K. N. Abazajian, M. A. Acero, S. K. Agarwalla, A. A. Aguilar-Arevalo, C. H. Albright, S. Antusch, C. A. Arguelles, A. B. Balantekin, G. Barenboim and V. Barger, et al. [arXiv:1204.5379 [hep-ph]].
  * (6) C. Giunti and M. Laveder, Phys. Rev. C 83 (2011), 065504 doi:10.1103/PhysRevC.83.065504 [arXiv:1006.3244 [hep-ph]].
  * (7) G. Mention, M. Fechner, T. Lasserre, T. A. Mueller, D. Lhuillier, M. Cribier and A. Letourneau, Phys. Rev. D 83 (2011), 073006 doi:10.1103/PhysRevD.83.073006 [arXiv:1101.2755 [hep-ex]].
  * (8) P. Huber, Phys. Rev. C 84 (2011), 024617 [erratum: Phys. Rev. C 85 (2012), 029901] doi:10.1103/PhysRevC.85.029901 [arXiv:1106.0687 [hep-ph]].
  * (9) D. Adey et al. [Daya Bay], Phys. Rev. Lett. 123 (2019) no.11, 111801 doi:10.1103/PhysRevLett.123.111801 [arXiv:1904.07812 [hep-ex]].
  * (10) A. C. Hayes, J. L. Friar, G. T. Garvey, G. Jungman and G. Jonkmans, Phys. Rev. Lett. 112 (2014), 202501 doi:10.1103/PhysRevLett.112.202501 [arXiv:1309.4146 [nucl-th]].
  * (11) P. Adamson et al. [MINOS+ and Daya Bay], Phys. Rev. Lett. 125 (2020) no.7, 071801 doi:10.1103/PhysRevLett.125.071801 [arXiv:2002.00301 [hep-ex]].
  * (12) A. Nicolaidis, G. Tsirigoti and J. Hansson, [arXiv:hep-ph/9904415 [hep-ph]].
  * (13) O. Yasuda, doi:10.1142/9789812777003_0017 [arXiv:hep-ph/0102166 [hep-ph]].
  * (14) H. Nunokawa, O. L. G. Peres and R. Zukanovich Funchal, Phys. Lett. B 562 (2003), 279-290 doi:10.1016/S0370-2693(03)00603-8 [arXiv:hep-ph/0302039 [hep-ph]].
  * (15) M. G. Aartsen et al. [IceCube], Phys. Rev. D 102 (2020) no.5, 052009 doi:10.1103/PhysRevD.102.052009 [arXiv:2005.12943 [hep-ex]].
  * (16) M. Maltoni, T. Schwetz, M. A. Tortola and J. W. F. Valle, New J. Phys. 6 (2004), 122 doi:10.1088/1367-2630/6/1/122 [arXiv:hep-ph/0405172 [hep-ph]].
  * (17) V. A. Naumov, Phys. Lett. B 529 (2002), 199-211 doi:10.1016/S0370-2693(02)01258-3 [arXiv:hep-ph/0112249 [hep-ph]].
  * (18) M. G. Aartsen et al. [IceCube], Phys. Rev. D 91 (2015), 122004 doi:10.1103/PhysRevD.91.122004 [arXiv:1504.03753 [astro-ph.HE]].
  * (19) G. L. Fogli, E. Lisi, A. Marrone and D. Montanino, Phys. Lett. B 425 (1998), 341-344 doi:10.1016/S0370-2693(98)00077-X [arXiv:hep-ph/9711421 [hep-ph]].
  * (20) R. Foot, R. R. Volkas and O. Yasuda, Phys. Lett. B 433 (1998), 82-87 doi:10.1016/S0370-2693(98)00706-0 [arXiv:hep-ph/9802287 [hep-ph]].
  * (21) J. G. Learned and S. Pakvasa, Astropart. Phys. 3 (1995), 267-274 doi:10.1016/0927-6505(94)00043-3 [arXiv:hep-ph/9405296 [hep-ph]].
  * (22) R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Phys. Rev. D 58 (1998), 093009 doi:10.1103/PhysRevD.58.093009 [arXiv:hep-ph/9807264 [hep-ph]].
  * (23) M. Honda, T. Kajita, K. Kasahara, S. Midorikawa and T. Sanuki, Phys. Rev. D 75 (2007), 043006 doi:10.1103/PhysRevD.75.043006 [arXiv:astro-ph/0611418 [astro-ph]].
  * (24) D. Chirkin, [arXiv:hep-ph/0407078 [hep-ph]].