# Revisiting series expansions of neutrino oscillation and decay probabilities
in matter

Jesper Grönroos  Tommy Ohlsson  Sampsa Vihonen  Department of Physics, School
of Engineering Sciences, KTH Royal Institute of Technology,
AlbaNova University Center, Roslagstullsbacken 21, SE-106 91 Stockholm, Sweden
The Oskar Klein Centre, AlbaNova University Center, Roslagstullsbacken 21,
SE-106 91 Stockholm, Sweden<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present analytic expressions for three-flavor neutrino oscillations in
presence of invisible neutrino decay and matter effects. Using the well-known
Cayley–Hamilton formalism, the leading-order terms are derived for oscillation
probabilities in all major channels assuming the neutrino mass eigenstate
$\nu_{3}$ to decay into sterile neutrinos. Our work extends and complements
previous studies utilizing the Cayley–Hamilton theorem, providing the series
expansions for $\nu_{e}\rightarrow\nu_{e}$, $\nu_{e}\rightarrow\nu_{\mu}$,
$\nu_{e}\rightarrow\nu_{\tau}$, $\nu_{\mu}\rightarrow\nu_{\mu}$,
$\nu_{\mu}\rightarrow\nu_{\tau}$ and $\nu_{\tau}\rightarrow\nu_{\tau}$. The
accuracy of the analytical formulas is investigated by comparing the results
with numerically calculated probabilities. We also comment on the implications
on unitarity violation.

##  I Introduction

Neutrino decay was one of the alternative solutions that were proposed to
describe the observed zenith angle dependence on the early atmospheric
neutrino data in Super-Kamiokande in 1998 [1]. According to the decay
paradigm, neutrinos can decay into sterile neutrinos and some other particles
to achieve the muon neutrino disappearance. If neutrinos were to be Dirac
particles, the decay product could consist of a right-handed $SU(2)$ singlet
and a complex scalar $\chi$ with lepton number $-2$ and weak isospin and
hypercharge both equal to zero [2, 3]. In the case neutrinos were Majorana
particles, the decay product would include a Majoron $J$ instead [4, 5]. In
both cases, the decay products are invisible to a detector. While the
attention on neutrino decay has declined as a description of the solar
neutrino deficit, there are still models that predict invisible neutrino
decay. In such cases, neutrino decay leads to sub-leading corrections that
could be probed in future neutrino experiments.

Experimental research on neutrino oscillation physics has witnessed more than
two decades of uninterrupted success. After the discovery of neutrino
oscillations, numerous experiments probing neutrino oscillations with reactor,
solar, atmospheric and accelerator origin have maintained that neutrino
oscillations can be described with the standard oscillation parameters, three
mixing angles $\theta_{12}$, $\theta_{13}$ and $\theta_{23}$, two mass-squared
differences $\Delta m_{21}^{2}$ and $\Delta m_{31}^{2}$, and the Dirac _CP_
phase $\delta$. The remaining unknowns of the standard parameters include the
sign of $\Delta m_{31}^{2}$, which can follow either normal neutrino mass
ordering, $\Delta m_{31}^{2}>$ 0, or inverted ordering, $\Delta m_{31}^{2}<$
0\. Also unknown is the true value of $\delta$, which may be _CP_ -violating,
$\sin\delta\neq$ 0, or _CP_ -conserving, $\sin\delta=$ 0\. The increased
precision on the standard neutrino mixing parameters has also led to stringent
bounds on beyond the Standard Model physics models that predict neutrino
decay. The most constraining bounds on the invisible decay of the neutrino
mass states $m_{1}$ and $m_{2}$ arise from the combined analysis of solar
neutrino data and KamLAND reactor neutrino data [6], $\tau_{1}/m_{1}>$
4$\times$10-3 s/eV and $\tau_{2}/m_{2}>$ 7$\times$10-4 s/eV (95% C.L.), where
$\tau_{i}$ denotes the lifetime of mass state $\nu_{i}$, $i=$ 1, 2, 3. Global
analysis on long-baseline and atmospheric neutrino data [7] has set
furthermore a lower bound on the heaviest mass state with $\tau_{3}/m_{3}>$
2.9$\times$10-10 s/eV. Regarding the constraints on complete neutrino decay,
supernova neutrino data from SN1987 sets lower bound on electron neutrino
lifetime $\tau_{\nu_{e}}>$ 5.0$\times$10${}^{5}(m_{\nu_{e}}/\text{eV})$ [8],
which translates into $\tau_{i}/m_{i}\geq$ 1.2$\times$105 s/eV for mass states
$m_{1}$ and $m_{2}$ at 90% CL [9]. With the coming of the next-generation
neutrino telescopes like IceCube Gen-2 [10] and KM3NeT [11] and laboratory
experiments like DUNE [12], ESSnuSB [13] and HyperKamiokande [14], the limits
on invisible neutrino decay are expected to improve.

Approximating neutrino oscillation probabilities with neutrino decay is a non-
trivial matter. In the most general formalization of neutrino decay, the
effective Hamiltonian consists of Hermitian and anti-Hermitian parts, the
latter of which accounts for neutrino decay. As the eigenvalues of the
Hermitian and anti-Hermitian parts do not necessarily match, the two matrices
must be diagonalized independently. Compact analytic expressions for three-
flavor oscillations were derived in Ref. [15] and reviewed in Refs. [16, 17].
Compact analytic forms were first obtained for non-commuting Hermitian and
anti-Hermitian operators in two neutrino flavors [18] and later for three
flavors [19] using Zassenhaus expansions [20]. The neutrino oscillation
probabilities are currently known for $P_{ee}$ up to order
$\mathcal{O}(\lambda^{2})$ and for $P_{e\mu}$ to $\mathcal{O}(\lambda^{3})$,
where $\lambda\equiv$ 0.2. The survival probability $P_{\mu\mu}$ has also been
derived to zeroth-order.

In the present work, we provide the full set of neutrino oscillation
probabilities in the scenario where the heaviest neutrino mass state is able
to decay. Employing the method introduced in Ref. [21], we present series
expansions for the neutrino oscillation probabilities $P_{ee}$, $P_{e\mu}$,
$P_{e\tau}$, $P_{\mu\mu}$, $P_{\mu\tau}$ and $P_{\tau\tau}$ in terms of the
small parameters $\alpha\equiv\Delta m_{21}^{2}/\Delta m_{31}^{2}$ and
$s_{13}\equiv\sin\theta_{13}$. Matter effects are taken into account
throughout this work. The analytical expressions presented are given for
$P_{ee}$ up to order $\mathcal{O}(\lambda^{2})$, and for $P_{e\mu}$,
$P_{e\tau}$, $P_{\mu\mu}$, $P_{\mu\tau}$ and $P_{\tau\tau}$ to order
$\mathcal{O}(\lambda^{3})$, respectively. While the results derived for
$P_{ee}$ and $P_{e\mu}$ are in agreement with those presented in Ref. [19],
the other compact formulas are novel results.

This article is organized as follows. We review the phenomenology of invisible
neutrino decay in section II, discussing both generalized decay where all
neutrino mass states can decay and also the special case where only the
heaviest neutrino mass state decays. The perturbative methods that are used in
this work are also briefly described. Compact analytical formulas for
oscillation probabilities with neutrino decay are presented in section III,
giving a full account on matter effects arising from neutrino propagation in
matter. The oscillation probabilities are also obtained in the vacuum limit in
this section. We moreover present the unitarity relations in presence of
neutrino decay. The applications and limitations of these formulas are
discussed in section IV and concluding remarks are presented in section V. The
Cayley–Hamilton formalism and perturbative diagonalization are described in
appendices A and B, respectively.

##  II Methodology

The aim of this work is to obtain compact analytical expressions for neutrino
oscillation probabilities in presence of neutrino decay and matter effects. To
achieve this, we compute the time-evolution operator
$S=\text{e}^{i\mathcal{H}L}$ for effective decay Hamiltonian $\mathcal{H}$ and
neutrino propagation length $L$ and obtain probabilities for
$\nu_{\alpha}\rightarrow\nu_{\beta}$ as
$P_{\alpha\beta}=|S_{\beta\alpha}|^{2}$. In this section, we review the
effective decay Hamiltonian for general decay matrix $\Gamma$ and also for the
case where one active neutrino decays into sterile neutrino. We also describe
the perturbative methods that are used to obtain the compact analytical
formulas.

###  II.1 Effective decay Hamiltonian

Neutrino decay can be described with the effective Hamiltonian that consists
of Hermitian and anti-Hermitian parts:

$\mathcal{H}_{m}=H_{m}-\frac{i}{2}\Gamma_{m}.$ (1)

The Hermitian $H_{m}$ is defined in the mass basis and contains the dynamics
of standard neutrino oscillations. The anti-Hermitian component $\Gamma_{m}$
on the other hand is a $3\times 3$ matrix that contains the elements leading
to neutrino decay. In this case, the decay of neutrino mass state $\nu_{i}$ is
best characterized by the term exp[$(-m_{i}/\tau_{i})/(L/E)$], where
$\tau_{i}$ represents the lifetime while $L$ and $E$, respectively, give the
propagation distance and energy of the decaying neutrino state $\nu_{i}$.

The anti-Hermitian component in equation (1) can be written in
Weisskopf–Wigner approximation:

$\Gamma_{ij}=2\pi\sum_{k}\langle\nu_{i}|\mathcal{H}^{\prime}|\phi_{k}\rangle\langle\phi_{k}|\mathcal{H}^{\prime}|\nu_{j}\rangle\,\delta(E_{k}-E),$
(2)

where $|\phi_{k}\rangle$ is the final state representing decay products with
energy $E_{k}$. The interaction term $\mathcal{H}^{\prime}$ is responsible for
neutrino decay and it arises from beyond the Standard Model physics.

