---
author:
- 'Chee Sheng Fong$^{1}$'
- 'Hisakazu Minakata$^{2,3}$'
- 'Hiroshi Nunokawa$^{4}$'
title: 'Non-unitary evolution of neutrinos in matter and the leptonic unitarity test\'
---

IFT-UAM/CSIC-17-117

Introduction {#sec:introduction}
============

Studies of neutrino oscillation entered into a “matured phase” after the structure of the three-flavor lepton mixing [@Maki:1962mu] is elucidated. The long-lasted discovery phase of neutrino oscillation has been unambiguously concluded by the Super-Kamiokande (Super-K) atmospheric neutrino observation which discovered neutrino oscillation and hence neutrino mass [@Fukuda:1998mi]. It was followed by the KamLAND reactor and the solar neutrino experiments which uncovered the three-flavor nature of the mixing by observing oscillation and/or nonadiabatic flavor conversion of neutrinos in matter [@Mikheev:1986gs; @Wolfenstein:1977ue] in the 1-2 sector [@Eguchi:2002dm; @Ahmad:2002jz].[^1] The last step of understanding the three-flavor structure of neutrino oscillation was carried out by the reactor [@An:2016ses; @RENO:2015ksa; @Schoppmann:2016iww] and the accelerator [@Abe:2017vif; @Adamson:2016tbq] measurement of $\theta_{13}$. It lefts only the two remaining unknowns in the standard three-flavor mixing paradigm, that is, measurement of CP violating phase and determination of neutrino mass ordering.[^2]

The completion of the theory of the three-flavor neutrino mixing, however, necessitates the paradigm test. A well-known example of such efforts is to verify unitarity of the quark CKM matrix [@Olive:2016xmw]. We have argued in ref. [@Fong:2016yyh] that we may need a different strategy to test leptonic unitarity. That is, first prepare a generic framework which describes unitarity violation at certain energy scale, and then confront it to experimental data. We contrasted the two typical alternatives, unitarity violation by new physics at high ($E \gg m_{W}$) and low ($E \ll m_{W}$) energy scales, which are dubbed as high-scale and low-scale unitarity violation, respectively.[^3] They differ in certain characteristic features, such as absence (low-scale) or presence (high-scale) of violation of flavor universality and zero-distance flavor transition.

In a previous paper, we have proposed a model-independent framework for testing leptonic unitarity assuming low-scale unitarity violation [@Fong:2016yyh]. Our framework is based on the three active plus $N$ sterile lepton (called neutrino) system, which is unitary in the whole $(3+N)$ dimensional state space but restriction to observables in the active neutrino subspace renders the theory non-unitary in that subspace. It is referred to as the “$(3+N)$ space unitary model”. We have shown that the restriction of sterile state masses to $0.1 \,\text{eV}^2 \lsim m^2_{J} \lsim 1 \,\text{MeV}^2$ is sufficient to make the observables sterile-sector model independent. That is, the neutrino oscillation probability can be written in such a way that it is independent of details of the sterile mass spectrum and mixing with active neutrinos. The model-independent nature of the framework will be translated into that of the constraints obtained, thereby making leptonic unitarity test more powerful. See ref. [@Parke:2015goa; @Blennow:2016jkn] for a comprehensive analysis of the currently available neutrino data with the active plus sterile framework, and [@Blennow:2016jkn; @Tang:2017khg] for analyses of the future experiments.

Can one distinguish between high-scale and low-scale unitarity violation? The answer is in principle yes, given the above difference between them, which can only be done by detecting violation of flavor universality or zero-distance flavor transition. In ref. [@Fong:2016yyh] we have pointed out another way of distinguishing them by observing the probability leaking term in the oscillation probability, which signals existence of energetically accessible sterile states. Unfortunately, this important trace of the hidden sterile sector has neither been mentioned in the context of unitarity test, nor included into the foregoing analyses before our analysis for JUNO [@JUNO]-like setting.

In this paper we will give a comprehensive treatment of the $(3+N)$ space unitary model. Our particular emphasis in this paper is twofold:

-   We formulate a novel perturbative framework for the $(3+N)$ space unitary model, which we call “[*small unitarity-violation perturbation theory*]{}”. It allows us to calculate the oscillation probability in the presence of matter effect comparable in size to the vacuum effect.

-   We show that the framework can be used in dual modes: It serves (1) as a suitable framework for leptonic unitarity test in neutrino oscillation experiments, and (2) as a hunting tool for unitarity violation effects, which could serve for another way of distinguishing low-scale unitarity violation from high-scale one.

It must be remarked that the sterile sector model-independent nature of the $(3+N)$ space unitary model is demonstrated in ref. [@Fong:2016yyh] only in vacuum and in matter to first order in matter perturbation theory. It is the first goal of this paper to formulate the appropriate framework to address this question, and show that the model-independence holds after inclusion of sizeable matter effect. In fact, we will observe that the same condition on the sterile neutrino masses guarantees this property. It may be useful in the application of our framework to some of the ongoing and next generation neutrino oscillation experiments [@Abe:2017vif; @Adamson:2016tbq; @Abe:2017aap; @Collaboration:2011ym; @Abe:2015zbg; @Acciarri:2015uup; @TheIceCube-Gen2:2016cap; @Adrian-Martinez:2016zzs]. The discussion of the second point above, our second goal, will follow once the formulation of perturbation theory is completed.

It may be worth to note that there exists difference between low- and high-scale unitarity violation with regard to need for and implication to future neutrino experiments. The scenario of high-scale unitarity violation has been studied extensively in the literature [@Antusch:2006vwa; @FernandezMartinez:2007ms; @Antusch:2009pm; @Antusch:2009gn; @Antusch:2014woa; @Escrihuela:2015wra; @Fernandez-Martinez:2016lgt; @Blennow:2016jkn; @Escrihuela:2016ube].[^4] In this case, due to preserved $SU(2) \times U(1)$ symmetry at high scales, it is conceivable that the constraints from measurements using probes in the charged lepton sector play a dominant role. On the other hand, in the case of low-scale unitarity violation, neutrino oscillation experiments will play key role in constraining unitarity violation.

Essence of the present and the previous papers {#sec:essence}
==============================================

In this section, we present essence of the present and the previous [@Fong:2016yyh] papers in which we try to construct an adequate formulation to describe unitarity violation at low energies, $E \ll m_{W}$.

Unitary 3 active + $N$ sterile neutrino system {#sec:3+N-system}
----------------------------------------------

We have shown in the $(3+N)$ space unitary model that the active neutrino oscillation probability can be written in a sterile sector model-independent way under the constraint on sterile state masses to $0.1\, \text{eV}^2 \lsim m^2_{J} \lsim 1\, \text{MeV}^2$ both in vacuum and to first order in matter effect [@Fong:2016yyh]. The lower limit of the sterile mass range comes from the requirement that fast oscillations due to active-sterile and sterile-sterile mass differences are averaged out due to decoherence. Whereas the upper limit is for sterile neutrinos to take part in reactor neutrino experiments, which may be relaxed to $m^2_{J} \lsim (1 - 10) \,\text{GeV}^2$ for accelerator neutrinos. Then, the obvious question is whether the conclusion remains the same when all order effect of matter is taken into account. We will answer the question in the positive in this paper.

Non-unitary evolution of neutrinos in vacuum {#sec:nonunitarity-vacuum}
--------------------------------------------

Despite large state space of $(3+N)$ dimensions the resultant expression of the oscillation probability has a simple form under the above stated restriction on the sterile neutrino mass spectrum. In vacuum it has the form $$\begin{aligned}
P(\nu_\beta \rightarrow \nu_\alpha) &=& 
\mathcal{C}_{\alpha \beta} + 
\left| \sum_{j=1}^{3} U_{\alpha j} U^{*}_{\beta j} \right|^2 - 
2 \sum_{j \neq k} 
\mbox{Re} 
\left( U_{\alpha j} U_{\beta j}^* U_{\alpha k}^* U_{\beta k} \right) 
\sin^2 \frac{ ( \Delta_{k} - \Delta_{j} ) x  }{ 2 }
\nonumber\\
&-&
\sum_{j \neq k} \mbox{Im} 
\left( U_{\alpha j} U_{\beta j}^* U_{\alpha k}^* U_{\beta k} \right) 
\sin ( \Delta_{k} - \Delta_{j} ) x, 
\label{P-beta-alpha-ave-vac}\end{aligned}$$ where $$\begin{aligned}
\mathcal{C}_{\alpha \beta} \equiv 
\sum_{J=4}^{3+N}
\vert W_{\alpha J} \vert^2 \vert W_{\beta J} \vert^2, 
\label{Cab} \end{aligned}$$ in the appearance ($\alpha \neq \beta$) as well as in the disappearance channels ($\alpha = \beta$).[^5] The $(3 \times N)$ $W$ matrix elements bridge between active and sterile state space, see eq. (\[U-parametrize\]) for definition.

The characteristic features of (\[P-beta-alpha-ave-vac\]) are:

1.  The non-unitary matrix $U$ replaces the standard unitary three-flavor mixing matrix often parametrized with Particle Data Group convention $U_{\text{\tiny PDG}}$ [@Olive:2016xmw].

2.  Probability leakage term $\mathcal{C}_{\alpha \beta}$ appears reflecting the nature of low-energy unitarity violation in which the probability can flow out into the sterile state space from active neutrino space.

3.  Due to non-unitarity, $\delta_{\alpha \beta}$ term in the unitary case is modified to $\left| \sum_{j=1}^{3} U_{\alpha j} U^{*}_{\beta j} \right|^2$.

The points 2 and 3 above are important ones and the clarifying remarks about them are in order:

-   Presence or absence of the probability leakage term $\mathcal{C}_{\alpha \beta}$ distinguishes between low-energy and high-energy unitarity violation [@Fong:2016yyh]. Unfortunately, $\mathcal{C}_{\alpha \beta}$ may be small because it is of fourth order in $W$.

-   Difference in normalization factor, the second term in (\[P-beta-alpha-ave-vac\]), between unitary and non-unitary cases is of order $\sim W^4$ ($\sim W^2$) in the appearance (disappearance) channels.

To understand the latter point, we notice that unitarity in the $(3+N)$ space unitary model can be written as $$\begin{aligned}
\delta_{\alpha \beta} = 
\sum_{j=1}^{3} U_{\alpha j} U^{*}_{\beta j} + \sum_{J=4}^{N+3} W_{\alpha J} W^{*}_{\beta J}. 
\label{unitarity0}\end{aligned}$$ Then, $\left| \sum_{j=1}^{3} U_{\alpha j} U^{*}_{\beta j} \right|^2 = \left| \sum_{J=4}^{N+3} W_{\alpha J} W^{*}_{\beta J} \right|^2$ in the appearance channel ($\alpha \neq \beta$), and $\left( \sum_{j=1}^{3} \vert U_{\alpha j} \vert^2 \right)^2 = \left( 1 -  \sum_{J=4}^{N+3} \vert W_{\alpha J} \vert^2 \right)^2 = 1 - \mathcal{O} (W^2)$ in the disappearance channel ($\alpha = \beta$), which justifies the above statement.

We emphasize, therefore, that the probability leaking term $\mathcal{C}_{\alpha \beta}$ and the another constant term $\left| \sum_{j=1}^{3} U_{\alpha j} U^{*}_{\beta j} \right|^2$ in the oscillation probabilities are the same order, $\mathcal{O} (W^4)$, in the appearance channels. Hence, we do not see any good reasons why the former can be ignored, as was done in the existing literatures. It is also worth to note that $\mathcal{O} (W^2)$ difference in normalization in the disappearance channel would make detection of unitarity violation more feasible. It is one of the reasons for high sensitivity to unitarity violation that could be reached in disappearance measurement in the JUNO-like setting [@Fong:2016yyh].

Non-unitary evolution of neutrinos in matter {#sec:nonunitarity-matter}
--------------------------------------------

We have also examined in ref. [@Fong:2016yyh] the question of whether inclusion of the matter effect alters the above features of non-unitary evolution of neutrino states in vacuum. We have found that as far as we remain in the region of unitarity violating element $\vert W \vert \simeq 0.1$[^6] or larger, the matter effect does not alter the above features of the oscillation probability in (\[P-beta-alpha-ave-vac\]) under the same restriction on sterile masses. Notice that $\vert W \vert \simeq 0.1$ implies that the unitarity violating effect in the probability is of the order of $\vert W \vert^4 \sim 10^{-4}$, except for the $\mathcal{O} (W^2)$ difference in normalization constant in the disappearance probability. It is practically the limit of order of magnitude that can be explored by the next generation neutrino oscillation experiments.

We must note, however, that this conclusion is based on our treatment to first order in matter perturbation theory. Hence, the statement is not a very convincing one. Given the fact that setup of some of the next generation long-baseline (LBL) experiments require consideration of matter effect of comparable size with the vacuum effect, it is clear that a better treatment is necessary to understand the influence of the matter effect in the $(3+N)$ model. In particular, the conditions that allows us to make the active space oscillation probabilities sterile-sector model independent have to be worked out. This task will be carried out in section \[sec:formulation\].

As an outcome of the honest computation using the small unitarity-violation perturbation theory in the $(3+N)$ space unitary model to order $W^4$, we postulate the following theorem:

[**Uniqueness theorem**]{}

-   All the matter dependent perturbative corrections in $W$ in the oscillation probability vanish or can be ignored after averaging over the fast oscillations and using the suppression due to the large energy denominators, leaving only the probability leakage term $\mathcal{C}_{\alpha \beta}$, the first term in (\[P-beta-alpha-ave-vac\]).

It must be remarked here that unitarity violation effects which are hidden in non-unitary active space mixing matrix $U$ produces zeroth- to higher order effects of $W$. The above theorem is only about the terms generated by explicit perturbative corrections in $W$.[^7]

As a result of the explicit computation, we show that the oscillation probability in matter between active flavor neutrinos in the $(3+N)$ space unitary model to fourth order in $W$ in our small unitarity-violation perturbation theory can be written as $$\begin{aligned}
P(\nu_\beta \rightarrow \nu_\alpha) &=& 
\mathcal{C}_{\alpha \beta} + 
\left| \sum_{j=1}^{3} U_{\alpha j} U^{*}_{\beta j} \right|^2 
\nonumber\\
&-&
2 \sum_{j \neq k} 
\mbox{Re} 
\left[ (UX)_{\alpha j} (UX)_{\beta j}^* (UX)_{\alpha k}^* (UX)_{\beta k} \right] 
\sin^2 \frac{ ( h_{k} - h_{j} ) x  }{ 2 }
\nonumber\\
&-&
\sum_{j \neq k} \mbox{Im} 
\left[ (UX)_{\alpha j} (UX)_{\beta j}^* (UX)_{\alpha k}^* (UX)_{\beta k} \right] 
\sin ( h_{k} - h_{j} ) x, 
\label{P-beta-alpha-final}\end{aligned}$$ where $h_{i}$ $(i=1,2,3)$ denote the energy eigenvalues of zeroth-order states of active neutrinos in matter, and $X$ is the unitary matrix which diagonalize the zeroth-order Hamiltonian used to formulate our perturbation theory. $\mathcal{C}_{\alpha \beta}$ is given in (\[Cab\]). The expression is valid under the restriction on sterile neutrino masses $0.1\, \text{eV}^2 \lsim m^2_{J} \lsim 1\, \text{MeV}^2$ for $\vert W \vert^4 \gsim 10^{-4}$. For the restriction needed on the sterile state masses for smaller $W$, see section \[sec:probability-2nd\].

Notice that the normalization term, the second term both in (\[P-beta-alpha-ave-vac\]) and (\[P-beta-alpha-final\]), comes out from the contribution of zeroth-order Hamiltonian which contains all orders effect of the matter potential. It is true despite the fact that it is a vacuum effect, which indeed exists in vacuum.[^8]

The expression (\[P-beta-alpha-final\]) is a very transparent result in the sense that (1) the vacuum non-unitary mixing matrix $U$ is “dressed” in a simple way by the matter effect represented by $X$, the unitary matrix which diagonalizes the zeroth-order Hamiltonian, and (2) the probability leaking term and the normalization term stay as they are in vacuum.

Enhanced higher-order $W$ corrections and nature of unitarity violation  {#sec:higher-order-corrections}
------------------------------------------------------------------------

We will point out that, despite our above theorem, there are regions of neutrino energy and baseline that condition for suppression due to the large energy denominators is not fully effective. In this region outside validity of the theorem we show that second order correction terms in $W$, together with the leaking term $\mathcal{C}_{\alpha \beta}$, may not be totally negligible, and it could be detectable. If it were the case, it could distinguish low-scale unitarity violation from high-scale one.

 Small unitarity-violation perturbation theory of neutrino oscillation in matter  {#sec:formulation}
=================================================================================

We formulate a perturbation theory of the $(3+N)$ state unitary model using an expansion parameter of matrix elements of $W$ signifying unitarity violation effects, assuming it small. In the main text we mostly confine ourselves to the formulas to second order in $W$, but include fourth order terms whenever it is necessary. Our formulation of the perturbative framework in this section will be done aiming at constructing a model-independent framework for leptonic unitarity test. Usage of the same probability formula as a hunting tool of unitarity violation and discriminator between low-scale and high-scale unitarity violation will be discussed in section \[sec:UV-low-high-E\].

For simplicity, we take the uniform number density approximation for electrons and neutrons in matter. However, extension to the varying density case is, in principle, straightforward as far as adiabaticity holds.

3 active plus $N$ sterile neutrino system in the flavor basis {#sec:flavor-basis}
-------------------------------------------------------------

The $S$ matrix describes possible flavor changes after traversing a distance $x$ $$\begin{aligned}
\nu_{\alpha} (x) = S_{\alpha \beta} \nu_{\beta} (0), 
\label{def-S}\end{aligned}$$ and the oscillation probability is given by $$\begin{aligned}
P(\nu_{\beta} \rightarrow \nu_{\alpha}; x)= 
\vert S_{\alpha \beta} \vert^2.  
\label{def-P}\end{aligned}$$ The neutrino evolution in flavor basis in the $(3+N)$ space unitary model is governed by the Schrödinger equation $$\begin{aligned}
i \frac{d}{dx} \nu = H \nu. 
\label{evolution}\end{aligned}$$ Given the flavor basis Hamiltonian $H$, the $S$ matrix is given by $$\begin{aligned}
S = T \text{exp} \left[ -i \int^{x}_{0} dx^{\prime} H(x^{\prime}) \right], 
\label{S-matrix-def}\end{aligned}$$ where $T$ symbol indicates the “time ordering” (in fact “space ordering” here). The right-hand side of (\[evolution\]) may be written as $e^{-i H x}$ for the case of constant matter density.

The flavor basis Hamiltonian $H$ is $(3+N) \times (3+N)$ matrix: $$\begin{aligned}
H = {\bf U}
\left[
\begin{array}{cccccc}
\Delta_{1} & 0 & 0 & 0 & 0 & 0 \\
0 & \Delta_{2} & 0 & 0 & 0 & 0 \\
0 & 0 & \Delta_{3} & 0 & 0 & 0 \\
0 & 0 & 0 & \Delta_{4} & 0 & 0 \\
0 & 0 & 0 & 0 & \cdot \cdot \cdot & 0 \\
0 & 0 & 0 & 0 & 0 & \Delta_{3+N} \\
\end{array}
\right] 
{\bf U}^{\dagger} 
+
\left[
\begin{array}{cccccc}
\Delta_{A} - \Delta_{B} & 0 & 0 & 0 & 0 & 0 \\
0 & - \Delta_{B} & 0 & 0 & 0 & 0 \\
0 & 0 & - \Delta_{B} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & \cdot \cdot \cdot & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}
\right],
\label{hamiltonian}\end{aligned}$$ where $$\begin{aligned}
\Delta_{i}  \equiv \frac{ m^2_{i} }{2E} 
\hspace{4mm}
(i = 1,2,3),
\hspace{8mm}
\Delta_{J}  \equiv \frac{ m^2_{J} }{2E}
\hspace{4mm}
(J = 4, \cdot \cdot \cdot, 3+N).
\label{Delta-def}\end{aligned}$$ Here, $m_{i}$ ($m_{J}$) denote the mass of mostly active (sterile) neutrinos and $E$ is the neutrino energy. $\Delta_{A}$ and $\Delta_{B}$ are related to Wolfenstein’s matter potential [@Wolfenstein:1977ue] due to CC and NC reactions, $a$ and $b$, as $$\begin{aligned}
\Delta_{A}  \equiv \frac{ a }{2E}, 
\hspace{10mm}
\Delta_{B}  \equiv \frac{ b }{2E}, 
\label{Delta-ab-def}\end{aligned}$$ where $$\begin{aligned}
a &=&  
2 \sqrt{2} G_F N_e E \approx 1.52 \times 10^{-4} \left( \frac{Y_e \rho}{\rm g\,cm^{-3}} \right) \left( \frac{E}{\rm GeV} \right) {\rm eV}^2, 
\nonumber \\
b &=& \sqrt{2} G_F N_n E = \frac{1}{2} \left( \frac{N_n}{N_e} \right) a. 
\label{matt-potential}\end{aligned}$$ Here, $G_F$ is the Fermi constant, $N_e$ and $N_n$ are the electron and neutron number densities in matter. $\rho$ and $Y_e$ denote, respectively, the matter density and number of electron per nucleon in matter. In (\[hamiltonian\]), ${\bf U}$ denotes the flavor mixing matrix which relates $(3+N)$ dimensional flavor neutrino states to the vacuum mass eigenstate basis as $\nu_{\zeta} = {\bf U}_{\zeta z} \tilde{\nu}_{z}$, where $\zeta$ runs over active flavor $\alpha = e,\mu,\tau$ and sterile flavor $s = s_1,\cdot \cdot \cdot,s_N$ indices, $z$ runs over mostly active $i=1,2,3$ and mostly sterile mass eigenstate indices $J=4,5,\cdot \cdot \cdot, N+3$.

For simplicity, we introduce a compact notation which writes $(3+N) \times (3+N)$ matrix in a form of $2 \times 2$ matrix. By defining the active $3 \times 3$ matter potential matrix $$\begin{aligned}
A =
\left[
\begin{array}{ccc}
\Delta_{A} - \Delta_{B} & 0 & 0 \\
0 & - \Delta_{B} & 0 \\
0 & 0 & - \Delta_{B} \\
\end{array}
\right] 
\label{matter-pot2}\end{aligned}$$ the flavor basis Hamiltonian is written as $$\begin{aligned}
H = 
{\bf U} 
\left[
\begin{array}{cc}
{\bf \Delta_{a} } & 0 \\
0 & {\bf \Delta_{s} } \\
\end{array}
\right] 
{\bf U}^{\dagger} 
+
\left[
\begin{array}{cc}
A & 0 \\
0 & 0 \\
\end{array}
\right] 
\equiv H_{ \text{vac} } + H_{ \text{matt} } 
\label{flavor-hamiltonian}\end{aligned}$$ where ${\bf \Delta_{a} } = \text{diag} ( \Delta_{1}, \Delta_{2}, \Delta_{3})$ and ${\bf \Delta_{s} } = \text{diag} (\Delta_{4}, \Delta_{5}, \cdot \cdot \cdot, \Delta_{N+3})$.

As an application of our framework, we anticipate leptonic unitarity test in the LBL accelerator neutrino experiments which utilize atmospheric-scale neutrino oscillations. Therefore, we assume that the system satisfies the following conditions in formulating our perturbation theory $$\begin{aligned}
\frac{ \Delta m^2_{31} L }{2E} \sim \frac{ \Delta m^2_{32} L }{2E} \sim \mathcal{O} (1), 
\hspace{10mm}
\frac{ a L }{2E} \sim \frac{ b L }{2E} \sim \mathcal{O} (1), 
\label{assumption}\end{aligned}$$ where $L$ denotes baseline, $\Delta m^2_{ji} \equiv m_j^2 - m_i^2$, and $$\begin{aligned}
\frac{ a L }{2E} 
&=& \sqrt{2} G_F N_e L
= 0.58  
\left(\frac{\rho}{3 \text{g/cm}^3}\right)
\left(\frac{L}{1000 \mbox{km}}\right). 
\label{aL-2E}\end{aligned}$$ They probably ensure that our oscillation probability formulas have applicability to the terrestrial LBL and atmospheric neutrino experiments with baseline up to $\sim 10^4$ km and energies from low to high, up to $E \sim 100$ GeV. More precise discussions on where our formulas are valid will be given in sections \[sec:energy-denominator\] and \[sec:higher-order\].

Vacuum mass eigenstate basis, or tilde basis {#sec:tilde-basis}
--------------------------------------------

To formulate perturbative treatment it is convenient to consider the vacuum mass eigenstate basis, the tilde basis, introduced in the previous section $$\begin{aligned}
\tilde{\nu}_{z} = ({\bf U}^{\dagger})_{z \zeta} \nu_{\zeta}.
\label{tilde-basis}\end{aligned}$$ The tilde basis Hamiltonian is related to the flavor basis one as $$\begin{aligned}
\tilde{H} = {\bf U}^{\dagger} H {\bf U}. 
\label{tilde-hamiltonian}\end{aligned}$$ The explicit form of $\tilde{H}$ is given by $$\begin{aligned}
\tilde{H} &=& \tilde{H}_{ \text{vac} } + \tilde{H}_{ \text{matt} } = 
\left[
\begin{array}{cc}
{\bf \Delta_{a} } & 0 \\
0 & {\bf \Delta_{s} } \\
\end{array}
\right] 
+ 
{\bf U}^{\dagger} 
\left[
\begin{array}{cc}
A & 0 \\
0 & 0 \\
\end{array}
\right] 
{\bf U}. 
\label{H-tilde-3+N}\end{aligned}$$ We parametrize the $(3+N) \times (3+N)$ dimensional flavor mixing matrix ${\bf U}$ as $$\begin{aligned}
{\bf U} = \left[
\begin{array}{cc}
U & W \\
Z & V \\
\end{array}
\right]. 
\label{U-parametrize}\end{aligned}$$ The matrix $U$ and $V$ are $3 \times 3$ and $N \times N$ matrices, respectively, and $W$ and $Z$ have sizes that just fill in the space. In our $(3+N)$ model, unitarity is obeyed in the whole $(3+N)$ state space: $$\begin{aligned}
{\bf U} {\bf U}^{\dagger} &=& 
\left[
\begin{array}{cc}
U U^{\dagger} + W W^{\dagger} & U Z^{\dagger} + W V^{\dagger} \\ 
Z U^{\dagger} + V W^{\dagger}  & Z Z^{\dagger} + V V^{\dagger} \\
\end{array}
\right] = 
\left[
\begin{array}{cc}
{\bf 1}_{3 \times 3} & 0 \\
0 & {\bf 1}_{N \times N} \\
\end{array}
\right], 
\nonumber \\
{\bf U}^{\dagger} {\bf U} &=&  
\left[
\begin{array}{cc}
U^{\dagger} U + Z^{\dagger} Z & U^{\dagger} W + Z^{\dagger} V  \\
W^{\dagger} U + V^{\dagger} Z & W^{\dagger} W + V^{\dagger} V \\
\end{array}
\right] =
\left[
\begin{array}{cc}
{\bf 1}_{3 \times 3} & 0 \\
0 & {\bf 1}_{N \times N} \\
\end{array}
\right].
\label{unitarity}\end{aligned}$$

Then, the Hamiltonian $\tilde{H}$ in vacuum mass eigenstate basis is given by $$\begin{aligned}
\tilde{H} &=& 
\left[
\begin{array}{cc}
{\bf \Delta_{a} } & 0 \\
0 & {\bf \Delta_{s} } \\
\end{array}
\right] 
+
\left[
\begin{array}{cc}
U^{\dagger} A U & U^{\dagger} A W \\
W^{\dagger} A U & W^{\dagger} A W \\
\end{array}
\right]. 
\label{tilde-H}\end{aligned}$$ As in vacuum, the neutrino oscillation is governed only by the $U$ and $W$ matrices, and is independent of $Z$ and $V$ matrices. It is natural that $V$ matrix does not show up in physical Hamiltonian matrix because the rotations inside sterile basis does not have any physical meaning, if we observe the system only by the Standard Model interactions.

Formulating small unitarity-violation perturbation theory {#sec:UV-perturb}
---------------------------------------------------------

We now construct the small unitarity-violation perturbation theory. It is natural to consider the framework in which the tilde-basis Hamiltonian $\tilde{H}$ is decomposed into the un-perturbed and perturbed parts, $\tilde{H}_{0} + \tilde{H}_{1}$, as follows: $$\begin{aligned}
\tilde{H}_{0} = 
\left[
\begin{array}{cc}
{\bf \Delta_{a} } + U^{\dagger} A U & 0 \\
0 & {\bf \Delta_{s} } \\
\end{array}
\right], 
\hspace{10mm}
\tilde{H}_{1} = 
\left[
\begin{array}{cc}
0 & U^{\dagger} A W \\
W^{\dagger} A U & W^{\dagger} A W \\
\end{array}
\right]. 
\label{tilde-H0+H1} \end{aligned}$$ Therefore, what we mean by “expansion by unitarity violation effect” is an expansion by the $W$ matrix elements.[^9] We assume, for simplicity, that all the $W$ matrix elements are small and have the same order $\epsilon_{s}$. Then, $3 \times N$ ($N \times 3$) sub-matrix elements in $\tilde{H}_{ \text{matt} }$ are of order $\epsilon_{s}$, while the pure sterile space $N \times N$ sub-matrix elements are of order $\epsilon_{s}^2$. For simplicity, we often use the expression “expanding to order $W^n$” which means to order $\epsilon_{s}^n$ in this paper.

### Hat basis

To formulate perturbation theory with $\tilde{H}_{0}$ and $\tilde{H}_{1}$ given above we transform to a basis in which the un-perturbed part of the Hamiltonian is diagonal, which we call the “hat basis”. Since the $3 \times 3$ sub-matrix ${\bf \Delta_{a} } + U^{\dagger} A U$ in $\tilde{H}_{0}$ is hermitian, it can be diagonalized by the unitary transformation $$\begin{aligned}
X^{\dagger} \left(  {\bf \Delta_{a} } + U^{\dagger} A U \right) X = 
\left[
\begin{array}{ccc}
h_{1} & 0 & 0 \\
0 & h_{2} & 0 \\
0 & 0 & h_{3} \\
\end{array}
\right] \equiv {\bf h} 
\label{H0-diag}\end{aligned}$$ with $X$ being the $3 \times 3$ unitary matrix. Then, $\tilde{H}_{0}$ can be diagonalized by using $$\begin{aligned}
{\bf X} \equiv 
\left[
\begin{array}{cc}
X & 0 \\
0 & 1 \\
\end{array}
\right] 
\label{bfX-def}\end{aligned}$$ as $$\begin{aligned}
&& 
{\bf X}^{\dagger} \tilde{H}_{0} {\bf X}
= \left[
\begin{array}{cc}
X^{\dagger} \left(  {\bf \Delta_{a} } + U^{\dagger} A U \right) X & 0 \\
0 & {\bf \Delta_{s} } \\
\end{array}
\right] 
= \left[
\begin{array}{cc}
{\bf h} & 0 \\
0 & {\bf \Delta_{s} } \\
\end{array}
\right] 
\equiv \hat{H}_{0}, 
\label{hat-H0}\end{aligned}$$ the zeroth-order Hamiltonian in the hat basis. Since $\hat{H}_{0}$ is diagonal it is easy to compute $e^{ \pm i \hat{H}_{0} x}$: $$\begin{aligned}
e^{ \pm i \hat{H}_{0} x} = 
\left[
\begin{array}{cc}
e^{ \pm i {\bf h} x } & 0 \\
0 & e^{ \pm i {\bf \Delta_{s} x }   } \\
\end{array}
\right]. 
\label{exp-H0}\end{aligned}$$ Then, the perturbed Hamiltonian is given by $$\begin{aligned}
\hat{H}_{1} &=& {\bf X}^{\dagger} \tilde{H}_{1} {\bf X} 
=
\left[
\begin{array}{cc}
0 & (UX)^{\dagger} A W \\
W^{\dagger} A (UX) & W^{\dagger} A W \\
\end{array}
\right]. 
\label{hat-H1}\end{aligned}$$ The eigenvalues of $\tilde{H}_{0}$ is therefore $h_{1}$, $h_{2}$, $h_{3}$, and $\Delta_{J}$ ($J=4, \cdot \cdot \cdot, 3+N$). Therefore, the sterile neutrino masses are affected neither by the active states nor the matter potential in our zeroth-order unperturbed basis. It must be a good approximation because we have assumed that the sterile masses are much heavier than the active ones, and we are interested in the energy region $a \sim \Delta m^2_{31}$.

To do real calculations of the $S$ matrix elements we must solve the zeroth order Hamiltonian $\tilde{H}_{0}$. This task will be carried out in section \[sec:exact-solution-zeroth\], in which we derive explicit expressions of the eigenvalues $h_{i}$ and the unitary matrix $X$.

Now, we formulate perturbation theory with the hat basis Hamiltonian, $\hat{H}_{0}$ in (\[hat-H0\]) and $\hat{H}_{1}$ in (\[hat-H1\]) after a clarifying note in the next subsection.

### The relationship between quantities in various bases

So far we have introduced the tilde- and the hat-basis: $$\begin{aligned}
\tilde{H} = {\bf U}^{\dagger} H {\bf U}, 
\hspace{10mm}
\hat{H} = {\bf X}^{\dagger} \tilde{H}  {\bf X}, 
\label{tilde-hat-relation}\end{aligned}$$ where ${\bf X}$ is given by eq. (\[bfX-def\]). Therefore, $$\begin{aligned}
\hat{H} 
= \left( {\bf U}  {\bf X} \right)^{\dagger} H \left( {\bf U}  {\bf X} \right).
\label{flavor-hat-relation}\end{aligned}$$ Or $$\begin{aligned}
H = \left( {\bf U}  {\bf X} \right) \hat{H} \left( {\bf U}  {\bf X} \right)^{\dagger}, 
\hspace{10mm}
S = \left( {\bf U}  {\bf X} \right) \hat{S} \left( {\bf U}  {\bf X} \right)^{\dagger}.
\label{flavor-hat-relation2}\end{aligned}$$ Notice that both ${\bf U}$ and ${\bf X}$ are unitary, and hence ${\bf U} {\bf X}$ is unitary too. The relationship between wave functions of various basis are given by $$\begin{aligned}
\hat{\nu}_{y} &=& 
{\bf X}^{\dagger}_{y z} \tilde{\nu}_{z} 
= \left( {\bf U} {\bf X} \right)^{\dagger}_{y \zeta} \nu_{\zeta},
\nonumber \\
\nu_{\zeta} &=& 
\left( {\bf U} {\bf X} \right)_{\zeta y} 
\hat{\nu}_{y}. 
\label{hat-basis}\end{aligned}$$ where $y$ denote the hat-basis indices. Using the explicit parametrization of the ${\bf U}$ matrix we have $$\begin{aligned}
{\bf U}  {\bf X} &=& 
\left[
\begin{array}{cc}
U & W \\
Z & V \\
\end{array}
\right] 
\left[
\begin{array}{cc}
X & 0 \\
0 & 1 \\
\end{array}
\right] 
= 
\left[
\begin{array}{cc}
U X & W \\
Z X & V \\
\end{array}
\right], 
\nonumber \\
\left( {\bf U}  {\bf X} \right)^{\dagger} &=& 
\left[
\begin{array}{cc}
(UX)^{\dagger} & \left( Z X \right)^{\dagger} \\
W^{\dagger} & V^{\dagger} \\
\end{array}
\right].
\label{UX-parametrize}\end{aligned}$$ It may be helpful for our discussions later to understand the relationship between $S$ and $\hat{S}$ matrix elements. For this purpose, we denote them in the block form $$\begin{aligned}
S = 
\left[
\begin{array}{cc}
S_{aa} & S_{aS}  \\
S_{Sa} & S_{SS} \\
\end{array}
\right], 
\hspace{10mm}
\hat{S} = 
\left[
\begin{array}{cc}
\hat{S}_{aa} & \hat{S}_{aS}  \\
\hat{S}_{Sa} & \hat{S}_{SS} \\
\end{array}
\right],
\label{S-hatS}\end{aligned}$$ where the subscripts $a$ and $S$ indicate that they act (for the right index) to the active or the sterile subspaces. Notice that $S_{aa}$ and $S_{aS}$, for example, are $3 \times 3$ and $3 \times N$ matrices, respectively. Then, the relationship between $S$ and $\hat{S}$ matrix elements can be written explicitly as $$\begin{aligned}
S_{aa} &=& 
(U X) \hat{S}_{aa} (UX)^{\dagger} + (U X) \hat{S}_{aS} W^{\dagger} + 
W \hat{S}_{Sa} (UX)^{\dagger} + W \hat{S}_{SS} W^{\dagger},
\nonumber \\
S_{aS} &=& 
(U X) \hat{S}_{aa} \left( Z X \right)^{\dagger} + (U X) \hat{S}_{aS} V^{\dagger} + 
W \hat{S}_{Sa} \left( Z X \right)^{\dagger} + W \hat{S}_{SS} V^{\dagger},
\nonumber \\
S_{Sa} &=& 
\left( Z X \right) \hat{S}_{aa} (U X)^{\dagger} + \left( Z X \right) \hat{S}_{aS} W^{\dagger}  + 
V \hat{S}_{Sa} (U X)^{\dagger} + V \hat{S}_{SS} W^{\dagger},
\nonumber \\
S_{SS} &=& 
\left( Z X \right) \hat{S}_{aa} \left( Z X \right)^{\dagger} + \left( Z X \right) \hat{S}_{aS} V^{\dagger}  + 
V \hat{S}_{Sa} \left( Z X \right)^{\dagger} + V \hat{S}_{SS} V^{\dagger}. 
\label{S-Shat-relation}\end{aligned}$$

### Computation of $\hat{S}$ matrix elements {#sec:hatS-matrix}

To calculate $\hat {S} (x) = \exp \left[ -i \int^{x}_{0} dx \hat{H} (x)  \right] $ we define $\Omega(x)$ as $$\begin{aligned}
\Omega(x) = e^{i \hat{H}_{0} x} \hat{S} (x).
\label{def-omega}\end{aligned}$$ $\Omega(x)$ obeys the evolution equation $$\begin{aligned}
i \frac{d}{dx} \Omega(x) = H_{1} \Omega(x),
\label{omega-evolution}\end{aligned}$$ where $$\begin{aligned}
H_{1} \equiv e^{i \hat{H}_{0} x} \hat{H}_{1} e^{-i \hat{H}_{0} x}. 
\label{def-H1}\end{aligned}$$ Then, $\Omega(x)$ can be computed perturbatively as $$\begin{aligned}
\Omega(x) &=& 1 + 
(-i) \int^{x}_{0} dx' H_{1} (x') + 
(-i)^2 \int^{x}_{0} dx' H_{1} (x') \int^{x'}_{0} dx'' H_{1} (x'') 
\nonumber \\
&+& 
(-i)^3 \int^{x}_{0} dx' H_{1} (x') \int^{x'}_{0} dx'' H_{1} (x'') \int^{x''}_{0} dx''' H_{1} (x''') 
\nonumber \\
&+& 
(-i)^4 \int^{x}_{0} dx' H_{1} (x') \int^{x'}_{0} dx'' H_{1} (x'') \int^{x''}_{0} dx''' H_{1} (x''') \int^{x'''}_{0} dx'''' H_{1} (x'''') +  \cdot \cdot \cdot, 
\nonumber \\
\label{Omega-expand}\end{aligned}$$ where the “space-ordered” form in (\[Omega-expand\]) is essential because of the highly nontrivial spatial dependence in $H_{1}$. Upon obtaining $\Omega(x)$, $\hat{S}$ matrix can be obtained as $$\begin{aligned}
\hat{S} (x) = e^{- i \hat{H}_{0} x} \Omega(x). 
\label{hatS-Omega}\end{aligned}$$ By knowing $\hat{S}$ matrix elements, the $S$ matrix is obtained by using (\[flavor-hat-relation2\]).

The perturbing Hamiltonian $H_{1}$ defined in (\[def-H1\]) has a structure $$\begin{aligned}
H_{1} = \left[
\begin{array}{cc}
0 & e^{ i {\bf h} x} (UX)^{\dagger} A W e^{ - i {\bf \Delta_{s} } x}  \\
e^{ i {\bf \Delta_{s} } x} W^{\dagger} A (UX) e^{ - i {\bf h} x}  & e^{ i {\bf \Delta_{s} } x} W^{\dagger} A W e^{ - i {\bf \Delta_{s} } x} \\
\end{array}
\right]. 
\label{H1-matrix}\end{aligned}$$ That is, $(H_{1})_{i j} = 0$ in the whole active neutrino subspace. The non-vanishing elements of $H_{1}$ are as follows: $$\begin{aligned}
(H_{1})_{i J} &=& 
e^{- i ( \Delta_{J} - h_{i} ) x} \left\{ (UX)^{\dagger} A W \right\}_{i J}, 
\nonumber \\
(H_{1})_{J i} &=& 
e^{ - i ( h_{i} - \Delta_{J} ) x} \left\{ W ^{\dagger} A (UX) \right\}_{J i}, 
\nonumber \\
(H_{1})_{J K} &=& 
e^{- i ( \Delta_{K} - \Delta_{J} ) x} \left\{ W ^{\dagger} A W \right\}_{J K}. 
\label{H1-elements}\end{aligned}$$ Inserting eq. (\[H1-elements\]) into (\[Omega-expand\]), we can compute all the $\Omega$ matrix elements. The simplest ones in first order in $H_{1}$, the second term in (\[Omega-expand\]), are given by $$\begin{aligned}
\Omega_{i j} [1] &=& 0, 
\nonumber \\
\Omega_{i J} [1] &=& 
\frac{e^{- i ( \Delta_{J} - h_{i} ) x} - 1 }{ ( \Delta_{J} - h_{i} )  }
\left\{ (UX)^{\dagger} A W \right\}_{i J}, 
\nonumber \\
\Omega_{J i} [1] &=& 
- \frac{e^{ i ( \Delta_{J} - h_{i} ) x} - 1 }{ ( \Delta_{J} - h_{i} ) }
\left\{ W ^{\dagger} A (UX) \right\}_{J i}, 
\nonumber \\
\Omega_{J K} \vert_{J \neq K} [1] &=& 
\frac{e^{- i ( \Delta_{K} - \Delta_{J} ) x} - 1 }{ ( \Delta_{K} - \Delta_{J} )  }
\left\{ W ^{\dagger} A W \right\}_{J K}, 
\nonumber \\
\Omega_{J J} [1] &=& 
(-i x) \left\{ W ^{\dagger} A W \right\}_{J J}, 
\label{Omega-1st-order}\end{aligned}$$ which serve as a building block of the perturbation series because of the structure in (\[Omega-expand\]). The notation “\[1\]” implies that the terms come from first order perturbation with $H_{1}$. For more about notations, see appendix \[sec:hatS-elements\].

We need to compute up to fourth order in $H_{1}$ because we want to keep all the order $W^4$ terms. The requirement arises because the probability leaking term, whose observation is crucial to distinguish between low-energy and high-energy unitarity violation, is of order $W^4$. The other normalization term, the second term in (\[P-beta-alpha-ave-vac\]), also deviates from the one in unitary case by a quantity of order $W^4$ in the appearance channels, but in an implicit way. The resulting expressions of $\hat{S}$ matrix elements to order $W^4$ are summarized in appendix \[sec:hatS-elements\].

There exists important consistency check in the calculation. That is, the identity relation between $\hat{S}$ matrix elements that follows from generalized T invariance:[^10] $$\begin{aligned}
\hat{S}_{i j} (U, W, X, A) &=& \hat{S}_{j i} (U^*, W^*, X^*, A^*),
\nonumber \\
\hat{S}_{i J} (U, W, X, A) &=& \hat{S}_{J i} (U^*, W^*, X^*, A^*),
\nonumber \\\hat{S}_{I J} (U, W, X, A) &=& \hat{S}_{J I} (U^*, W^*, X^*, A^*), 
\label{T-invariance}\end{aligned}$$ where $\hat{S}_{J i}$ is obtained by performing the exchange $h_{i} \leftrightarrow \Delta_{J}$ in $\hat{S}_{i J}$. The generalized T invariance relation will be explicitly verified by the computed results of $\hat{S}$ matrix elements in appendix \[sec:hatS-elements\].[^11]

To carry out a complete consistency check, we have obtained all the $\hat{S}$ matrix elements (including $\hat{S}_{I J}$) to fourth order in $W$, and verified their generalized T invariance. However, only the ones which are required to obtain $S$ matrix elements to fourth order are exhibited in appendix \[sec:hatS-elements\].

### Computation of $S$ matrix elements {#sec:S-matrix}

Given the results of $\hat{S}$ matrix elements it is straightforward to calculate $S$ matrix elements by using the formulas in eq. (\[flavor-hat-relation2\]). The active neutrino space $S$ matrix elements can be written in perturbative forms, $S_{\alpha \beta} = S_{\alpha \beta}^{(0)} + S_{\alpha \beta}^{(2)} + S_{\alpha \beta}^{(4)}$, where $$\begin{aligned}
S_{\alpha \beta}^{(0)} &=& 
\sum_{k l} (UX)_{\alpha k} (UX)^*_{\beta l} 
\hat{S}_{kl}^{(0)}, 
\nonumber \\
S_{\alpha \beta}^{(2)} &=&
\sum_{k l} (UX)_{\alpha k} (UX)^*_{\beta l} 
\hat{S}_{kl}^{(2)} 
+
\sum_{k L} (UX)_{\alpha k} W^*_{\beta L} 
\hat{S}_{kL}^{(1)} 
\nonumber \\
&+&
\sum_{K l} W_{\alpha K} (UX)^*_{\beta l} 
\hat{S}_{K l}^{(1)}
+
\sum_{K L} W_{\alpha K} W^*_{\beta L}
\hat{S}_{KL}^{(0)}, 
\nonumber \\
S_{\alpha \beta}^{(4)} &=&
\sum_{k l} (UX)_{\alpha k} (UX)^*_{\beta l} 
\hat{S}_{kl}^{(4)} 
+
\sum_{k L} (UX)_{\alpha k} W^*_{\beta L} 
\hat{S}_{kL}^{(3)} 
\nonumber \\
&+&
\sum_{K l} W_{\alpha K} (UX)^*_{\beta l} 
\hat{S}_{K l}^{(3)}
+
\sum_{K L} W_{\alpha K} W^*_{\beta L}
\hat{S}_{KL}^{(2)}.
\label{Sab-hatSab}\end{aligned}$$ Using (\[Sab-hatSab\]) the explicit expressions of $S$ matrix elements can be easily obtained with use of $\hat{S}$ matrix elements given in appendix \[sec:hatS-elements\]. For example, $S_{\alpha \beta}$ in zeroth and second orders in $W$ are given, respectively, by $$\begin{aligned}
S_{\alpha \beta}^{(0)} &=& 
\sum_{k} (UX)_{\alpha k} (UX)^*_{\beta k} 
e^{-i h_{k} x}, 
\label{S-alpha-beta-0th}\end{aligned}$$ and $$\begin{aligned}
&& S_{\alpha \beta}^{(2)} =  
\sum_{k, K} 
\frac{ 1 }{ \Delta_{K} - h_{k} } 
\left[
(ix) e^{- i h_{k} x} + \frac{e^{- i \Delta_{K} x} - e^{- i h_{k} x} }{ ( \Delta_{K} - h_{k} )  } 
\right]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\
&-& 
\sum_{k \neq l} \sum_{K}  
\frac{ 1 }{ ( h_{l}  - h_{k} ) (\Delta_{K} - h_{k}) (\Delta_{K} - h_{l}) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{k} \right) e^{- i h_{l} x} 
- \left( \Delta_{K} - h_{l} \right) e^{- i h_{k} x} 
- ( h_{l}  - h_{k} )  e^{- i \Delta_{K} x} 
\biggr]
\nonumber \\
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\nonumber \\
&+& 
\sum_{k, K} 
\frac{e^{- i \Delta_{K} x} - e^{- i h_{k} x} }{ ( \Delta_{K} - h_{k} ) } 
\biggl[ 
(UX)_{\alpha k} W^*_{\beta K} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
+
W_{\alpha K} (UX)^*_{\beta k} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\biggr]
\nonumber \\
&+&
\sum_{K} 
e^{- i \Delta_{K} x} 
W_{\alpha K} W^*_{\beta K}. 
\label{S-alpha-beta-2nd}\end{aligned}$$

The oscillation probability to second order in $W$ {#sec:probability-2nd}
--------------------------------------------------

In this section, we discuss the oscillation probability to second order in $W$. It is to illuminate the principle of calculation, how averaging over the fast oscillation works, and to show which constraints are obtained on the sterile state masses by the requirement of suppression by the energy denominator to make these sterile-sector model dependent terms negligible.

Of course, we will calculate in this paper all the oscillation probabilities $P(\nu_\beta \rightarrow \nu_\alpha)$ in matter to fourth order in $W$ to keep the necessary term, the probability leaking term $\mathcal{C}_{\alpha \beta}$, as mentioned earlier. The key features of the fourth-order terms will be described in the next section \[sec:probability-4th\].

The following formulas include the cases of both disappearance ($\alpha=\beta$) and appearance ($\alpha \neq \beta$) channels. The oscillation probability $P(\nu_\beta \rightarrow \nu_\alpha)$ is given to second order in $W$ as $$\begin{aligned}
&&P(\nu_\beta \rightarrow \nu_\alpha)^{(0+2)} = 
\left| S^{(0)}_{\alpha \beta} \right|^2 
+ 2 \mbox{Re} \left[ \left( S^{(0)}_{\alpha \beta} \right)^{*} S^{(2)}_{\alpha \beta} \right] 
\nonumber \\
&=&
\sum_{k} (UX)_{\alpha k} (UX)^*_{\beta k} 
(UX)^*_{\alpha k} (UX)_{\beta k} 
+ \sum_{k \neq l} (UX)_{\alpha k} (UX)^*_{\beta k} 
(UX)^*_{\alpha l} (UX)_{\beta l} 
e^{-i ( h_{k} - h_{l} ) x}
\nonumber \\ 
&+& 
2 \mbox{Re} 
\biggl\{
\sum_{m}
\sum_{k, K} 
\frac{ 1 }{ \Delta_{K} - h_{k} } 
\left[
(ix) e^{- i ( h_{k} - h_{m} ) x} + \frac{e^{- i ( \Delta_{K} - h_{m} ) x} - e^{- i ( h_{k} - h_{m} ) x} }{ ( \Delta_{K} - h_{k} )  } 
\right]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta k} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\
&-& 
\sum_{m}
\sum_{k \neq l} \sum_{K} 
\frac{ 1 }{ ( h_{l}  - h_{k} ) (\Delta_{K} - h_{k}) (\Delta_{K} - h_{l}) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{k} \right) e^{- i ( h_{l} - h_{m} ) x} 
- \left( \Delta_{K} - h_{l} \right) e^{- i ( h_{k} - h_{m} ) x} 
- ( h_{l}  - h_{k} )  e^{- i ( \Delta_{K} - h_{m} ) x} 
\biggr]
\nonumber \\
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\nonumber \\
&+& 
\sum_{m}
\sum_{k, K} 
\frac{e^{- i (\Delta_{K} - h_{m} ) x} - e^{- i ( h_{k} - h_{m} ) x} }{ ( \Delta_{K} - h_{k} ) } 
\biggl[ 
(UX)_{\alpha k} W^*_{\beta K} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\nonumber \\
&+& 
W_{\alpha K} (UX)^*_{\beta k} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\biggr]
\nonumber \\ 
&+& 
\sum_{m}
\sum_{K} 
e^{- i (\Delta_{K} - h_{m} ) x} 
W_{\alpha K} W^*_{\beta K} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\biggr\}. 
\label{P-beta-alpha-0th+2nd}\end{aligned}$$ Notice that there is no matter dependent terms without suppression either by high-frequency oscillations $\propto \cos (\Delta_{K} - h_{m}) x$ (or $\sin$), or by large energy denominators $\propto \frac{ 1 }{ \Delta_{K} - h_{k} }$.

We take averaging over fast oscillations due to active-sterile and sterile-sterile mass squared differences which leads to $$\begin{aligned}
\left\langle \sin \Delta_{J i} x \right\rangle 
&\approx& 
\left\langle \sin \Delta_{J K} x \right\rangle \approx 0, 
\nonumber \\ 
\left\langle \cos \Delta_{J i} x 
\right\rangle &\approx& 
\left\langle \sin \Delta_{J K} x
\right\rangle \approx 0,
\label{average-out}\end{aligned}$$ where $\langle ... \rangle$ stands for averaging over neutrino energy within the uncertainty of energy resolution, as well as averaging over uncertainty of distance between production and detection points of neutrinos.[^12] The second approximate equalities in assume that there is no accidental degeneracy among the sterile state masses. That is, we assume that the relation $|\Delta m^2_{JK}| \gg |\Delta m^2_{31}|$ always holds.

After averaging out the fast oscillations, $P(\nu_\beta \rightarrow \nu_\alpha)$ is given to second order in $W$ by $$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(0+2)} 
\nonumber \\
&=& 
\left| \sum_{j=1}^{3} U_{\alpha j} U^{*}_{\beta j} \right|^2 - 
2 \sum_{j \neq k} 
\mbox{Re} 
\left[ (UX)_{\alpha j} (UX)_{\beta j}^* (UX)_{\alpha k}^* (UX)_{\beta k} \right] 
\sin^2 \frac{ ( h_{k} - h_{j} ) x  }{ 2 }
\nonumber\\
&-&
\sum_{j \neq k} \mbox{Im} 
\left[ (UX)_{\alpha j} (UX)_{\beta j}^* (UX)_{\alpha k}^* (UX)_{\beta k} \right] 
\sin ( h_{k} - h_{j} ) x 
\nonumber \\
&+&
2 \mbox{Re} 
\biggl\{
\sum_{m}
\sum_{k, K} 
\frac{ 1 }{ \Delta_{K} - h_{k} } 
\left[
(ix) e^{- i ( h_{k} - h_{m} ) x} 
- \frac{ e^{- i ( h_{k} - h_{m} ) x} }{ ( \Delta_{K} - h_{k} )  } 
\right]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta k} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\
&-& 
\sum_{m}
\sum_{k \neq l} \sum_{K} 
\frac{ 1 }{ ( h_{l}  - h_{k} ) (\Delta_{K} - h_{k}) (\Delta_{K} - h_{l}) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{k} \right) e^{- i ( h_{l} - h_{m} ) x} 
- \left( \Delta_{K} - h_{l} \right) e^{- i ( h_{k} - h_{m} ) x} 
\biggr]
\nonumber \\
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\nonumber \\
&-& 
\sum_{m}
\sum_{k, K} 
\frac{ e^{- i ( h_{k} - h_{m} ) x} }{ ( \Delta_{K} - h_{k} ) } 
\biggl[ 
(UX)_{\alpha k} W^*_{\beta K} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\nonumber \\
&+& 
W_{\alpha K} (UX)^*_{\beta k} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\biggr]
\biggr\},
\label{P-beta-alpha-2nd-averaged}\end{aligned}$$ where the zeroth-order terms are rewritten in a familiar way. However, the expression in (\[P-beta-alpha-2nd-averaged\]) does not mean a great simplification. That is, the matter dependent terms with sterile sector model-dependence remain.

 Suppression by the energy denominator {#sec:energy-denominator}
--------------------------------------

We, therefore, have to examine the effect of suppression by the large energy denominator which characterizes transition between active-sterile sates, $1/ ( \Delta_{K} - h_{k} )$. We demand that the matter dependent terms in (\[P-beta-alpha-2nd-averaged\]) be smaller than the probability leaking and the normalization terms of order $\sim W^4$. It leads to $$\begin{aligned}
\biggl | \frac{ AA L }{ ( \Delta_{J} - h_{i} ) } \biggr | 
< W^2, 
\hspace{4mm}
\biggl | \frac{ AA }{ ( h_{k} - h_{j} ) ( \Delta_{J} - h_{i} ) } \biggr | 
< W^2, 
\hspace{4mm}
\text{and}
\hspace{4mm}
\biggl | \frac{ A }{ ( \Delta_{J} - h_{i} ) } \biggr | 
< W^2, 
\label{denominator-size}\end{aligned}$$ where $L$ is the baseline distance and $i$ and $J$ denote, respectively, generic indices for active and sterile states. For notational convenience, we define $\lambda_{i}$ $(i=1,2,3)$ to be the eigenvalues of $3 \times 3$ submatrix $2E \tilde{H}_0$ in (\[tilde-H0+H1\]) corresponding to the active neutrino mass squared in matter. Then, $h_{i} = \frac{ \lambda_{i} }{2E}$. As one can easily see[^13] the last one in (\[denominator-size\]) gives the severest constraint (taking the matter potential due to CC in $A$ and removing the factor $\frac{1}{2E}$) $$\begin{aligned}
\biggl | \frac{ a }{ m^2_{J} - \lambda_{i} } \biggr | 
\approx 
\biggl | \frac{ a }{ \Delta m^2_{J i} } \biggr | 
< W^2.
\label{suppression-cond}\end{aligned}$$ Notice that, in order for the first inequality in (\[suppression-cond\]) to be valid, we have restricted the energy region for a given matter density such that $\lambda_{i}$ remain in the order of active neutrino masses. Roughly speaking, it corresponds to $- 50 \,\text{ (g/cm}^3) \text{GeV} \lsim Y_{e} \rho E \lsim 50 \, \text{ (g/cm}^3) \text{GeV}$. See e.g., figure 3 of ref. [@Minakata:2015gra]. Clearly, it excludes the interesting region of “IceCube resonance” due to sterile neutrino mass of eV scales [@Nunokawa:2003ep], for which an entirely different theoretical framework would be necessary.

Then, we notice that in a regime $W^2 \sim 10^{-2}$, the condition in (\[suppression-cond\]) is valid given the estimation (assuming $Y_{e} = 0.5$) $$\begin{aligned}
\biggl | \frac{ a }{ \Delta m^2_{J i} } \biggr | = 2.13 \times 10^{-3} 
\left(\frac{ \Delta m^2_{J i} }{ 0.1~\mbox{eV}^2}\right)^{-1}
\left(\frac{\rho}{2.8 \,\text{g/cm}^3}\right) \left(\frac{E}{1~\mbox{GeV}}\right), 
\label{rA-def-value}\end{aligned}$$ unless $\rho E \gsim 10 \, \text{ (g/cm}^3) \text{GeV}$. That is, the second-order matter dependent correction terms can be ignored in comparison with $\mathcal{O} (W^4)$ terms if $\Delta m^2_{J k} \gsim 0.1$ eV$^2$, which is already required in vacuum. If we want to treat the regime $W^2 \gsim 10^{-n}$, we need to limit the sterile masses to $\Delta m^2_{J k} \simeq m^2_{J} \gsim 10^{(n-3)}$ eV$^2$ to keep our $(3+N)$ space unitary model insensitive to details of the sterile sector [@Fong:2016yyh]. We note, however, that terms of order $W^4 \sim 10^{-4}$ may be the limit of exploration for near future neutrino oscillation experiments.

The condition (\[suppression-cond\]) is identical with the one obtained using the first order matter perturbation theory [@Fong:2016yyh], which may look strange to the readers. Let us understand the reason why taking care of all order matter effect does not alter the condition obtained by first-order treatment in matter perturbation theory. The matter-dependent term in the zeroth-order Hamiltonian $\tilde{H}_{0}$ only involves $U$ matrix, but no $W$ matrix. Since we treat $\tilde{H}_{0}$ in an unperturbed fashion it produces all-order effect of the matter potential which is however independent of $W$ matrix elements. On the other hand, perturbative effects that come from single or double powers in $W$ in $\hat{H}_{1}$ are always accompanied by the matter potential in the form of $WA$ or $W^{\dagger} A$, as in eq. (\[H1-matrix\]). That is, perturbative effect of $W$ is always accompanied by matter potential, and hence can always be dealt with matter perturbation theory. It is the reason why the matter perturbation theory is able to yield the same condition on sterile masses as obtained in a fuller treatment of matter effect done in this paper.[^14]

The oscillation probability in fourth order in $W$ {#sec:probability-4th}
--------------------------------------------------

The oscillation probability in fourth order in $W$ contains the two terms $$\begin{aligned}
&&P(\nu_\beta \rightarrow \nu_\alpha)^{(4)} = 
\left| S^{(2)}_{\alpha \beta} \right|^2 
+ 2 \mbox{Re} \left[ \left( S^{(0)}_{\alpha \beta} \right)^{*} S^{(4)}_{\alpha \beta} \right]. 
\label{P-beta-alpha-4th-def}\end{aligned}$$ We will show in appendix \[sec:second-order-square\] that the first term in (\[P-beta-alpha-4th-def\]), after averaging over the fast oscillations and using the suppression by energy denominator as discussed in the previous section, leaves the unique term, the probability leaking term $\mathcal{C}_{\alpha \beta}$ in (\[P-beta-alpha-ave-vac\]). We will also show in appendix \[sec:interference\] that the second term in (\[P-beta-alpha-4th-def\]), under the same treatment for the first term, gives vanishing contribution. Therefore, no matter-dependent term remains after using energy denominator suppression and averaging over the fast oscillations.

In conclusion, the oscillation probability in matter between active flavor neutrinos in the $(3+N)$ space unitary model to fourth order in $W$ in our small unitarity-violation perturbation theory can be written as in eq. (\[P-beta-alpha-final\]) in section \[sec:essence\]. We hope that it serves as a useful tool to test leptonic unitarity in various ongoing and future neutrino oscillation experiments.

 Analytical and numerical methods for solving non-unitary evolution in matter  {#sec:analytical-numerical}
==============================================================================

In this section, we describe the numerical and analytical methods for calculating the neutrino oscillation probability by solving non-unitary evolution in matter.

 Numerical method for calculating neutrino oscillation probability  {#sec:numerical}
-------------------------------------------------------------------

We describe a numerical method for computing the oscillation probability in matter. This method can be used, assuming adiabaticity, in cases with varying matter density. We show that in zeroth order in $W$ the system simplifies to an evolution equation in the $3 \times 3$ active subspace.

We solve the Schrödinger equation in the vacuum mass eigenstate basis (“tilde basis”), $\tilde{\nu}_{z} = ({\bf U}^{\dagger})_{z \zeta} \nu_{\zeta}$ with Hamiltonian $\tilde{H}$ in (\[tilde-H\]): $$\begin{aligned}
i \frac{d}{dx} 
\left[
\begin{array}{c}
\tilde{\nu}_{i} \\
\tilde{\nu}_{J} \\
\end{array}
\right] = 
\left[
\begin{array}{cc}
{\bf \Delta_{a} } + U^{\dagger} A U & U^{\dagger} A W \\
W^{\dagger} A U & {\bf \Delta_{s} } + W^{\dagger} A W \\
\end{array}
\right] \left[
\begin{array}{c}
\tilde{\nu}_{i} \\
\tilde{\nu}_{J} \\
\end{array}
\right],
\label{Schroedinger-eq}\end{aligned}$$ where $i = 1,2,3$ and $J= 4,5,\cdot \cdot \cdot,3+N$ denote mostly active and mostly sterile neutrino mass eigenstate labels, respectively. The initial condition with only active component implies $$\begin{aligned}
\tilde{\nu}_{i} (0) &=& \sum_{\alpha} (U^{\dagger})_{i \alpha} \nu_{\alpha} (0), \nonumber \\
\tilde{\nu}_{J} (0) &=& \sum_{\alpha} (W^{\dagger})_{J \alpha} \nu_{\alpha} (0). 
\label{initial-condition}\end{aligned}$$ Using the solution of equation (\[Schroedinger-eq\]), we need the wave function of active flavor component to calculate the probability at baseline $x=L$. $$\begin{aligned}
\nu_{\alpha} (L) &=& 
\sum_{i} U_{\alpha i} \tilde{\nu}_{i} (L) 
+ \sum_{J} W_{\alpha J} \tilde{\nu}_{J} (L).
\label{final-condition}\end{aligned}$$ Therefore, in the mass-basis formulation only $U$ and $W$ are involved, which is consistent with our experience in $W$ perturbation theory. An apparent contradiction to this property that one faces in the evolution equation in the flavor basis is resolved in appendix \[sec:flavor-basis-evolution\].

A drawback of this method is that we have to solve explicitly the evolution of the sterile states which are coupled to the active states. Then, we need to specify the sterile sector model, and have to know how to deal with averaging over the fast modes.

We notice, however, that in the zeroth-order in $W$ the system simplifies. Since the Hamiltonian $\tilde{H}$ is block-diagonal it suffices to solve the equation only in the $3 \times 3$ active neutrino subspace: $$\begin{aligned}
i \frac{d}{dx} \nu_{i} = 
\sum_{j}
\left( {\bf \Delta_{a} } + U^{\dagger} A U \right)_{ij} 
\nu_{j}. 
\label{Schroedinger-eq-0th}\end{aligned}$$ The initial condition (\[initial-condition\]) and final reverse-back formula (\[final-condition\]) involve only $U$ matrix elements. Therefore, the oscillation probability in the zeroth-order in $W$ can be calculable in a manner independent of sterile sector models.[^15]

 An exact solution of zeroth-order oscillation probability  {#sec:exact-solution-zeroth}
-----------------------------------------------------------

Here, we describe a method for obtaining the analytical solution of the zeroth-order Hamiltonian. The exact solution, as well as the numerical one described in the previous section, provides the basis for computing the higher order corrections in $W$.

We calculate an exact form of the oscillation probability $P(\nu_\beta \rightarrow \nu_\alpha)$ in leading order in our perturbative framework, the one in (\[P-beta-alpha-final\]) except for $\mathcal{C}_{\alpha \beta}$, in the case of uniform matter density.

The zeroth-order $S$ matrix element $S_{\alpha \beta}^{(0)}$ in (\[S-alpha-beta-0th\]) can be written as $$\begin{aligned}
S_{\alpha \beta}^{(0)} &=& 
\sum_{i, j} U_{\alpha i} U^*_{\beta j} 
\left( 
\sum_{k} X_{i k} X^*_{j k} e^{-i h_{k} x}
\right),
\label{S-alpha-beta-0th-KTY}\end{aligned}$$ and the factor in parenthesis can be calculated by the KTY technique [@Kimura:2002wd]. We want to diagonalize the Hamiltonian $$\begin{aligned}
H_{0} \equiv 
\frac{1}{2E} 
\left\{  
\left[
\begin{array}{ccc}
m^2_{1} & 0 & 0 \\
0 & m^2_{2} & 0 \\
0 & 0 & m^2_{3} \\
\end{array}
\right] + 
U^{\dagger} \left[
\begin{array}{ccc}
a - b & 0 & 0 \\
0 & -b & 0 \\
0 & 0 & -b \\
\end{array}
\right] U 
\right\} ,
\label{H0-def}\end{aligned}$$ the active $3 \times 3$ block of $\tilde{H}_0$ in (\[tilde-H0+H1\]). We have defined in eq. (\[H0-diag\]) the unitary matrix $X$ which diagonalize $H_{0}$ as $$\begin{aligned}
H_{0} = 
\frac{1}{2E} 
X \left[
\begin{array}{ccc}
\lambda_{1} & 0 & 0 \\
0 & \lambda_{2} & 0 \\
0 & 0 & \lambda_{3} \\
\end{array}
\right] X^{\dagger}
\equiv 
H_{d}. 
\label{H0-diagonal}\end{aligned}$$ For our notational convenience we call this form of $H_{0}$ as $H_{d}$. Note that $h_{i} = \frac{ \lambda_{i} }{2E}$ where $$\begin{aligned}
\lambda_{1} & = & \frac{T}{3}-\frac{1}{3}\sqrt{T^{2}-3A}\cos\left[\frac{1}{3}\arccos\left(\frac{2T^{3}-9AT+27D}{2\left(T^{2}-3A\right)^{3/2}}\right)\right]\nonumber \\
 &  & -\frac{1}{\sqrt{3}}\sqrt{T^{2}-3A}\sin\left[\frac{1}{3}\arccos\left(\frac{2T^{3}-9AT+27D}{2\left(T^{2}-3A\right)^{3/2}}\right)\right],\\
\lambda_{2} & = & \frac{T}{3}-\frac{1}{3}\sqrt{T^{2}-3A}\cos\left[\frac{1}{3}\arccos\left(\frac{2T^{3}-9AT+27D}{2\left(T^{2}-3A\right)^{3/2}}\right)\right]\nonumber \\
 &  & +\frac{1}{\sqrt{3}}\sqrt{T^{2}-3A}\sin\left[\frac{1}{3}\arccos\left(\frac{2T^{3}-9AT+27D}{2\left(T^{2}-3A\right)^{3/2}}\right)\right],\\
\lambda_{3} & = & \frac{T}{3}+\frac{2}{3}\sqrt{T^{2}-3A}\cos\left[\frac{1}{3}\arccos\left(\frac{2T^{3}-9AT+27D}{2\left(T^{2}-3A\right)^{3/2}}\right)\right],\end{aligned}$$ with $$\begin{aligned}
T & = & (2E)\, {\rm Tr} H_0,\\
A & = & (2E)^2\, {\rm Tr}\left({\rm Adj} H_0 \right),\\
D & = & (2E)^3 \det H_0.\end{aligned}$$ The adjugate of $H_{0}$ is defined as $\text{Adj} H_{0} \equiv (H_{0})^{-1} \text{det} H_{0}$. Notice that $T$, $A$ and $D$ are invariant under unitary transformation of $ H_0\to K H_0 K^{\dagger}$ with $K$ any unitary matrix and so are $\lambda_i$.

Following the notation in [@Kimura:2002wd] we define $p_{ij}$ and $q_{ij}$ as ($i,j=1,2,3$) $$\begin{aligned}
\frac{ p_{ij} }{ 2E } &\equiv& \left( H_{0} \right)_{ij}, 
\nonumber \\
\frac{ q_{ij} }{ (2E)^2 } &\equiv& \left( \text{Adj} H_{0} \right)_{ij}. 
\label{pq-def}\end{aligned}$$ Notice that $p_{ij}$ and $q_{ij}$ are written only by the known (or given) quantities. Then, the equations $$\begin{aligned}
\left( H_{d} \right)_{ij} &=& \frac{ p_{ij} }{ 2E }, 
\nonumber \\
\left( \text{Adj} H_{d} \right)_{ij} &=& \frac{ q_{ij} }{ (2E)^2 }, 
\label{KTY-eq}\end{aligned}$$ together with unitarity of $X$, become the equations to determine $X X^{\dagger}$: $$\begin{aligned}
X_{i 1} X_{j 1}^* + X_{i 2} X_{j 2}^* + X_{i 3} X_{j 3}^* &=& \delta_{ij}, 
\nonumber \\
\lambda_{1} X_{i 1} X_{j 1}^* + \lambda_{2} X_{i 2} X_{j 2}^* + \lambda_{3}X_{i 3} X_{j 3}^* &=& p_{ij}, 
\nonumber \\
\lambda_{2} \lambda_{3} X_{i 1} X_{j 1}^* + \lambda_{3} \lambda_{1} X_{i 2} X_{j 2}^* + \lambda_{1} \lambda_{2} X_{i 3} X_{j 3}^* &=& q_{ij}. 
\label{KTY-eq-explicit}\end{aligned}$$ They lead to the solution ($k=1,2,3$) $$\begin{aligned}
X_{i k} X_{j k}^* = 
\frac{ q_{ij} + p_{ij} \lambda_{k} - \delta_{ij} \lambda_{k} ( \lambda_{l} + \lambda_{m} ) }{ (\lambda_{l} -\lambda_{k} ) (\lambda_{m} -\lambda_{k} ) },
\label{KTY-eq-solution}\end{aligned}$$ where $k,l,m$ is cyclic, and sum over $k$ is not implied in (\[KTY-eq-solution\]).

Therefore, to zeroth-order in $W$ expansion, the $S$ matrix elements are given by $$\begin{aligned}
S_{\alpha \beta}^{(0)} &=& 
\sum_{k} 
\left(
\sum_{i, j} 
U_{\alpha i} 
\left[ q_{ij} + p_{ij} \lambda_{k} - \delta_{ij} \lambda_{k} ( \lambda_{l} + \lambda_{m} ) \right] 
U^*_{\beta j} 
\right)
\frac{ e^{-i h_{k} x} }{ (\lambda_{l} -\lambda_{k} ) (\lambda_{m} -\lambda_{k} ) }, 
\label{S-alpha-beta-0th-final}\end{aligned}$$ and the oscillation probability by $P(\nu_\beta \rightarrow \nu_\alpha) = \vert S_{\alpha \beta}^{(0)} \vert^2$.

Finally, armed with the solution , we can also calculate all higher order terms in oscillation probability for e.g. those in eq.  since only such combination $X_{ik} X_{jk}^*$ (no sum over $k$ implied) can appear.

 Low-scale vs. high-scale unitarity violation and $W$ corrections {#sec:UV-low-high-E}
=================================================================

 Low-scale versus high-scale unitarity violation  {#sec:low-vs-high}
-------------------------------------------------

In leptonic unitarity test, a clear understanding of the relationship between low-scale and high-scale unitarity violation may be one of the key issues. While the presence of probability leaking term $\mathcal{C}_{\alpha \beta}$, if detected, clearly testifies for low-scale unitarity violation [@Fong:2016yyh], it would be better to provide a global view of the difference between them, looking for alternative ways to reach the goal. In this section, we would like to give a preliminary discussion on this topics.

One would guess, on intuitive ground, that at zeroth order in $W$ our system describes high-scale unitarity violation. There is no “$W$ corrections” in high-scale unitarity violation because the energy scale is so high that the high-mass sector is truncated. It is in agreement with the formulations in ref. [@Blennow:2016jkn] with which we share the same evolution equation (\[Schroedinger-eq-0th\]) in the vacuum mass eigenstate basis. See also [@Antusch:2006vwa].[^16] Based on this perspective, we just conclude this subsection by stating the generic characteristic differences between high-scale and low-scale unitarity violation:

-   Our system calculated with evolution equations (\[Schroedinger-eq-0th\]) in leading (zeroth) order in $W$ expansion, which is common to both high- and low-scale unitarity violation, would serve as a first indicator to the existence of leptonic unitarity violation. Existence of second and higher order corrections in $W$, if detected, uniquely identifies the case for low-scale unitarity violation.

-   Presence (low-$E$) or absence (high-$E$) of the probability leaking term $\mathcal{C}_{\alpha \beta}$ in (\[Cab\]) remains to be the best possible way to distinguish between high- and low-scale unitarity violation.

 Oscillation probabilities with and without unitarity violation  {#sec:probabilities}
----------------------------------------------------------------

Under the hope to uncover in which region of parameter space the effect of non-unitarity is most prominent, we first compare the oscillation probabilities with and without unitarity violation. We examine the three channels $\nu_{\mu} \rightarrow \nu_{e}$, $\nu_{\mu} \rightarrow \nu_{\tau}$, and $\nu_{\mu} \rightarrow \nu_{\mu}$. We do not intend to do quantitative analyses, nor attempt to cover the whole parameter space. Yet, we try to give the readers a feeling on how large and where are the effects of unitarity violation. Therefore, we present the results just for a specific choice of the parameters. Also, the use of uniform matter density $\rho = 3.2~{\rm g\,cm}^{-3} $ over the entire baseline is not realistic.[^17]

Notice that the results given in section \[sec:probabilities\] apply to both high-scale unitarity violation as well as low-scale one in its leading order in $W$. The probability leaking term is not included in the analysis in this section \[sec:probabilities\], but the effect is discussed in section \[sec:correction-terms\].

###  A brief note on the parameter choice  {#sec:parameter-choice}

Let us start by explaining briefly how the unitarity violation parameters are chosen. We go to the $(3+1)$ model in which the constraints on the parameters are best understood [@Kopp:2013vaa; @deGouvea:2015euy; @Fernandez-Martinez:2016lgt; @TheIceCube:2016oqi]. In consistent with the current ones we have chosen: $\sin^2\theta_{14} = 0.02$, $\sin^2\theta_{24} = 0.01$, and $\sin^2\theta_{34} = 0.1$ for $\Delta m^2_{41} = 0.1$ eV$^2$, and set all the CP phases to zero. Then, we cut out the $3\times3$ active neutrino mixing matrix, which is non-unitary.[^18] For the Standard Model mixing parameters in $U_{\text{\tiny PDG}}$, we take $\sin^2\theta_{12} = 0.3$, $\sin^2\theta_{23} = 0.5$, $\sin^2 (2\theta_{13}) = 0.09$, and the mass squared differences $\Delta m_{21}^2 = 7.4 \times 10^{-5}$ eV$^2$ and $\Delta m_{31}^2 = 2.4 \times 10^{-3}$ eV$^2$, and set the CP phase $\delta_{\text{CP}}$ to zero.

![ In the upper panel (a), presented is the iso-contour of $P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{non-unitary} }^{(0)}$ in space spanned by neutrino energy $E$ and baseline $L$. In the lower panel (b), the iso-contour of the difference $\Delta P (\nu_{\mu} \rightarrow \nu_{e}) \equiv P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{standard} }- P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{non-unitary} }^{(0)}$ is presented. For the values of unitarity-violating as well as the standard mixing parameters taken, see the text. []{data-label="fig:Pmue_energy_dist"}](Pmue_energy_dist_non_unitary_small_size.jpeg "fig:"){width="85.00000%"} ![ In the upper panel (a), presented is the iso-contour of $P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{non-unitary} }^{(0)}$ in space spanned by neutrino energy $E$ and baseline $L$. In the lower panel (b), the iso-contour of the difference $\Delta P (\nu_{\mu} \rightarrow \nu_{e}) \equiv P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{standard} }- P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{non-unitary} }^{(0)}$ is presented. For the values of unitarity-violating as well as the standard mixing parameters taken, see the text. []{data-label="fig:Pmue_energy_dist"}](Pmue_energy_dist_difference_small_size.jpeg "fig:"){width="85.00000%"}

![ In the upper panel (a), presented is the iso-contour of $P(\nu_{\mu} \rightarrow \nu_{\tau})_{ \text{non-unitary} }^{(0)}$ in $E-L$ space. In the lower panel (b), the iso-contour of the difference $\Delta P (\nu_{\mu} \rightarrow \nu_{\tau}) \equiv P(\nu_{\mu} \rightarrow \nu_{\tau})_{ \text{standard} }- P(\nu_{\mu} \rightarrow \nu_{\tau})_{ \text{non-unitary} }^{(0)}$ is presented. The parameters used are the same as in figure \[fig:Pmue\_energy\_dist\]. []{data-label="fig:Pmutau_energy_dist"}](Pmutau_energy_dist_non_unitary_small_size.jpeg "fig:"){width="85.00000%"} ![ In the upper panel (a), presented is the iso-contour of $P(\nu_{\mu} \rightarrow \nu_{\tau})_{ \text{non-unitary} }^{(0)}$ in $E-L$ space. In the lower panel (b), the iso-contour of the difference $\Delta P (\nu_{\mu} \rightarrow \nu_{\tau}) \equiv P(\nu_{\mu} \rightarrow \nu_{\tau})_{ \text{standard} }- P(\nu_{\mu} \rightarrow \nu_{\tau})_{ \text{non-unitary} }^{(0)}$ is presented. The parameters used are the same as in figure \[fig:Pmue\_energy\_dist\]. []{data-label="fig:Pmutau_energy_dist"}](Pmutau_energy_dist_difference_small_size.jpeg "fig:"){width="85.00000%"}

![ In the upper panel (a), presented is the iso-contour of $P(\nu_{\mu} \rightarrow \nu_{\mu})_{ \text{non-unitary} }^{(0)}$ in $E-L$ space. In the lower panel (b), the iso-contour of the difference $\Delta P (\nu_{\mu} \rightarrow \nu_{\mu}) \equiv P(\nu_{\mu} \rightarrow \nu_{\mu})_{ \text{standard} }- P(\nu_{\mu} \rightarrow \nu_{\mu})_{ \text{non-unitary} }^{(0)}$ is presented. The parameters used are the same as in figure \[fig:Pmue\_energy\_dist\]. []{data-label="fig:Pmumu_energy_dist"}](Pmumu_energy_dist_non_unitary_small_size.jpeg "fig:"){width="85.00000%"} ![ In the upper panel (a), presented is the iso-contour of $P(\nu_{\mu} \rightarrow \nu_{\mu})_{ \text{non-unitary} }^{(0)}$ in $E-L$ space. In the lower panel (b), the iso-contour of the difference $\Delta P (\nu_{\mu} \rightarrow \nu_{\mu}) \equiv P(\nu_{\mu} \rightarrow \nu_{\mu})_{ \text{standard} }- P(\nu_{\mu} \rightarrow \nu_{\mu})_{ \text{non-unitary} }^{(0)}$ is presented. The parameters used are the same as in figure \[fig:Pmue\_energy\_dist\]. []{data-label="fig:Pmumu_energy_dist"}](Pmumu_energy_dist_difference_small_size.jpeg "fig:"){width="85.00000%"}

###  $P(\nu_{\mu} \rightarrow \nu_{e})$ 

In figure \[fig:Pmue\_energy\_dist\]-(a) (upper panel) and (b) (lower panel), presented are the iso-contours of $P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{non-unitary} }^{(0)}$ and $\Delta P (\nu_{\mu} \rightarrow \nu_{e}) \equiv P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{standard} }- P(\nu_{\mu} \rightarrow \nu_{e})_{ \text{non-unitary} }^{(0)}$ in $E - L$ space. Here, the superscript $(0)$ implies that it is calculated in zeroth-order in $W$ using (\[Schroedinger-eq-0th\]) with appropriate initial condition and final projection to flavor eigenstate. In most of the $E - L$ space $P(\nu_{\mu} \rightarrow \nu_{e})$ is small. However, we identify the two regions where $P(\nu_{\mu} \rightarrow \nu_{e})$ is relatively large, $\gsim 0.3$. One of them is at low energy, $E \lsim \text{a few hundred}$ MeV, and baseline $L \gsim$1000 km. The other one is a region $E \sim 10$ GeV and $L \sim$ 10000 km. The former may be understood as due to the solar MSW enhancement, and the latter as the atmospheric MSW enhancement [@Mikheev:1986gs; @Wolfenstein:1977ue]. Roughly speaking, the regions with relatively large $| \Delta P (\nu_{\mu} \rightarrow \nu_{e}) |$ overlap with these regions.

###  $P(\nu_{\mu} \rightarrow \nu_{\tau})$ and $P(\nu_{\mu} \rightarrow \nu_{\mu})$ 

In figures \[fig:Pmutau\_energy\_dist\] and \[fig:Pmumu\_energy\_dist\], the same quantities (in each upper (a) and lower (b) panel) are presented but in $\nu_{\mu} \rightarrow \nu_{\tau}$ and $\nu_{\mu} \rightarrow \nu_{\mu}$ channels, respectively. In contrast to $\nu_{\mu} \rightarrow \nu_{e}$ channel, $P(\nu_{\mu} \rightarrow \nu_{\tau})$ and $P(\nu_{\mu} \rightarrow \nu_{\mu})$ contours are globally “vacuum effect dominated”, apart from the solar MSW region, both in the standard (not shown) and the non-unitary cases. The first oscillation peak of $P(\nu_{\mu} \rightarrow \nu_{\tau})$ scales roughly as the vacuum oscillation peak does, $L / 10^3 \,\text{km} = 0.33 E / 1 \, \text{GeV}$. This feature is more or less seen in $P(\nu_{\mu} \rightarrow \nu_{e})$, but $P(\nu_{\mu} \rightarrow \nu_{\tau})$ has a higher peak height $\simeq 0.7-0.8$, and the effect of atmospheric MSW enhancement is less prominent.

For $P(\nu_{\mu} \rightarrow \nu_{\mu})_{ \text{non-unitary} }^{(0)}$, roughly speaking, the relation $P(\nu_{\mu} \rightarrow \nu_{\mu}) \approx 1 - P(\nu_{\mu} \rightarrow \nu_{\tau}) $ holds in region where $P(\nu_{\mu} \rightarrow \nu_{e})$ is small. It must be the case in the unitary case, but even in non-unitary case the relation holds approximately because unitarity violation is small in our choice of the parameters. Therefore, $P(\nu_{\mu} \rightarrow \nu_{\mu})$ is large in region where $P(\nu_{\mu} \rightarrow \nu_{\tau})$ is small, and vice versa, as seen in figure \[fig:Pmumu\_energy\_dist\]. It appears that the anticorrelation is inherited to the relationship between $\Delta P (\nu_{\mu} \rightarrow \nu_{\mu})$ and $\Delta P (\nu_{\mu} \rightarrow \nu_{\tau})$.

Due to relatively large $\Delta P (\nu_{\mu} \rightarrow \nu_{\tau})$ in the first and second oscillation maxima over the entire baseline, or similarly because of relatively large depletion of $\Delta P (\nu_{\mu} \rightarrow \nu_{\mu})$ there, it should be possible to detect the effect of non-unitarity if the detector has a good $\tau$ (in the former channel) and $\mu$ (in the latter channel) detection capabilities. It would be ideal that the detector has sensitivities to the both channels because then one can take anticorrelation between $\mu$ and $\tau$ yields to enhance the sensitivity to unitarity violation.

Some comments on observational aspects: In the two regions where $| \Delta P (\nu_{\mu} \rightarrow \nu_{e}) |$ is large, and $\Delta P (\nu_{\mu} \rightarrow \nu_{\mu})$ in energy region $E \lsim 10$ GeV may be explored by high-statistics atmospheric neutrino observation by Super-K/Hyper-K, or DUNE [@Abe:2017aap; @Abe:2015zbg; @Acciarri:2015uup]. The atmospheric MSW enhanced region of $P(\nu_{\mu} \rightarrow \nu_{e})$ would be a good target for PINGU extensions of IceCube and KM3NeT-ORCA [@TheIceCube-Gen2:2016cap; @Adrian-Martinez:2016zzs]. $P(\nu_{\mu} \rightarrow \nu_{\tau})$ and $P(\nu_{\mu} \rightarrow \nu_{\mu})$ would be explored by them, with possibility of seeing anticorrelation between $\mu$ and $\tau$ yields. Although it is very interesting to investigate these experimental prospects, a detailed examination of these questions is beyond the scope of this paper.

 The probability leaking and $W$ correction terms  {#sec:correction-terms}
--------------------------------------------------

As we stated in section \[sec:low-vs-high\], uncovering the effect of $W$ correction or the probability leaking term $\mathcal{C}_{\alpha \beta}$ in eq. (\[Cab\]) would offer another way of distinguishing low-scale unitarity violation from high-scale one. Let us give a brief sketch of how and where we might see visible effect.

###  A further note on the parameter choice  {#sec:parameter-choice2}

To discuss $W$ correction and the probability leaking term we have to determine the $W$ matrix. Given the non-unitary $U$ matrix there is a way to construct the $W$ matrix. In general, it is given by $$W = S \sqrt{w} R, 
\label{W-construction}$$ where $S$ is a $3\times 3$ matrix which diagonalizes ${\bf 1}_{3\times 3} - UU^\dagger$, $w$ is diagonal matrix which consists of eigenvalues of ${\bf 1}_{3\times 3} - UU^\dagger$, and $R$ is an arbitrary $3\times N$ complex matrix obeying $RR^\dagger = \bf{1}_{3 \times 3}$. The construction makes sense for $N \geq 3$. Therefore, for a given $N$ there is a large arbitrariness on the choice of the $W$ matrix, and hence on the sizes of the $W$ corrections and $\mathcal{C}_{\alpha \beta}$.

Lacking a guiding principle of how to choose the $R$ matrix in (\[W-construction\]), we examine the cases with largest and smallest possible values of $\mathcal{C}_{\alpha \beta}$ for given values of unitarity violation $1 - \sum_{j=1}^3 |U_{\alpha j}|^2$ ($\alpha=e,\mu,\tau$). It is shown that in the $(3+N)$ model $\mathcal{C}_{\alpha \beta}$ is bounded from above and below as [@Fong:2016yyh] $$\frac{1}{N}
\biggl( 1 - \sum_{j=1}^3 |U_{\alpha j}|^2 \biggr) 
\biggl( 1 - \sum_{j=1}^3 |U_{\beta j}|^2 \biggr) 
\leq
\mathcal{C}_{\alpha \beta}
\leq
\biggl( 1 - \sum_{j=1}^3 |U_{\alpha j}|^2 \biggr)
\biggl( 1 - \sum_{j=1}^3 |U_{\beta j}|^2 \biggr).
\label{Cab-bound}$$ In the $(3+1)$ model, the $W$ matrix elements are unique, with the upper and lower bound being equal. For the numbers given in section \[sec:parameter-choice\], we have $W_{e4} = 0.141$, $W_{\mu 4} = 0.099$, and $W_{\tau 4} = 0.141$ assuming that they are real. Then, the leaking constants have the unique values, $\mathcal{C}_{e \mu}^{(N=1)} = 2 \times 10^{-4}$, $\mathcal{C}_{\mu \mu}^{(N=1)} = 9.6 \times 10^{-5}$, and $\mathcal{C}_{\tau \mu}^{(N=1)} = 9.5 \times 10^{-4}$. The lower bound is realized in the “universal scaling” model described in appendix \[sec:scaling-model\], which predicts $W_{\alpha J} = \frac{ 1 }{ \sqrt{N} } W_{\alpha 4}^{(N=1)}$ ($J=4,5, \cdot \cdot \cdot, 3+N$).[^19] It is shown in appendix \[sec:scaling-model\] that under the assumption of equal sterile state masses the universal scaling model predicts the same $W^2$ correction terms as those of the $(3+1)$ model.

### How large are the $W$ corrections and $\mathcal{C}_{\alpha \beta}$?  {#sec:how-large}

Let us go back to the expression of the oscillation probability to second order in $W$, eq. (\[P-beta-alpha-2nd-averaged\]), in section \[sec:probability-2nd\] to know where we might see visible effects. If we enter into the region $\rho E \gg 10 \, \text{ (g/cm}^3) \text{GeV}$ at around the first oscillation maximum, the first two terms in $2 \mbox{Re} \{ \cdot \cdot \cdot \}$ in (\[P-beta-alpha-2nd-averaged\]) can become large apart from $W^2$ suppression, $$\begin{aligned}
\biggl | \frac{ AA L }{ ( \Delta_{J} - h_{i} ) } \biggr | 
&\sim& 
\biggl | \frac{ AA }{ ( h_{k} - h_{j} ) ( \Delta_{J} - h_{i} ) } \biggr | 
= 0.27 
\left(\frac{ \Delta m^2_{J i} }{ 0.1 \mbox{eV}^2}\right)^{-1}
\left( \frac{ \rho E }{ 100 (\text{g/cm}^3) \mbox{GeV} }\right)^2. 
\label{enhanced-case}\end{aligned}$$ It should be compared to (\[denominator-size\]). After taking account of $W^2$ suppression of $\sim 0.01$ (assuming $W \simeq 0.1$), $| \frac{ AA L }{ ( \Delta_{J} - h_{i} ) } W^2 | \sim 3 \times 10^{-2}$ at $E \sim 100$ GeV, assuming $\Delta m^2_{J i} =0.1$ eV$^2$.

![ The sum of the order $W^2$ correction terms (see eq. (\[P-beta-alpha-2nd-averaged\])) plus the probability leaking term $\mathcal{C}_{\mu \alpha}$ (see eq. (\[Cab\]) for definition) in $P(\nu_{\mu} \rightarrow \nu_{\alpha})$, namely, $\delta P(\nu_{\mu} \rightarrow \nu_{\alpha})
\equiv P(\nu_{\mu} \rightarrow \nu_{\alpha}) - 
P(\nu_{\mu} \rightarrow \nu_{\alpha})^{(0)}$, are plotted assuming a common $m_J^2 = 0.1$ eV$^2$. The top, middle and bottom panels are for $\alpha = e, \tau$, and $\mu$, respectively. In each panel the three cases are shown: $N=1$ case with maximal $\mathcal{C}_{\mu \alpha}$ (solid line), the universal scaling model with $N=3$ (dotted line), and the order $W^2$ correction only (dashed line). The last case corresponds to the universal scaling model with $N=\infty$. The blue lines are for $E=10$ GeV, and the red for $E=100$ GeV. The leaking constants in the $N=1$ model (shown without superscript $(N=1)$ in the legend) have values $\mathcal{C}_{e \mu} = 2 \times 10^{-4}$, $\mathcal{C}_{\tau \mu} = 9.5 \times 10^{-4}$, and $\mathcal{C}_{\mu \mu} = 9.6 \times 10^{-5}$. []{data-label="fig:W-correction"}](Delta_Prob_W_corrections.pdf){width="120.00000%"}

To know more quantitatively their sizes, we fix the parameters as discussed in section \[sec:parameter-choice2\] with a common $m_J^2 = 0.1$ eV$^2$,[^20] and plot in figure \[fig:W-correction\] the order $W^2$ correction terms in $P(\nu_{\mu} \rightarrow \nu_{\alpha})$ in eq. (\[P-beta-alpha-2nd-averaged\]) plus the probability leaking term $\mathcal{C}_{\mu \alpha}$, $\alpha = e$ (top panel), $\alpha = \tau$ (middle panel), and $\alpha = \mu$ (bottom panel). Under the approximation of ignoring the fourth (and higher-order) corrections, they are identical to $\delta P(\nu_{\mu} \rightarrow \nu_{\alpha})  \equiv P(\nu_{\mu} \rightarrow \nu_{\alpha}) - P(\nu_{\mu} \rightarrow \nu_{\alpha})^{(0)}$. In each panel the three cases are examined. $N=1$ case with maximal $\mathcal{C}_{\mu \alpha}$ (solid line), the universal scaling model with $N=3$ (dotted line), and the order $W^2$ correction only (dashed line). The last case corresponds to the universal scaling model with $N=\infty$. The blue lines are for $E=10$ GeV, and the red for $E=100$ GeV.

We will first focus on the appearance channels $\nu_\mu \to \nu_e$ and $\nu_\mu \to \nu_\tau$. At $E=10$ GeV (100 GeV) $\delta P$ depends very much on the above three cases, $N=1$, $N=3$, and $N=\infty$ for baseline $L$ of several 100 km ($L \gsim 1000$ km). The maximum value of $| \delta P |$ is always given by the case of maximal (minimal) $\mathcal{C}_{\mu \alpha}$ for positive (negative) $\delta P(\nu_{\mu} \rightarrow \nu_{\alpha})$ shown by the solid (dashed) lines. These maximal values of $| \delta P(\nu_{\mu} \rightarrow \nu_{\alpha}) |$ are, roughly speaking, $\simeq 10^{-3}$ for $\nu_{\mu} \rightarrow \nu_{\tau}$, and $\simeq \text{a few} \times 10^{-4}$ for $\nu_{\mu} \rightarrow \nu_{e}$. The effect might be visible for the former, though it might be challenging for the latter channel.[^21] At longer distance, however, we see enhancement. At $E=10$ GeV, we observe a factor of several enhancement both in $| \delta P (\nu_{\mu} \rightarrow \nu_{\alpha})|$ for $\nu_{\mu} \rightarrow \nu_{e}$ and $\nu_{\mu} \rightarrow \nu_{\tau}$ in region $L \gsim 3000$ km, which may provide a clearer signature. The similar tendency exists at $E=100$ GeV, but in a less pronounced way. On the other hand, $\delta P(\nu_{\mu} \rightarrow \nu_{\alpha})$ flips sign at around $1000-3000$ km for $\nu_{\mu} \rightarrow \nu_{e}$, and $3000-6000$ km for $\nu_{\mu} \rightarrow \nu_{\tau}$ channels. It produces, assuming detector’s sensitivity, a peculiar zenith angle dependence. The relevant energy region of $\rho E = 50 - 1000 \text{ (g/cm}^3) \text{GeV}$ may be explored, for example, by atmospheric neutrino observation by Deep Core, PINGU, or KM3NeT-ORCA [@Collaboration:2011ym; @TheIceCube-Gen2:2016cap; @Adrian-Martinez:2016zzs].

For the disappearance channel $\nu_\mu \to \nu_\mu$, $|\delta P(\nu_{\mu} \rightarrow \nu_{\mu})| \sim 10^{-3} (10^{-2})$ for $E = 10$ GeV (100 GeV). In this case, the contribution from $\mathcal{C}_{\mu \mu}$ is subdominant compared to $W^2$ correction terms.

Since $\mathcal{C}_{\alpha \beta}$ is a constant term in the oscillation probability, it can in principle be distinguished from the other normalization term which shares $U$ matrix element dependences with the oscillation terms. In particular, they can dominate for large $m_J^2$ since the $W^2$ correction terms are suppressed by at least $\sim 1/m_J^2$. In this case, they will be the sole indicator of low-scale unitarity violation. In general (though not in the $N=1$ model), the order $W^2$ terms depend upon details of the sterile sector, e.g., matrix structure of $W$. Therefore, once the effect is seen it would give us useful information on the structure of low-scale leptonic unitarity violation.

 Some remaining theoretical issues  {#sec:theoretical}
===================================

In this section, we will give some remarks on the theoretical basis in our framework, basic one as well as on its perturbative aspects. They include our treatment of decoherence, generic structure of higher-order corrections and its relation to the “Uniqueness theorem” (see section \[sec:nonunitarity-matter\]), and absence of enhancement due to small solar mass splitting denominator.

 Decoherence imposed onto coherent evolution system  {#sec:decoherence}
----------------------------------------------------

We have started with the Schrödinger equation (\[evolution\]) with Hamiltonian (\[hamiltonian\]) assuming that all the neutrino states remain coherent. We have shown in this and the previous papers that the coherence between active and sterile, and sterile and sterile states are not maintained for sterile mass differences larger than $0.1$ eV$^2$. The effect of decoherence is taken into account by making average over the fast oscillations. We feel it desirable for the current treatment be replaced by the real quantum mechanical one using wave packets, in which the effect of decoherence would automatically come in. But, we do believe that our present framework is able to describe effectively the right physics derived from such improved treatment.

 Smallness of expansion parameters and higher order corrections  {#sec:higher-order}
----------------------------------------------------------------

Here, we discuss general structure of the perturbation series without recourse to averaging out the fast oscillations. The effective expansion parameters in our perturbative framework are the following four, $$\begin{aligned}
\frac{A W }{ \Delta_{J} - h_{i} }, 
\hspace{10mm} 
\frac{A W }{ h_{j} - h_{i} }, 
\hspace{10mm} 
A L W, 
\hspace{6mm} 
\text{and}
\hspace{6mm} 
W. 
\label{expansion-parameters}\end{aligned}$$ We already saw them, except for the last one, in the discussion in section \[sec:energy-denominator\], and it can be seen by inspecting the expressions of the oscillation probabilities up to the fourth orders given in section \[sec:probability-2nd\] and appendix \[sec:expression-probability-4th\]. Formally, the expansion parameter is the first one in (\[expansion-parameters\]) in view of (\[Omega-expand\]) with $\Omega [1]$, the kernel, in (\[Omega-1st-order\]). But, the spacial integration in (\[Omega-expand\]) produces different effective expansion parameters, the second and the third ones in (\[expansion-parameters\]). The extra factor of $W$’s without the kinematical factors is provided when transforming from the $\hat{S}$ to $S$ matrices, as seen in section \[sec:S-matrix\].

For simplicity of the discussion in this section, we limit ourselves to the case of $|W| \sim 0.1$. Under the same conditions we have imposed in section \[sec:energy-denominator\] the first one in (\[expansion-parameters\]) is $\simeq 7.6 \times 10^{-4}$ for $\Delta m^2_{J i} = 0.1$ eV$^2$ and $\rho E = 10 \text{ (g/cm}^3) \text{GeV}$, while the second and the third, which are comparable at around the first oscillation maximum, are estimated as $2.3 \times 10^{-2}$ for the same condition. Therefore, the smallness of the expansion parameter is ensured unless $\rho E \gg 10 \text{ (g/cm}^3) \text{GeV}$. In fact, a close examination of the order $W^4$ terms in the oscillation probability (see appendix \[sec:expression-probability-4th\]) shows that all the formally $W^4$ terms are actually further suppressed. The largest term in the fourth-order oscillation probabilities is of the one suppressed by a factor $\left\vert \left( \frac{A W }{ \Delta_{J} - h_{i} } \right) \left( A L W \right) W^2 \right\vert \lsim 1.7 \times 10^{-7}$, which is as small as $\sim10^{-4}$ even in the case $|W|=0.5$. Therefore, we expect that the formula for the oscillation probability in (\[P-beta-alpha-final\]) works under much relaxed conditions than the one in (\[suppression-cond\]).

 On Uniqueness theorem and matter-dependent dynamical phase  {#sec:U-theorem}
------------------------------------------------------------

We have shown in sections \[sec:probability-2nd\] and \[sec:probability-4th\] that there is no surviving matter dependent correction term in the oscillation probability up to order $W^4$ after averaging out fast oscillations and using the suppression by energy denominators. Should we expect that this feature is stable against higher order corrections beyond order $W^4$, as postulated by the Uniqueness theorem, in our perturbative framework? We argue that the answer is [*Yes*]{}.

We first note that higher-order corrections in terms of $W$ are computed by using $\Omega [1]$ as the kernel, as indicated in eq. (\[Omega-expand\]). Notice also that all the elements of $\Omega [1]$, except for $\Omega [1]_{JJ}$, carry the energy denominator, as shown in (\[Omega-1st-order\]). Then, higher order correction terms are always accompanied by the energy denominators which are composed of some of the first three in (\[expansion-parameters\]), and therefore they are suppressed. The unique exception for it is the terms generated only by $\Omega [1]_{JJ}$ which lacks the energy denominator. Therefore, apart from this special case, we have shown that higher-order corrections in $W$ does not produce the surviving terms after averaging over fast oscillation and using the energy denominator suppression. It is consistent with what we saw in our explicit computation to order $W^4$. This concludes our justification of the Uniqueness theorem.

We need to clear up the issue of special type of perturbative correction terms which involve only $\Omega [1]_{JJ}$ as the kernel in (\[Omega-expand\]). It produces the unique form of $\hat{S}_{JJ}$ as $$\begin{aligned}
\hat{S}_{JJ} = 
e^{ - i \Delta_{J} x } 
\sum_{n} 
\frac{ (-ix)^n }{n !} \left\{ (W^{\dagger} A W)_{JJ} \right\}^n, 
\label{hatS-JJ-term}\end{aligned}$$ a collection of terms of matter-dependent higher order renormalization to $\sum_{J} W_{\alpha J} W_{\beta J}^*$, the probability leaking term at the amplitude level. However, it exponentiates and has contribution to the $S$ matrix element as[^22] $$\begin{aligned}
S_{\alpha \beta} = 
\sum_{J} 
W_{\alpha J} W_{\beta J}^* 
\exp { \left[ -i \left\{ \Delta_{J} + (W^{\dagger} A W )_{JJ} \right\} x \right] }. 
\label{S-alpha-beta-WW-term}\end{aligned}$$ The unique form of $S$ matrix, in principle, raises an interesting issue of dynamically generated phase produced jointly by unitarity violation and matter effect.[^23] In our setting, however, it either disappears from the amplitude squared, or has vanishing effect when the high frequency oscillation is averaged out.

 Absence of enhancement due to small solar mass splitting denominator  {#sec:no-enhancement}
----------------------------------------------------------------------

In perturbation theory one has to sum up intermediate states including off mass shell states. Therefore, even though we sit in the kinematic region where atmospheric-scale oscillations are large the energy denominator can become small, to the order of solar $\Delta m^2$ mass splitting. Then, one might question whether the correction terms blow up at the small denominator, which would invalidate our perturbative treatment.

Fortunately, one can show that the “singularity” which could be produced in the limit of small solar mass splitting always cancels against the small numerator of the similar size. This problem exists already in the second-order expression of the oscillation probability (\[P-beta-alpha-0th+2nd\]). See the second term in second order (in $W$) term. If we denote $h_{l} - h_{k} \equiv \epsilon$ the term would have $1/\epsilon$ singularity in the limit of $\epsilon \rightarrow 0$. However, one can see by inspection by eye that the expression inside the square parenthesis is antisymmetric under $l \leftrightarrow k$, and hence it is of order $\epsilon$ or higher. Therefore, the singularity cancels. Notice that the antisymmetry under $l \leftrightarrow k$ is not required for the whole expression including the matrix element factor.

The situation is a little bit more complicated in the fourth-order expression of the oscillation probability given in appendix \[sec:expression-probability-4th\]. In addition to $1/\epsilon$ singularity similar to the one we already saw, there exist apparent singularity of $1/\epsilon^2$ type. See, for example, the second term in (\[P-beta-alpha-W4-H4-diag\]) and the last term in (\[P-beta-alpha-W4-H4-single\]). But, an explicit calculation shows that the $1/\epsilon^2$ singularity always cancels against order $\epsilon^2$ numerator in the limit of small solar splitting.

This phenomenon is reminiscent of the finiteness of the oscillation probability at the small solar mass splitting limit in helio-perturbation theory with the unique expansion parameter $\frac{\Delta m^2_{21}}{\Delta m^2_{31}}$ (or a renormalized one), see e.g., [@Minakata:2015gra] and the references therein. Possible interpretation of applicability of the perturbative framework to the region of solar level crossing has been discussed [@Xu:2015kma; @Ge:2016dlx]. Another example for the similar phenomena is the one at the small atmospheric mass splitting limit with additional expansion parameter $\sin \theta_{13}$. In this case it is observed that near the atmospheric resonance region not only the oscillation probability is finite but also its accuracy improves when the higher order terms to fourth order in $\sin \theta_{13}$ is added [@Asano:2011nj].

Then, one might ask if our small unitarity violation perturbation theory gives quantitatively accurate result at around the denominator with small solar mass splitting. However, we note that this problem is not relevant in our case because all these terms with apparent singularities vanish after averaging over the high-frequency oscillations and using the suppression by the sterile mass denominators. Yet, we must remark that if we investigate possible enhancement of the correction terms outside the condition (\[suppression-cond\]), as done in section \[sec:correction-terms\], the quantitative accuracy of the expression may become an issue.

Concluding remarks {#sec:conclusion}
==================

In this paper, we have presented a comprehensive treatment of the three active plus $N$ sterile neutrino model in the context of leptonic unitarity test. We have formulated an appropriate perturbative framework with expansion in small unitarity violating $W$ matrix elements, while keeping (non-$W$ suppressed) matter effect to all orders.

What we have done in this paper is mainly threefold:

-   We have shown that the same condition on sterile state masses $0.1\, \text{eV}^2 \lsim m^2_{J} \lsim 1\, \text{MeV}^2$, as we imposed in vacuum, is sufficient to make the $(3+N)$ model sterile-sector model independent, apart from $N$ dependence in the lower bound on probability leaking constant $\mathcal{C}_{\alpha \beta}$. We have shown in an explicit computation to fourth order in $W$ that these higher order correction terms either vanish or are negligible after averaging over fast oscillations and using energy denominator suppression, leaving only the leaking term $\mathcal{C}_{\alpha \beta}$. We have argued by postulating the “Uniqueness theorem” that this feature prevails to all orders in $W$ perturbation theory.

-   As an outcome of our framework, we have derived very simple formulas for oscillation probabilities between active neutrinos in the $(3+N)$ model. We have analyzed them in $\nu_{\mu} \rightarrow \nu_{\alpha}$ channel, $\alpha=e,\mu,\tau$, to know in which region of energy and baseline the effects of unitarity violation are large. We have observed relatively large effects in $\nu_{\mu} \rightarrow \nu_{\mu}$ and $\nu_{\mu} \rightarrow \nu_{\tau}$ channels and pointed out, though qualitatively, that anticorrelation of signals between them would enhance the sensitivity to unitarity violating effects.

-   We have discussed the question of how to distinguish low-scale unitarity violation from high-scale one. We have pointed out that the second order $W$ correction to the leading order, if detected, would signal low-scale unitarity violation. This is to add to the way we have discussed in [@Fong:2016yyh], i.e., to detect the probability leaking term $\mathcal{C}_{\alpha \beta}$ which testifies for the existence of “sterile”, or undetectable but communicable, sector at low energy scales.

Notice that in “constraining mode” of unitarity violation, the model-independence of the framework translates into a universal nature of the bounds, thereby making them more powerful. Whereas in “discovery mode” of unitarity violation, the model dependence, in particular through the $W$ dependent correction terms, is welcome because it serves for identifying the structure of the sterile sector.

During the course of this work, we have obtained the new results and had some interesting observations including:

-   We have obtained an exact solution of the Hamiltonian (\[H0-def\]) for uniform matter density. It describes neutrino evolution in low-scale unitarity violation in zeroth-order in $W$, which applies also to the case of high-scale unitarity violation. The key parts of the solution can also be utilized to calculate the higher order corrections in $W$. When applied to each shell inside the earth, it could provide a semi-quantitative way of simulating non-unitary neutrino evolution for the terrestrial experiments.

-   We have observed outside of the region of validity of our Uniqueness theorem there is a limited phase space, $\rho E \gg 10 \text{ (g/cm}^3) \text{GeV}$, in which the second order $W$ correction terms can be sizeable. If detected, it may allow us to probe structure of $W$ matrix elements which bridges between active and sterile sectors.

-   We have examined a toy model of equally distributed $W$ matrix elements within each flavor, in which the probability leaking term scales as $1/N$. Likewise, possible detection of $\mathcal{C}_{\alpha \beta}$ may provide information of the hidden sterile sector.

We emphasize that, unlike in the case of high-scale unitarity violation, neutrino experiment is the most powerful way to execute leptonic unitarity test in the scenarios of low-scale unitarity violation. Nonetheless, we have to admit that our observations on what we could do for experimental detection of possible non-unitarity effects are too qualitative to make any definitive claim for possible detection in the future. Clearly, the better analyses are called for.

While we worked exclusively on the $(3+N)$ state unitary model as a model of low-scale unitarity violation, we do not know if it is the unique choice, or it merely reflects our ignorance. Even in the case there exist more generic class of models for low-scale unitarity violation, the phenomenon of probability leaking is likely to survive. It is because the probability leaking must take place whenever the extra light sector exists and communicates with the three active neutrinos.

One of the authors (H.M.) thanks Enrique Fernandez-Martinez for interesting discussions about the relationship between high-scale and low-scale unitarity violation. He expresses a deep gratitude to Instituto Física Teórica, UAM/CSIC in Madrid, for its support via “Theoretical challenges of new high energy, astro and cosmo experimental data” project, Ref: 201650E082. He had been a member of Yachay Tech for 14 months, at that time the first research-oriented university in Ecuador [@Yachay-story], during which he was warmly supported by Ecuadorian people. He thanks kind supports by ICTP-SAIFR (FAPESP grant 2016/01343-7), UNICAMP (FAPESP grant 2014/19164-6), and PUC-Rio (CNPq) which enabled him to visit these institutions where part of this work was done. C.S.F. is supported by FAPESP under grants 2013/01792-8 and 2012/10995-7. H.N. is supported by CNPq. This work was supported in part by the Fermilab Neutrino Physics Center.

$\hat{S}$ matrix elements {#sec:hatS-elements}
=========================

The method of computing $\hat{S}$ matrix elements is outlined in section \[sec:hatS-matrix\]. In this appendix \[sec:hatS-elements\] we carry out this task. With the expression of $H_{1}$ in (\[H1-matrix\]) we compute perturbatively the matrix elements of $\Omega (x)$. Then, $\hat{S} (x) = e^{- i \hat{H}_{0} x} \Omega(x)$.

We denote computated results of the matrix elements of $\hat{S}$ as $\hat{S} [n]$ to indicate that it is the one that comes from $n$-th order contribution in $H_{1}$. Since the elements of $H_{1}$ are of order either $W$ or $W^2$, $\hat{S} [n]$ generally has order $W^n$ or higher. To show that a particular contribution is of order $W^m$ we use the superscript $``(m)''$. That is, $\hat{S}^{(m)} [n]$ denotes contribution to $\hat{S}$ that arizes from $n$-th order perturbative contribution in $H_{1}$ and is of order $W^m$.

Contribution to $\hat{S}$ matrix elements from zeroth and first order in $H_{1}$ {#sec:hatS-0th-1st}
--------------------------------------------------------------------------------

The zeroth and first order $\hat{S}$ matrix elements can be calculated as follows: $$\begin{aligned}
\hat{S}_{ii}^{(0+1)} &=& 
\left( e^{-i \hat{H}_{0} x} \right)_{i k} (\Omega_{k i}) +
\left( e^{-i \hat{H}_{0} x} \right)_{i K} (\Omega_{K i}) = 
e^{-i h_{i} x} (\Omega_{i i}) = e^{-i h_{i} x}, 
\nonumber \\
\hat{S}_{i j}^{(1)} \vert_{i \neq j} &=& 
\left( e^{-i \hat{H}_{0} x} \right)_{i k} (\Omega_{k j}) +
\left( e^{-i \hat{H}_{0} x} \right)_{i K} (\Omega_{K j}) = 
e^{-i h_{i} x} (\Omega_{i j}) = 0, 
\nonumber \\
\hat{S}_{i J}^{(1)} &=& 
\left( e^{-i \hat{H}_{0} x} \right)_{i j} (\Omega_{j J}) +
\left( e^{-i \hat{H}_{0} x} \right)_{i K} (\Omega_{K J}) = 
e^{- i h_{i} x} (\Omega_{i J}) = 
\frac{e^{- i \Delta_{J} x} - e^{- i h_{i} x} }{ ( \Delta_{J} - h_{i} )  }
\left\{ (UX)^{\dagger} A W \right\}_{i J}, 
\nonumber \\
\hat{S}_{J i}^{(1)} &=& 
\left( e^{-i \hat{H}_{0} x} \right)_{J k} (\Omega_{k i}) +
\left( e^{-i \hat{H}_{0} x} \right)_{J K} (\Omega_{K i}) = 
\frac{e^{- i \Delta_{J} x} - e^{ - i h_{i} x} }{ ( \Delta_{J} - h_{i} ) }
\left\{ W ^{\dagger} A (UX) \right\}_{J i}, 
\nonumber \\ 
\hat{S}_{J K}^{(1)} \vert_{J \neq K} &=& 
\left( e^{-i \hat{H}_{0} x} \right)_{J i} (\Omega_{i K}) +
\sum_{I} \left( e^{-i \hat{H}_{0} x} \right)_{J I} (\Omega_{I K}) 
= \frac{e^{- i \Delta_{K} x } - e^{- i \Delta_{J} x} }{ ( \Delta_{K} - \Delta_{J} )  }
\left( W ^{\dagger} A W \right)_{J K}, 
\nonumber \\ 
\hat{S}_{J J}^{(0+1)} &=& 
\sum_{i} \left( e^{-i \hat{H}_{0} x} \right)_{J i} (\Omega_{i J}) +
\sum_{I} \left( e^{-i \hat{H}_{0} x} \right)_{J I} (\Omega_{I J}) = 
e^{-i \Delta_{J} x} 
\left[
1 - (i x) \left( W ^{\dagger} A W \right)_{J J} 
\right].
\label{hat-S-elements-1st}\end{aligned}$$

Contribution to $\hat{S}$ matrix elements from second order in $H_{1}$ {#sec:hatS-2nd}
----------------------------------------------------------------------

Likewise, $\hat{S}$ matrix elements can be calculated in second order in $\hat{H}_{1}$ by using the formula for $\Omega$ in (\[Omega-expand\]) and $\hat{S}$-$\Omega$ relation in (\[hatS-Omega\]) as $$\begin{aligned}
\hat{S}_{i i}^{(2)} [2] &=& 
\sum_{K} \left[
(ix) e^{- i h_{i} x} + \frac{e^{- i \Delta_{K} x} - e^{- i h_{i} x} }{ ( \Delta_{K} - h_{i} )  } 
\right]
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\frac{ 1 }{ \Delta_{K} - h_{i} } 
\left\{ W ^{\dagger} A (UX) \right\}_{K i}, 
\nonumber \\ 
\hat{S}_{i j}^{(2)} \vert_{i \neq j} [2] &=& 
- \sum_{K} \left[ 
\frac{e^{- i h_{j} x} - e^{- i h_{i} x}  }{ ( h_{j}  - h_{i} ) }
- \frac{e^{- i \Delta_{K} x} - e^{- i h_{i} x}  }{ ( \Delta_{K} - h_{i} ) } 
\right]
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\frac{ 1 }{ \Delta_{K} - h_{j} } 
\left\{ W ^{\dagger} A (UX) \right\}_{K j}, 
\nonumber \\
\hat{S}_{i J}^{(3)} [2] &=& 
- \left[
(ix) e^{- i \Delta_{J} x} + 
\frac{ e^{- i \Delta_{J} x} - e^{- i h_{i} x} }{ ( \Delta_{J} - h_{i} ) } 
\right]
\left\{ (UX)^{\dagger} A W \right\}_{i J} 
\frac{ 1 }{ \Delta_{J} - h_{i} } 
\left\{ W^{\dagger} A W \right\}_{J J} 
\nonumber \\
&+& 
\sum_{K \neq J}
\left[
\frac{e^{- i \Delta_{J} x} - e^{- i h_{i} x} }{ ( \Delta_{J} - h_{i} ) } - 
\frac{e^{- i \Delta_{K} x} - e^{- i h_{i} x} }{ ( \Delta_{K} - h_{i} ) } 
\right]
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\frac{ 1 }{ \Delta_{J} - \Delta_{K} } 
\left\{ W^{\dagger} A W \right\}_{K J}, 
\nonumber \\ 
\hat{S}_{J i}^{(3)} [2] &=& 
- \left[
(ix) e^{- i \Delta_{J} x} + 
\frac{ e^{- i \Delta_{J} x} - e^{ - i h_{i} x} }{ ( \Delta_{J} - h_{i} ) } 
\right]
\left\{ W^{\dagger} A W \right\}_{J J} 
\frac{ 1 }{ \Delta_{J} - h_{i} } 
\left\{ W^{\dagger} A (UX) \right\}_{J i} 
\nonumber \\ 
&-& 
\sum_{K \neq J}
\left[
\frac{e^{- i h_{i} x} - e^{- i \Delta_{J} x} }{ h_{i} - \Delta_{J} } - 
\frac{e^{- i \Delta_{K} x} - e^{- i \Delta_{J} x} }{ ( \Delta_{K} - \Delta_{J} ) } 
\right]
\left\{ W^{\dagger} A W \right\}_{J K} 
\frac{ 1 }{ \Delta_{K} - h_{i} } 
\left\{ W^{\dagger} A (UX) \right\}_{K i}, 
\nonumber \\ 
\hat{S}_{I J}^{(2+4)} \vert_{I \neq J} [2] &=& 
\sum_{k}
\left[\frac{e^{- i \Delta_{J} x} - e^{- i \Delta_{I} x} }{ ( \Delta_{J} - \Delta_{I} ) } - 
\frac{e^{- i h_{k} x} - e^{- i \Delta_{I} x} }{ ( h_{k} - \Delta_{I} ) } 
\right]
\left\{ W^{\dagger} A (UX) \right\}_{I k} 
\frac{ 1 }{ \Delta_{J} - h_{k} } 
\left\{ (UX)^{\dagger} A W \right\}_{k J} 
\nonumber \\ 
&+& 
\left[ (ix) e^{- i \Delta_{I} x} + 
\frac{e^{- i \Delta_{J} x} - e^{- i \Delta_{I} x} }{ ( \Delta_{J} - \Delta_{I} ) } 
\right]
\left\{ W^{\dagger} A W \right\}_{I I} 
\frac{ 1 }{ \Delta_{J} - \Delta_{I} } 
\left\{ W^{\dagger} A W \right\}_{I J} 
\nonumber \\ 
&-& 
\left[ (ix) e^{- i \Delta_{J} x} + 
\frac{e^{- i \Delta_{J} x} - e^{- i \Delta_{I} x} }{ ( \Delta_{J} - \Delta_{I} ) } 
\right]
\left\{ W^{\dagger} A W \right\}_{I J} 
\frac{ 1 }{ \Delta_{J} - \Delta_{I} } 
\left\{ W^{\dagger} A W \right\}_{J J} 
\nonumber \\ 
&+& 
\sum_{K \neq I, J} 
\left[ \frac{e^{- i \Delta_{J} x} - e^{- i \Delta_{I} x} }{ ( \Delta_{J} - \Delta_{I} ) } - 
\frac{e^{- i \Delta_{K} x} - e^{- i \Delta_{I} x} }{ ( \Delta_{K} - \Delta_{I} ) }
\right]
\left\{ W^{\dagger} A W \right\}_{I K} 
\frac{ 1 }{ \Delta_{J} - \Delta_{K} } 
\left\{ W^{\dagger} A W \right\}_{K J}, 
\nonumber \\
\hat{S}_{I I}^{(2+4)} [2] &=& 
- \sum_{k} 
\left[ (ix) e^{- i \Delta_{I} x} + 
\frac{e^{- i h_{k} x} - e^{- i \Delta_{I} x} }{ h_{k} - \Delta_{I} } 
\right]
\left\{ W^{\dagger} A (UX) \right\}_{I k} 
\frac{ 1 }{ \Delta_{I} - h_{k} } 
\left\{ (UX)^{\dagger} A W \right\}_{k I} 
\nonumber \\ 
&+& \frac{ (- ix)^2 }{ 2 } e^{- i \Delta_{I} x} 
\left( \left\{ W^{\dagger} A W \right\}_{I I} \right)^2 
\nonumber \\
&-&
\sum_{K \neq I} 
\left[ (ix) e^{- i \Delta_{I} x} + 
\frac{e^{- i \Delta_{K} x} - e^{- i \Delta_{I} x}  }{ \Delta_{K} - \Delta_{I} } 
\right]
\left\{ W^{\dagger} A W \right\}_{I K} 
\frac{ 1 }{ \Delta_{I} - \Delta_{K} } 
\left\{ W^{\dagger} A W \right\}_{K I}.
\label{hatS-2nd-order}\end{aligned}$$ We now discuss generalized T invariance of second order $\hat{S}$ matrix elements. Let us first examine $\hat{S}_{i j} \vert_{i \neq j} [2]$. We first note that the matrix element transforms under generalized T transformation $i \leftrightarrow j$, $U \rightarrow U^*$, $W \rightarrow W^*$, etc. as $$\begin{aligned}
\left\{ (UX)^{\dagger} A W \right\}_{i K} &\rightarrow & 
\left\{ W^{\dagger} A (UX) \right\}_{K j}, 
\nonumber \\
\left\{ W^{\dagger} A (UX) \right\}_{K j} &\rightarrow & 
\left\{ (UX)^{\dagger} A W \right\}_{i K}.
\label{matrix-T-transf}\end{aligned}$$ Therefore, the matrix element part is consistent with generalized T invariance. However, the rest of $\hat{S}_{i j} [2] \vert_{i \neq j}$ does not appear to be manifestly generalized T invariant. But, one notices that $\hat{S}_{i j} \vert_{i \neq j} [2]$ can be written as $$\begin{aligned}
&& \hat{S}_{i j}^{(2)} \vert_{i \neq j} [2] = 
- \sum_{K} 
\frac{ 1 }{ ( h_{j}  - h_{i} ) (\Delta_{K} - h_{i}) (\Delta_{K} - h_{j}) } 
\nonumber \\
&\times& 
\left[ 
\left( e^{- i h_{j} x} - e^{- i h_{i} x} \right) \Delta_{K} 
- ( h_{j}  - h_{i} )  e^{- i \Delta_{K} x} 
+ h_{j} e^{- i h_{i} x} - h_{i} e^{- i h_{j} x} 
\right]
\nonumber \\
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K j}.
\label{hatS-ij[2]}\end{aligned}$$ Then, $\hat{S}_{i j} \vert_{i \neq j} [2] (U, W, \text{etc.} )= \hat{S}_{j i} \vert_{i \neq j} [2] (U^*, W^*, \text{etc.} )$ holds, that is, generalized T invariance is maintained.

Similarly, the first and the fourth terms in $\hat{S}_{I J} \vert_{I \neq J} [2]$ are not manifestly generalized T invariant. But, they can also be written as manifestly generalized T invariant forms as $$\begin{aligned}
&& \hat{S}_{I J}^{(2+4)} \vert_{I \neq J} [2] = 
\sum_{k} 
\frac{ 1 }{ ( \Delta_{J} - \Delta_{I} ) (\Delta_{J} - h_{k}) (\Delta_{I} - h_{k}) } 
\nonumber \\
&\times& 
\left[ 
\Delta_{J} e^{- i \Delta_{I} x} - \Delta_{I} e^{- i \Delta_{J} x} 
+ \left( e^{- i \Delta_{J} x} - e^{- i \Delta_{I} x} \right) h_{k} 
-  ( \Delta_{J} - \Delta_{I} ) e^{- i h_{k} x} 
\right]
\nonumber \\
&\times&
\left\{ W^{\dagger} A (UX) \right\}_{I k} 
\left\{ (UX)^{\dagger} A W \right\}_{k J} 
\nonumber \\ 
&+& 
(ix) e^{- i \Delta_{I} x} \frac{ 1 }{ \Delta_{J} - \Delta_{I} } 
\left\{ W^{\dagger} A W \right\}_{I I} 
\left\{ W^{\dagger} A W \right\}_{I J} 
\nonumber \\ 
&-& 
(ix) e^{- i \Delta_{J} x} \frac{ 1 }{ \Delta_{J} - \Delta_{I} }
\left\{ W^{\dagger} A W \right\}_{I J} 
\left\{ W^{\dagger} A W \right\}_{J J} 
\nonumber \\ 
&+&
\frac{e^{- i \Delta_{J} x} - e^{- i \Delta_{I} x} }{ ( \Delta_{J} - \Delta_{I} )^2 } \left\{ W^{\dagger} A W \right\}_{I I} 
\left\{ W^{\dagger} A W \right\}_{I J} 
\nonumber \\ 
&-&
\frac{e^{- i \Delta_{J} x} - e^{- i \Delta_{I} x} }{ ( \Delta_{J} - \Delta_{I} )^2 } \left\{ W^{\dagger} A W \right\}_{I J} 
\left\{ W^{\dagger} A W \right\}_{J J} 
\nonumber \\ 
&-& 
\sum_{K \neq I, J} 
\frac{ 1 }{ ( \Delta_{J} - \Delta_{I} ) ( \Delta_{K} - \Delta_{I} ) ( \Delta_{K} - \Delta_{J} ) } 
\nonumber \\
&\times& \left[ 
\Delta_{J} e^{- i \Delta_{I} x} - \Delta_{I} e^{- i \Delta_{J} x} 
+ \left( e^{- i \Delta_{J} x} - e^{- i \Delta_{I} x} \right) \Delta_{K}
-  ( \Delta_{J} - \Delta_{I} ) e^{- i \Delta_{K} x} 
\right]
\nonumber \\
&\times&
\left\{ W^{\dagger} A W \right\}_{I K} 
\left\{ W^{\dagger} A W \right\}_{K J}. 
\label{hatS-IJ(2+4)}\end{aligned}$$ Under generalized T transformation, the second term goes to the third term (both $ix$ terms), and vice versa, and hence they are invariant. The situation is completely the same as in the fourth and the fifth terms. The first and the last terms are invariant in an analogous way as $\hat{S}_{i j} \vert_{i \neq j} [2]$.

For $\hat{S}_{i J} [2]$ and $\hat{S}_{J i} [2]$ generalized T invariance holds as they are. For $\hat{S}_{i i} [2]$ and $\hat{S}_{I I} [2]$ the matrix factor is self-invariant, and therefore generalized T invariance holds.

What we should do in the rest of appendix is to compute $\hat{S}$ matrix elements perturbatively to fourth order in $H_{1}$. In presenting the results of computation of $\hat{S}$ matrix elements, however, we change strategy of our description. That is,

-   We present $\hat{S}$ matrix elements at fixed order in $W$. To do this we, of course, have to sum up contributions that come from different order in $H_{1}$. For example, $\hat{S}_{i J}^{(3)} = \hat{S}_{i J}^{(3)} [3] + \hat{S}_{i J}^{(3)} [2]$.

-   We present only the terms which are required to compute $S$ matrix elements to order $W^4$. In view of the relations between $\hat{S}$ and $S$ matrix elements given in eq. (\[Sab-hatSab\]), $\hat{S}_{I J}^{(3)}$, $\hat{S}_{I J}^{(4)}$, and $\hat{S}_{i J}^{(4)}$ (and $\hat{S}_{J i}^{(4)}$) are all unnecessary.

-   We only give the results of manifestly generalized T invariant form of $\hat{S}$ matrix elements with which it must be straightforward to prove generalized T invariance.

Order $W^2$ $\hat{S}$ matrix elements from Second-order in $H_{1}$  {#sec:hatSIJ-2nd}
-------------------------------------------------------------------

The second order in $H_{1}$ contribution $\hat{S}_{I J} [2]$ (and $\hat{S}_{I I} [2]$) has both order $W^2$ and $W^4$ terms. To compute $S$ matrix elements to order $W^4$ we need only the former term because the latter produces order $W^6$ terms, as one can see in (\[Sab-hatSab\]).

Then, it is necessary to remind the readers that the expressions of order $W^2$ terms in $\hat{S}_{I J} [2]$ and $\hat{S}_{I I} [2]$ already exist as the first terms in (\[hatS-IJ(2+4)\]) and (\[hatS-2nd-order\]), respectively.

Order $W^3$ $\hat{S}$ matrix elements $\hat{S}_{i J}^{(3)}$ and $\hat{S}_{J i}^{(3)}$  {#sec:hatSiJ-2nd-]}
--------------------------------------------------------------------------------------

By adding the contribution of the third and second order in $H_{1}$ we obtain for $\hat{S}$ matrix elements $\hat{S}_{i J}^{(3)}$ at order $W^3$ $$\begin{aligned}
&& \hat{S}_{i J}^{(3)} [2+3]
\nonumber \\
&=&
- \frac{ 1 }{ \Delta_{J} - h_{i} } 
\left[ (ix) e^{- i \Delta_{J} x} + \frac{ e^{- i \Delta_{J} x} - e^{- i h_{i} x} }{ ( \Delta_{J} - h_{i} ) } 
\right]
\left\{ (UX)^{\dagger} A W \right\}_{i J} 
\left\{ W^{\dagger} A W \right\}_{J J} 
\nonumber \\
&+& 
\sum_{K \neq J} 
\frac{ 1 }{ ( \Delta_{J} - \Delta_{K} ) ( \Delta_{J} - h_{i} ) ( \Delta_{K} - h_{i} ) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{i} \right) e^{- i \Delta_{J} x} 
- \left( \Delta_{J} - h_{i} \right) e^{- i \Delta_{K} x} 
- \left( \Delta_{K} - \Delta_{J} \right) e^{- i h_{i} x} 
\biggr]
\nonumber \\
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W^{\dagger} A W \right\}_{K J} 
\nonumber \\ 
&-& 
\frac{ 1 }{ ( \Delta_{J} - h_{i} )^2 } 
\biggl[
(ix) \left( e^{- i h_{i} x} + e^{- i \Delta_{J} x} \right)
+ 2 \frac{e^{- i \Delta_{J} x} - e^{- i h_{i} x} }{ ( \Delta_{J} - h_{i} ) } 
\biggr] 
\nonumber \\ 
&\times& 
\left\{ (UX)^{\dagger} A W \right\}_{i J} 
\left\{ W^{\dagger} A (UX) \right\}_{J i} 
\left\{ (UX)^{\dagger} A W \right\}_{i J} 
\nonumber \\ 
&+& 
\sum_{k \neq i} 
\biggl[
- \frac{ (ix) e^{- i \Delta_{J} x} }{ ( \Delta_{J} - h_{i} )( \Delta_{J} - h_{k} ) } 
+ \frac{ 1 }{ ( h_{i} - h_{k} ) (\Delta_{J} - h_{i} )^2 (\Delta_{J} - h_{k} )^2 } 
\nonumber \\ 
&\times& 
\biggl\{
(\Delta_{J} - h_{k} )^2 e^{- i h_{i} x} 
- (\Delta_{J} - h_{i} )^2 e^{- i h_{k} x} 
+ ( h_{i} - h_{k} )( h_{i} + h_{k} - 2 \Delta_{J} ) e^{- i \Delta_{J} x} 
\biggr\} 
\biggr] 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i J} 
\left\{ W^{\dagger} A (UX) \right\}_{J k} 
\left\{ (UX)^{\dagger} A W \right\}_{k J} 
\nonumber \\ 
&+& 
\sum_{K \neq J} 
\biggl[
- \frac{ (ix) e^{- i h_{i} x} }{ (\Delta_{J} - h_{i}) (\Delta_{K} - h_{i}) }   
+ \frac{ 1 }{ (\Delta_{J} - \Delta_{K}) (\Delta_{J} - h_{i})^2 ( \Delta_{K} - h_{i} )^2 }
\nonumber \\ 
&\times&
\biggl\{
( \Delta_{K} - h_{i} )^2 e^{- i \Delta_{J} x} 
- (\Delta_{J} - h_{i})^2 e^{- i \Delta_{K} x} 
+ e^{- i h_{i} x} (\Delta_{J} - \Delta_{K}) (\Delta_{J} + \Delta_{K} - 2 h_{i} )
\biggr\}
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W^{\dagger} A (UX) \right\}_{K i} 
\left\{ (UX)^{\dagger} A W \right\}_{i J} 
\nonumber \\ 
&+& 
\sum_{K \neq J} 
\sum_{k \neq i} 
\frac{ 1 }{ (\Delta_{J} - \Delta_{K}) ( h_{k}  - h_{i} ) (\Delta_{J} - h_{i}) (\Delta_{J} - h_{k}) (\Delta_{K} - h_{i}) (\Delta_{K} - h_{k}) } 
\nonumber \\ 
&\times&
\biggl[
( h_{k}  - h_{i} ) 
\biggl\{(\Delta_{K} - h_{i}) (\Delta_{K} - h_{k}) e^{- i \Delta_{J} x} 
- (\Delta_{J} - h_{i}) (\Delta_{J} - h_{k}) e^{- i \Delta_{K} x}
\biggr\}
\nonumber \\ 
&-&
(\Delta_{J} - \Delta_{K}) 
\biggl\{ (\Delta_{J} - h_{k}) (\Delta_{K} - h_{k}) e^{- i h_{i} x} 
- (\Delta_{J} - h_{i}) (\Delta_{K} - h_{i}) e^{- i h_{k} x}
\biggr\}
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k J}, 
\label{hatS-iJ-W3-H2+3}\end{aligned}$$ and similarly, for $\hat{S}_{J i}^{(3)}$ $$\begin{aligned}
&& \hat{S}_{J i}^{(3)} [2+3]
\nonumber \\ 
&=& 
- \frac{ 1 }{ \Delta_{J} - h_{i} } 
\left[
(ix) e^{- i \Delta_{J} x} + 
\frac{ e^{- i \Delta_{J} x} - e^{ - i h_{i} x} }{ ( \Delta_{J} - h_{i} ) } 
\right]
\left\{ W^{\dagger} A W \right\}_{J J} 
\left\{ W^{\dagger} A (UX) \right\}_{J i} 
\nonumber \\
&+& 
\sum_{K \neq J} 
\frac{ 1 }{ ( \Delta_{J} - \Delta_{K} ) ( \Delta_{J} - h_{i} ) ( \Delta_{K} - h_{i} ) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{i} \right) e^{- i \Delta_{J} x} 
- \left( \Delta_{J} - h_{i} \right) e^{- i \Delta_{K} x} 
- \left( \Delta_{K} - \Delta_{J} \right) e^{- i h_{i} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ W^{\dagger} A W \right\}_{J K}  
\left\{ W^{\dagger} A (UX) \right\}_{K i} 
\nonumber \\ 
&-& 
\frac{ 1 }{ (\Delta_{J} - h_{i})^2 } 
\biggl[
(ix) \left( e^{- i h_{i} x} + e^{- i \Delta_{J} x} \right)
+ 2 \frac{e^{- i \Delta_{J} x} - e^{- i h_{i} x} }{ ( \Delta_{J} - h_{i} ) } 
\biggr]
\nonumber \\ 
&\times& 
\left\{ W^{\dagger} A (UX) \right\}_{J i} 
\left\{ (UX)^{\dagger} A W \right\}_{i J} 
\left\{ W^{\dagger} A (UX) \right\}_{J i} 
\nonumber \\ 
&+& 
\sum_{k \neq i} 
\biggl[
- \frac{ (ix) e^{- i \Delta_{J} x} }{ ( \Delta_{J} - h_{i} )( \Delta_{J} - h_{k} ) } 
+ \frac{ 1 }{ ( h_{i} - h_{k} ) (\Delta_{J} - h_{i} )^2 (\Delta_{J} - h_{k} )^2 } 
\nonumber \\ 
&\times& 
\biggl\{
(\Delta_{J} - h_{k} )^2 e^{- i h_{i} x} 
- (\Delta_{J} - h_{i} )^2 e^{- i h_{k} x} 
+ ( h_{i} - h_{k} )( h_{i} + h_{k} - 2 \Delta_{J} ) e^{- i \Delta_{J} x} 
\biggr\}
\biggr] 
\nonumber \\ 
&\times& 
\left\{ W^{\dagger} A (UX) \right\}_{J k} 
\left\{ (UX)^{\dagger} A W \right\}_{k J} 
\left\{ W^{\dagger} A (UX) \right\}_{J i} 
\nonumber \\ 
&+& 
\sum_{K \neq J} 
\biggl[
- \frac{ (ix) e^{- i h_{i} x} }{ (\Delta_{J} - h_{i}) (\Delta_{K} - h_{i}) }  
+ \frac{ 1 }{ (\Delta_{J} - \Delta_{K}) (\Delta_{J} - h_{i})^2 ( \Delta_{K} - h_{i} )^2 }
\nonumber \\ 
&\times&
\biggl\{
( \Delta_{K} - h_{i} )^2 e^{- i \Delta_{J} x} 
- (\Delta_{J} - h_{i})^2 e^{- i \Delta_{K} x} 
+ e^{- i h_{i} x} (\Delta_{J} - \Delta_{K}) (\Delta_{J} + \Delta_{K} - 2 h_{i} )
\biggr\}
\biggr]
\nonumber \\ 
&\times&
\left\{ W^{\dagger} A (UX) \right\}_{J i} 
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W^{\dagger} A (UX) \right\}_{K i} 
\nonumber \\ 
&+& 
\sum_{K \neq J} 
\sum_{k \neq i} 
\frac{ 1 }{ (\Delta_{J} - \Delta_{K}) ( h_{k}  - h_{i} ) (\Delta_{J} - h_{i}) (\Delta_{J} - h_{k}) (\Delta_{K} - h_{i}) (\Delta_{K} - h_{k}) } 
\nonumber \\ 
&\times&
\biggl[
( h_{k}  - h_{i} ) 
\biggl\{ (\Delta_{K} - h_{i}) (\Delta_{K} - h_{k}) e^{- i \Delta_{J} x} 
- (\Delta_{J} - h_{i}) (\Delta_{J} - h_{k}) e^{- i \Delta_{K} x}
\biggr\}
\nonumber \\ 
&-&
(\Delta_{J} - \Delta_{K}) 
\biggl\{ (\Delta_{J} - h_{k}) (\Delta_{K} - h_{k}) e^{- i h_{i} x} 
- (\Delta_{J} - h_{i}) (\Delta_{K} - h_{i}) e^{- i h_{k} x}
\biggr\}
\biggr]
\nonumber \\ 
&\times&
\left\{ W^{\dagger} A (UX) \right\}_{J k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W^{\dagger} A (UX) \right\}_{K i}. 
\label{hatS-Ji-W3-H2+3}\end{aligned}$$ Using (\[hatS-iJ-W3-H2+3\]) and (\[hatS-Ji-W3-H2+3\]), one can easily confirm generalized T invariance $\hat{S}_{i J}^{(3)} (U, W, \text{etc}) = \hat{S}_{J i}^{(3)} (U^*, W^*, \text{etc})$.

Order $W^4$ $\hat{S}$ matrix elements $\hat{S}_{i j}^{(4)} \vert_{i \neq j} [3]$ and $\hat{S}_{i j}^{(4)} \vert_{i \neq j} [4]$  {#sec:hatSij-3nd-4th}
--------------------------------------------------------------------------------------------------------------------------------

The third order in $H_{1}$ contribution to order $W^4$ $\hat{S}$ matrix elements $\hat{S}_{i j}^{(4)} \vert_{i \neq j}$ is given by $$\begin{aligned}
&& \hat{S}_{ij}^{(4)} \vert_{i \neq j} [3] 
\nonumber \\
&=& 
- \sum_{K} 
(ix) e^{- i \Delta_{K} x} 
\frac{ 1 }{ \Delta_{K} - h_{i} } \frac{ 1 }{ \Delta_{K} - h_{j} } 
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A W \right\}_{K K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K j} 
\nonumber \\ 
&+& 
\sum_{K} 
\frac{ 1 }{  ( h_{j}  - h_{i} ) ( \Delta_{K} - h_{i} )^2  ( \Delta_{K} - h_{j} )^2 } 
\nonumber \\ 
&\times& 
\biggl[
( h_{j} - h_{i} ) ( h_{j} + h_{i} - 2  \Delta_{K} ) e^{- i \Delta_{K} x} + 
( \Delta_{K} - h_{i} )^2 e^{ - i h_{j} x} - 
( \Delta_{K} - h_{j} )^2 e^{ - i h_{i} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A W \right\}_{K K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K j} 
\nonumber \\ 
&+& 
\sum_{K \neq L} 
\frac{ 1 }{ ( h_{j}  - h_{i} ) ( \Delta_{L} - \Delta_{K} ) ( \Delta_{K} - h_{i} ) ( \Delta_{K} - h_{j} ) ( \Delta_{L} - h_{i} ) ( \Delta_{L} - h_{j} )   } 
\nonumber \\ 
&\times& 
\biggl[ 
( h_{j}  - h_{i} ) 
\biggl\{ ( \Delta_{K} - h_{i} ) ( \Delta_{K} - h_{j} ) e^{- i \Delta_{L} x} - 
( \Delta_{L} - h_{i} ) ( \Delta_{L} - h_{j} ) e^{- i \Delta_{K} x} \biggr\} 
\nonumber \\ 
&+&
( \Delta_{L} - \Delta_{K} ) 
\biggl\{ ( \Delta_{K} - h_{i} ) ( \Delta_{L} - h_{i} ) e^{- i h_{j} x} - 
( \Delta_{K} - h_{j} ) ( \Delta_{L} - h_{j} ) e^{- i h_{i} x} \biggr\} 
\biggl]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A W \right\}_{K L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L j}. 
\label{hatS-3rd-order-ij-3}\end{aligned}$$ For bookkeeping purpose we decompose the fourth order in $H_{1}$ contribution to order $W^4$ $\hat{S}$ matrix elements $\hat{S}_{i j}^{(4)} \vert_{i \neq j}$ into the two terms $$\begin{aligned}
\hat{S}_{ij}^{(4)} [4] \vert_{i \neq j} = 
\hat{S}_{ij}^{(4)} [4] \vert_{i \neq j}~(\text{single sum}) + 
\hat{S}_{ij}^{(4)} [4] \vert_{i \neq j}~(\text{double sum})
\label{hatS-4th-order-ij-T-transf}\end{aligned}$$ where $$\begin{aligned}
&& \hat{S}_{ij}^{(4)} [4] \vert_{i \neq j}~(\text{single sum})
\nonumber \\
&=& 
\sum_{K}
\biggl[
(ix) e^{ - i h_{i} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( h_{j} - h_{i} ) } 
- (ix) e^{ - i \Delta_{K} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{K} - h_{j} ) } 
\nonumber \\
&+&
\frac{ 1 }{ ( \Delta_{K} - h_{j} )^2 ( h_{j} - h_{i} )^2 } 
e^{ - i h_{j} x}
\nonumber \\
&-&
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^3 ( h_{j} - h_{i} )^2 } 
\left( \Delta_{K} + 2 h_{j} - 3 h_{i} \right) 
e^{ - i h_{i} x} 
\nonumber \\
&+&
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^3 ( \Delta_{K} - h_{j} )^2 } 
\left( h_{i} + 2 h_{j}  - 3 \Delta_{K} \right) 
e^{ - i \Delta_{K} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K i} 
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K j} 
\nonumber \\
&+& 
\sum_{K}
\biggl[
- (ix) e^{ - i h_{j} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{j} )^2 ( h_{j} - h_{i} ) } 
- (ix) e^{ - i \Delta_{K} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{j} )^2 ( \Delta_{K} - h_{i} ) } 
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{K} - h_{j} )^3 } 
\left( h_{j} + 2 h_{i} - 3 \Delta_{K} \right) 
e^{ - i \Delta_{K} x} 
\nonumber \\
&-& 
\frac{ 1 }{ ( \Delta_{K} - h_{j} )^3 ( h_{j} - h_{i} )^2 } 
\left( \Delta_{K} + 2 h_{i} - 3 h_{j}  \right) 
e^{ - i h_{j} x} 
+ \frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( h_{j} - h_{i} )^2 } 
e^{ - i h_{i} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K j} 
\left\{ (UX)^{\dagger} A W \right\}_{j K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K j} 
\nonumber \\ 
&+& 
\sum_{K} \sum_{k \neq i, j} 
\biggl[
- (ix) e^{- i \Delta_{K} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} ) ( \Delta_{K} - h_{j} ) ( \Delta_{K} - h_{k} ) } 
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{K} - h_{j} )^2 ( \Delta_{K} - h_{k} ) } 
\left( h_{i} + h_{j} - 2 \Delta_{K}  \right)
e^{- i \Delta_{K} x} 
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{j} )^2 ( \Delta_{K} - h_{k} ) ( h_{j} - h_{i} ) } 
e^{ - i h_{j} x} 
- 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{K} - h_{k} ) ( h_{j} - h_{i} ) } 
e^{ - i h_{i} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K j}, 
\label{hatS-ij-W4-H4-single}\end{aligned}$$ and $$\begin{aligned}
&& \hat{S}_{ij}^{(4)} [4] \vert_{i \neq j}~(\text{double sum})
\nonumber \\
&=& 
\sum_{K \neq L} 
\biggl[
(ix) e^{- i h_{i} x} 
\frac{ 1 }{ (\Delta_{K} - h_{i} ) (\Delta_{L} - h_{i} ) ( h_{j}  - h_{i} ) } 
\nonumber \\
&-& 
\frac{ 1 }{ (\Delta_{K} - \Delta_{L} ) (\Delta_{L} - h_{i} )^2 (\Delta_{L} - h_{j} ) } 
e^{- i \Delta_{L} x} 
+ \frac{ 1 }{ (\Delta_{K} - \Delta_{L} ) (\Delta_{K} - h_{i} )^2 (\Delta_{K} - h_{j} ) } 
e^{- i \Delta_{K} x}
\nonumber \\
&+& 
\frac{ 1 }{ (\Delta_{K} - h_{j} ) (\Delta_{L} - h_{j} ) ( h_{j}  - h_{i} )^2 } e^{- i h_{j} x} 
\nonumber \\
&-& 
\frac{ 1 }{ (\Delta_{K} - h_{i} )^2 (\Delta_{L} - h_{i} )^2 ( h_{j}  - h_{i} )^2 } 
\biggl\{
3 h_{i}^2 - 2 h_{i} h_{j} + \left( h_{j} - 2 h_{i} \right)  (\Delta_{K} + \Delta_{L} ) + \Delta_{K} \Delta_{L} 
\biggr\}
e^{- i h_{i} x} 
\biggr] 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K i} 
\left\{ (UX)^{\dagger} A W \right\}_{i L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L j} 
\nonumber \\
&+&
\sum_{K \neq L} 
\biggl[
- (ix) e^{- i h_{j} x} 
\frac{ 1 }{ (\Delta_{K} - h_{j} ) (\Delta_{L} - h_{j} ) ( h_{j}  - h_{i} ) } 
\nonumber \\
&+& 
\frac{ 1 }{ (\Delta_{K} - \Delta_{L} ) (\Delta_{K} - h_{i} ) (\Delta_{K} - h_{j} )^2 } 
e^{- i \Delta_{K} x}
- \frac{ 1 }{ (\Delta_{K} - \Delta_{L} ) (\Delta_{L} - h_{i} ) (\Delta_{L} - h_{j} )^2 } 
e^{- i \Delta_{L} x} 
\nonumber \\
&+&
\frac{ 1 }{ (\Delta_{K} - h_{i} ) (\Delta_{L} - h_{i} ) ( h_{j}  - h_{i} )^2 } e^{- i h_{i} x} 
\nonumber \\
&-& 
\frac{ 1 }{ (\Delta_{K} - h_{j} )^2 (\Delta_{L} - h_{j} )^2 ( h_{j}  - h_{i} )^2 } 
\biggl\{
3 h_{j}^2 - 2 h_{i} h_{j} - \left( 2 h_{j} - h_{i} \right)  (\Delta_{K} + \Delta_{L} ) + \Delta_{K} \Delta_{L} 
\biggr\}
e^{- i h_{j} x} 
\biggr] 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K j} 
\left\{ (UX)^{\dagger} A W \right\}_{j L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L j} 
\nonumber \\
&+&
\sum_{K \neq L} \sum_{k \neq i, j} 
\biggl[ 
\frac{1}{ ( h_{k} - h_{i} ) ( h_{k} - h_{j} ) ( h_{j} - h_{i} ) ( \Delta_{K} - h_{i} ) ( \Delta_{K} - h_{j} ) ( \Delta_{L} - h_{i} ) ( \Delta_{L} - h_{j} ) } 
\nonumber \\
&+& 
\biggl\{ 
( h_{k} - h_{j} ) (\Delta_{K} - h_{j} ) ( \Delta_{L} - h_{j} ) e^{ - i h_{i} x} 
-
( h_{k} - h_{i} ) (\Delta_{K} - h_{i} ) ( \Delta_{L} - h_{i} ) e^{ - i h_{j} x} 
\biggr\} 
\nonumber \\
&+& 
\frac{1}{ ( h_{k} - h_{i} ) ( h_{k} - h_{j} ) ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{k} ) } 
e^{ - i h_{k} x}
\nonumber \\
&+& 
\frac{1}{ ( \Delta_{K} - h_{i} ) ( \Delta_{K} - h_{j} ) ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{i} ) ( \Delta_{L} - h_{j} ) ( \Delta_{L} - h_{k} ) ( \Delta_{L} - \Delta_{K} ) } 
\nonumber \\ 
&\times&
\biggl\{  ( \Delta_{K} - h_{i} ) ( \Delta_{K} - h_{j} ) ( \Delta_{K} - h_{k} ) e^{- i \Delta_{L} x} 
- ( \Delta_{L} - h_{i} ) ( \Delta_{L} - h_{j} ) ( \Delta_{L} - h_{k} ) e^{- i \Delta_{K} x}
\biggr\}
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L j}. 
\label{hatS-ij-W4-H4-double}\end{aligned}$$

Order $W^4$ $\hat{S}$ matrix elements $\hat{S}_{i i}^{(4)} [3]$ and $\hat{S}_{i i}^{(4)} [4]$  {#sec:hatSii-3nd-4th}
----------------------------------------------------------------------------------------------

$\hat{S} _{i i}^{(4)} [3]$ is given by $$\begin{aligned}
&& \hat{S}_{ii}^{(4)} [3]
\nonumber \\
&=& 
- \sum_{K} 
\biggl[ 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 } 
\{ (ix) e^{- i \Delta_{K} x} + (ix) e^{- i h_{i} x} \}
+ 
\frac{ 2 }{ ( \Delta_{K} - h_{i} )^3 } 
\left( e^{- i \Delta_{K} x} - e^{- i h_{i} x} \right)
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A W \right\}_{K K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K i} 
\nonumber \\ 
&+& 
\sum_{K \neq L} 
\biggl[ 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{K} - \Delta_{L} )} 
e^{- i \Delta_{K} x} 
-
\frac{ 1 }{ ( \Delta_{L} - h_{i} )^2 ( \Delta_{K} - \Delta_{L} )} 
e^{- i \Delta_{L} x} 
\nonumber \\ 
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{L} - h_{i} )^2 }  
\left( \Delta_{K} + \Delta_{L} - 2 h_{i} \right) 
e^{- i h_{i} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A W \right\}_{K L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L i}. 
\label{hatS-3rd-order-ii-T-transf}\end{aligned}$$ $\hat{S} _{i i}^{(4)} [4]$ is given by $$\begin{aligned}
&& \hat{S} _{i i}^{(4)} [4] 
\nonumber \\
&=& 
\sum_{K} 
\biggl[
- \frac{x^2}{2} e^{- i h_{i} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 } 
- 
(ix) e^{- i h_{i} x} 
\frac{ 2 }{ ( \Delta_{K} - h_{i} )^3 } 
- 
(ix) e^{- i \Delta_{K} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^3 } 
\nonumber \\
&-&
\frac{ 3 }{ ( \Delta_{K} - h_{i} )^4 } 
\left( e^{- i \Delta_{K} x} - e^{- i h_{i} x} \right)
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K i} 
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K i} 
\nonumber \\
&+& 
\sum_{K} \sum_{k \neq i} 
\biggl[
(ix) e^{- i h_{i} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( h_{k}  - h_{i} ) } 
- 
(ix) e^{- i \Delta_{K} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{K} - h_{k} ) } 
\nonumber \\
&+& 
\frac{ ( h_{i} + 2 h_{k} - 3 \Delta_{K} ) }{ ( \Delta_{K} - h_{i} )^3 ( \Delta_{K} - h_{k} )^2 } 
e^{- i \Delta_{K} x} 
+ 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 ( h_{k}  - h_{i} )^2 } 
e^{- i h_{k} x} 
\nonumber \\
&-&
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^3 ( h_{k}  - h_{i} )^2 } 
\left( \Delta_{K} + 2 h_{k} - 3 h_{i} \right)
e^{- i h_{i} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K i} 
\nonumber \\
&+& 
\sum_{K \neq L} 
\biggl[
- \frac{x^2}{2} e^{- i h_{i} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} ) ( \Delta_{L} - h_{i} ) } 
\nonumber \\
&-& (ix) e^{- i h_{i} x} 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{L} - h_{i} )^2 }
\left( \Delta_{K} + \Delta_{L} - 2 h_{i} \right) 
\nonumber \\
&-& 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^3 ( \Delta_{L} - \Delta_{K} ) } 
e^{- i \Delta_{K} x} 
+ \frac{ 1 }{ ( \Delta_{L} - h_{i} )^3 ( \Delta_{L} - \Delta_{K} ) } 
e^{- i \Delta_{L} x} 
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{i} )^3 ( \Delta_{L} - h_{i} )^3 } 
\biggl\{
\Delta_{L}^2 + \Delta_{L} \Delta_{K} + \Delta_{K}^2 - 3 h_{i} ( \Delta_{L} + \Delta_{K} ) + 3 h_{i}^2  
\biggr\}
e^{- i h_{i} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K i} 
\left\{ (UX)^{\dagger} A W \right\}_{i L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L i} 
\nonumber \\
&+& 
\sum_{K \neq L} \sum_{k \neq i} 
\biggl[
\frac{1}{ ( \Delta_{K} - h_{i} ) ( \Delta_{L} - h_{i} ) ( h_{k} - h_{i} ) }
(ix) e^{- i h_{i} x} 
\nonumber \\
&-& 
\frac{1}{ ( \Delta_{K} - h_{i} )^2 ( \Delta_{L} - h_{i} )^2 ( h_{k} - h_{i} )^2 } 
\biggl\{ 
\Delta_{K} \Delta_{L} + ( h_{k} - 2 h_{i} ) ( \Delta_{K} + \Delta_{L} ) + 3 
h_{i}^2 - 2 h_{k} h_{i}
\biggr\}
e^{- i h_{i} x} 
\nonumber \\
&+&
\frac{1}{ ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{k} ) ( h_{k} - h_{i} )^2 } e^{- i h_{k} x} 
\nonumber \\
&+& 
\frac{1}{ ( \Delta_{K} - \Delta_{L} ) ( \Delta_{K} - h_{i} )^2 ( \Delta_{K} - h_{k} ) } e^{- i \Delta_{K} x}
- \frac{1}{ ( \Delta_{K} - \Delta_{L} ) ( \Delta_{L} - h_{i} )^2 ( \Delta_{L} - h_{k} ) } e^{- i \Delta_{L} x} 
\biggr]
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{i K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L i}. 
\label{hatS-4th-order-ii-4th}\end{aligned}$$