The difficulty of obtaining compact formulas with the effective Hamiltonian
$\mathcal{H}_{m}$ arises from the anti-Hermitian decay component $\Gamma_{m}$
in equation (1). Since neutrino mass and decay eigenstates do not necessarily
need to coincide, the Hermitian and anti-Hermitian components $H_{m}$ and
$\Gamma_{m}$ do not always commute and neutrino oscillation probabilities must
therefore be expressed in terms of $[H_{m},\Gamma_{m}]\neq$ 0\. Neutrino
interactions with matter further complicate the calculation, as vacuum
eigenstates of the Hamiltonian are generally different from the eigenstates in
matter.

When neutrinos are propagating in matter, the generalized Hamiltonian can be
written as follows

$\mathcal{H}_{f}(\Gamma)=U\left[\frac{\Delta
m_{31}^{2}}{2E}\begin{pmatrix}0&0&0\\\ 0&\Delta m_{21}^{2}&0\\\ 0&0&\Delta
m_{31}^{2}\end{pmatrix}-\frac{i}{2}\Gamma\right]U^{\dagger}+\begin{pmatrix}V_{CC}&0&0\\\
0&0&0\\\ 0&0&0\end{pmatrix}.$ (3)

Here $\mathcal{H}_{f}$ is given in flavor basis, $U$ is the
Pontecorvo–Maki–Nakagawa–Sakata matrix and $V_{CC}=\sqrt{2}G_{F}N_{e}$ is the
standard matter potential consisting of Fermi coupling constant $G_{f}$ and
electron density $N_{e}$, respectively.

It is often more convenient to focus on a special case where only one neutrino
mass state undergoes invisible decay. In the case of normal neutrino mass
ordering, the decaying mass state would be $\nu_{3}$. In such case it is
enough to take one decay parameter $\gamma$ and give effective Hamiltonian as
follows

$\mathcal{H}_{f}(\gamma)=\frac{\Delta
m_{31}^{2}}{2E}U\left[\begin{pmatrix}0&0&0\\\ 0&\alpha&0\\\
0&0&1\end{pmatrix}-i\begin{pmatrix}0&0&0\\\ 0&0&0\\\
0&0&\gamma\end{pmatrix}\right]U^{\dagger}+\begin{pmatrix}V_{CC}&0&0\\\
0&0&0\\\ 0&0&0\end{pmatrix},$ (4)

where $\alpha\equiv\Delta m_{21}^{2}/\Delta m_{31}^{2}$ is a small parameter.
The decay parameter $\gamma$ can then be written as $\gamma\equiv
m_{3}/(\tau_{3}\Delta m_{31}^{2})$ and it is always real.

###  II.2 Time-evolution operator

Neutrino oscillation probabilities are generally calculated with the time-
evolution operator $S\equiv\text{e}^{i\mathcal{H}_{f}L}$, where
$\mathcal{H}_{f}=\mathcal{H}_{f}(\gamma)$ is the effective Hamiltonian
depending on the decay parameter $\gamma$ and $L$ is the propagation distance
of neutrinos. Oscillation probabilities for channels
$\nu_{\alpha}\rightarrow\nu_{\beta}$ ($\alpha$, $\beta=e$, $\mu$, $\tau$) are
then obtained as the square of the time-evolution operator, that is,
$P_{\alpha\beta}=\left|S_{\beta\alpha}\right|^{2}$. We derive analytical
probabilities by calculating series expansions in terms of $\alpha$ and
$\sin\theta_{13}$. Such expansion can be carried out with $N$ terms with the
use of the Cayley–Hamilton formalism [22, 23, 24]. We provide more information
on the Cayley–Hamilton formalism in appendix A.

Working in the flavor basis, the time-evolution operator for neutrino
propagation can be written as

$\begin{split}S=\text{e}^{i{\mathcal{H}_{f}L}}&=\frac{\text{e}^{i{E_{1}}L}}{(E_{1}-E_{2})(E_{1}-E_{3})}\Big{[}E_{2}E_{3}\bm{I}-(E_{2}+E_{3})\mathcal{H}+\mathcal{H}^{2}\Big{]}\\\
&+\frac{\text{e}^{iE_{2}L}}{(E_{2}-E_{3})(E_{2}-E_{1})}\Big{[}E_{3}E_{1}\bm{I}-(E_{3}+E_{1})\mathcal{H}+\mathcal{H}^{2}\Big{]}\\\
&+\frac{\text{e}^{iE_{3}L}}{(E_{3}-E_{1})(E_{3}-E_{2})}\Big{[}E_{1}E_{2}\bm{I}-(E_{1}+E_{2})\mathcal{H}+\mathcal{H}^{2}\Big{]},\end{split}$
(5)

where $E_{i}$ are the eigenvalues of the Hamiltonian and $\bm{I}$ is the
identity matrix. For three-flavor oscillations, $S$ is a complex $3\times
3$-matrix. We derive the energy eigenvalues of equation (5) up to first order
in $\alpha$ and $s_{13}$:

$\displaystyle E_{1}$ $\displaystyle\simeq\frac{\Delta
m_{31}^{2}}{2E}\bigg{(}A+\alpha
s_{12}^{2}+s_{13}^{2}\frac{A(1-i\gamma)}{A-1+i\gamma}\bigg{)},$ (6)
$\displaystyle E_{2}$ $\displaystyle\simeq\frac{\Delta m_{31}^{2}}{2E}\alpha
c_{12}^{2},$ (7) $\displaystyle E_{3}$ $\displaystyle\simeq\frac{\Delta
m_{31}^{2}}{2E}\bigg{(}1-i\gamma-
s_{13}^{2}\frac{A(1-i\gamma)}{A-1+i\gamma}\bigg{)},$ (8)

where $A\equiv V_{CC}E/\Delta m_{31}^{2}$ is the standard matter potential
arising from neutrino interactions in matter and $\gamma$ is the decay
parameter in the simplest scenario where the heaviest neutrino decays into a
sterile neutrino. We find the eigenvalues presented in equations (6)–(8) to be
consistent with those those shown in Ref. [19]. We are also able to recover
the eigenvalues derived for the standard oscillations in Ref. [21] by applying
$\gamma\rightarrow 0$ and $A\rightarrow 0$.

In order to obtain compact analytical formulas for neutrino oscillations in
matter with invisible neutrino decay, we employ the Cayley–Hamilton formalism
to write series expansions using the expressions provided in equations (5) and
(6)–(8). For convenience, we apply the following simplifications

$\displaystyle B$ $\displaystyle\equiv A-1-\gamma^{2},$ (9) $\displaystyle C$
$\displaystyle\equiv A\gamma,$ (10) $\displaystyle D$
$\displaystyle\equiv(A-1)^{2}+\gamma^{2}.$ (11)

Throughout this work, we employ perturbative diagonalization to treat the
leptonic mixing matrix $U$. This method is used to obtain the eigenvalues
$E_{i}$ of the effective Hamiltonian. More details on this process can be
found in appendix B. We also perform multivariate series expansions using
small parameters $\alpha$ and $\sin\theta_{13}$. We carry out the multivariate
series expansions in `Mathematica`. For convenience, the results are expressed
in terms of the book-keeping parameter $\lambda\equiv 0.2$, which allows to
express the expansion parameters as $s_{13}\simeq
0.14\sim\mathcal{O}(\lambda)$ and $\alpha\simeq
0.03\sim\mathcal{O}(\lambda^{2})$.

Probabilities involving only $\nu_{\mu}$ and $\nu_{\tau}$ are more complicated
to derive than those involving $\nu_{e}$. In those cases, the expressions grow
so large that it is more beneficial to compute them as sum of terms to all
relevant orders of parameter $\lambda$. As such, we write probabilities
involving only $\nu_{\mu}$ and $\nu_{\tau}$ as

$P_{\alpha\beta}=P_{\alpha\beta}^{(0)}+P_{\alpha\beta}^{(1)}+P_{\alpha\beta}^{(2)}+P_{\alpha\beta}^{(3)}+\mathcal{O}(\lambda^{4}).$
(12)

##  III Oscillation probabilities with neutrino decay and matter effects

In this section, we present a complete set of series expansions for
oscillation probabilities with neutrino decay. The probabilities are obtained
with matter effects in section III.1. The expansions are also presented in the
vacuum limit $A\rightarrow 0$ and in unitarity relations, which are discussed
in sections III.2 and III.3, respectively.

###  III.1 Neutrino oscillation and decay probabilities in matter

We first present the full set of neutrino oscillation probabilities in
presence of neutrino decay and matter effects. The expansions are obtained in
the scenario where the heaviest light neutrino mass state $\nu_{3}$ decays
into the sterile state $\nu_{s}$ and neutrino masses follow the normal
ordering where $m_{1}<m_{2}\ll m_{3}$.

For $\nu_{e}\rightarrow\nu_{e}$ oscillations, we acquire the formula

$P_{ee}=1-\frac{2s_{13}^{2}}{D^{2}}\Bigg{\\{}\Big{(}B^{2}-C^{2}+2ACD\Delta\Big{)}-\bigg{[}\big{(}B^{2}-C^{2}\big{)}\cos[2(A-1)\Delta]+2BC\sin\big{[}2(A-1)\Delta\big{]}\bigg{]}\text{e}^{-2\gamma\Delta}\Bigg{\\}},$
(13)

where $\Delta\equiv\Delta m_{31}^{2}L/(4E)$ and $B$, $C$ and $D$ are functions
of the decay parameter $\gamma$ and the matter potential $A$ as given in
equations (9)–(11). As one can see, the expansion for $P_{ee}$ exists up to
$\mathcal{O}(\lambda^{2})$ and it is fully provided in equation (13). We note
that the result is in agreement with that of Ref. [19], where the probability
$P_{ee}$ was originally derived with the Cayley–Hamilton formalism.

In case of $\nu_{e}\rightarrow\nu_{\mu}$ oscillations, we obtain the compact
formula

$\begin{split}P_{e\mu}&=s_{13}^{2}s_{23}^{2}\frac{1+\gamma^{2}}{D}\Big{\\{}1-2\cos\big{[}2(A-1)\Delta\big{]}\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{\\}}\\\
&+\frac{4\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}}{D}\frac{\sin(A\Delta)}{A}\Bigg{\\{}B\bigg{[}\sin(A\Delta-\delta)+\sin\big{[}(A-2)\Delta+\delta\big{]}\text{e}^{-2\gamma\Delta}\bigg{]}\\\
&+C\bigg{[}\cos(A\Delta-\delta)-\cos\big{[}(A-2)\Delta+\delta\big{]}\text{e}^{-2\gamma\Delta}\bigg{]}\Bigg{\\}},\end{split}$
(14)