Structure of $S$ matrix elements {#sec:structure-S-matrix}
================================

Here we present details in computation of the $S$ matrix elements.

The $S$ matrix elements $S_{\alpha \beta}^{(4)}$
------------------------------------------------

We decompose $S_{\alpha \beta}^{(4)}$ into the following three pieces (include both $\alpha \neq \beta$ and $\alpha = \beta$) $$\begin{aligned}
S_{\alpha \beta}^{(4)} &=& S_{\alpha \beta}^{(4)} [3+4] + S_{\alpha \beta}^{(4)} [3] + S_{\alpha \beta}^{(4)} [2] 
\label{Sab-4th}\end{aligned}$$ where $[n]$ implies that the term comes from $n$-th order perturbation of $H_{1}$. To prevent too long expression, we decompose the first term in (\[Sab-4th\]) as $$\begin{aligned}
S_{\alpha \beta}^{(4)} [3+4] 
&=& 
S_{\alpha \beta}^{(4)} [3]_{ \text{ diag } } + S_{\alpha \beta}^{(4)} [4]_{ \text{ diag } } + 
S_{\alpha \beta}^{(4)} [3]_{ \text{ offdiag } } + S_{\alpha \beta}^{(4)} [4]_{ \text{ offdiag } } 
\label{Sab-4th-3+4}\end{aligned}$$ where ($n = 3, 4$) “diag” and “offdiag”, respectively, implies $$\begin{aligned}
S_{\alpha \beta}^{(4)} [n]_{ \text{ diag } } 
&=& 
\sum_{k} (UX)_{ik} 
\left( \hat{S}_{kk}^{(4)} [n] \right)
\left\{ (UX)^{\dagger} \right\}_{kj}, 
\nonumber \\
S_{\alpha \beta}^{(4)} [n]_{ \text{ offdiag } }
&=& 
\sum_{k \neq l} (UX)_{ik} 
\left( \hat{S}_{kl}^{(4)} [n] \right)
\left\{ (UX)^{\dagger} \right\}_{lj}. \end{aligned}$$ The latter two terms in (\[Sab-4th\]) are given, respectively, by $$\begin{aligned}
S_{\alpha \beta}^{(4)} [3] 
&=& 
\sum_{k L} (UX)_{ik} \hat{S}_{kL}^{(3)} \left\{ (W^{\dagger}) \right\}_{L j} 
+ \sum_{K l} W_{iK} \hat{S}_{Kl}^{(3)} \left\{ (UX)^{\dagger} \right\}_{lj}, 
\nonumber \\
S_{\alpha \beta}^{(4)} [2] &=& 
\sum_{K} W_{iK} \hat{S}_{KK}^{(2)} \left\{ (W^{\dagger}) \right\}_{K j} + 
\sum_{K \neq L} W_{iK} \hat{S}_{KL}^{(2)} \left\{ (W^{\dagger}) \right\}_{L j}. 
\label{S-alpha-beta-4th-[3]}\end{aligned}$$ We do not display explicitly the expression of each term in (\[Sab-4th\]). But, the notation of $S_{\alpha \beta}^{(4)} [n]_{ \text{ diag } }$ and $S_{\alpha \beta}^{(4)} [n]_{ \text{ offdiag } }$ will be transported to the notation for the oscillation probability such that $2 \mbox{Re} \left[ \left( S^{(0)}_{\alpha \beta} \right)^{*} S_{\alpha \beta}^{(4)} [n]_{ \text{ diag } } \right]$. Similarly, to make the equation fit to a single page we present the first and the second terms of $S_{\alpha \beta}^{(4)} [3]$ in (\[S-alpha-beta-4th-\[3\]\]) separately, $S_{\alpha \beta}^{(4)} [3]_\text{First} = \sum_{k L} (UX)_{\alpha k} W^*_{\beta L} \hat{S}_{kL}^{(3)}$ and $S_{\alpha \beta}^{(4)} [3]_\text{Second} = \sum_{L k} W_{\alpha L} (UX)^*_{\beta k} \hat{S}_{L k}^{(3)}$, whose notations are also transported to the oscillation probability.

Expression of the oscillation probability in fourth order in $W$ {#sec:expression-probability-4th}
================================================================

The oscillation probability to second order in $W$ is given in eq. (\[P-beta-alpha-0th+2nd\]) in section \[sec:probability-2nd\]. What is left is, therefore, the expressions of the oscillation probability in fourth order in $W$, the explicit form of the two terms in (\[P-beta-alpha-4th-def\]), $P(\nu_\beta \rightarrow \nu_\alpha) = 
\left| S^{(2)}_{\alpha \beta} \right|^2 
+ 2 \mbox{Re} \left[ \left( S^{(0)}_{\alpha \beta} \right)^{*} S^{(4)}_{\alpha \beta} \right] $.

Second order $S$ matrix squared term: $\left| S^{(2)}_{\alpha \beta} \right|^2$ {#sec:second-order-square}
-------------------------------------------------------------------------------

The $S$ matrix element $S^{(2)}_{\alpha \beta}$ in eq. (\[S-alpha-beta-2nd\]) contains four terms. To prevent too long expressions, we divide $\left| S^{(2)}_{\alpha \beta} \right|^2$ into the two terms, one sum of each term squared and the other one composed of cross terms. The first one is given by $$\begin{aligned}
&& \left| S^{(2)}_{\alpha \beta} \right|^2_{\text{1st}} = 
\sum_{k, K} \sum_{l, L} 
\frac{ 1 }{ ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{l} ) } 
\nonumber \\ 
&\times&
\biggl[
x^2 e^{- i ( h_{k} - h_{l} ) x} 
- (ix) \frac{e^{- i ( \Delta_{K} - h_{l} ) x} - e^{- i ( h_{k} - h_{l} ) x} }{ ( \Delta_{K} - h_{k} )  } 
+ (ix) \frac{e^{- i ( h_{k} - \Delta_{L} ) x} - e^{- i ( h_{k} - h_{l} ) x} }{ ( \Delta_{L} - h_{l} )  } 
\nonumber \\ 
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{l} ) } 
\biggl\{
e^{- i ( \Delta_{K} - \Delta_{L} ) x} + e^{- i ( h_{k} - h_{l} ) x} 
- e^{- i ( \Delta_{K} - h_{l} ) x} - e^{- i ( h_{k} - \Delta_{L} ) x} 
\biggr\}
\biggr] 
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\
&\times& 
(UX)_{\alpha l}^* (UX)_{\beta l} 
\left\{ (UX)^{\dagger} A W \right\}_{l L} 
\left\{ W^{\dagger} A (UX) \right\}_{L l} 
\nonumber \\ 
&+& 
\sum_{k \neq m} \sum_{K}  \sum_{l \neq n} \sum_{L}  
\frac{ 1 }{ ( h_{m}  - h_{k} ) (\Delta_{K} - h_{k}) (\Delta_{K} - h_{m}) } 
\frac{ 1 }{ ( h_{n}  - h_{l} ) (\Delta_{L} - h_{l}) (\Delta_{L} - h_{n}) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{k} \right) e^{- i h_{m} x} 
- \left( \Delta_{K} - h_{m} \right) e^{- i h_{k} x} 
- ( h_{m}  - h_{k} )  e^{- i \Delta_{K} x} 
\biggr]
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{L} - h_{l} \right) e^{+ i h_{n} x} 
- \left( \Delta_{L} - h_{n} \right) e^{+ i h_{l} x} 
- ( h_{n}  - h_{l} )  e^{+ i \Delta_{L} x} 
\biggr]
\nonumber \\
&\times&
(UX)_{\alpha k} (UX)^*_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K m} 
\nonumber \\
&\times& 
(UX)_{\alpha l}^* (UX)_{\beta n} 
\left\{ (UX)^{\dagger} A W \right\}_{n L} 
\left\{ W^{\dagger} A (UX) \right\}_{L l} 
\nonumber \\
&+& 
\sum_{k, K} \sum_{l, L} 
\frac{1}{ ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{l} ) } 
\left( e^{- i \Delta_{K} x} - e^{- i h_{k} x} \right) 
\left(  e^{+ i \Delta_{L} x} - e^{+ i h_{l} x} \right)
\nonumber \\
&\times&
\biggl[ 
(UX)_{\alpha k} W^*_{\beta K} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
+
W_{\alpha K} (UX)^*_{\beta k} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\biggr]
\nonumber \\
&\times& 
\biggl[ 
(UX)_{\alpha l}^* W_{\beta L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L l}
+
W_{\alpha L}^* (UX)_{\beta l} 
\left\{ (UX)^{\dagger} A W \right\}_{l L}
\biggr]
\nonumber \\
&+& 
\sum_{K} 
\vert W_{\alpha K} \vert^2 \vert W_{\beta K} \vert^2
+ \sum_{K \neq L} 
e^{- i ( \Delta_{K} - \Delta_{L} ) x} 
W_{\alpha K} W^*_{\beta K} W_{\alpha L}^* W_{\beta L}.  
\label{S(2)-squared-1}\end{aligned}$$ Except for the first term in (\[S(2)-squared-1\]) we did not try to unify the two exponential factors because the expressions become cumbersome. Apart from the last line in (\[S(2)-squared-1\]) all the terms are suppressed by the two energy denominators with $\Delta m^2_{J k}$ which doubly suppress the active-sterile state transition. The first term in the last line is the probability leaking term mentioned in section \[sec:nonunitarity-vacuum\].

The second term of $\left| S^{(2)}_{\alpha \beta} \right|^2$ (interference terms) is given by $$\begin{aligned}
&& \left| S^{(2)}_{\alpha \beta} \right|^2_{\text{2nd}} = 
- 2 \mbox{Re} 
\biggl\{ 
\sum_{k, K} 
\sum_{l \neq m} \sum_{L} 
\frac{ 1 }{ ( \Delta_{K} - h_{k} ) (\Delta_{L} - h_{l}) (\Delta_{L} - h_{m}) ( h_{m}  - h_{l} ) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{L} - h_{l} \right) e^{- i h_{m} x} 
- \left( \Delta_{L} - h_{m} \right) e^{- i h_{l} x} 
- ( h_{m}  - h_{l} )  e^{- i \Delta_{L} x} 
\biggr]
\left[
- (ix) e^{+ i h_{k} x} + \frac{e^{+ i \Delta_{K} x} - e^{+ i h_{k} x} }{ ( \Delta_{K} - h_{k} )  } 
\right]
\nonumber \\
&\times&
(UX)_{\alpha l} (UX)^*_{\beta m} 
(UX)_{\alpha k}^* (UX)_{\beta k} 
\nonumber \\
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{l L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\nonumber \\
&+& 
2 \mbox{Re} 
\biggl\{ 
\sum_{k, K} \sum_{l, L} 
\left[
- (ix) e^{+ i h_{k} x} + \frac{e^{+ i \Delta_{K} x} - e^{+ i h_{k} x} }{ ( \Delta_{K} - h_{k} )  } 
\right]
\frac{e^{- i \Delta_{L} x} - e^{- i h_{l} x} }{ ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{l} ) } 
\nonumber \\ 
&\times&
(UX)_{\alpha k}^* (UX)_{\beta k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\nonumber \\ 
&\times&
\biggl[ 
(UX)_{\alpha l} W^*_{\beta L} 
\left\{ (UX)^{\dagger} A W \right\}_{l L} 
+
W_{\alpha L} (UX)^*_{\beta l} 
\left\{ W ^{\dagger} A (UX) \right\}_{L l}
\biggr]
\biggr\}
\nonumber \\
&+& 
2 \mbox{Re} 
\biggl\{
\sum_{k, K} \sum_{L} 
\left[
- (ix) e^{- i ( \Delta_{L} - h_{k} ) x} + \frac{ e^{- i ( \Delta_{L} - \Delta_{K} ) x} - e^{- i ( \Delta_{L} - h_{k} ) x} }{ ( \Delta_{K} - h_{k} )  } 
\right] 
\frac{ 1 }{ \Delta_{K} - h_{k} } 
\nonumber \\ 
&\times&
(UX)_{\alpha k}^* (UX)_{\beta k} 
W_{\alpha L} W^*_{\beta L} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\biggr\}
\nonumber \\
&-& 
2 \mbox{Re} 
\biggl\{
\sum_{k \neq m} \sum_{K} \sum_{l, L} 
\frac{ 1 }{ ( h_{m}  - h_{k} ) (\Delta_{K} - h_{k}) (\Delta_{K} - h_{m}) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{k} \right) e^{+ i h_{m} x} 
- \left( \Delta_{K} - h_{m} \right) e^{+ i h_{k} x} 
- ( h_{m}  - h_{k} )  e^{+ i \Delta_{K} x} 
\biggr]
\frac{e^{- i \Delta_{L} x} - e^{- i h_{l} x} }{ ( \Delta_{L} - h_{l} ) } 
\nonumber \\
&\times& 
(UX)^*_{\alpha k} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{m K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\nonumber \\
&\times& 
\biggl[ 
(UX)_{\alpha l} W^*_{\beta L} 
\left\{ (UX)^{\dagger} A W \right\}_{l L} 
+
W_{\alpha L} (UX)^*_{\beta l} 
\left\{ W ^{\dagger} A (UX) \right\}_{L l}
\biggr]
\nonumber \\
&-& 
2 \mbox{Re} 
\biggl\{
\sum_{k \neq m} \sum_{K} \sum_{L} 
\frac{ 1 }{ ( h_{m}  - h_{k} ) (\Delta_{K} - h_{k}) (\Delta_{K} - h_{m}) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{k} \right) e^{- i ( \Delta_{L} - h_{m} ) x} 
- \left( \Delta_{K} - h_{m} \right) e^{- i ( \Delta_{L} - h_{k} ) x}
- ( h_{m}  - h_{k} ) e^{- i ( \Delta_{L} - \Delta_{K} ) x}
\biggr]
\nonumber \\
&\times& 
(UX)^*_{\alpha k} (UX)_{\beta m} 
W_{\alpha L} W^*_{\beta L} 
\left\{ (UX)^{\dagger} A W \right\}_{m K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k}
\biggr\}
\nonumber \\
&+& 
2 \mbox{Re} 
\biggl\{
\sum_{K} \sum_{l, L} 
\frac{e^{- i ( \Delta_{L} - \Delta_{K} ) x} - e^{- i ( h_{l} - \Delta_{K} ) x} }{ ( \Delta_{L} - h_{l} ) } 
\nonumber \\
&\times&
\biggl[ 
W^*_{\alpha K} W_{\beta K} 
(UX)_{\alpha l} W^*_{\beta L} 
\left\{ (UX)^{\dagger} A W \right\}_{l L} 
+
W^*_{\alpha K} W_{\beta K} 
W_{\alpha L} (UX)^*_{\beta l} 
\left\{ W ^{\dagger} A (UX) \right\}_{L l}
\biggr]
\biggr\}.
\label{S(2)-squared-2}\end{aligned}$$

Interference terms of the type $2 \mbox{Re} \left[ \left( S^{(0)}_{\alpha \beta} \right)^{*} S^{(4)}_{\alpha \beta} \right]$ {#sec:interference}
----------------------------------------------------------------------------------------------------------------------------

We classify the fourth order in $W$ contribution of the interference terms into 8 terms: $$\begin{aligned}
P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{ \text{interference} }
&=& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{1st} 
+ P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{2nd} 
+ P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{3rd} 
\nonumber \\
&+&
P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{4th-s} 
+ P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{4th-d} 
\nonumber \\
&+& 
P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{5th-1st}
+ P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{5th-2nd} 
+ P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{6th}. 
\nonumber \\\end{aligned}$$ The nature of each term is explicitly indicated as follows: $$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{1st} 
\equiv 
2 \mbox{Re} \left[ \left( 
S^{(0)}_{\alpha \beta} \right)^{*} 
S_{\alpha \beta}^{(4)} [3]_{ \text{ diag } } 
\right]
\nonumber \\
&=& 
2 \mbox{Re} 
\biggl\{
- \sum_{K} \sum_{k} 
\sum_{m} 
\biggl[ 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 } 
\left\{ (ix) + \frac{ 2 }{ ( \Delta_{K} - h_{k} ) } \right\}
\left( e^{- i ( \Delta_{K} - h_{m} ) x} + e^{- i ( h_{k} - h_{m} ) x} \right) 
\biggr]
\nonumber \\ 
&\times& 
(UX)^*_{\alpha m} (UX)_{\beta m} 
(UX)_{\alpha k} (UX)^*_{\beta k}
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A W \right\}_{K K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\ 
&+& 
\sum_{K \neq L} \sum_{k} 
\sum_{m} 
\biggl[ 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 ( \Delta_{K} - \Delta_{L} )} 
e^{- i ( \Delta_{K} - h_{m} )  x} 
-
\frac{ 1 }{ ( \Delta_{L} - h_{k} )^2 ( \Delta_{K} - \Delta_{L} )} 
e^{- i ( \Delta_{L} - h_{m} ) x} 
\nonumber \\ 
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 ( \Delta_{L} - h_{k} )^2 }  
\left( \Delta_{K} + \Delta_{L} - 2 h_{k} \right) 
e^{- i ( h_{k} - h_{m} ) x} 
\biggr]
\nonumber \\ 
&\times& 
(UX)^*_{\alpha m} (UX)_{\beta m} 
(UX)_{\alpha k} (UX)^*_{\beta k}
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A W \right\}_{K L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L k} 
\biggr\}. 
\nonumber \\ 
\label{P-beta-alpha-W4-H3-diag}\end{aligned}$$ $$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{2nd} 
\equiv 
2 \mbox{Re} \left[ \left( 
S^{(0)}_{\alpha \beta} \right)^{*} 
S_{\alpha \beta}^{(4)} [4]_{ \text{ diag } } 
\right] 
\nonumber \\
&=& 
2 \mbox{Re} 
\biggl\{ 
\sum_{n} \sum_{k} \sum_{K} 
\biggl[
- \frac{x^2}{2} \frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 } 
e^{- i ( h_{k} - h_{n} ) x} 
- 
\frac{ 2 (ix) }{ ( \Delta_{K} - h_{k} )^3 } 
e^{- i ( h_{k} - h_{n} ) x} 
\nonumber \\
&-& 
\frac{ (ix) }{ ( \Delta_{K} - h_{k} )^3 } 
e^{- i ( \Delta_{K} - h_{n} ) x} 
-
\frac{ 3 }{ ( \Delta_{K} - h_{k} )^4 } 
\left( e^{- i ( \Delta_{K} - h_{n} ) x} - e^{- i ( h_{k} - h_{n} ) x} \right)
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\
&+& 
\sum_{n} \sum_{k} \sum_{K} \sum_{m \neq k} 
\biggl[
\frac{ (ix) }{ ( \Delta_{K} - h_{k} )^2 ( h_{m}  - h_{k} ) } 
e^{- i ( h_{k} - h_{n} ) x} 
\nonumber \\
&-& 
\frac{ (ix) }{ ( \Delta_{K} - h_{k} )^2 ( \Delta_{K} - h_{m} ) } 
e^{- i ( \Delta_{K} - h_{n} ) x} 
+
\frac{ ( h_{k} + 2 h_{m} - 3 \Delta_{K} ) }{ ( \Delta_{K} - h_{k} )^3 ( \Delta_{K} - h_{m} )^2 } 
e^{- i ( \Delta_{K} - h_{n} ) x} 
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{m} )^2 ( h_{m}  - h_{k} )^2 } 
e^{- i ( h_{m} - h_{n} ) x}
- 
\frac{ \left( \Delta_{K} + 2 h_{m} - 3 h_{k} \right) }{ ( \Delta_{K} - h_{k} )^3 ( h_{m}  - h_{k} )^2 } 
e^{- i ( h_{k} - h_{n} ) x} 
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K m} 
\left\{ (UX)^{\dagger} A W \right\}_{m K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\
&+& 
\sum_{n} \sum_{k} \sum_{K \neq L} 
\biggl[
- \frac{x^2}{2} 
\frac{ 1 }{ ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{k} ) } 
e^{- i ( h_{k} - h_{n} ) x} 
- (ix) 
\frac{ \left( \Delta_{K} + \Delta_{L} - 2 h_{k} \right) }{ ( \Delta_{K} - h_{k} )^2 ( \Delta_{L} - h_{k} )^2 } 
e^{- i ( h_{k} - h_{n} ) x} 
\nonumber \\
&-& 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^3 ( \Delta_{L} - \Delta_{K} ) } 
e^{- i ( \Delta_{K} - h_{n} ) x} 
+ \frac{ 1 }{ ( \Delta_{L} - h_{k} )^3 ( \Delta_{L} - \Delta_{K} ) } 
e^{- i ( \Delta_{L} - h_{n} ) x} 
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^3 ( \Delta_{L} - h_{k} )^3 } 
\biggl\{
\Delta_{L}^2 + \Delta_{L} \Delta_{K} + \Delta_{K}^2 - 3 h_{k} ( \Delta_{L} + \Delta_{K} ) + 3 h_{k}^2  
\biggr\}
e^{- i ( h_{k} - h_{n} ) x} 
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L k} 
\nonumber \\
&+& 
\sum_{n} \sum_{k} \sum_{K \neq L} \sum_{m \neq k} 
\biggl[
\frac{ (ix) }{ ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{k} ) ( h_{m} - h_{k} ) }
e^{- i ( h_{k} - h_{n} ) x} 
\nonumber \\
&-& 
\frac{1}{ ( \Delta_{K} - h_{k} )^2 ( \Delta_{L} - h_{k} )^2 ( h_{m} - h_{k} )^2 } 
\biggl\{ 
\Delta_{K} \Delta_{L} + ( h_{m} - 2 h_{k} ) ( \Delta_{K} + \Delta_{L} ) + 3 
h_{k}^2 - 2 h_{m} h_{k}
\biggr\}
e^{- i ( h_{k} - h_{n} ) x} 
\nonumber \\
&+&
\frac{1}{ ( \Delta_{K} - h_{m} ) ( \Delta_{L} - h_{m} ) ( h_{m} - h_{k} )^2 } 
e^{- i ( h_{m} - h_{n} ) x} 
\nonumber \\
&+& 
\frac{1}{ ( \Delta_{K} - \Delta_{L} ) ( \Delta_{K} - h_{k} )^2 ( \Delta_{K} - h_{m} ) } e^{- i ( \Delta_{K} - h_{n} ) x} 
- \frac{1}{ ( \Delta_{K} - \Delta_{L} ) ( \Delta_{L} - h_{k} )^2 ( \Delta_{L} - h_{m} ) } e^{- i ( \Delta_{L} - h_{n} ) x} 
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K m} 
\left\{ (UX)^{\dagger} A W \right\}_{m L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L k} 
\biggr\}. 
\label{P-beta-alpha-W4-H4-diag}\end{aligned}$$ $$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{3rd} 
\equiv 
2 \mbox{Re} \left[ \left( 
S^{(0)}_{\alpha \beta} \right)^{*} 
S_{\alpha \beta}^{(4)} [4]_{ \text{ offdiag } } 
\right]
\nonumber \\
&=& 
2 \mbox{Re} 
\biggl\{ 
\sum_{m} 
\sum_{k \neq l } \sum_{K} 
\biggl[ 
- \frac{ (ix) }{ ( \Delta_{K} - h_{k} ) (\Delta_{K} - h_{l} ) } 
e^{- i ( \Delta_{K} - h_{m} ) x} 
+ \frac{ 1 }{  ( h_{l}  - h_{k} ) ( \Delta_{K} - h_{k} )^2  ( \Delta_{K} - h_{l} )^2 } 
\nonumber \\ 
&\times& 
\biggl\{
( h_{l} - h_{k} ) ( h_{l} + h_{k} - 2  \Delta_{K} ) 
e^{- i ( \Delta_{K} - h_{m} ) x} 
+ ( \Delta_{K} - h_{k} )^2 e^{ - i ( h_{l} - h_{m} ) x} 
- ( \Delta_{K} - h_{l} )^2 e^{ - i ( h_{k} - h_{m} ) x} 
\biggr\}
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A W \right\}_{K K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\nonumber \\ 
&+& 
\sum_{m} 
\sum_{k \neq l } \sum_{K \neq L } 
\frac{ 1 }{ ( h_{l}  - h_{k} ) ( \Delta_{L} - \Delta_{K} ) ( \Delta_{K} - h_{k} ) ( \Delta_{K} - h_{l} ) ( \Delta_{L} - h_{k} ) ( \Delta_{L} - h_{l} )   } 
\nonumber \\ 
&\times& 
\biggl[ 
( h_{l}  - h_{k} ) 
\biggl\{ ( \Delta_{K} - h_{k} ) ( \Delta_{K} - h_{l} ) 
e^{- i ( \Delta_{L} - h_{m} ) x} - 
( \Delta_{L} - h_{k} ) ( \Delta_{L} - h_{l} ) 
e^{- i ( \Delta_{K} - h_{m} ) x}
\biggr\} 
\nonumber \\ 
&+& 
( \Delta_{L} - \Delta_{K} ) 
\biggl\{ ( \Delta_{K} - h_{k} ) ( \Delta_{L} - h_{k} ) e^{- i ( h_{l} - h_{m} ) x} - 
( \Delta_{K} - h_{l} ) ( \Delta_{L} - h_{l} ) e^{- i ( h_{k} - h_{m} ) x} \biggr\} 
\biggl]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha m} (UX)_{\beta m} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A W \right\}_{K L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L l} 
\biggr\}.
\label{P-beta-alpha-W4-H3-offdiag}\end{aligned}$$ $$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{4th-s} 
\equiv 
2 \mbox{Re} \left[ \left( 
S^{(0)}_{\alpha \beta} \right)^{*} 
S_{\alpha \beta}^{(4)} [4]_{ \text{ offdiag } } (\text{single})
\right]
\nonumber \\
&=& 
2 \mbox{Re} 
\biggl\{
\sum_{n} 
\sum_{k \neq l } 
\sum_{K}
\biggl[
\frac{ (ix) }{ ( \Delta_{K} - h_{k} )^2 ( h_{l} - h_{k} ) } 
e^{ - i ( h_{k} -  h_{n} ) x} 
- 
\frac{ (ix) }{ ( \Delta_{K} - h_{k} )^2 ( \Delta_{K} - h_{l} ) } 
e^{ - i ( \Delta_{K} -  h_{n} ) x} 
\nonumber \\
&+&
\frac{ 1 }{ ( \Delta_{K} - h_{l} )^2 ( h_{l} - h_{k} )^2 } 
e^{ - i ( h_{l} -  h_{n} ) x}
- \frac{ 1 }{ ( \Delta_{K} - h_{k} )^3 ( h_{l} - h_{k} )^2 } 
\left( \Delta_{K} + 2 h_{l} - 3 h_{k} \right) 
e^{ - i ( h_{k} -  h_{n} ) x} 
\nonumber \\
&+&
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^3 ( \Delta_{K} - h_{l} )^2 } 
\left( h_{k} + 2 h_{l}  - 3 \Delta_{K} \right) 
e^{ - i ( \Delta_{K} -  h_{n} ) x} 
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l}
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\nonumber \\
&+& 
\sum_{n} 
\sum_{k \neq l }  
\sum_{K}
\biggl[
- 
\frac{ (ix) }{ ( \Delta_{K} - h_{l} )^2 ( h_{l} - h_{k} ) } 
e^{ - i ( h_{l} -  h_{n} ) x} 
- 
\frac{ (ix) }{ ( \Delta_{K} - h_{l} )^2 ( \Delta_{K} - h_{k} ) } 
e^{ - i ( \Delta_{K} -  h_{n} ) x}
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 ( \Delta_{K} - h_{l} )^3 } 
\left( h_{l} + 2 h_{k} - 3 \Delta_{K} \right) 
e^{ - i ( \Delta_{K} -  h_{n} ) x}
\nonumber \\
&-& 
\frac{ 1 }{ ( \Delta_{K} - h_{l} )^3 ( h_{l} - h_{k} )^2 } 
\left( \Delta_{K} + 2 h_{k} - 3 h_{l}  \right) 
e^{ - i ( h_{l} -  h_{n} ) x} 
+ \frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 ( h_{l} - h_{k} )^2 } 
e^{ - i ( h_{k} -  h_{n} ) x} 
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\left\{ (UX)^{\dagger} A W \right\}_{l K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\nonumber \\ 
&+& 
\sum_{n} 
\sum_{k \neq l } 
\sum_{K} \sum_{m \neq k, l} 
\biggl[
\frac{ (ix) }{ ( \Delta_{K} - h_{k} ) ( \Delta_{K} - h_{l} ) ( \Delta_{K} - h_{m} ) } 
e^{ - i ( \Delta_{K} -  h_{n} ) x} 
\nonumber \\
&-& 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 ( \Delta_{K} - h_{l} )^2 ( \Delta_{K} - h_{m} )^2 } 
\nonumber \\
&\times&
\left\{ 
3 \Delta_{K}^2 - 2 \Delta_{K} \left( h_{k} + h_{l} + h_{m} \right) + \left( h_{k} h_{l}+ h_{l} h_{m} + h_{m} h_{k} \right)
\right\}
e^{ - i ( \Delta_{K} -  h_{n} ) x}
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{m} )^2 ( h_{m} - h_{k} ) ( h_{m} - h_{l} ) } 
e^{ - i ( h_{m} -  h_{n} ) x} 
\nonumber \\
&+& 
\frac{ 1 }{ ( \Delta_{K} - h_{l} )^2 ( h_{l} - h_{m} ) ( h_{l} - h_{k} ) } 
e^{ - i ( h_{l} -  h_{n} ) x} 
- 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^2 ( h_{k} - h_{m} ) ( h_{l} - h_{k} ) } 
e^{ - i ( h_{k} -  h_{n} ) x} 
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K m} 
\left\{ (UX)^{\dagger} A W \right\}_{m K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\biggr\}. 
\label{P-beta-alpha-W4-H4-single}\end{aligned}$$ $$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{4th-d} 
\equiv 
2 \mbox{Re} \left[ \left( 
S^{(0)}_{\alpha \beta} \right)^{*} 
S_{\alpha \beta}^{(4)} [4]_{ \text{ offdiag } } (\text{double})  
\right] 
\nonumber \\
&=& 
2 \mbox{Re} 
\biggl\{
\sum_{n} 
\sum_{k \neq l } 
\sum_{K \neq L} 
\biggl[ 
\frac{ (ix) }{ (\Delta_{K} - h_{k} ) (\Delta_{L} - h_{k} ) ( h_{l}  - h_{k} ) } 
e^{- i ( h_{k} - h_{n} ) x} 
\nonumber \\
&-& 
\frac{ 1 }{ (\Delta_{K} - \Delta_{L} ) (\Delta_{L} - h_{k} )^2 (\Delta_{L} - h_{l} ) } 
e^{- i ( \Delta_{L} - h_{n} ) x} 
+ \frac{ 1 }{ (\Delta_{K} - \Delta_{L} ) (\Delta_{K} - h_{k} )^2 (\Delta_{K} - h_{l} ) } 
e^{- i ( \Delta_{K} - h_{n} ) x}
\nonumber \\
&+& 
\frac{ 1 }{ (\Delta_{K} - h_{l} ) (\Delta_{L} - h_{l} ) ( h_{l}  - h_{k} )^2 } 
e^{- i ( h_{l} - h_{n} ) x} 
\nonumber \\
&-& 
\frac{ 1 }{ (\Delta_{K} - h_{k} )^2 (\Delta_{L} - h_{k} )^2 ( h_{l}  - h_{k} )^2 } 
\biggl\{
3 h_{k}^2 - 2 h_{k} h_{l} + \left( h_{l} - 2 h_{k} \right)  (\Delta_{K} + \Delta_{L} ) + \Delta_{K} \Delta_{L} 
\biggr\}
e^{- i ( h_{k} - h_{n} ) x} 
\biggr] 
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L l} 
\nonumber \\
&+& 
\sum_{n} 
\sum_{k \neq l } 
\sum_{K \neq L} 
\biggl[
- 
\frac{ (ix) }{ (\Delta_{K} - h_{l} ) (\Delta_{L} - h_{l} ) ( h_{l}  - h_{k} ) } 
e^{- i ( h_{l} - h_{n} ) x} 
\nonumber \\
&+& 
\frac{ 1 }{ (\Delta_{K} - \Delta_{L} ) (\Delta_{K} - h_{k} ) (\Delta_{K} - h_{l} )^2 } 
e^{- i ( \Delta_{K} - h_{n} ) x}
- \frac{ 1 }{ (\Delta_{K} - \Delta_{L} ) (\Delta_{L} - h_{k} ) (\Delta_{L} - h_{l} )^2 } 
e^{- i ( \Delta_{L} - h_{n} ) x} 
\nonumber \\
&+&
\frac{ 1 }{ (\Delta_{K} - h_{k} ) (\Delta_{L} - h_{k} ) ( h_{l}  - h_{k} )^2 } 
e^{- i ( h_{k} - h_{n} ) x} 
\nonumber \\
&-& 
\frac{ 1 }{ (\Delta_{K} - h_{l} )^2 (\Delta_{L} - h_{l} )^2 ( h_{l}  - h_{k} )^2 } 
\biggl\{
3 h_{l}^2 - 2 h_{k} h_{l} - \left( 2 h_{l} - h_{k} \right)  (\Delta_{K} + \Delta_{L} ) + \Delta_{K} \Delta_{L} 
\biggr\}
e^{- i ( h_{l} - h_{n} ) x} 
\biggr] 
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K l} 
\left\{ (UX)^{\dagger} A W \right\}_{l L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L l} 
\nonumber \\
&+& 
\sum_{n} 
\sum_{k \neq l } 
\sum_{K \neq L} \sum_{m \neq k, l} 
\biggl[ 
\frac{1}{ ( h_{m} - h_{k} ) ( h_{m} - h_{l} ) ( h_{l} - h_{k} ) ( \Delta_{K} - h_{k} ) ( \Delta_{K} - h_{l} ) ( \Delta_{L} - h_{k} ) ( \Delta_{L} - h_{l} ) } 
\nonumber \\
&+& 
\biggl\{ 
( h_{m} - h_{l} ) (\Delta_{K} - h_{l} ) ( \Delta_{L} - h_{l} ) 
e^{- i ( h_{k} - h_{n} ) x} 
-
( h_{m} - h_{k} ) (\Delta_{K} - h_{k} ) ( \Delta_{L} - h_{k} ) 
e^{- i ( h_{l} - h_{n} ) x} 
\biggr\} 
\nonumber \\
&+& 
\frac{1}{ ( h_{m} - h_{k} ) ( h_{m} - h_{l} ) ( \Delta_{K} - h_{m} ) ( \Delta_{L} - h_{m} ) } 
e^{- i ( h_{m} - h_{n} ) x} 
\nonumber \\
&+& 
\frac{1}{ ( \Delta_{K} - h_{k} ) ( \Delta_{K} - h_{l} ) ( \Delta_{K} - h_{m} ) ( \Delta_{L} - h_{k} ) ( \Delta_{L} - h_{l} ) ( \Delta_{L} - h_{m} ) ( \Delta_{L} - \Delta_{K} ) } 
\nonumber \\ 
&\times&
\biggl\{  ( \Delta_{K} - h_{k} ) ( \Delta_{K} - h_{l} ) ( \Delta_{K} - h_{m} ) 
e^{- i ( \Delta_{L} - h_{n} ) x} 
- ( \Delta_{L} - h_{k} ) ( \Delta_{L} - h_{l} ) ( \Delta_{L} - h_{m} ) 
e^{- i ( \Delta_{K} - h_{n} ) x}
\biggr\}
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} (UX)^*_{\beta l} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\nonumber \\ 
&\times&
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W ^{\dagger} A (UX) \right\}_{K m} 
\left\{ (UX)^{\dagger} A W \right\}_{m L} 
\left\{ W ^{\dagger} A (UX) \right\}_{L l} 
\biggr\}.
\label{P-beta-alpha-W4-H4-double}\end{aligned}$$