where the probability is provided up to $\mathcal{O}(\lambda^{3})$. Again, we
find the result to be in agreement with the one derived in Ref. [19], where
the expansion for $P_{e\mu}$ was originally derived using the Cayley–Hamilton
formalism.

As the first novel result of this work, we provide the series expansion for
$\nu_{e}\rightarrow\nu_{\tau}$ oscillations up to $\mathcal{O}(\lambda^{3})$
as

$\displaystyle\begin{split}P_{e\tau}&=s_{13}^{2}c_{23}^{2}\frac{1+\gamma^{2}}{D}\Big{\\{}1-2\cos\big{[}2(A-1)\Delta\big{]}\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{\\}}\\\
&-\frac{4\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}}{D}\frac{\sin(A\Delta)}{A}\Bigg{\\{}B\bigg{[}\sin(A\Delta-\delta)+\sin\big{[}(A-2)\Delta+\delta\big{]}\text{e}^{-2\gamma\Delta}\bigg{]}\\\
&+C\bigg{[}\cos(A\Delta-\delta)-\cos\big{[}(A-2)\Delta+\delta\big{]}\text{e}^{-2\gamma\Delta}\bigg{]}\Bigg{\\}}.\end{split}$
(15)

The result presented in equation (15) has not been presented elsewhere.

In the case of $\nu_{\mu}\rightarrow\nu_{\mu}$, the probability has been
previously computed to zeroth order in Ref. [19]. We reproduce their result by
obtaining the following expression

$P_{\mu\mu}^{(0)}=1-s_{23}^{2}\Big{(}1-\text{e}^{-4\gamma\Delta}\Big{)}-c_{23}^{2}s_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}.$
(16)

We additionally present the higher-order terms for
$\nu_{\mu}\rightarrow\nu_{\mu}$ up to $\mathcal{O}(\lambda^{3})$. They are
given by

$\begin{split}P_{\mu\mu}^{(2)}&=\frac{2s_{13}^{2}s_{23}^{2}}{D^{2}}\bigg{\\{}\bigg{[}\left(B^{2}-C^{2}\right)\cos(2A\Delta)-2BC\sin(2A\Delta)\bigg{]}c_{23}^{2}\\\
&+\bigg{[}\left(B^{2}-C^{2}\right)\cos\big{[}2(A-1)\Delta\big{]}-2BC\sin\big{[}2(A-1)\Delta\big{]}\bigg{]}s_{23}^{2}\text{e}^{-2\gamma\Delta}\\\
&-\bigg{[}\cos(2\Delta)c_{23}^{2}+s_{23}^{2}\text{e}^{-2\gamma\Delta}\bigg{]}\bigg{(}B^{2}-C^{2}-2ACD\Delta\bigg{)}\text{e}^{-2\gamma\Delta}\\\
&+2AB\big{(}\gamma+D\Delta\big{)}\sin(2\Delta)c_{23}^{2}\text{e}^{-2\gamma\Delta}\bigg{\\}}+4\alpha
c_{12}^{2}c_{23}^{2}s_{23}^{2}\Delta\sin(2\Delta)\text{e}^{-2\gamma\Delta}\end{split}$
(17)

and

$\begin{split}P_{\mu\mu}^{(3)}&=\frac{4\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}\cos\delta}{AD}\Bigg{\\{}A\big{(}1-A\big{)}s_{23}^{2}\text{e}^{-4\gamma\Delta}+D\bigg{[}c_{23}^{2}-\cos(2\theta_{23})\cos(2\Delta)\text{e}^{-2\gamma\Delta}\bigg{]}\\\
&+B\bigg{[}\cos(2A\Delta)c_{23}^{2}-\bigg{(}\cos(2\Delta)c_{23}^{2}-\cos\big{[}2(A-1)\Delta\big{]}s_{23}^{2}\bigg{)}\text{e}^{-2\gamma\Delta}\bigg{]}\\\
&-C\bigg{[}\sin(2A\Delta)c_{23}^{2}-\bigg{(}\sin(2\Delta)c_{23}^{2}+\sin\big{[}2(A-1)\Delta\big{]}s_{23}^{2}\bigg{)}\text{e}^{-2\gamma\Delta}\bigg{]}\Bigg{\\}}.\end{split}$
(18)

The first-order term $P_{\mu\mu}^{(1)}$ is found to be zero and is therefore
omitted. The higher-order terms shown in equations (17) and (18) are presented
for the first time in this work.

For $\nu_{\mu}\rightarrow\nu_{\tau}$ oscillations, the zeroth-order term is

$P_{\mu\tau}^{(0)}=c_{23}^{2}s_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}.$
(19)

The first-order term is again zero and not shown here. To second-order, we
reach

$\begin{split}P_{\mu\tau}^{(2)}&=-\frac{2s_{13}^{2}c_{23}^{2}s_{23}^{2}}{D^{2}}\bigg{\\{}\bigg{[}\left(B^{2}-C^{2}\right)\cos(2A\Delta)-2BC\sin(2A\Delta)\bigg{]}\\\
&-\bigg{[}\left(B^{2}-C^{2}\right)\cos\big{[}2(A-1)\Delta\big{]}-2BC\sin\big{[}2(A-1)\Delta\big{]}\bigg{]}\text{e}^{-2\gamma\Delta}\\\
&-\bigg{[}\cos(2\Delta)-\text{e}^{-2\gamma\Delta}\bigg{]}\bigg{(}B^{2}-C^{2}-2ACD\Delta\bigg{)}\text{e}^{-2\gamma\Delta}\\\
&+2AB\big{(}\gamma+D\Delta\big{)}\sin(2\Delta)\text{e}^{-2\gamma\Delta}\bigg{\\}}-4\alpha
c_{12}^{2}c_{23}^{2}s_{23}^{2}\Delta\sin(2\Delta)\text{e}^{-2\gamma\Delta}.\end{split}$
(20)

The third-order term in $P_{\mu\tau}$ is more cumbersome. It can be written as

$\begin{split}P_{\mu\tau}^{(3)}&=\frac{2\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}}{AD}\bigg{\\{}\bigg{[}A\big{(}1-A\big{)}\text{e}^{-4\gamma\Delta}-D\Big{(}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}\Big{)}\bigg{]}\cos(2\theta_{23})\cos\delta-C\text{e}^{-4\gamma\Delta}\sin\delta\\\
&-B\bigg{[}\bigg{(}\cos(2A\Delta)-\Big{[}\cos(2\Delta)+\cos\big{[}2(A-1)\Delta\big{]}\Big{]}\text{e}^{-2\gamma\Delta}\bigg{)}\cos(2\theta_{23})\cos\delta\\\
&+\bigg{(}\sin(2A\Delta)-\Big{[}\sin(2\Delta)+\sin\big{[}2(A-1)\Delta\big{]}\Big{]}\text{e}^{-2\gamma\Delta}\bigg{)}\sin\delta\bigg{]}\\\
&+C\bigg{[}\bigg{(}\sin(2A\Delta)-\Big{[}\sin(2\Delta)+\sin\big{[}2(A-1)\Delta\big{]}\Big{]}\text{e}^{-2\gamma\Delta}\bigg{)}\cos(2\theta_{23})\cos\delta\\\
&-\bigg{(}\cos(2A\Delta)-\Big{[}\cos(2\Delta)+\cos\big{[}2(A-1)\Delta\big{]}\Big{]}\text{e}^{-2\gamma\Delta}\bigg{)}\sin\delta\bigg{]}\bigg{\\}}.\end{split}$
(21)

Equations (19)–(21) show the probability for $\nu_{\mu}\rightarrow\nu_{\tau}$
transition up to $\mathcal{O}(\lambda^{3})$. The expansions are novel results
and they have not been presented elsewhere.

We finally obtain the expansions for the $\nu_{\tau}\rightarrow\nu_{\tau}$
channel. This probability has not been presented elsewhere, and we show the
result for the first time. The zeroth-order term can be written as

$\begin{split}P_{\tau\tau}^{(0)}&=1-c_{23}^{2}\Big{(}1-\text{e}^{-4\gamma\Delta}\Big{)}-c_{23}^{2}s_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]},\end{split}$
(22)

whereas the second-order term is

$\begin{split}P_{\tau\tau}^{(2)}&=\frac{2s_{13}^{2}c_{23}^{2}}{D^{2}}\bigg{\\{}\bigg{[}\left(B^{2}-C^{2}\right)\cos(2A\Delta)-2BC\sin(2A\Delta)\bigg{]}s_{23}^{2}\\\
&+\bigg{[}\left(B^{2}-C^{2}\right)\cos\big{[}2(A-1)\Delta\big{]}-2BC\sin\big{[}2(A-1)\Delta\big{]}\bigg{]}c_{23}^{2}\text{e}^{-2\gamma\Delta}\\\
&-\bigg{[}\cos(2\Delta)s_{23}^{2}+c_{23}^{2}\text{e}^{-2\gamma\Delta}\bigg{]}\bigg{(}B^{2}-C^{2}-2ACD\Delta\bigg{)}\text{e}^{-2\gamma\Delta}\\\
&+2AB\big{(}\gamma+D\Delta\big{)}\sin(2\Delta)s_{23}^{2}\text{e}^{-2\gamma\Delta}\bigg{\\}}+4\alpha
c_{12}^{2}c_{23}^{2}s_{23}^{2}\Delta\sin(2\Delta)\text{e}^{-2\gamma\Delta}\end{split}$
(23)

and the third-order one is