$$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{5th-1st} 
\equiv 
2 \mbox{Re} \left[ \left( 
S^{(0)}_{\alpha \beta} \right)^{*} 
S_{\alpha \beta}^{(4)} [3]_\text{First} 
\right] 
\nonumber \\
&=& 
2 \mbox{Re} 
\biggl\{
- \sum_{n} 
\sum_{k L} 
\frac{ 1 }{ \Delta_{L} - h_{k} } 
\left[ (ix) e^{- i ( \Delta_{L} - h_{n} ) x} + \frac{ e^{- i ( \Delta_{L} - h_{n} ) x} - 
e^{- i ( h_{k} - h_{n} ) x} }{ ( \Delta_{K} - h_{k} ) } 
\right]
\nonumber \\
&\times&
(UX)_{\alpha k} W^*_{\beta L} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\left\{ W^{\dagger} A W \right\}_{L L} 
\nonumber \\
&+& 
\sum_{n} 
\sum_{k L} 
\sum_{K \neq L} 
\frac{ 1 }{ ( \Delta_{L} - \Delta_{K} ) ( \Delta_{L} - h_{k} ) ( \Delta_{K} - h_{k} ) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{k} \right) 
e^{- i ( \Delta_{L} - h_{n} ) x}
- \left( \Delta_{L} - h_{k} \right) 
e^{- i ( \Delta_{K} - h_{n} ) x}
- \left( \Delta_{K} - \Delta_{L} \right) 
e^{- i ( h_{k} - h_{n} ) x}
\biggr]
\nonumber \\
&\times&
(UX)_{\alpha k} W^*_{\beta L} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W^{\dagger} A W \right\}_{K L} 
\nonumber \\ 
&-& 
\sum_{n} 
\sum_{k L} 
\frac{ 1 }{ ( \Delta_{L} - h_{k} )^2 } 
\biggl[
(ix) \left( e^{- i ( h_{k} - h_{n} ) x} + e^{- i ( \Delta_{L} - h_{n} ) x} \right)
+ 2 \frac{e^{- i ( \Delta_{L} - h_{n} ) x} - e^{- i ( h_{k} - h_{n} ) x} }{ ( \Delta_{L} - h_{k} ) } 
\biggr] 
\nonumber \\ 
&\times& 
(UX)_{\alpha k} W^*_{\beta L} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\left\{ W^{\dagger} A (UX) \right\}_{L k} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\nonumber \\ 
&+& 
\sum_{n} 
\sum_{k L} 
\sum_{m \neq k} 
\biggl[
- \frac{ (ix) e^{- i ( \Delta_{L} - h_{n} ) x} }{ ( \Delta_{L} - h_{k} )( \Delta_{L} - h_{m} ) } 
+ \frac{ 1 }{ ( h_{k} - h_{m} ) (\Delta_{L} - h_{k} )^2 (\Delta_{L} - h_{m} )^2 } 
\nonumber \\ 
&\times& 
\biggl\{
(\Delta_{L} - h_{m} )^2 
e^{- i ( h_{k} - h_{n} ) x} 
- (\Delta_{L} - h_{k} )^2 
e^{- i ( h_{m} - h_{n} ) x} 
+ ( h_{k} - h_{m} )( h_{k} + h_{m} - 2 \Delta_{L} ) e^{- i ( \Delta_{L} - h_{n} ) x} 
\biggr\} 
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} W^*_{\beta L} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\left\{ W^{\dagger} A (UX) \right\}_{L m} 
\left\{ (UX)^{\dagger} A W \right\}_{m L} 
\nonumber \\ 
&+& 
\sum_{n} 
\sum_{k L} 
\sum_{K \neq L} 
\biggl[
- \frac{ (ix) e^{- i ( h_{k} - h_{n} ) x} }{ (\Delta_{L} - h_{k}) (\Delta_{K} - h_{k}) }   
+ \frac{ 1 }{ (\Delta_{L} - \Delta_{K}) (\Delta_{L} - h_{k})^2 ( \Delta_{K} - h_{k} )^2 }
\nonumber \\ 
&\times&
\biggl\{
( \Delta_{K} - h_{k} )^2 
e^{- i ( \Delta_{L} - h_{n} ) x} 
- (\Delta_{L} - h_{k})^2 
e^{- i ( \Delta_{K} - h_{n} ) x}
+ (\Delta_{L} - \Delta_{K}) (\Delta_{L} + \Delta_{K} - 2 h_{k} ) 
e^{- i ( h_{k} - h_{n} ) x} 
\biggr\}
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} W^*_{\beta L} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\nonumber \\ 
&+& 
\sum_{n} 
\sum_{k L} 
\sum_{K \neq L} 
\sum_{m \neq k} 
\frac{ 1 }{ (\Delta_{L} - \Delta_{K}) ( h_{m}  - h_{k} ) (\Delta_{L} - h_{k}) (\Delta_{L} - h_{m}) (\Delta_{K} - h_{k}) (\Delta_{K} - h_{m}) } 
\nonumber \\ 
&\times&
\biggl[
( h_{m}  - h_{k} ) 
\biggl\{(\Delta_{K} - h_{k}) (\Delta_{K} - h_{m}) 
e^{- i ( \Delta_{L} - h_{n} ) x} 
- (\Delta_{L} - h_{k}) (\Delta_{L} - h_{m}) 
e^{- i ( \Delta_{K} - h_{n} ) x}
\biggr\}
\nonumber \\ 
&-&
(\Delta_{L} - \Delta_{K}) 
\biggl\{ (\Delta_{L} - h_{m}) (\Delta_{K} - h_{m}) 
e^{- i ( h_{k} - h_{n} ) x} 
- (\Delta_{L} - h_{k}) (\Delta_{K} - h_{k}) 
e^{- i ( h_{m} - h_{n} ) x} 
\biggr\}
\biggr]
\nonumber \\ 
&\times&
(UX)_{\alpha k} W^*_{\beta L} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W^{\dagger} A (UX) \right\}_{K m} 
\left\{ (UX)^{\dagger} A W \right\}_{m L} 
\biggr\}. 
\label{P-beta-alpha-W4-H3-First} \end{aligned}$$

$$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{5th-2nd} 
\equiv 
2 \mbox{Re} \left[ \left( 
S^{(0)}_{\alpha \beta} \right)^{*} 
S_{\alpha \beta}^{(4)} [3]_\text{Second} 
\right] 
\nonumber \\
&=& 
2 \mbox{Re} 
\biggl\{
- \sum_{n} 
\sum_{L k} 
\frac{ 1 }{ \Delta_{L} - h_{k} } 
\left[
(ix) e^{- i ( \Delta_{L} - h_{n} ) x} 
+ 
\frac{ e^{- i ( \Delta_{L} - h_{n} ) x} - e^{- i ( h_{k} - h_{n} ) x} }{ ( \Delta_{L} - h_{k} ) } 
\right]
\nonumber \\
&\times& 
W_{\alpha L} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ W^{\dagger} A W \right\}_{L L} 
\left\{ W^{\dagger} A (UX) \right\}_{L k} 
\nonumber \\
&+& 
\sum_{n} 
\sum_{L k} 
\sum_{K \neq L} 
\frac{ 1 }{ ( \Delta_{L} - \Delta_{K} ) ( \Delta_{L} - h_{k} ) ( \Delta_{K} - h_{k} ) } 
\nonumber \\
&\times& 
\biggl[
\left( \Delta_{K} - h_{k} \right) 
e^{- i ( \Delta_{L} - h_{n} ) x} 
- \left( \Delta_{L} - h_{k} \right) 
e^{- i ( \Delta_{K} - h_{n} ) x} 
- \left( \Delta_{K} - \Delta_{L} \right) 
e^{- i ( h_{k} - h_{n} ) x} 
\biggr]
\nonumber \\ 
&\times&
W_{\alpha L} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ W^{\dagger} A W \right\}_{L K}  
\left\{ W^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\ 
&-& 
\sum_{n} 
\sum_{L k} 
\frac{ 1 }{ (\Delta_{L} - h_{k})^2 } 
\biggl[
(ix) \left( e^{- i ( h_{k} - h_{n} ) x} + e^{- i ( \Delta_{L} - h_{n} ) x}  \right)
+ 2 \frac{ e^{- i ( \Delta_{L} - h_{n} ) x} - e^{- i ( h_{k} - h_{n} ) x}  }{ ( \Delta_{L} - h_{k} ) } 
\biggr]
\nonumber \\ 
&\times& 
W_{\alpha L} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ W^{\dagger} A (UX) \right\}_{L k} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\left\{ W^{\dagger} A (UX) \right\}_{L k} 
\nonumber \\ 
&+& 
\sum_{n} 
\sum_{L k} 
\sum_{m \neq k} 
\biggl[
- \frac{ (ix) e^{- i ( \Delta_{L} - h_{n} ) x} }{ ( \Delta_{L} - h_{k} )( \Delta_{L} - h_{m} ) } 
+ \frac{ 1 }{ ( h_{k} - h_{m} ) (\Delta_{L} - h_{k} )^2 (\Delta_{L} - h_{m} )^2 } 
\nonumber \\ 
&\times& 
\biggl\{
(\Delta_{L} - h_{m} )^2 
e^{- i ( h_{k} - h_{n} ) x} 
- (\Delta_{L} - h_{k} )^2 
e^{- i ( h_{m} - h_{n} ) x} 
+ ( h_{k} - h_{m} )( h_{k} + h_{m} - 2 \Delta_{L} ) e^{- i ( \Delta_{L} - h_{n} ) x} 
\biggr\}
\biggr]
\nonumber \\ 
&\times& 
W_{\alpha L} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ W^{\dagger} A (UX) \right\}_{L m} 
\left\{ (UX)^{\dagger} A W \right\}_{m L} 
\left\{ W^{\dagger} A (UX) \right\}_{L k} 
\nonumber \\ 
&+& 
\sum_{n} 
\sum_{L k} 
\sum_{K \neq L} 
\biggl[
- \frac{ (ix) e^{- i ( h_{k} - h_{n} ) x} }{ (\Delta_{L} - h_{k}) (\Delta_{K} - h_{k}) }  
+ \frac{ 1 }{ (\Delta_{L} - \Delta_{K}) (\Delta_{L} - h_{k})^2 ( \Delta_{K} - h_{k} )^2 }
\nonumber \\ 
&\times&
\biggl\{
( \Delta_{K} - h_{k} )^2 
e^{- i ( \Delta_{L} - h_{n} ) x} 
- (\Delta_{L} - h_{k})^2 
e^{- i ( \Delta_{K} - h_{n} ) x} 
+ 
(\Delta_{L} - \Delta_{K}) (\Delta_{L} + \Delta_{K} - 2 h_{k} ) e^{- i ( h_{k} - h_{n} ) x} 
\biggr\}
\biggr]
\nonumber \\ 
&\times&
W_{\alpha L} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ W^{\dagger} A (UX) \right\}_{L k} 
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\left\{ W^{\dagger} A (UX) \right\}_{K k} 
\nonumber \\ 
&+& 
\sum_{n} 
\sum_{L k} 
\sum_{K \neq L} 
\sum_{m \neq k} 
\frac{ 1 }{ (\Delta_{L} - \Delta_{K}) ( h_{m}  - h_{k} ) (\Delta_{L} - h_{k}) (\Delta_{L} - h_{m}) (\Delta_{K} - h_{k}) (\Delta_{K} - h_{m}) } 
\nonumber \\ 
&\times&
\biggl[
( h_{m}  - h_{k} ) 
\biggl\{ (\Delta_{K} - h_{k}) (\Delta_{K} - h_{m}) 
e^{- i ( \Delta_{L} - h_{n} ) x} 
- (\Delta_{L} - h_{k}) (\Delta_{L} - h_{m}) 
e^{- i ( \Delta_{K} - h_{n} ) x} 
\biggr\}
\nonumber \\ 
&-&
(\Delta_{L} - \Delta_{K}) 
\biggl\{ (\Delta_{L} - h_{m}) (\Delta_{K} - h_{m}) 
e^{- i ( h_{k} - h_{n} ) x}
- (\Delta_{L} - h_{k}) (\Delta_{K} - h_{k}) 
e^{- i ( h_{m} - h_{n} ) x}
\biggr\}
\biggr]
\nonumber \\ 
&\times&
W_{\alpha L} (UX)^*_{\beta k} 
(UX)^*_{\alpha n} (UX)_{\beta n} 
\left\{ W^{\dagger} A (UX) \right\}_{L m} 
\left\{ (UX)^{\dagger} A W \right\}_{m K} 
\left\{ W^{\dagger} A (UX) \right\}_{K k} 
\biggr\}. 
\label{P-beta-alpha-W4-H3-Second} \end{aligned}$$