$\begin{split}P_{\tau\tau}^{(3)}&=\frac{-4\alpha
s_{13}c_{12}^{2}s_{12}^{2}c_{23}^{2}s_{23}^{2}\cos\delta}{AD}\bigg{\\{}A\big{(}1-A\big{)}c_{23}^{2}\text{e}^{-4\gamma\Delta}+D\Big{(}s_{23}^{2}+\cos(2\theta_{23})\cos(2\Delta)\text{e}^{-2\gamma\Delta}\Big{)}\\\
&+B\bigg{[}\cos(2A\Delta)s_{23}^{2}-\bigg{(}\cos(2\Delta)s_{23}^{2}-\cos\big{[}2(A-1)\Delta\big{]}c_{23}^{2}\bigg{)}\text{e}^{-2\gamma\Delta}\bigg{]}\\\
&-C\bigg{[}\sin(2A\Delta)s_{23}^{2}-\bigg{(}\sin(2\Delta)s_{23}^{2}-\sin\big{[}2(A-1)\Delta\big{]}c_{23}^{2}\bigg{)}\text{e}^{-2\gamma\Delta}\bigg{]}\bigg{\\}},\end{split}$
(24)

respectively. The first-order term is absent in $P_{\tau\tau}$ and is not
presented here.

###  III.2 Vacuum limits of the neutrino oscillation and decay probabilities

We now obtain the vacuum limits for the neutrino oscillation probabilities
with invisible neutrino decay. The vacuum limits are acquired by applying
$A\rightarrow 0$ in the series expansions presented in equations (14)–(24).

For the $\nu_{e}\rightarrow\nu_{e}$ case, the vacuum probability reads as

$P_{ee}^{\text{vac}}=1-2s_{13}^{2}\Big{[}1-\cos(2\Delta)\text{e}^{-2\gamma\Delta}\Big{]}.$
(25)

The probabilities for $\nu_{e}\rightarrow\nu_{\mu}$ and
$\nu_{e}\rightarrow\nu_{\tau}$ channels are similarly given by the following
equations

$P_{e\mu}^{\text{vac}}=s_{13}^{2}s_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}+4\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}\Delta\Big{[}\sin(2\Delta-\delta)+\text{e}^{-2\gamma\Delta}\sin\delta\Big{]},$
(26)
$P_{e\tau}^{\text{vac}}=s_{13}^{2}c_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}-4\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}\Delta\Big{[}\sin(2\Delta-\delta)+\text{e}^{-2\gamma\Delta}\sin\delta\Big{]}.$
(27)

For $\nu_{\mu}\rightarrow\nu_{\mu}$ and $\nu_{\mu}\rightarrow\nu_{\tau}$
channels, we arrive at

$\begin{split}P_{\mu\mu}^{\text{vac}}&=1-s_{23}^{2}\Big{(}1-\text{e}^{-4\gamma\Delta}\Big{)}-c_{23}^{2}s_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}+4\alpha
c_{12}^{2}c_{23}^{2}s_{23}^{2}\Delta\sin(2\Delta)\text{e}^{-2\gamma\Delta}\\\
&+s_{13}^{2}s_{23}^{2}\Big{\\{}1-\text{e}^{-4\gamma\Delta}+\cos(2\theta_{23})\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}\Big{\\}}\\\
&-\frac{2\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}\cos\delta}{(1+\gamma^{2})}\bigg{\\{}1-\text{e}^{-4\gamma\Delta}+\cos(2\theta_{23})\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}\\\
&-2\sin(2\Delta)\text{e}^{-2\gamma\Delta}\Big{[}\cos(2\theta_{23})\Big{(}\big{(}1+\gamma^{2}\big{)}\Delta+\gamma\Big{)}-1\Big{]}\bigg{\\}},\end{split}$
(28)
$\begin{split}P_{\mu\tau}^{\text{vac}}&=c_{23}^{2}s_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}-4\alpha
c_{12}^{2}c_{23}^{2}s_{23}^{2}\Delta\sin(2\Delta)\text{e}^{-2\gamma\Delta}\\\
&-2s_{13}^{2}c_{23}^{2}s_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}\\\
&+\frac{2\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}}{(1+\gamma^{2})}\bigg{\\{}\sin\delta\Big{(}1-\text{e}^{-2\gamma\Delta}\Big{)}\Big{[}2\Delta(1+\gamma^{2})-\gamma\Big{(}1-\text{e}^{-2\gamma\Delta}\Big{)}\Big{]}\\\
&-\cos(2\theta_{23})\cos\delta\Big{[}\Big{(}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{)}+2\big{(}1+\gamma^{2}\big{)}\Delta\sin(2\Delta)\text{e}^{-2\gamma\Delta}\Big{]}\bigg{\\}}.\end{split}$
(29)

We finally present the vacuum limit for the $\nu_{\tau}\rightarrow\nu_{\tau}$
channel

$\begin{split}P_{\tau\tau}^{\text{vac}}&=1-c_{23}^{2}\Big{(}1-\text{e}^{-4\gamma\Delta}\Big{)}-c_{23}^{2}s_{23}^{2}\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}+4\alpha
c_{12}^{2}c_{23}^{2}s_{23}^{2}\Delta\sin(2\Delta)\text{e}^{-2\gamma\Delta}\\\
&+s_{13}^{2}c_{23}^{2}\Big{\\{}1-\text{e}^{-4\gamma\Delta}+\cos(2\theta_{23})\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}\Big{\\}}\\\
&+\frac{2\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}\cos\delta}{(1+\gamma^{2})}\bigg{\\{}1-\text{e}^{-4\gamma\Delta}+\cos(2\theta_{23})\Big{[}1-2\cos(2\Delta)\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}\\\
&+2\sin(2\Delta)\text{e}^{-2\gamma\Delta}\Big{[}2c_{23}^{2}\big{(}1+\gamma^{2}\big{)}\Delta-\gamma\Big{]}\bigg{\\}}.\end{split}$
(30)

###  III.3 Unitarity relations with neutrino decay

Including invisible decay into neutrino oscillation probabilities lead to
violation of unitarity. In this section, we derive the unitarity relations for
the $e$, $\mu$ and $\tau$ channels. The results are presented in equations
(31)–(33):

$\sum_{\beta}P_{e\beta}=1+\frac{s_{13}^{2}}{D^{2}}\bigg{\\{}3C^{2}-B^{2}-4AC(D-2)\Delta+(1+\gamma^{2})D\text{e}^{-4\gamma\Delta}+4C\text{e}^{-2\gamma\Delta}\Big{[}B\sin\big{[}2(A-1)\Delta]-C\cos\big{[}2(A-1)\Delta]\Big{]}\bigg{\\}},$
(31)

$\begin{split}\sum_{\beta}P_{\mu\beta}&=1-\big{(}1-\text{e}^{-4\gamma\Delta}\big{)}s_{23}^{2}-\frac{s_{13}^{2}s_{23}^{2}}{D^{2}}\bigg{\\{}\big{(}1-\text{e}^{-4\gamma\Delta}\big{)}(1+\gamma^{2})(A+B)-A^{2}\Big{[}(1+\gamma^{2})-(1-3\gamma^{2})\text{e}^{-4\gamma\Delta}\Big{]}\\\
&-4ACD\Delta\text{e}^{-4\gamma\Delta}+4C\text{e}^{-2\gamma\Delta}\Big{[}B\sin\big{[}2(A-1)\Delta\big{]}+C\cos\big{[}2(A-1)\Delta\big{]}\Big{]}\bigg{\\}}\\\
&-\frac{2\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}}{AD}\bigg{\\{}C\Big{[}1-2\cos\big{[}2(A-1)\Delta\big{]}\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}\sin\delta\\\
&-\Big{[}A(A-1)\big{(}1-\text{e}^{-4\gamma\Delta}\big{)}-2C\cos(2\theta_{23})\sin\big{[}2(A-1)\Delta\big{]}\text{e}^{-2\gamma\Delta}\Big{]}\cos\delta\bigg{\\}}\end{split}$
(32)

and

$\begin{split}\sum_{\beta}P_{\tau\beta}&=1-\big{(}1-\text{e}^{-4\gamma\Delta}\big{)}c_{23}^{2}-\frac{s_{13}^{2}c_{23}^{2}}{D^{2}}\bigg{\\{}\big{(}1-\text{e}^{-4\gamma\Delta}\big{)}(1+\gamma^{2})(A+B)-A^{2}\Big{[}(1+\gamma^{2})-(1-3\gamma^{2})\text{e}^{-4\gamma\Delta}\Big{]}\\\
&-4ACD\Delta\text{e}^{-4\gamma\Delta}+4C\text{e}^{-2\gamma\Delta}\Big{[}B\sin\big{[}2(A-1)\Delta\big{]}+C\cos\big{[}2(A-1)\Delta\big{]}\Big{]}\bigg{\\}}\\\
&+\frac{2\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}}{AD}\bigg{\\{}C\Big{[}1-2\cos\big{[}2(A-1)\Delta\big{]}\text{e}^{-2\gamma\Delta}+\text{e}^{-4\gamma\Delta}\Big{]}\sin\delta\\\
&-\Big{[}A(A-1)\big{(}1-\text{e}^{-4\gamma\Delta}\big{)}-2C\sin\big{[}2(A-1)\Delta\big{]}\text{e}^{-2\gamma\Delta}\Big{]}\cos\delta\bigg{\\}}.\end{split}$
(33)

The transposed relations $\sum_{\alpha}P_{\alpha e}$,
$\sum_{\alpha}P_{\alpha\mu}$ and $\sum_{\alpha}P_{\alpha\tau}$ can be obtained
by applying the transformation $\delta\rightarrow-\delta$ in equations
(31)–(33), respectively. We can immediately see that the unitarity relations
reduce to unity in both cases when $\gamma\rightarrow 0$.

##  IV Discussion

In this section, we discuss the qualitative effects that occur for the
neutrino oscillation probabilities in presence of neutrino decay. We begin by
obtaining the probability formulas presented in section III at the no-decay
limit and comparing the results with those derived in Ref. [21]. We also
investigate the violation of unitarity and estimate the accuracy of the
expansion formulas presented in this work.

###  IV.1 No-decay limit

At the no-decay limit $\gamma\rightarrow 0$, the analytical probabilities
presented in this work are expected to follow the standard three-flavor
paradigm. The auxiliary parameters $B$, $C$ and $D$ introduced in equations
(9)–(11) then read as

$\displaystyle B_{\gamma\to 0}$ $\displaystyle=A-1,$ (34) $\displaystyle
C_{\gamma\to 0}$ $\displaystyle=0,$ (35) $\displaystyle D_{\gamma\to 0}$
$\displaystyle=\left(A-1\right)^{2}.$ (36)

One can readily see that $B$, $C$ and $D$ are finite at the no-decay limit.
The vanishing parameter $C$ is furthermore never present in denominators and
does not induce divergences in the probability formulas.

We take a moment to compare our results with the analytical probabilities that
were derived in Ref. [21]. As the authors of Ref. [21] used the same
methodology, that is, obtaining analytical expressions for the neutrino
oscillation probabilities via perturbative diagonalization in Cayley–Hamilton
formalism, it is expected that the results derived in our work agree with
those presented in Ref. [21]. Indeed, we find that the expansions we presented
for the probabilities $P_{ee}$, $P_{e\mu}$ and $P_{e\tau}$ in equations (13),
(14) and (15), respectively, are consistent with the probability formulas
discussed in Ref. [21]. Our findings for the probabilities $P_{\mu\mu}$,
$P_{\mu\tau}$ and $P_{\tau\tau}$ in equations (16)–(18), (19)–(21) and
(22)–(24) on the other hand are not in perfect agreement with the expansions
presented in Ref. [21]. When the results for $P_{\mu\mu}$, $P_{\mu\tau}$ and
$P_{\tau\tau}$ are compared to Ref. [21], we find that a third-order term is
present in all three probabilities. It is given by

$\Delta P_{\alpha\beta}^{\text{no decay}}=\alpha
s_{13}\sin(2\theta_{12})\sin(4\theta_{23})\sin^{2}\Delta\cos\delta=8\alpha
s_{13}c_{12}s_{12}c_{23}s_{23}\cos(2\theta_{23})\sin^{2}\Delta\cos\delta,$
(37)

where $\alpha,\beta=\mu,\tau$. This additional term is not present in Ref.
[21] and it originates from expanding to lower order in $\alpha$ while
applying perturbative diagonalization to the effective Hamiltonian. When we
take $\alpha\sim\mathcal{O}(\lambda^{2})$, the terms involving $\alpha^{2}$
are disregarded in the energy eigenvalues (6)–(8). The first energy eigenvalue
for example reads

$E_{1}=\frac{\Delta m_{31}^{2}}{2E}\bigg{[}A+\alpha
s_{12}^{2}+s_{13}^{2}\frac{A(1-i\gamma)}{A-1+i\gamma}\bigg{]}+\mathcal{O}(\lambda^{4}),$
(38)

whereas the corresponding one in Ref. [21] is given by

$E_{1}\simeq\frac{\Delta m_{31}^{2}}{2E}\left[A+\alpha
s_{12}^{2}+s_{13}^{2}\frac{A}{A-1}+\alpha^{2}\frac{\sin^{2}(2\theta_{12})}{4A}\right].$
(39)

The fact that the term proportional to $\alpha^{2}$ is present in equation
(39), but not in equation (38) is not surprising. It is a trade-off that
results from simplifying the complicated expressions that arise from
calculating neutrino oscillation probabilities with neutrino decay. By
expanding the exponential as a function of $\alpha^{2}$ in energy eigenvalues
$E_{1}$, $E_{2}$ and $E_{3}$, we obtain additional terms that are linear in
$\alpha$ from the chain rule. In some instances these extra terms multiply
$s_{13}$ to generate third-order terms. This scenario was not taken into
account in the derivation of the probability formulas in the present work and
we consider it the most probable source to the mixed term $\alpha s_{13}$ that
was obtained for the energy eigenvalues in Ref. [21]. This is also likely the
case in Ref. [19], where authors reached the same energy eigenvalues as those
we presented in equations (6)–(8). We estimate the magnitude of equation (37)
to be $10^{-3}$ using the experimentally determined values for the neutrino
oscillation parameters [25].

Some expansions of neutrino oscillation probabilities with neutrino decay have
been performed in vacuum. An example of such publication is Ref. [26], where
authors present vacuum formulas for $P_{\mu e}$ and $P_{\mu\mu}$ using the so-
called one-mass-squared-difference (OMSD) approximation. In this scenario,
expansions are carried out assuming $\Delta m_{12}^{2}=0$. We find that our
results for $P_{\mu e}$ and $P_{\mu\mu}$ are in agreement with those derived
in Ref. [26] in the vacuum limit $A\rightarrow 0$.

###  IV.2 Violation of unitarity

In the absence of neutrino decay, the oscillation probabilities for neutrinos
satisfy the unitarity relations

$\sum_{\alpha}P_{\alpha\beta}=\sum_{\beta}P_{\alpha\beta}=1,$ (40)

where $\alpha$ and $\beta$ run over the flavor states $e,\mu$ and $\tau$. This
ensures that each neutrino must either remain in its original state or undergo
a transition to another flavor state.

In perturbative expansion in small parameters, the unitarity relations shown
in equation (40) are satisfied up to a given order. It has been shown that the
number of independent probabilities can be reduced to five [21]. When neutrino
decay is allowed, however, the evolution is no longer unitary. Equations
(41)–(43) show the corresponding unitarity relations when probabilities
derived under neutrino decay are included:

$\displaystyle\sum_{\beta}P_{e\beta}$
$\displaystyle=1+\mathcal{O}(\lambda^{2}),$ (41)
$\displaystyle\sum_{\beta}P_{\mu\beta}$
$\displaystyle=1-\big{(}1-\text{e}^{-4\gamma\Delta}\big{)}s_{23}^{2}+\mathcal{O}(\lambda^{2}),$
(42) $\displaystyle\sum_{\beta}P_{\tau\beta}$
$\displaystyle=1-\big{(}1-\text{e}^{-4\gamma\Delta}\big{)}c_{23}^{2}+\mathcal{O}(\lambda^{2}).$
(43)

It is observed that all three equations display corrections of order
$\mathcal{O}(\lambda^{2})$. The sums in the second and third equations (42)
and (43) additionally contain the zeroth-order corrections
$(1-\text{e}^{-4\gamma\Delta})s_{23}^{2}$ and
$(1-\text{e}^{-4\gamma\Delta})c_{23}^{2}$, respectively.

In the Cayley–Hamilton formalism, the reversed transition probabilities are
obtained by simply replacing $\delta\rightarrow-\delta$, whereby
$P_{\alpha\beta}(-\delta)=P_{\beta\alpha}(\delta)$. This relationship prevails
in spite of unitarity violation. As the survival probability formulas
$P_{ee}$, $P_{\mu\mu}$ and $P_{\tau\tau}$ derived in this work only depend on
$\delta$ through $\cos\delta$, they are unchanged under the transformation
$\delta\rightarrow-\delta$. The relations corresponding to the sums
$\sum_{\alpha}P_{\alpha\beta}$ are therefore omitted.

We calculate the neutrino oscillation probabilities using the probability
formulas (13)–(24) in section III. We use the standard oscillation parameters
$\theta_{12}$, $\theta_{13}$, $\theta_{23}$, $\delta$, $\Delta m_{21}^{2}$ and
$\Delta m_{31}^{2}$ and assume the values that are listed in Table 1. These
values are adopted from Ref. [25] assuming normal neutrino mass ordering. The
allowed parameter values are also shown at 1$\sigma$, 2$\sigma$ and 3$\sigma$
C.L. The resulting unitarity relations are presented in Figure 1.

As we can already tell from equations (41)–(43), the first equation has
corrections of second order and higher. This indicates that unitarity is only
minutely violated and the probability sum $P_{ee}+P_{e\mu}+P_{e\tau}$ is close
to unity. This behavior can also be confirmed in Figure 1. We note that the
deviation from unitarity is the largest for the lowest energies, where the sum
increases monotonically. In equations (42) and (43), we observe a more
significant departure from unitarity at low energies. The comparatively higher
impact on $P_{\mu e}+P_{\mu\mu}+P_{\mu\tau}$ and $P_{\tau
e}+P_{\tau\mu}+P_{\tau\tau}$ is also evident in equations (42) and (43). In
this regard, the zeroth-order corrections proportional to
$1-\text{e}^{-4\gamma\Delta}$ become notable at low energies where $\Delta$ is
large. Moreover, it is observed that violation of unitarity is most profound
in equation (42), where muon neutrino survival and conversion probabilities
are summed over. The qualitative behavior is similar for the probabilities in
equation (43). Examining the full expressions in section III.3, we note that
these probabilities are indeed nearly identical. Since we have assumed
$\theta_{23}$ to be 49∘, we have $\cos\theta_{23}<\sin\theta_{23}$ and
equation (43) becomes larger than equation (42). We also see this behavior in
Figure 1.

Parameter | Unit | Best fit $\bm{\pm 1\sigma}$ |  $\bm{2\sigma}$ range |  $\bm{3\sigma}$ range
---|---|---|---|---
$\Delta m_{21}^{2}$ | $10^{-5}\text{eV}^{2}$ | $7.50_{-0.20}^{+0.22}$ | $7.12-7.93$ | $6.94-8.14$
$|\Delta m_{31}^{2}|$ | $10^{-3}\text{eV}^{2}$ | $2.55_{-0.03}^{+0.02}$ | $2.49-2.60$ | $2.47-2.63$
$\theta_{12}$ | ∘ | $34.3\pm 1.0$ | $32.3-36.4$ | $31.4-37.4$
$\theta_{13}$ | ∘ | $8.53_{-0.12}^{+0.13}$ | $8.27-8.79$ | $8.13-8.92$
$\theta_{23}$ | ∘ | $49.26\pm 0.79$ | $47.37-50.71$ | $41.20-51.33$
$\delta$ | $\pi$ | $1.08_{-0.12}^{+0.13}$ | $0.84-1.42$ | $0.71-1.99$
Table 1: The values of the neutrino oscillation parameters used in this work.
This table is adopted from Ref. [25]. Figure 1: Unitarity relations
$\sum_{\beta}P_{\alpha\beta}$ plotted in the energy range $0.5$ GeV–$10$ GeV
and baseline length $1300$ km, where $\alpha,\beta=e,\mu,\tau$. The neutrino
oscillation probabilities are computed for decay parameter $\gamma=0.1$ and
matter density $\rho=3$ g/cm3.