$$\begin{aligned}
&& P(\nu_\beta \rightarrow \nu_\alpha)^{(4)}_{6th} 
\equiv 
2 \mbox{Re} \left[ \left( 
S^{(0)}_{\alpha \beta} \right)^{*} 
S_{\alpha \beta}^{(4)} [2] 
\right]
\nonumber \\
&=& 
2 \mbox{Re} 
\biggl\{
- \sum_{n} 
\sum_{K} \sum_{k} 
\frac{ 1 }{ \Delta_{K} - h_{k} } 
\left[ (ix) e^{- i ( \Delta_{K} - h_{n} ) x} 
+ 
\frac{e^{- i ( h_{k} - h_{n} ) x} - e^{- i ( \Delta_{K} - h_{n} ) x} }{ h_{k} - \Delta_{K} } 
\right]
\nonumber \\
&\times& 
(UX)^*_{\alpha n} (UX)_{\beta n} 
W_{\alpha K} W^*_{\beta K} 
\left\{ W^{\dagger} A (UX) \right\}_{K k}  
\left\{ (UX)^{\dagger} A W \right\}_{k K} 
\nonumber \\
&+& 
\sum_{n} 
\sum_{K \neq L} 
\sum_{k} 
\frac{ 1 }{ ( \Delta_{L} - \Delta_{K} ) (\Delta_{L} - h_{k}) (\Delta_{K} - h_{k}) } 
\nonumber \\
&\times& 
\biggl[
\Delta_{L} e^{- i ( \Delta_{K} - h_{n} ) x} 
- \Delta_{K} e^{- i ( \Delta_{L} - h_{n} ) x} 
+ \left( e^{- i ( \Delta_{L} - h_{n} ) x} - e^{- i ( \Delta_{K} - h_{n} ) x} \right) h_{k} 
-  ( \Delta_{L} - \Delta_{K} ) 
e^{- i ( h_{k} - h_{n} ) x} 
\biggr]
\nonumber \\
&\times& 
(UX)^*_{\alpha n} (UX)_{\beta n} 
W_{\alpha K} W^*_{\beta L} 
\left\{ W^{\dagger} A (UX) \right\}_{K k} 
\left\{ (UX)^{\dagger} A W \right\}_{k L} 
\biggr\}. 
\label{P-beta-alpha-W4-H2} \end{aligned}$$