###  IV.3 Accuracy of the series expansions

We now investigate the accuracy of the series expansions that we presented in
section III. The accuracy of the series expansion method is discussed in great
detail in Ref. [21]. In the present work, we use the double expansion method
where both $\alpha$ and $s_{13}$ are used as expansion parameters. It is
argued in Ref. [21] that this method is more accurate than the ones where
either $\alpha$ or $s_{13}$ is considered as the expansion parameter. The
expansions involving both $\alpha$ and $s_{13}$ are valid when

$\alpha\Delta=\frac{\Delta m_{21}^{2}L}{4E}\ll 1\quad\iff\quad\frac{L}{E}\ll
10^{4}\text{ km/GeV}.$ (44)

Conversely, the probability formulas derived in this work become inaccurate
when either the baseline length $L$ is very long or the neutrino energy $E$ is
very low.

Figure 2 illustrates the accuracy of the probabilities that are obtained from
the analytical probabilities derived in this work. The figure shows the
approximate probabilities for decay parameter $\gamma=0.1$, baseline length
$L=1300$ km and neutrino energies $E=0.5$ GeV–$10$ GeV. The configuration
corresponds to that of the DUNE setup. The probabilities are calculated from
equations (13)–(24) and they are depicted with blue solid lines. The figure
also shows probabilities that have been obtained numerically from the time-
evolution operator as described in section II. The numerically computed
probabilities are shown with yellow dashed lines.

The accuracy condition shown in equation (44) is in agreement with the
probabilities shown in Figure 2. The numerically computed probabilities are
notably different from the approximate ones only at low energies where
$E\lesssim 1.5$ GeV. We note that the discrepancy appears more profound in the
electron neutrino probabilities $P_{ee}$, $P_{e\mu}$ and $P_{e\tau}$, where
the probabilities are constrained to a narrow interval. In general, we find
the analytically computed probabilities to be in a good agreement with the
numerically calculated ones in the DUNE setup.

The absolute errors $\left|\Delta
P_{\alpha\beta}\right|\equiv\big{|}P_{\alpha\beta}^{\text{expansion}}-P_{\alpha\beta}^{\text{numerical}}\big{|}$
for $\alpha,\beta=e,\mu,\tau$ are further illustrated with contour plots in
Figure 3. The figure shows the differences between the approximate and
numerically calculated probabilities in logarithmic scale for neutrino
energies $E=0.1$ GeV–$50$ GeV and baseline lengths $L=100$ km–$10~{}000$ km.
As is shown in the legend bar, regions with darker shades correspond to
smaller absolute errors and those with lighter shades to larger errors,
respectively. We observe that the absolute errors are the smallest in the
lower-right triangular half-plane, where the errors are of order $10^{-3}$.
The series expansion method appears to be the most accurate in regions where
neutrino energies are high and baseline lengths are short. The absolute errors
then grow as neutrino energies become lower and baseline lengths get longer.
The absolute errors also display oscillatory features, as there are regions
where the approximate probabilities are relatively close to the numerical ones
despite the accuracy condition (44) no longer being valid.

(a) Survival probability $P_{ee}$.

(b) Transition probability $P_{e\mu}$.

(c) Transition probability $P_{e\tau}$.

(d) Survival probability $P_{\mu\mu}$.

(e) Transition probability $P_{\mu\tau}$.

(f) Survival probability $P_{\tau\tau}$.

Figure 2: Neutrino oscillation probabilities with invisible decay neutrino
decay of $\gamma=$ 0.1. The probabilities have been plotted for the series
expansions (blue solid line) for neutrino energies $E=0.5$ GeV–$10$ GeV and
baseline length $L=1300$ km. For comparison, numerically computed
probabilities are also plotted (yellow dashed line).

(a) Absolute error $|\Delta P_{ee}|$.

(b) Absolute error $|\Delta P_{e\mu}|$.

(c) Absolute error $|\Delta P_{e\tau}|$.

(d) Absolute error $|\Delta P_{\mu\mu}|$.

(e) Absolute error $|\Delta P_{\mu\tau}|$.

(f) Absolute error $|\Delta P_{\tau\tau}|$.

Figure 3:  Absolute differences for numerically computed probabilities and
probabilities obtained from the series expansions. The differences are shown
for probabilities (a) $P_{ee}$, (b) $P_{e\mu}$, (c) $P_{e\tau}$, (d)
$P_{\mu\mu}$, (e) $P_{\mu\tau}$ and (f) $P_{\tau\tau}$ assuming neutrino
energies $E=0.1$ GeV–$50$ GeV and baseline lengths $L=100$ km–$10~{}000$ km.

Having analytical expressions for the neutrino oscillation probabilities is
useful even at an approximate level, as they demonstrate the functional
dependence on the neutrino oscillation parameters. It is also important to
know the limitations of the approximate probabilities for different baseline
lengths and neutrino energies. In the rest of this section, we examine the
relevance of the series expansions in greater detail.

In Figure 4, we show the transition probability $P_{e\mu}$ and the absolute
error $\Delta P_{e\mu}$ for neutrino energies $E=0.15$ GeV–$3$ GeV and the
baseline length $L=295$ km. This configuration corresponds to the experiment
setup of T2K [27], which is a long-baseline neutrino oscillation experiment
currently taking data in Japan. The series expansion for the conversion
probability $P_{\mu e}$ can be obtained from equation (14) through the
transformation $\delta\rightarrow-\delta$. As before, we take the decay
parameter to be $\gamma=0.1$ and the average matter density $\rho=3$ g/cm3.
The absolute error is computed as the difference of the numerically computed
probability $P_{e\mu}$ and the matching series expansion formula. We observe
in Figure 4 that the analytical and numerically calculated probabilities are
in good agreement with neutrino energies above $0.4$ GeV. Figure 4(a) shows
that the locations of maxima and minima are predicted correctly by the series
expansion formulas, with the first oscillation maximum taking place at about
$0.65$ GeV and the first minimum above $3$ GeV. The absolute errors in Figure
4(b) indicate that the absolute error is the most significant at $0.15$ GeV,
where the error stays below $0.02$. We remark that the absolute error has its
its minimum at $0.6$ GeV, where the error is of order $10^{-5}$. At higher
energies, the absolute error grows to about $10^{-3}$ and starts falling
again. The approximate probabilities are therefore most suitable to study
probabilities at $0.6$ GeV, which also coincides with the peak energy of the
T2K setup.

We next investigate the accuracy of the series expansion that has been derived
for the muon neutrino survival probability $P_{\mu\mu}$. The performance of
the series expansion is shown in Figure 5, where the probability is plotted
with both the approximate method and numerical calculation. The approximate
and numerically acquired probabilities are shown with the blue solid and
yellow dashed lines, respectively. The approximate probabilities are computed
using the series expansion presented in equations (16)–(18). The probabilities
are plotted for neutrino energies $E=0.5$ GeV–$10$ GeV and the baseline length
$L=1300$ km. The setup therefore matches with the DUNE setup [12], which is a
neutrino oscillation experiment currently under construction in USA. Figure
5(a) shows that the series expansion approximates the numerically computed
probabilities well, while Figure 5(b) confirms that the absolute error
$|\Delta P_{\mu\mu}|$ stays below $10^{-2}$ for all considered neutrino
energies and below $3\times 10^{-3}$ for neutrino energies higher than $0.75$
GeV. We note that in the absence of neutrino decay, the probability
corresponding to $\nu_{\mu}\rightarrow\nu_{\mu}$ channel varies between zero
and unity in the energies considered in Figure 5. The effect of neutrino decay
is therefore significant for $P_{\mu\mu}$ for the chosen value $\gamma=0.1$.
In comparison, the effect appears to be less noticeable in the conversion
probability $P_{e\mu}$. Many high-energy accelerator neutrino experiments like
DUNE use muon neutrino beams and are therefore suitable for the search of
invisible neutrino decay.

The accuracy of the series expansions of $P_{\mu e}$ and $P_{\mu\mu}$ have
previously been studied in Ref. [19], where the probabilities have been
computed for the baseline length of the DUNE setup, $L=1300$ km, and the magic
baseline, $L=7000$ km. We note that our results in Figure 4 are mostly in
agreement with the probabilities that were provided in Ref. [19]. The main
difference comes from the absolute error, which is found to be an order of
magnitude smaller for the analytical formulas that have been obtained in the
present work. In the energy range $0.75$ GeV–$3$ GeV, in particular, the
magnitude of the absolute error in Figure 5(b) is of order $10^{-3}$, whereas
the magnitude of the corresponding error in Ref. [19] is of order $10^{-2}$.
This difference in the accuracy arises from the order of the series expansion
used in both works, which is $\mathcal{O}(\lambda^{3})$ in our work and
$\mathcal{O}(\lambda^{2})$ in Ref. [19].

(a) Transition probability $P_{e\mu}$.

(b) Approximation error $|\Delta P_{e\mu}|$.

Figure 4: The transition probability $P_{e\mu}$ and the absolute error
$|\Delta P_{e\mu}|$, plotted in the energy range $E=0.15$ GeV–$3$ GeV with
baseline length $L=295$ km, corresponding to the conditions at the T2K setup.

We finally consider the tau neutrino appearance probability $P_{\mu\tau}$.
Figure 6 shows the probability and its absolute error in the conditions of the
OPERA setup [28], which entails neutrino energies $E=0.3$ GeV–$8$ GeV and the
baseline length $L=730$ km. This energy range includes the first three
oscillation maxima and two minima of $P_{\mu\tau}$. As before, we have plotted
the probability assuming the decay parameter $\gamma=0.1$ and matter density
$\rho=3$ g/cm3. The probabilities and their errors look qualitatively similar
to those plotted for $P_{\mu\mu}$ in Figure 5. The maximum value of the
absolute error $\Delta P_{\mu\tau}=0.02$ is found at the lowest considered
energy $E=0.3$ GeV, while neutrino energies higher than $0.5$ GeV yield errors
that are of order $10^{-3}$. There are three local minima in the absolute
error curve, with the lowest values being of order 10-5.