Neutrino evolution equation in flavor basis {#sec:flavor-basis-evolution}
===========================================

The Schrödinger equation takes the form with flavor basis Hamiltonian $H$ in (\[flavor-hamiltonian\]) $$\begin{aligned}
i \frac{d}{dx} 
\left[
\begin{array}{c}
\nu_{a} \\
\nu_{s} \\
\end{array}
\right] = 
\left[
\begin{array}{cc}
U {\bf \Delta_{a} } U^{\dagger} + W {\bf \Delta_{s} } W^{\dagger} + A & 
U {\bf \Delta_{a} } Z^{\dagger} + W {\bf \Delta_{s} } V^{\dagger}  \\
Z {\bf \Delta_{a} } U^{\dagger} + V {\bf \Delta_{s} } W^{\dagger} & 
Z {\bf \Delta_{a} } Z^{\dagger} + V {\bf \Delta_{s} } V^{\dagger} \\
\end{array}
\right]
\left[
\begin{array}{c}
\nu_{a} \\
\nu_{s} \\
\end{array}
\right] 
\label{Schroedinger-eq-flavor-basis}\end{aligned}$$ where $\nu_{a}$ ($\nu_{s}$) denotes $3$ ($N$) component vector in active (sterile) space. Apparently, the system depends not only on $U$ and $W$, but also on $Z$ and $V$ matrix elements, which is not the case in our treatment using the mass eigenstate basis.

Here, we show that the dependence on $Z$ and $V$ is superficial. Since there is no physical meaning of the particular basis for the sterile sector fermions we can redefine it by doing the transformation $$\begin{aligned}
\left[
\begin{array}{c}
\nu_{a} \\
\nu_{s} \\
\end{array}
\right] 
\rightarrow 
\left[
\begin{array}{cc}
1 & 0 \\
0 & Y \\
\end{array}
\right] 
\left[
\begin{array}{c}
\nu_{a} \\
\nu_{s} \\
\end{array}
\right] 
\equiv 
\left[
\begin{array}{c}
\nu_{a}' \\
\nu_{s}' \\
\end{array}
\right] 
$$ where $Y$ is a $N \times N$ unitary matrix. In the primed basis the Hamiltonian becomes $$\begin{aligned}
H' \equiv 
\left[
\begin{array}{cc}
1 & 0 \\
0 & Y \\
\end{array}
\right] 
H
\left[
\begin{array}{cc}
1 & 0 \\
0 & Y^{\dagger} \\
\end{array}
\right] =
\left[
\begin{array}{cc}
( U {\bf \Delta_{a} } U^{\dagger} + W {\bf \Delta_{s} } W^{\dagger} + A ) & 
( U {\bf \Delta_{a} } Z^{\dagger} + W {\bf \Delta_{s} } V^{\dagger} ) Y^{\dagger} \\
Y ( Z {\bf \Delta_{a} } U^{\dagger} + V {\bf \Delta_{s} } W^{\dagger} ) & 
Y ( Z {\bf \Delta_{a} } Z^{\dagger} + V {\bf \Delta_{s} } V^{\dagger} ) Y^{\dagger} \\
\end{array}
\right]. 
$$ We can arbitrarily choose $Y=V^{\dagger}$. Then, one can show by using unitarity (\[unitarity\]) that $$\begin{aligned}
H' = 
\left[
\begin{array}{cc}
( U {\bf \Delta_{a} } U^{\dagger} + W {\bf \Delta_{s} } W^{\dagger} + A ) & 
- U {\bf \Delta_{a} } U^{\dagger} W + W {\bf \Delta_{s} } ( {\bf 1} - W^{\dagger} W )\\
- W^{\dagger} U {\bf \Delta_{a} } U^{\dagger} + ( {\bf 1} - W^{\dagger} W ) {\bf \Delta_{s} } W^{\dagger} & 
W^{\dagger} U {\bf \Delta_{a} } 
U^{\dagger} W  + ( {\bf 1} - W^{\dagger} W ) {\bf \Delta_{s} } 
( {\bf 1} - W^{\dagger} W ) \end{array}
\right]. 
\nonumber \\ \end{aligned}$$ Therefore, our system depends only on $U$ and $W$, and the dependence on $Z$ and $V$ is superficial.

Universal scaling model of $N$ sterile sector  {#sec:scaling-model}
==============================================

Suppose that we obtain a particular parametrization of $U$ matrix by taking $N=1$ sterile sector, as we did in section \[sec:probabilities\]. In this $(3+1)$ model, the $W$ matrix elements are completely determined, up to phase, by unitarity $$\begin{aligned}
\vert W_{\alpha 4} \vert^2 = 1 - \sum_{j=1}^3 |U_{\alpha j}|^2. \end{aligned}$$ Now, we attempt to create a toy model of $N$ sterile sector by “universal scaling”. We postulate that all the $W$ matrix elements are real and equal: $$\begin{aligned}
W_{\alpha 4} = W_{\alpha 5} = \cdot \cdot \cdot  W_{\alpha N+3} =
\frac{ 1 }{ \sqrt{N} }
\biggl( 1 - \sum_{j=1}^3 |U_{\alpha j}|^2 \biggr)^{1/2}\end{aligned}$$ which is consistent with $(3+N)$ space unitarity. In this universal scaling model, the order $W^2$ correction terms in (\[P-beta-alpha-2nd-averaged\]) remains unchanged provided that we further assume that all the sterile masses are equal.[^24] It is because the $W$ matrix elements enter into the $W^2$ terms in the form $$\begin{aligned}
\sum_{K} W_{\alpha K} W^{\dagger}_{K \beta} 
\frac{ 1 }{ ( \Delta_{K} - h_{k} )^n  }, 
\label{W2-scaling} \end{aligned}$$ where $n=1$ or 2.

However, the leaking constant $\mathcal{C}_{\alpha \beta}$ becomes smaller by a factor of $N$ in the universal scaling model. In the $(3+1)$ model, $\mathcal{C}_{\alpha \beta}$ takes the largest value, the upper limit in eq. (\[Cab-bound\]). Because $\mathcal{C}_{\alpha \beta}$ is fourth order in $W$ it is evident that in the universal scaling model, $$\mathcal{C}_{\alpha \beta} = 
\frac{1}{N}
\biggl( 1 - \sum_{j=1}^3 |U_{\alpha j}|^2 \biggr) 
\biggl( 1 - \sum_{j=1}^3 |U_{\beta j}|^2 \biggr), 
\label{Cab-USM}$$ which is the lower limit of (\[Cab-bound\]).

[99]{}

Z. Maki, M. Nakagawa and S. Sakata, “Remarks on the unified model of elementary particles,” Prog. Theor. Phys.  [**28**]{} (1962) 870. doi:10.1143/PTP.28.870 Y. Fukuda [*et al.*]{} \[Super-Kamiokande Collaboration\], “Evidence for oscillation of atmospheric neutrinos,” Phys. Rev. Lett.  [**81**]{} (1998) 1562 doi:10.1103/PhysRevLett.81.1562 \[hep-ex/9807003\]. S. P. Mikheev and A. Y. Smirnov, “Resonance Amplification of Oscillations in Matter and Spectroscopy of Solar Neutrinos,” Sov. J. Nucl. Phys.  [**42**]{} (1985) 913 \[Yad. Fiz.  [**42**]{} (1985) 1441\]. L. Wolfenstein, “Neutrino Oscillations in Matter,” Phys. Rev. D [**17**]{} (1978) 2369. doi:10.1103/PhysRevD.17.2369 K. Eguchi [*et al.*]{} \[KamLAND Collaboration\], “First results from KamLAND: Evidence for reactor anti-neutrino disappearance,” Phys. Rev. Lett.  [**90**]{} (2003) 021802 doi:10.1103/PhysRevLett.90.021802 \[hep-ex/0212021\]. Q. R. Ahmad [*et al.*]{} \[SNO Collaboration\], “Direct evidence for neutrino flavor transformation from neutral current interactions in the Sudbury Neutrino Observatory,” Phys. Rev. Lett.  [**89**]{} (2002) 011301 doi:10.1103/PhysRevLett.89.011301 \[nucl-ex/0204008\]. B. T. Cleveland, T. Daily, R. Davis, Jr., J. R. Distel, K. Lande, C. K. Lee, P. S. Wildenhain and J. Ullman, “Measurement of the solar electron neutrino flux with the Homestake chlorine detector,” Astrophys. J.  [**496**]{} (1998) 505. doi:10.1086/305343 K. S. Hirata [*et al.*]{} \[Kamiokande-II Collaboration\], “Real time, directional measurement of B-8 solar neutrinos in the Kamiokande-II detector,” Phys. Rev. D [**44**]{} (1991) 2241 Erratum: \[Phys. Rev. D [**45**]{} (1992) 2170\]. doi:10.1103/PhysRevD.44.2241, 10.1103/PhysRevD.45.2170 W. Hampel [*et al.*]{} \[GALLEX Collaboration\], “GALLEX solar neutrino observations: Results for GALLEX IV,” Phys. Lett. B [**447**]{} (1999) 127. doi:10.1016/S0370-2693(98)01579-2 J. N. Abdurashitov [*et al.*]{} \[SAGE Collaboration\], “Solar neutrino flux measurements by the Soviet-American Gallium Experiment (SAGE) for half the 22 year solar cycle,” J. Exp. Theor. Phys.  [**95**]{} (2002) 181 \[Zh. Eksp. Teor. Fiz.  [**122**]{} (2002) 211\] doi:10.1134/1.1506424 \[astro-ph/0204245\]. S. Fukuda [*et al.*]{} \[Super-Kamiokande Collaboration\], “Solar B-8 and hep neutrino measurements from 1258 days of Super-Kamiokande data,” Phys. Rev. Lett.  [**86**]{} (2001) 5651 doi:10.1103/PhysRevLett.86.5651 \[hep-ex/0103032\]. B. Aharmim [*et al.*]{} \[SNO Collaboration\], “Combined Analysis of all Three Phases of Solar Neutrino Data from the Sudbury Neutrino Observatory,” Phys. Rev. C [**88**]{} (2013) 025501 doi:10.1103/PhysRevC.88.025501 \[arXiv:1109.0763 \[nucl-ex\]\]. F. P. An [*et al.*]{} \[Daya Bay Collaboration\], “Measurement of electron antineutrino oscillation based on 1230 days of operation of the Daya Bay experiment,” Phys. Rev. D [**95**]{} (2017) no.7, 072006 doi:10.1103/PhysRevD.95.072006 \[arXiv:1610.04802 \[hep-ex\]\]. J. H. Choi [*et al.*]{} \[RENO Collaboration\], “Observation of Energy and Baseline Dependent Reactor Antineutrino Disappearance in the RENO Experiment,” Phys. Rev. Lett.  [**116**]{} (2016) no.21, 211801 doi:10.1103/PhysRevLett.116.211801 \[arXiv:1511.05849 \[hep-ex\]\]. S. Schoppmann \[Double Chooz Collaboration\], “Latest results of Double Chooz,” PoS HQL [**2016**]{} (2017) 010.

K. Abe [*et al.*]{} \[T2K Collaboration\], “Measurement of neutrino and antineutrino oscillations by the T2K experiment including a new additional sample of $\nu_e$ interactions at the far detector,” arXiv:1707.01048 \[hep-ex\]. P. Adamson [*et al.*]{} \[NOvA Collaboration\], “First measurement of electron neutrino appearance in NOvA,” Phys. Rev. Lett.  [**116**]{} (2016) no.15, 151806 doi:10.1103/PhysRevLett.116.151806 \[arXiv:1601.05022 \[hep-ex\]\]. M. Hartz (for the T2K Collaboration), “T2K Neutrino Oscillation Results with Data up to 2017 Summer”, KEK Colloquium, August 4, 2017.

C. Patrignani [*et al.*]{} \[Particle Data Group\], “Review of Particle Physics,” Chin. Phys. C [**40**]{} (2016) no.10, 100001 and 2017 update. doi:10.1088/1674-1137/40/10/100001 C. S. Fong, H. Minakata and H. Nunokawa, “A framework for testing leptonic unitarity by neutrino oscillation experiments,” JHEP [**1702**]{} (2017) 114 doi:10.1007/JHEP02(2017)114 \[arXiv:1609.08623 \[hep-ph\]\].

A. E. Nelson and J. Walsh, “Short Baseline Neutrino Oscillations and a New Light Gauge Boson,” Phys. Rev. D [**77**]{} (2008) 033001 doi:10.1103/PhysRevD.77.033001 \[arXiv:0711.1363 \[hep-ph\]\]. M. Pospelov and J. Pradler, “Elastic scattering signals of solar neutrinos with enhanced baryonic currents,” Phys. Rev. D [**85**]{} (2012) 113016 Erratum: \[Phys. Rev. D [**88**]{} (2013) no.3, 039904\] doi:10.1103/PhysRevD.85.113016, 10.1103/PhysRevD.88.039904 \[arXiv:1203.0545 \[hep-ph\]\]. R. Harnik, J. Kopp and P. A. N. Machado, “Exploring nu Signals in Dark Matter Detectors,” JCAP [**1207**]{} (2012) 026 doi:10.1088/1475-7516/2012/07/026 \[arXiv:1202.6073 \[hep-ph\]\]. S. Parke and M. Ross-Lonergan, “Unitarity and the three flavor neutrino mixing matrix,” Phys. Rev. D [**93**]{} (2016) no.11, 113009 doi:10.1103/PhysRevD.93.113009 \[arXiv:1508.05095 \[hep-ph\]\]. M. Blennow, P. Coloma, E. Fernandez-Martinez, J. Hernandez-Garcia and J. Lopez-Pavon, “Non-Unitarity, sterile neutrinos, and Non-Standard neutrino Interactions,” JHEP [**1704**]{} (2017) 153 doi:10.1007/JHEP04(2017)153 \[arXiv:1609.08637 \[hep-ph\]\]. J. Tang, Y. Zhang and Y. F. Li, “Probing Direct and Indirect Unitarity Violation in Future Accelerator Neutrino Facilities,” Phys. Lett. B [**774**]{} (2017) 217 doi:10.1016/j.physletb.2017.09.055 \[arXiv:1708.04909 \[hep-ph\]\].

F. An [*et al.*]{} \[JUNO Collaboration\], “Neutrino Physics with JUNO,” J. Phys. G [**43**]{}, no. 3, 030401 (2016) doi:10.1088/0954-3899/43/3/030401 \[arXiv:1507.05613 \[physics.ins-det\]\]. K. Abe [*et al.*]{} \[Super-Kamiokande Collaboration\], “Atmospheric neutrino oscillation analysis with external constraints in Super-Kamiokande I-IV,” arXiv:1710.09126 \[hep-ex\]. R. Abbasi [*et al.*]{} \[IceCube Collaboration\], “The Design and Performance of IceCube DeepCore,” Astropart. Phys.  [**35**]{} (2012) 615 doi:10.1016/j.astropartphys.2012.01.004 \[arXiv:1109.6096 \[astro-ph.IM\]\]. K. Abe [*et al.*]{} \[Hyper-Kamiokande Proto- Collaboration\], “Physics potential of a long-baseline neutrino oscillation experiment using a J-PARC neutrino beam and Hyper-Kamiokande,” PTEP [**2015**]{} (2015) 053C02 doi:10.1093/ptep/ptv061 \[arXiv:1502.05199 \[hep-ex\]\]. R. Acciarri [*et al.*]{} \[DUNE Collaboration\], “Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE) : Volume 2: The Physics Program for DUNE at LBNF,” arXiv:1512.06148 \[physics.ins-det\]. M. G. Aartsen [*et al.*]{} \[IceCube Collaboration\], “PINGU: A Vision for Neutrino and Particle Physics at the South Pole,” J. Phys. G [**44**]{} (2017) no.5, 054006 doi:10.1088/1361-6471/44/5/054006 \[arXiv:1607.02671 \[hep-ex\]\]. S. Adrián-Martínez [*et al.*]{}, “Intrinsic limits on resolutions in muon- and electron-neutrino charged-current events in the KM3NeT/ORCA detector,” JHEP [**1705**]{} (2017) 008 doi:10.1007/JHEP05(2017)008 \[arXiv:1612.05621 \[physics.ins-det\]\]. S. Antusch, C. Biggio, E. Fernandez-Martinez, M. B. Gavela and J. Lopez-Pavon, “Unitarity of the Leptonic Mixing Matrix,” JHEP [**0610**]{} (2006) 084 doi:10.1088/1126-6708/2006/10/084 \[hep-ph/0607020\]. E. Fernandez-Martinez, M. B. Gavela, J. Lopez-Pavon and O. Yasuda, “CP-violation from non-unitary leptonic mixing,” Phys. Lett. B [**649**]{} (2007) 427 doi:10.1016/j.physletb.2007.03.069 \[hep-ph/0703098\]. S. Antusch, M. Blennow, E. Fernandez-Martinez and J. Lopez-Pavon, “Probing non-unitary mixing and CP-violation at a Neutrino Factory,” Phys. Rev. D [**80**]{} (2009) 033002 doi:10.1103/PhysRevD.80.033002 \[arXiv:0903.3986 \[hep-ph\]\]. S. Antusch, S. Blanchet, M. Blennow and E. Fernandez-Martinez, “Non-unitary Leptonic Mixing and Leptogenesis,” JHEP [**1001**]{} (2010) 017 doi:10.1007/JHEP01(2010)017 \[arXiv:0910.5957 \[hep-ph\]\]. S. Antusch and O. Fischer, “Non-unitarity of the leptonic mixing matrix: Present bounds and future sensitivities,” JHEP [**1410**]{} (2014) 094 doi:10.1007/JHEP10(2014)094 \[arXiv:1407.6607 \[hep-ph\]\]. F. J. Escrihuela, D. V. Forero, O. G. Miranda, M. Tórtola and J. W. F. Valle, “On the description of non-unitary neutrino mixing,” Phys. Rev. D [**92**]{} (2015) no.5, 053009 doi:10.1103/PhysRevD.92.053009 \[arXiv:1503.08879 \[hep-ph\]\]. E. Fernandez-Martinez, J. Hernandez-Garcia and J. Lopez-Pavon, “Global constraints on heavy neutrino mixing,” JHEP [**1608**]{} (2016) 033 doi:10.1007/JHEP08(2016)033 \[arXiv:1605.08774 \[hep-ph\]\]. F. J. Escrihuela, D. V. Forero, O. G. Miranda, M. Tórtola and J. W. F. Valle, “Probing CP violation with non-unitary mixing in long-baseline neutrino oscillation experiments: DUNE as a case study,” New J. Phys.  [**19**]{} (2017) no.9, 093005 doi:10.1088/1367-2630/aa79ec \[arXiv:1612.07377 \[hep-ph\]\].

N. Klop and A. Palazzo, “Imprints of CP violation induced by sterile neutrinos in T2K data,” Phys. Rev. D [**91**]{} (2015) no.7, 073017 doi:10.1103/PhysRevD.91.073017 \[arXiv:1412.7524 \[hep-ph\]\]. R. Gandhi, B. Kayser, M. Masud and S. Prakash, “The impact of sterile neutrinos on CP measurements at long baselines,” JHEP [**1511**]{} (2015) 039 doi:10.1007/JHEP11(2015)039 \[arXiv:1508.06275 \[hep-ph\]\]. S. K. Agarwalla, S. S. Chatterjee, A. Dasgupta and A. Palazzo, “Discovery Potential of T2K and NOvA in the Presence of a Light Sterile Neutrino,” JHEP [**1602**]{} (2016) 111 doi:10.1007/JHEP02(2016)111 \[arXiv:1601.05995 \[hep-ph\]\]. O. G. Miranda, M. Tortola and J. W. F. Valle, “New ambiguity in probing CP violation in neutrino oscillations,” Phys. Rev. Lett.  [**117**]{} (2016) no.6, 061804 doi:10.1103/PhysRevLett.117.061804 \[arXiv:1604.05690 \[hep-ph\]\]. S. F. Ge, P. Pasquini, M. Tortola and J. W. F. Valle, “Measuring the leptonic CP phase in neutrino oscillations with nonunitary mixing,” Phys. Rev. D [**95**]{} (2017) no.3, 033005 doi:10.1103/PhysRevD.95.033005 \[arXiv:1605.01670 \[hep-ph\]\].

D. Dutta and P. Ghoshal, “Probing CP violation with T2K, NO$\nu$A and DUNE in the presence of non-unitarity,” JHEP [**1609**]{} (2016) 110 doi:10.1007/JHEP09(2016)110 \[arXiv:1607.02500 \[hep-ph\]\]. D. Dutta, P. Ghoshal and S. Roy, “Effect of Non Unitarity on Neutrino Mass Hierarchy determination at DUNE, NO$\nu$A and T2K,” Nucl. Phys. B [**920**]{} (2017) 385 doi:10.1016/j.nuclphysb.2017.04.018 \[arXiv:1609.07094 \[hep-ph\]\]. J. Rout, M. Masud and P. Mehta, Phys. Rev. D [**95**]{}, no. 7, 075035 (2017) doi:10.1103/PhysRevD.95.075035 \[arXiv:1702.02163 \[hep-ph\]\].

H. Minakata and S. J. Parke, “Simple and Compact Expressions for Neutrino Oscillation Probabilities in Matter,” JHEP [**1601**]{} (2016) 180 doi:10.1007/JHEP01(2016)180 \[arXiv:1505.01826 \[hep-ph\]\]. H. Nunokawa, O. L. G. Peres and R. Zukanovich Funchal, “Probing the LSND mass scale and four neutrino scenarios with a neutrino telescope,” Phys. Lett. B [**562**]{} (2003) 279 doi:10.1016/S0370-2693(03)00603-8 \[hep-ph/0302039\]. K. Kimura, A. Takamura and H. Yokomakura, “Exact formulas and simple CP dependence of neutrino oscillation probabilities in matter with constant density,” Phys. Rev. D [**66**]{} (2002) 073005 doi:10.1103/PhysRevD.66.073005 \[hep-ph/0205295\]. J. Kopp, P. A. N. Machado, M. Maltoni and T. Schwetz, “Sterile Neutrino Oscillations: The Global Picture,” JHEP [**1305**]{} (2013) 050 doi:10.1007/JHEP05(2013)050 \[arXiv:1303.3011 \[hep-ph\]\]. A. de Gouvêa and A. Kobach, “Global Constraints on a Heavy Neutrino,” Phys. Rev. D [**93**]{}, no. 3, 033005 (2016) doi:10.1103/PhysRevD.93.033005 \[arXiv:1511.00683 \[hep-ph\]\]. M. G. Aartsen [*et al.*]{} \[IceCube Collaboration\], “Searches for Sterile Neutrinos with the IceCube Detector,” Phys. Rev. Lett.  [**117**]{} (2016) no.7, 071801 doi:10.1103/PhysRevLett.117.071801 \[arXiv:1605.01990 \[hep-ex\]\].

X. J. Xu, “Why is the neutrino oscillation formula expanded in $\Delta m^2_{21} / \Delta m^2_{31}$ still accurate near the solar resonance in matter?,” JHEP [**1510**]{} (2015) 090 doi:10.1007/JHEP10(2015)090 \[arXiv:1502.02503 \[hep-ph\]\]. S. F. Ge and A. Y. Smirnov, “Non-standard interactions and the CP phase measurements in neutrino oscillations at low energies,” JHEP [**1610**]{} (2016) 138 doi:10.1007/JHEP10(2016)138 \[arXiv:1607.08513 \[hep-ph\]\]. K. Asano and H. Minakata, “Large-Theta(13) Perturbation Theory of Neutrino Oscillation for Long-Baseline Experiments,” JHEP [**1106**]{} (2011) 022 doi:10.1007/JHEP06(2011)022 \[arXiv:1103.4387 \[hep-ph\]\].

Emiliano Rodríguez Mega, “Plans for a research powerhouse in the Andes begin to unravel”, Science (2017) doi:10.1126/science.aan7140.

[^1]: The unique citation of the solar neutrino measurement here must be understood as the representative of all the foregoing solar neutrino experiments [@Cleveland:1998nv; @Hirata:1991ub; @Hampel:1998xg; @Abdurashitov:2002nt; @Fukuda:2001nj; @Aharmim:2011vm].

[^2]: Recently, however, there exists accumulating indication that CP phase $\delta$ takes value around $\sim\frac{3\pi}{2}$ [@Hartz-KEK-colloquium].

[^3]: In recent years there exist a good amount of activities addressing the possibilities of new physics at low energies. For scenarios which involve light sterile neutrino(s), see, e.g., [@Nelson:2007yq; @Pospelov:2012gm; @Harnik:2012ni], and the references therein.

[^4]: Works have also been done on unitarity violation by sterile sector from somewhat different point of view, e.g., if it exists, how it could disturb measurement of lepton Kobayashi-Maskawa phase $\delta_{\rm CP}$, or mass ordering. See for example, [@Klop:2014ima; @Gandhi:2015xza; @Agarwalla:2016mrc; @Miranda:2016wdr; @Ge:2016xya; @Dutta:2016vcc; @Dutta:2016czj; @Rout:2017udo].

[^5]: To understand the formulas presented in this section, minimal explanation of definitions of the quantities may be helpful. Let the flavor mixing matrix $\bf{U}$ in $(3+N) \times (3+N)$ space that connect the flavor and mass eigenstates as $\nu_{\zeta} = {\bf U}_{\zeta z} \nu_{z}$, where $\zeta$ runs over active $\alpha = e,\mu,\tau$ and sterile flavor $s=s_1,\cdot \cdot \cdot,s_N$ indices, and $z$ over mostly active $i=1,2,3$ and mostly sterile mass eigenstate $J=4,5,\cdot \cdot \cdot, N+3$ indices. Let $U$ and $W$ be a part of $\bf{U}$ so that they connect the mass eigenstates to active neutrino flavor states, $\nu_{\alpha} = \sum_{i} (U)_{\alpha i} \nu_{i} + \sum_{J} (W)_{\alpha J} \nu_{J}$. The kinematical phase factors $\Delta_{j}$ and $\Delta_{J}$ are defined as $\Delta_{j} \equiv \frac{m^2_{j}}{2E}$ and $\Delta_{J} \equiv \frac{m^2_{J}}{2E}$, respectively. $m_{j}$ and $m_{J}$ denote the active and sterile neutrino masses, respectively.

[^6]: Speaking more precisely, we mean that all the $W$ matrix elements are assumed to be small, of the order of $\simeq 0.1$.

[^7]: Our result also means that taking care of all order matter effect does not change the feature obtained by first-order treatment in matter perturbation theory. That is, the same condition on sterile states masses derived by using first-order matter perturbation theory suffices to guarantee the absence of matter-dependent higher order correction terms in $W$. The reasons for this feature will be explained at the end of section \[sec:energy-denominator\].

[^8]: In matter, we have $\sum_{j=1}^{3} (UX)_{\alpha j} (UX)^{*}_{\beta j} 
    = \sum_{j,k,l=1}^{3} (U)_{\alpha k} (U)^{*}_{\beta l} X_{kj} X^*_{lj}
    = \sum_{k=1}^{3} (U)_{\alpha k} (U)^{*}_{\beta k}$ where in the last step, we have used the unitarity relation $\sum_{j=1}^3 X_{kj} X^*_{lj} = \delta_{kl}$.

[^9]: Through unitarity (\[unitarity\]), $U$ matrix elements have some dependence on $W$ matrix elements. We choose not to expand $U$ matrix elements by this $W$ dependence. In this sense, we use a “renormalized basis” (in the same sense as in ref. [@Minakata:2015gra]) in which some higher order effects are absorbed into the zeroth-order state.

[^10]: As in the Standard Model in particle physics T invariance is broken in our system only by complex numbers in the mixing matrix.

[^11]: Since $\hat{H}$ system is a consistent dynamical system it is legitimate and easier to verify generalized T invariance in the $\hat{S}$ level, though it can be done in the $S$ matrix level as well.

[^12]: To check the point of how the “averaging out the fast oscillation” procedure works, we numerically solved the $3+1$ system explicitly and confirmed that it does, as it should be.

[^13]: One can show that $L \sim \frac{ 2E }{ \Delta m^2_{31} } \sim \frac{1}{( h_{k} - h_{j} ) }$ holds at around the oscillation maximum of atmospheric-scale oscillations under the restriction we make for (\[suppression-cond\]). Then, the left-hand side of the first two inequalities receive an extra factor $LA \sim \frac{A}{( h_{k} - h_{j} ) } \sim \frac{ a }{ \Delta m^2_{31} } \simeq 0.1 \left(\frac{\rho}{2.8 \,\text{g/cm}^3}\right) \left(\frac{E}{1~\mbox{GeV}}\right)$, which further suppress it unless $\rho E \gsim 10\, \text{ (g/cm}^3) \text{GeV}$.

[^14]: We must remark, however, that this reasoning does not prove that the first order in matter perturbation theory is sufficient to obtain all the necessary conditions on the sterile state masses.

[^15]: As we remarked in footnote 7 the non-unitary mixing matrix $U$ has some $W$ dependence through unitarity of the ${\bf U}$ matrix in the whole $(3+N)$ space. Therefore, the nature of the eq. (\[Schroedinger-eq-0th\]) as the zeroth-order in $W$ is ambiguous. However, following [@Fong:2016yyh], we remain in the treatment with this “$W$ effect renormalized basis” in this paper.

[^16]: However, it appears that the flavor basis formulation of neutrino evolution in matter in high-scale unitarity violation poses highly nontrivial features such as non-Hermitian Hamiltonian [@Antusch:2006vwa], or the evolution equation $i \frac{d}{dx} \nu_{\alpha} = \sum_{j} \left[ U \left( {\bf \Delta_{a} } + U^{\dagger} A U \right) U^{\dagger} \right]_{\alpha \beta} \nu_{\beta}$ [@Escrihuela:2016ube]. The latter is not equivalent to (\[Schroedinger-eq-0th\]) in the vacuum mass eigenstate basis due to non-unitarity of the $U$ matrix.

[^17]: One can apply our formulas of $S$ matrix obtained under the constant matter density approximation to semi-realistic calculation for earth crossing neutrinos by using them in each shell (core, mantle, and crust regions, etc.) with proper connecting conditions at the boundaries.

[^18]: It can be re-parametrized in terms of the “$\alpha$ matrix parametrization” defined in ref. [@Escrihuela:2015wra]. The resultant values of $\alpha$ parameters are given as follows: $\alpha_{11} = 0.990$, $\alpha_{21} = - 0.0141$, $\alpha_{22} = 0.995$, $\alpha_{31} = -0.0445$, $\alpha_{32} = -0.0316$, $\alpha_{33} = 0.949$.

[^19]: This feature must be obvious if one goes back to the derivation of bound on $\mathcal{C}_{\alpha \beta}$ in [@Fong:2016yyh].

[^20]: We are aware that the assumption of equal sterile neutrino masses is contradictory to the assumption of no accidental degeneracy in the sterile mass spectrum we made in section \[sec:probability-2nd\]. It was done not to complicate term by term evaluation of the perturbative series, and to avoid using degenerate perturbation theory. Fortunately, we can remove this assumption to second order in $W$ in which no purely sterile sector energy denominator is involved.

[^21]: Of course, there is an issues of how to separate effects of $W^2$ correction terms from unitarity violation through $U$ matrix in leading order.

[^22]: It might be easier to obtain the phase factor if we use a different decomposition of $\tilde{H}$ from (\[tilde-H0+H1\]) by absorbing $W^{\dagger} A W$ into $\tilde{H}_{0}$.

[^23]: The phase itself needs not be small. Taking the matter potential of CC reaction and the earth diameter, $AL = 6.2 \left(\frac{\rho}{5 \text{g/cm}^3}\right) \left(\frac{L}{6,400 \mbox{km}}\right)$. Therefore, $ALW^2$ can be order unity for $|W| \simeq 0.4$.

[^24]: This statement applies also to the original expression (\[P-beta-alpha-0th+2nd\]).