(a) Transition probability $P_{\mu\mu}$.

(b) Approximation error $|\Delta P_{\mu\mu}|$.

Figure 5: The survival probability $P_{\mu\mu}$ and the absolute error
$|\Delta P_{\mu\mu}|$, plotted in the energy range $E=0.5$ GeV–$10$ GeV with
baseline length $L=1300$ km, corresponding to the conditions of the DUNE
setup.

(a) Transition probability $P_{\mu\tau}$.

(b) Approximation error $|\Delta P_{\mu\tau}|$.

Figure 6: The transition probability $P_{\mu\tau}$ and the absolute error
$|\Delta P_{\mu\tau}|$, plotted in the energy range $E=0.3$ GeV–$8$ GeV with
baseline length $730$ km, corresponding to the conditions of the OPERA setup.

In addition to the scenarios that are studied in this section, the probability
formulas that have been derived in this work are applicable to a wide range of
neutrino experiments. Long-baseline experiments including DUNE and T2HK
benefit not only from having compact analytical probability formulas for the
muon neutrino survival and conversion probabilities $P_{\mu\mu}$ and $P_{\mu
e}$, but understanding their data also requires knowledge on the electron
neutrino probability $P_{ee}$. Experiments with high-energy neutrino beams
like DUNE are furthermore expected to have limited sensitivity to tau neutrino
appearance and therefore benefit from the series expansions of $P_{\mu\tau}$
and $P_{e\tau}$. The electron neutrino survival probability $P_{ee}$ is also
relevant for reactor neutrino experiments like Daya Bay [29] and JUNO [30],
whereas tau neutrino appearance probabilities $P_{e\tau}$ and $P_{\mu\tau}$
may additionally be studied by observing atmospheric neutrinos in experiments
like IceCube [31]. The tau neutrino survival probability $P_{\tau\tau}$ is not
directly testable in neutrino experiments, but it can be indirectly measured
through unitarity relations [32], which have to be modified to account for
neutrino decay. We note that our results are not applicable to conditions
where the accuracy relation (44) is not valid or the matter density profile
cannot be approximated with a constant density. The probability formulas
presented in this work are therefore cannot be used to describe neutrino
oscillations in solar neutrino experiments like SNO [33] and BOREXINO [34].

##  V Summary and conclusions

In this paper, we have presented complete set of approximate formulas for
neutrino oscillation probabilities in presence of invisible neutrino decay. We
have derived compact analytical formulas for the oscillation channels
$\nu_{e}\rightarrow\nu_{e}$, $\nu_{e}\rightarrow\nu_{\mu}$,
$\nu_{e}\rightarrow\nu_{\tau}$, $\nu_{\mu}\rightarrow\nu_{\mu}$,
$\nu_{\mu}\rightarrow\nu_{\tau}$ and $\nu_{\tau}\rightarrow\nu_{\tau}$ by
taking series expansions in small parameters $\alpha$ and $s_{13}$. Our
results include matter effects with constant density. The probability formulas
have been obtained using the Cayley–Hamilton formalism and perturbative
diagonalization. Expressing results in terms of $\lambda=$ 0.2, the analytical
probabilities have been derived up to $\mathcal{O}(\lambda^{2})$ for the
probability $P_{ee}$ and $\mathcal{O}(\lambda^{3})$ for the probabilities
$P_{e\mu}$, $P_{e\tau}$, $P_{\mu\mu}$, $P_{\mu\tau}$ and $P_{\tau\tau}$,
respectively. We have additionally studied the accuracy of the analytical
probability formulas for a wide range of neutrino energies $E$ and baseline
lengths $L$ and examined the impact on unitarity relations. Below we summarize
the main findings of this work:

  * •

The probability formulas derived are valid for the accuracy condition
$\alpha\Delta=\Delta m_{21}^{2}L/(4E)\ll 1$, which holds up to $L/E\ll 10^{4}$
km/GeV. The absolute errors in this region are generally of order $10^{-3}$,
indicating a good agreement between numerically computed probabilities and the
analytical probability formulas.

  * •

The accuracy of the derived probability formulas is particularly good for the
conditions of T2K and T2HK, where the absolute errors for $P_{e\mu}$ are of
order $10^{-5}$ for neutrino energies $E\simeq 0.6$ GeV. We also found a good
agreement for the DUNE setup, where the absolute error pertaining to
$P_{\mu\mu}$ is of order $10^{-2}$ at the lowest neutrino energies and stays
under $3\times 10^{-3}$ for energies higher than $0.75$ GeV.

  * •

Neutrino oscillation probabilities with neutrino decay have previously been
studied with the series expansion method in Ref. [19]. We find that our
results for the probabilities $P_{e\mu}$ and $P_{\mu\mu}$ are in agreement
with the equivalent relations in Ref. [19]. Our work reproduces the formula
for $P_{e\mu}$ completely and extends the formula for $P_{\mu\mu}$ to include
terms up to $\mathcal{O}(\lambda^{3})$, respectively.

  * •

Probabilities without neutrino decay have been studied using the series
expansion method in Ref. [21]. Comparing the probability formulas derived at
the no-decay limit $\gamma\rightarrow 0$ in this work and the corresponding
formulas in Ref. [21], we find our formulas for $P_{ee}$, $P_{e\mu}$ and
$P_{e\tau}$ to be consistent with those in Ref. [21]. We also find a third-
order term in the formulas for $P_{\mu\mu}$, $P_{\mu\tau}$ and $P_{\tau\tau}$
that is present in Ref. [21] but not in this work. We estimate the difference
to be of order $10^{-3}$.

  * •

Probability formulas for calculating neutrino oscillation probabilities in
vacuum have also been studied in the OMSD limit [26]. We find the probability
formulas derived in this work to be in agreement with those in Ref. [26] after
applying the vacuum limit $A\rightarrow 0$.

  * •

Unitarity is generally violated in presence of neutrino decay. In the muon and
tau neutrino probabilities, we find the violation to be proportional to
$1-\text{e}^{-4\gamma\Delta}$ at the leading order, whereas in the electron
neutrino probabilities the violation is of the order
$\mathcal{O}(\lambda^{2})$.

In conclusion, we have presented series expansions to approximate the
probabilities for neutrino oscillations in constant matter density and one
decaying neutrino mass eigenstate. The formulas derived in this work are in
agreement with those provided in Ref. [19] for $P_{e\mu}$ and $P_{\mu\mu}$,
whereas the rest of our results are novel. We have provided a detailed study
on the numerical accuracy of the formulas and discussed their relevance for
various neutrino experiments.

##  Appendix A Matrix exponential in the Cayley–Hamilton formalism

In this section, we briefly describe the Cayley–Hamilton formalism. In order
to obtain the exponent form of an $N\times N$ matrix $\mathcal{M}$, one can
write the series

$\text{e}^{\mathcal{M}}=\sum_{n=0}^{\infty}\frac{1}{n!}\mathcal{M}^{n}.$ (45)

Following from linear algebra, the characteristic equation of $\mathcal{M}$ is

$\chi_{\mathcal{M}}(\lambda)\equiv\det(\mathcal{M}-\lambda\bm{I})=0,$ (46)

where $\chi_{\mathcal{M}}$ is the characteristic polynomial, $\bm{I}$ is the
$N\times N$ identity matrix and $\lambda$ is an eigenvalue of $\mathcal{M}$.
The characteristic polynomial of $\mathcal{M}$ can be written in terms of its
eigenvalues as

$\chi_{\mathcal{M}}(\lambda)=\lambda^{N}+a_{N-1}\lambda^{N-1}+\ldots+a_{1}\lambda+a_{0}=0,$
(47)

where $a_{0},a_{1},\ldots,a_{N-1}$ are some coefficients.

The Cayley–Hamilton theorem states that every matrix $\mathcal{M}$ satisfies
the characteristic equation $\chi_{\mathcal{M}}=0$, that is,

$\mathcal{M}^{N}=-a_{N-1}\mathcal{M}^{N-1}-\ldots-
a_{1}\mathcal{M}-a_{0}\bm{I},$ (48)

where the highest power of $\mathcal{M}$ is written in terms of its lower
powers. Therefore, for any $p\geq N$, there exist coefficients
$b_{0},b_{1},\ldots,b_{N-1}$ that satisfy the equation

$M^{p}=b_{N-1}\mathcal{M}^{N-1}+\ldots+b_{1}\mathcal{M}+b_{0}\bm{I}.$ (49)

The Cayley–Hamilton formalism makes it possible to reduce the infinite series
of exponential matrix $\text{e}^{\mathcal{M}}$ to a finite one, using a set of
coefficients $c_{0},c_{1},\ldots,c_{N-1}$ as

$\text{e}^{\mathcal{M}}=c_{0}\bm{I}+c_{1}\mathcal{M}+\ldots+c_{N-1}\mathcal{M}^{N-1}=\sum_{n=0}^{N-1}c_{n}\mathcal{M}^{n}.$
(50)

In neutrino oscillations, the formalism is useful in expressing the time-
evolution operator $S=\text{e}^{iHL}$ with a reduced set of operators. For
three active flavor states, the reduction can be performed in three terms.

##  Appendix B Perturbative diagonalization

In this work, we adopt the approach in Ref. [21] for the perturbative
diagonalization of the neutrino mixing matrix. Let us consider a $3\times 3$
matrix $M$ that can be diagonalized by another matrix $W$:

$M=W\hat{M}W^{\dagger},$ (51)

where $\hat{M}={\rm diag}(\lambda_{1},\lambda_{2},\lambda_{3})$ is the
diagonalized matrix and $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are
eigenvalues of $M$. Correspondingly, one may obtain the eigenvectors $v_{1}$,
$v_{2}$ and $v_{3}$ from the columns of the diagonalizing matrix $W$. In the
perturbative approach, the quantities can be decomposed into contributions
consisting of terms in various orders of perturbation parameters, such as
$M=M^{(0)}+M^{(1)}+M^{(2)}+\ldots$ . This method can be applied to both
eigenvalues and eigenvectors. In three-flavor oscillations, a convenient
choice for $M$ is

$M\equiv U_{13}U_{12}{\rm
diag}(0,\alpha,1-i\gamma)U_{12}^{\dagger}U_{13}^{\dagger}+{\rm diag}(A,0,0),$
(52)

where $U_{ij}$ are elements of the Pontecorvo–Maki–Nakagawa–Sakata matrix and
$A$ represents the matter potential. The Hamiltonian then takes the well-known
form

$\mathcal{H}=\frac{\Delta
m_{31}^{2}}{2E}U_{23}U_{\delta}MU_{\delta}^{\dagger}U_{23}^{\dagger}.$ (53)

This turns out to be a convenient choice, as the matter potential matrix is
unaffected by the action of $U_{23}$ and $U_{\delta}$. The matrix product
$U_{23}U_{\delta}W\equiv\tilde{U}$ therefore represents the leptonic mixing
matrix in matter. Writing equation (52) in full form and and grouping its
terms in orders of the expansion parameters $\alpha$ and $s_{13}$, we obtain

$\displaystyle M^{(0)}$ $\displaystyle=\begin{pmatrix}A&0&0\\\ 0&0&0\\\
0&0&1-i\gamma\end{pmatrix},$ (54) $\displaystyle M^{(1)}$
$\displaystyle=\begin{pmatrix}0&0&c_{13}s_{13}(1-i\gamma)\\\ 0&0&0\\\
c_{13}s_{13}(1-i\gamma)&0&0\end{pmatrix},$ (55) $\displaystyle M^{(2)}$
$\displaystyle=\begin{pmatrix}\alpha s_{12}^{2}+s_{13}^{2}(1-i\gamma)&\alpha
c_{13}c_{12}s_{12}&0\\\ \alpha c_{13}c_{12}s_{12}&\alpha c_{12}^{2}&0\\\
0&0&-s_{13}^{2}(1-i\gamma)\end{pmatrix}.$ (56)

Since $M^{(0)}$ is diagonal, the zeroth-order contributions to eigenvalues are
diagonal elements, whereas the corresponding corrections to the eigenvectors
are the standard Euclidean basis vectors $e_{1}$, $e_{2}$ and $e_{3}$. The
first and second order corrections to the eigenvalues on the other hand are
given as

$\displaystyle\lambda_{i}^{(1)}$ $\displaystyle=M_{ii}^{(1)},$ (57)
$\displaystyle\lambda_{i}^{(2)}$ $\displaystyle=M_{ii}^{(2)}+\sum_{j\neq
i}\frac{\left(M_{ij}^{(1)}\right)^{2}}{\lambda_{i}^{(0)}-\lambda_{j}^{(0)}},$
(58)

and to the eigenvectors as

$\displaystyle v_{i}^{(1)}$ $\displaystyle=\sum_{j\neq
i}\frac{M_{ij}^{(1)}}{\lambda_{i}^{(0)}-\lambda_{j}^{(0)}}e_{j},$ (59)
$\displaystyle v_{i}^{(2)}$ $\displaystyle=\sum_{j\neq
i}\frac{1}{\lambda_{i}^{(0)}-\lambda_{j}^{(0)}}\left[M_{ij}^{(2)}+\left(M^{(1)}v_{i}^{(1)}\right)_{j}-\lambda_{i}\left(v_{i}^{(1)}\right)_{j}\right]e_{j},$
(60)

respectively. The energy eigenvalues are therefore written for the three-
flavor case as

$E_{i}=\frac{\Delta m_{31}^{2}}{2E}\lambda_{i},$ (61)

where $i=$ 1, 2 and 3.

## References

  * Barger  _et al._ [1999] V. D. Barger, J. G. Learned, S. Pakvasa, and T. J. Weiler, Phys. Rev. Lett. 82, 2640 (1999), arXiv:astro-ph/9810121 .
  * Acker  _et al._ [1992] A. Acker, S. Pakvasa, and J. T. Pantaleone, Phys. Rev. D 45, 1 (1992).
  * Acker and Pakvasa [1994] A. Acker and S. Pakvasa, Phys. Lett. B 320, 320 (1994), arXiv:hep-ph/9310207 .
  * Gelmini and Roncadelli [1981] G. B. Gelmini and M. Roncadelli, Phys. Lett. B 99, 411 (1981).
  * Chikashige  _et al._ [1981] Y. Chikashige, R. N. Mohapatra, and R. D. Peccei, Phys. Lett. B 98, 265 (1981).
  * Berryman  _et al._ [2015] J. M. Berryman, A. de Gouvea, and D. Hernandez, Phys. Rev. D 92, 073003 (2015), arXiv:1411.0308 [hep-ph] .
  * Gonzalez-Garcia and Maltoni [2008] M. C. Gonzalez-Garcia and M. Maltoni, Phys. Lett. B 663, 405 (2008), arXiv:0802.3699 [hep-ph] .
  * Frieman  _et al._ [1988] J. A. Frieman, H. E. Haber, and K. Freese, Phys. Lett. B 200, 115 (1988).
  * Iváñez Ballesteros and Volpe [2023] P. Iváñez Ballesteros and M. C. Volpe, Phys. Lett. B 847, 138252 (2023), arXiv:2307.03549 [hep-ph] .
  * Aartsen  _et al._ [2021] M. G. Aartsen _et al._ (IceCube-Gen2), J. Phys. G 48, 060501 (2021), arXiv:2008.04323 [astro-ph.HE] .
  * Adrian-Martinez  _et al._ [2016] S. Adrian-Martinez _et al._ (KM3Net), J. Phys. G 43, 084001 (2016), arXiv:1601.07459 [astro-ph.IM] .
  * Acciarri  _et al._ [2015] R. Acciarri _et al._ (DUNE), arXiv:1512.06148 [physics.ins-det] (2015).
  * Alekou  _et al._ [2023] A. Alekou _et al._ (ESSnuSB), Universe 9, 347 (2023), arXiv:2303.17356 [hep-ex] .
  * Abe  _et al._ [2015] K. Abe _et al._ (Hyper-Kamiokande), PTEP 2015, 053C02 (2015), arXiv:1502.05199 [hep-ex] .
  * Lindner  _et al._ [2001] M. Lindner, T. Ohlsson, and W. Winter, Nucl. Phys. B 607, 326 (2001), arXiv:hep-ph/0103170 .
  * Abrahão  _et al._ [2015] T. Abrahão, H. Minakata, H. Nunokawa, and A. A. Quiroga, JHEP 11, 001, arXiv:1506.02314 [hep-ph] .
  * Ghoshal  _et al._ [2021] A. Ghoshal, A. Giarnetti, and D. Meloni, J. Phys. G 48, 055004 (2021), arXiv:2003.09012 [hep-ph] .
  * Chattopadhyay  _et al._ [2022] D. S. Chattopadhyay, K. Chakraborty, A. Dighe, S. Goswami, and S. M. Lakshmi, Phys. Rev. Lett. 129, 011802 (2022), arXiv:2111.13128 [hep-ph] .
  * Chattopadhyay  _et al._ [2023] D. S. Chattopadhyay, K. Chakraborty, A. Dighe, and S. Goswami, JHEP 01, 051, arXiv:2204.05803 [hep-ph] .
  * Kimura [2017] T. Kimura, PTEP 2017, 041A03 (2017), arXiv:1702.04681 [math-ph] .
  * Akhmedov  _et al._ [2004] E. K. Akhmedov, R. Johansson, M. Lindner, T. Ohlsson, and T. Schwetz, JHEP 04, 078, arXiv:hep-ph/0402175 .
  * Ohlsson and Snellman [2000a] T. Ohlsson and H. Snellman, J. Math. Phys. 41, 2768 (2000a), [Erratum: J. Math. Phys. 42, 2345 (2001)], arXiv:hep-ph/9910546 .
  * Ohlsson and Snellman [2000b] T. Ohlsson and H. Snellman, Phys. Lett. B 474, 153 (2000b), [Erratum: Phys. Lett. B 480, 419–419 (2000)], arXiv:hep-ph/9912295 .
  * Ohlsson [2001] T. Ohlsson, Phys. Scripta T 93, 18 (2001).
  * de Salas  _et al._ [2021] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola, and J. W. F. Valle, JHEP 02, 071, arXiv:2006.11237 [hep-ph] .
  * Choubey  _et al._ [2021] S. Choubey, M. Ghosh, D. Kempe, and T. Ohlsson, JHEP 05, 133, arXiv:2010.16334 [hep-ph] .
  * Abe  _et al._ [2014] K. Abe _et al._ (T2K), Phys. Rev. Lett. 112, 061802 (2014), arXiv:1311.4750 [hep-ex] .
  * Acquafredda  _et al._ [2009] R. Acquafredda _et al._ , JINST 4, P04018.
  * An  _et al._ [2017] F. P. An _et al._ (Daya Bay), Phys. Rev. D 95, 072006 (2017), arXiv:1610.04802 [hep-ex] .
  * An  _et al._ [2016] F. An _et al._ (JUNO), J. Phys. G 43, 030401 (2016), arXiv:1507.05613 [physics.ins-det] .
  * Aartsen  _et al._ [2019] M. G. Aartsen _et al._ (IceCube), Phys. Rev. D 99, 032007 (2019), arXiv:1901.05366 [hep-ex] .
  * Martinez-Soler and Minakata [2021] I. Martinez-Soler and H. Minakata, Phys. Rev. D 104, 093006 (2021), arXiv:2109.06933 [hep-ph] .
  * Berezinsky and Lissia [2001] V. Berezinsky and M. Lissia, Phys. Lett. B 521, 287 (2001), arXiv:hep-ph/0108108 .
  * Pallavicini  _et al._ [2023] M. Pallavicini _et al._ (BOREXINO), in  _16th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics and Relativistic Field Theories_ (2023).