#  Neutrino Oscillations in Matter using the
Adjugate of the Hamiltonian

Asli Abdullahi  Stephen J. Parke<EMAIL_ADDRESS><EMAIL_ADDRESS>Particle Theory
Dept., Fermi National Accelerator Laboratory, Batavia, IL, USA

(August 11, 2023)

###### Abstract

We revisit neutrino oscillations in constant matter density for a number of
different scenarios: three flavors with the standard Wolfenstein matter
potential, four flavors with standard matter potential and three flavors with
non-standard matter potentials. To calculate the oscillation probabilities for
these scenarios one must determine the eigenvalues and eigenvectors of the
Hamiltonians. We use a method for calculating the eigenvalues that is well
known, determination of the zeros of determinant of matrix $(\lambda I-H)$,
where H is the Hamiltonian, I the identity matrix and $\lambda$ is a scalar.
To calculate the associated eigenvectors we use a method that is little known
in the particle physics community, the calculation of the adjugate (transpose
of the cofactor matrix) of the same matrix, $(\lambda I-H)$. This method can
be applied to any Hamiltonian, but provides a very simple way to determine the
eigenvectors for neutrino oscillation in matter, independent of the complexity
of the matter potential. This method can be trivially automated using the
Faddeev-LeVerrier algorithm for numerical calculations. For the above
scenarios we derive a number of quantities that are invariant of the matter
potential, many are new such as the generalization of the Naumov-Harrison-
Scott identity for four or more flavors of neutrinos. We also show how these
matter potential independent quantities become matter potential dependent when
off-diagonal non-standard matter effects are included.

††preprint: FERMILAB-Pub-22-922-T

orcid AA: 0000-0002-6122-4986 orcid SJP: 0000-0003-2028-6782

##  I Introduction

The Wolfenstein matter effect Wolfenstein (1978) has been and continues to be
very important for exploring neutrino mixing and masses. It is the solution to
the solar neutrino anomalyMikheyev and Smirnov (1985); Parke (1986),
discovered by Davis et. al. Cleveland et al. (1995) and in particular the SNO
experiment Ahmad et al. (2002) determined $|U_{e2}|^{2}$
($\sin^{2}\theta_{12}\approx 0.3$) and that the sign111Our labelling of the
mass eigenstates is determined by $|U_{e1}|^{2}>|U_{e2}|^{2}>|U_{e3}|^{2}$ and
the form of the Pontecorvo, Maki, Nakagawa, Sakata (PMNS) matrix, U, as in the
PDG Workman et al. (2022).  of $\Delta m^{2}_{21}$ was positive , i.e. that
the state with most $\nu_{e}$, $\nu_{1}$, is lighter than state with the next
most $\nu_{e}$, $\nu_{2}$. The current long baseline (LBL) experiments, T2K
Abe et al. (2011) and NOvA Ayres et al. (2004), as well as the future
experiments, Hyper-Kamiokande (HK) Abe et al. (2018), Korean Neutrino
observatory (KNO) Abe et al. (2018) and DUNE Acciarri et al. (2015), will
determine the parameters $|U_{\mu 3}|^{2}$ ($\approx\sin^{2}\theta_{23}$) and
$\Delta m^{2}_{32}$ including its sign. If $|U_{\mu 3}|^{2}>1/2$ ($|U_{\mu
3}|^{2}<1/2$) then $\nu_{\mu}$ ($\nu_{\tau}$) dominates the neutrino mass
eigenstate with the least $\nu_{e}$, $\nu_{3}$. The sign of $\Delta
m^{2}_{32}$ determines the atmospheric mass ordering, whether $\nu_{3}$ is the
heaviest or lightest mass state in the neutrino spectrum. Most importantly,
these experiments will determine whether there is significant CP violation in
the neutrino sector via the determination of the CP violating Jarlskog
invariant Jarlskog (1985) (or the phase $\delta_{CP}$).

None of these LBL experiments are performed in vacuum as the neutrino beams
passes through 295 to 1300 km of the earth’s crust. Therefore the Wolfenstein
matter effect must be included in the calculation of the oscillation
probability and can have significant effects on the oscillation probability
and the determination of the neutrino parameters for all of the LBL
experiments. This is especially true for DUNE that will use the matter effects
to determine the atmospheric mass ordering but it is also true for all the
experiments for the determination of the size of CP violation, Mena and Parke
(2004). Although the density of the earth’s crust varies somewhat over the
baseline of all the above LBL experiments it has been shown in Kelly and Parke
(2018); King et al. (2020) that a constant density approximation gives
accurate enough oscillation probabilities given the statistics of these
experiments. For a review of CP violation in neutrino oscillations see
Nunokawa et al. (2008).

The neutrino oscillation probabilities for three neutrino flavors in constant
medium have been exactly calculated222Use of efficient numerical packages,
while excellent for generating figures and doing analysis, they do not provide
detailed physics insight or understanding. by many authors: Barger et al
Barger et al. (1980), Zaglauer & Schwarzer Zaglauer and Schwarzer (1988),
Harrison & Scott Harrison and Scott (2000), Ohlsson & Snellman Ohlsson and
Snellman (2000, 1999), Kimura et al Kimura et al. (2002a, b) as well as many
approximation schemes reviewed in Barenboim et al. (2019). All of the exact
methods first determine the squared masses of the neutrinos in the medium by
solving the characteristic equation, $\text{\bf Det}(\lambda I-H)=0$, where
$\lambda$’s are the eigenvalues of the Hamiltonian, $H$, and $I$ is the
identity matrix. Where these exact methods differ, is how one then determine
the eigenvectors associated with the eigenvalues. Using the methods implored
in the above papers the calculations of the eigenvectors are far from
transparent. The eigenvectors make up the elements of the mixing matrix for
the neutrino states in the medium and with these one can easily calculate the
oscillation probabilities. In this paper, we present a different, efficient
and simple way to determine the eigenvectors which is known in the mathematics
community but so far has not appeared in the neutrino physics literature. We
call this method the “adjugate method” as it requires the calculation of the
adjugate, the transpose of the co-factor matrix, of the Hamiltonian,
$\text{\bf Adj}(H)$. Because of the form of the neutrino Hamiltonian in the
flavor basis this adjugate method is particularly suited to neutrino
oscillation calculations as we will show in Sec. III. We use this method to
identify a number of quantities that are invariant, that is, independent of
the matter potential, in the various scenarios addressed. For the standard
three flavor case these identities are known whereas for the other scenarios
these identities are new such as the generalization of the Naumov-Harrison-
Scott identity for four or more flavors of neutrinos.

The outline of this paper is as follows: in Sec. II we revisit what is
required to calculate the oscillation probabilities in matter with particular
attention to the intrinsic CP violating term. Sec. III we discuss the various
methods for diagonalizing a Hamiltonian and present the “adjugate method” and
why it is particularly suited to neutrino oscillation calculations. Sec. IV is
revisiting the well known three flavor case to demonstrate how simple and
straight forward this method is. Sec. V is application of this method to four
neutrino flavors, three active and one sterile, with particularly attention to
the matter potential invariants such as the generalization of various
identities, including the Naumov-Harrison-Scott identity. Sec. VI addresses
the case of three flavors plus non-standard interactions (NSI). Next is the
conclusions, Sec. VII, followed by a number of very useful appendices.

##  II Calculation of the Oscillation Probabilities

For neutrinos propagating through a constant medium, the oscillation
probabilities can be easily calculated if you know the eigenvalues,
$\lambda_{i}$, of the Hamiltonian in the flavor bases and $V_{\alpha
i}V^{*}_{\beta i}$ of the normalized eigenvectors, $V_{\alpha i}$, no matter
how complicated the Hamiltonian. The $\lambda_{i}$’s and the $V_{\alpha i}$’s
are the medium equivalents of $m^{2}_{i}$ and elements of the PMNS matrix
$U_{\alpha i}$ in vacuum. In general the Hamiltonian is not the same for
neutrinos and anti-neutrinos, so the $\lambda_{i}$ and $V_{\alpha
i}V^{*}_{\beta i}$ need to be calculated for both. With these qualities the
neutrino oscillation probabilities are given by

$\displaystyle P(\nu_{\alpha}\to\nu_{\beta})$ $\displaystyle=$
$\displaystyle\left|~{}\sum_{j}V^{*}_{\alpha j}~{}V_{\beta
j}\text{e}^{-i\lambda_{j}L/(2E)}~{}\right|^{2}$ (1) $\displaystyle=$
$\displaystyle\delta_{\alpha\beta}-4\sum_{i>j}{\Re}(V_{\alpha i}V^{*}_{\beta
i}\,V^{*}_{\alpha j}V_{\beta j})\sin^{2}\left(\Delta_{ij}\right)$
$\displaystyle\quad\quad-8\sum_{i>j}{\Im}(V_{\alpha i}V^{*}_{\beta
i}\,V^{*}_{\alpha j}V_{\beta
j})\,\sin\Delta_{ij}\sin\Delta_{ik}\sin\Delta_{jk}\,,$

with $\Delta_{ij}\equiv(\lambda_{i}-\lambda_{j})L/(4E)$. There is no sum over
k and it is fixed for the calculation, typically k=1 is chosen. The last term
contains the intrinsic CP violating piece and is only non-zero when
$\alpha\neq\beta$. The usual way of writing this term, as in the PDG,

$2\sum_{i>j}{\Im}(V_{\alpha i}V^{*}_{\beta i}V^{*}_{\alpha j}V_{\beta
j})\sin\left(2\Delta_{ij}\right)$

appears to have a linear dependence in $(L/E)$, but this is illusionary, as
this term must be of order $(L/E)^{3}$ in the small (L/E) limit as can be seen
directly from the first line of this expression there can be no terms linear
in $(L/E)$. Therefore, it is more informative to rewrite the CP violating term
as in eq. 1, since

$\displaystyle 2\sum_{i>j}{\Im}(V_{\alpha i}V^{*}_{\beta i}V^{*}_{\alpha
j}V_{\beta j})\sin\left(2\Delta_{ij}\right)=8\sum_{i>j}{\Im}(V_{\alpha
i}V^{*}_{\beta i}\,V^{*}_{\alpha j}V_{\beta
j})\,\sin\Delta_{ij}\sin\Delta_{ik}\sin\Delta_{jk}$ (2)

for any fixed k. So for n-flavors there are $n$ ways to chose k, they are all
equivalent. Typically k=1 is chosen. There are only (n-1)(n-2)/2 non-zero
terms on the RHS, as the terms when i=k or j=k are zero. Therefore the number
of such terms is the same as the number of Dirac type phases in the n-flavor
PMNS matrix. A derivation of eq. 2 is given in appendix A. Details for the
case when n=4 are presented in Section V.

In the next section we will show how to calculate both $\lambda_{i}$ and
$V_{\alpha i}V^{*}_{\beta i}$ in a systematic, straightforward and simple way,
independent of the number of neutrino flavors and independent of the
complexity of the Hamiltonian. Although the method for calculating
$\lambda_{i}$ is well known, the closely associated method for calculating
$V_{\alpha i}V^{*}_{\beta i}$ has not appeared in the neutrino physics
literature before and to our surprise is not well known in the particle
physics community. Once the $\lambda_{i}$’s are known, there is a simple
matrix expression using the adjugate of a matrix that relates the
$\lambda_{i}$’s and the Hamiltonian in flavor basis to the $V_{\alpha
i}V^{*}_{\beta i}$’s.

##  III Diagonalizing a Hermitian Hamiltonian

Let H be an n $\times$ n Hermitian matrix that we want to diagonalize as
follows:

$\displaystyle H=V\Lambda V^{\dagger}$ (3)

where $\Lambda=\text{Diag}(\lambda_{1},\lambda_{2},\cdots,\lambda_{n})$ and
$V_{\alpha i}$ is a unitary matrix such that i-th column of $V$ is a
normalized eigenvector of H with eigenvalue $\lambda_{i}$, since
$HV=V\Lambda$.

It is well known that the eigenvalues are obtained be solving the
characteristic equation

$\displaystyle\text{Det}[\lambda_{i}I-H]=0\,.$ (4)

For $n\leq 4$ this can be done analytically, but for $n$ $>$ 4 this cannot be
done analytically in general.

It is well known in the mathematics community, but rarely known in the
particle physics community, that the components of the normalized,
i${}^{\text{th}}$ eigenvector are easily calculated using

$\displaystyle V_{\alpha i}V^{*}_{\beta i}=\frac{\text{\bf
Adj}[\lambda_{i}I-H]_{\alpha\beta}}{\Pi_{k\neq
i}(\lambda_{i}-\lambda_{k})}\,,$ (5)

where Adj is the adjugate of the matrix333See wikipedia Adjugate Matrix. The
adjugate of a matrix (the transpose of the co-factor matrix) is well defined
independent of whether or not the matrix is invertible. For invertible
matrices, the inverse of the matrix can be calculated using $H^{-1}=\text{\bf
Adj}[H]/\text{\bf Det}[H]$. The matrix $(\lambda_{i}I-H)$ is not invertible
since $\text{\bf Det}[\lambda_{i}I-H]=0$. , see Denton et al. (2022), i.e. the
transpose of the co-factor matrix and references therein444Note, adding $xI$
to the Hamiltonian leaves the RHS of this equation unchanged, as x is added to
all eigenvalues. The physics here is that, neutrino oscillations cannot
determine the absolute neutrino mass scale.. For neutrino oscillations one
needs exactly the quantities $V_{\alpha i}V^{*}_{\beta i}$ which are given
directly and simply from eq.5. If needed, this equation can also be used to
calculate the eigenvectors in a simple, straight forward manner, up to an
overall phase which is arbitrary. The quantity $V_{\alpha i}V^{*}_{\beta i}$
given by eq. 5 is independent of this arbitrary choice but can be used to
calculate the magnitude, $|V_{\alpha i}|$, and relative phase of each
component of an eigenvector, since this phase is given by $V_{\alpha
i}V^{*}_{\beta i}/(|V_{\alpha i}||V_{\beta i}|)$. A physicist’s style sketch
of why eq.5 is correct is given in appendix B.

With $\lambda_{i}$’s and the $V_{\alpha i}V^{*}_{\beta i}$’s constructed this
way, eq. 3 is satisfied and 555While the first of these is easily checked from
eq. 5, the latter is most easily checked in numerical calculations using
$\sum_{\beta}V_{\alpha i}V^{*}_{\beta i}V_{\beta j}V^{*}_{\gamma
j}=\delta_{ij}V_{\alpha i}V^{*}_{\gamma i}$.

$\displaystyle\sum_{i}V_{\alpha i}V^{*}_{\beta
i}=\delta_{\alpha\beta}\quad\text{and}\quad\sum_{\alpha}V_{\alpha
i}V^{*}_{\alpha j}=\delta_{ij}\,.$ (6)

Therefore the eigenvectors constructed this way are normalized and orthogonal.

Even less well known until recently, even in the mathematics community, is
that the magnitudes of the elements of V can be simple, obtained using the
eigenvector-eigenvalue identity: Denton et al. (2020) & Denton et al. (2022),

$\displaystyle|V_{\alpha
i}|^{2}=\frac{\Pi_{l}(\lambda_{i}-\lambda_{l}(h_{\alpha}))}{\Pi_{k\neq
i}(\lambda_{i}-\lambda_{k})}$ (7)

where $h_{\alpha}$ is the principal minor of H with the $\alpha$-th row and
column deleted. $\lambda_{l}(h_{\alpha})$ are the eigenvalues of $h_{\alpha}$,
i.e. the solutions to $\text{Det}[\lambda_{l}I-h_{\alpha}]=0$. Constructed
this way $\sum_{i}|V_{\alpha i}|^{2}=1=\sum_{\alpha}|V_{\alpha i}|^{2}$, i.e.
they are normalized appropriately. Eq. 7 is mathematically equivalent to eq. 5
when $\alpha=\beta$, however the square of the diagonal elements are more
simply obtained using eq. 7. For standard three flavor oscillations in matter,
eq.7 is sufficient, see Denton et al. (2022), however for more complicated
scenarios like 3+1 or NSIs then eq. 5 is required. Eq. 5 gives the most
straight forward and systematic way to evaluate the magnitude and relative
phase of all components of the eigenvectors and gives directly the
quantities,$V_{\alpha i}V^{*}_{\beta i}$, for neutrino oscillation
calculations.

Therefore, for the evaluation of the eigenvalues we need, $\text{\bf
Det}(\lambda I-H)$, and for the eigenvectors, $\text{\bf Adj}(\lambda I-H)$ is
also needed, see eq. 4 & eq.5. We first note that both are polynomials in
$\lambda$ of degree $n$ and (n-1) respectively, i.e.

$\displaystyle\text{\bf Det}[\lambda I-H]$
$\displaystyle\equiv\lambda^{n}+\lambda^{n-1}d_{1}+\cdots+d_{n}$ (8)
$\displaystyle\text{\bf Adj}[\lambda I-H]$
$\displaystyle\equiv\lambda^{n-1}A_{1}+\lambda^{n-2}A_{2}+\cdots+A_{n}\,.$ (9)

where the $d_{i}$’s are scalars and the $A_{i}$’s are n $\times$ n Hermitian
matrices.

There are at least five of different algebraic methods that can be used to
calculate $\text{\bf Det}[\lambda I-H]$ and $\text{\bf Adj}[\lambda I-H]$, we
review them here:

  1. 1.

Using the Leibniz formulae666 See wikipedia Leibniz for n $\times$ n to
calculate the determinant and elements of the adjugate of $(\lambda I-H)$.
However this standard method does not, in general, give the simplest results
without significant additional manipulations.

  2. 2.

Le Verrier-Faddeev algorithm: a more general method is the very simple and
recursive Le Verrier-Faddeev algorithm777See wikipedia Faddeev–LeVerrier
algorithm. , see also Hou (1998), given by

$\displaystyle A_{1}\equiv I_{n}\quad\quad$ $\displaystyle
d_{1}=-\text{Tr}(HA_{1})$ (10)
$\displaystyle\,\quad=-\text{Tr}(H)=-\sum_{j}\lambda_{j}$ $\displaystyle
A_{2}=HA_{1}+d_{1}I_{n}.\quad$ $\displaystyle\quad
d_{2}=-(1/2)\text{Tr}(HA_{2})$
$\displaystyle~{}\quad=H-\text{Tr}(H)I_{n}\quad$
$\displaystyle~{}\quad=(1/2)(\text{Tr}^{2}(H)-\text{Tr}(H^{2}))=\sum_{j1>j2}\lambda_{j1}\lambda_{j2}$
$\displaystyle\vdots\quad$ $\displaystyle\quad\quad\vdots$ $\displaystyle
A_{n}=HA_{n-1}+d_{n-1}I_{n}\quad$ $\displaystyle\quad
d_{n}=-(1/n)\text{Tr}(HA_{n})$ $\displaystyle\,=(-1)^{n-1}\text{\bf Adj}[H]$
$\displaystyle\,\quad\quad=(-1)^{n}\text{\bf
Det}(H)=(-1)^{n}\Pi^{n}_{j=1}\lambda_{j}\,.$

The next iteration gives

$A_{n+1}=(-1)^{n-1}[H\text{\bf Adj}(H)-\text{\bf Det}(H)I_{n}]=0,$

which is related to the Caley-Hamiltonian identity. Then $d_{n+1}=0$ follows,
terminating the series. As expected, from eq. 8, the d’s are the combination
of traces of $H^{m}$ (m=1,…, n) such that
$d_{m}=(-1)^{m}\sum_{j1>...>jm}\lambda_{j1}\cdots\lambda_{jm}$. In appendix C
the d’s and A’s are explicitly given for n$\leq$6.

It is worth noting here the following relationship between the trace of the
A’s and the d’s: $\text{\bf Tr}[A_{m+1}]=(n-m)d_{m}$, valid for m=0, $\cdots$,
n. (using $d_{0}\equiv 1$), and therefore giving

$\displaystyle\text{\bf Tr}[\text{\bf Adj}[\lambda
I-H]]=\frac{\partial}{\partial\lambda}\,{\bf Det}[\lambda I-H]$ (11)

for any $\lambda$.

  3. 3.

Lagrange method: the results of the Le Verrier-Faddeev algorithm can be re-
arranged, see appendix D, as

$\displaystyle\text{\bf Adj}(\lambda_{i}I-H)$ $\displaystyle={\bf\Pi}_{j\neq
i}(H-\lambda_{j}I)\,.$ (12)

This is the Lagrange method that has been used in many papers: Harrison and
Scott (2002), Fong (2022). The LHS of this equation only depends on the
eigenvalue $\lambda_{i}$ whereas the RHS depends on all the other eigenvalues,
$\lambda_{j}$’s with $j\neq i$. The LHS and RHS of eq. 12 are the two extremes
of the many ways to write this expression given the many relationships between
the traces of the powers of the Hamiltonian and products of the eigenvalues,
e.g. $\sum_{i}\lambda_{i}=\text{tr}(H)$,
$\sum_{i>j}\lambda_{i}\lambda_{j}=\frac{1}{2}[\text{tr}^{2}(H)-\text{tr}(H^{2})]$,
etc. This implies that the denominator of eq. 5 can also be written as
$\Pi_{j\neq i}(\lambda_{i}-\lambda_{j})=\text{\bf Tr}[{\bf\Pi}_{j\neq
i}(H-\lambda_{j}I)]$.

For numerical calculations, the Le Verrier-Faddeev algorithm and using the
Lagrange expression require a similar numbers of complex matrix operations.
However, the Le Verrier-Faddeev algorithm gives coefficients of $\lambda^{k}$
for both the determinant and the adjugate of $(\lambda I-H)$ in one very
simple iterative procedure. For purely numerical calculations, especially if
the matrices are large, there exist many very efficient numerical packages.
However for analytic neutrino oscillation calculations where the number of
neutrino flavors is $\leq$6, we will now argue that using the adjugate of
$(\lambda I-H)$ is simpler, due to the form of the Hamiltonian, as we will see
in the next item.

  4. 4.

Adjugate method: for neutrino oscillation calculations one can use the special
form of the Hamiltonian in the flavor basis as follows:

$\displaystyle H$ $\displaystyle=H_{0}+H_{1}$ (13)

where $H_{0}$ is the vacuum Hamiltonian in the flavor basis, so that,

$\displaystyle H_{0}$ $\displaystyle=UXU^{\dagger}$
$\displaystyle\text{with}\quad X$ $\displaystyle=\text{\bf
diag}(x_{1},x_{2},\cdots,x_{n})$ (14)

where $x_{i}=m^{2}_{i}/(2E)$. The adjugate of $H_{0}$ is simply given by

$\displaystyle\text{\bf Adj}[H_{0}]_{\alpha\beta}$
$\displaystyle=\sum_{i}U_{\alpha i}U^{*}_{\beta i}\,\Pi_{j\neq i}x_{j}\,,$
(15)

which follows from $\text{\bf Adj}[UXU^{\dagger}]=U\text{\bf
Adj}[X]U^{\dagger}$. $H_{1}$ contains the matter potential which for most
cases is a sparse matrix in the flavor basis. For the standard 3 flavor case
only one diagonal element of $H_{1}$ is non-zero and for 3+1 flavors only 2
diagonal elements are non-zero. If $H_{1}$ is sparse, then calculation of the
$\text{\bf Adj}[H]$ is simple and straight forward. A similar argument applies
to $\text{\bf Det}[H]$.

However, what is actually needed is $\text{\bf Det}[\lambda I-H]$ and
$\text{\bf Adj}[\lambda I-H]$, but these can be trivially calculated from
$\text{\bf Det}[H]$ and $\text{\bf Adj}[H]$, since

$\displaystyle\lambda I-H=U(\lambda I-X)U^{\dagger}-H_{1}\,.$ (16)

Therefore, replacing $x_{i}\rightarrow(x_{i}-\lambda)$ in $\text{\bf Det}[H]$
and $\text{\bf Adj}[H]$ one obtains the required results, since

$\displaystyle\text{\bf Det}[\lambda I-H(x_{i})]$
$\displaystyle=(-1)^{n}~{}\text{\bf Det}[H(x_{i}-\lambda)]\,,$
$\displaystyle\text{\bf Adj}[\lambda I-H(x_{i})]$
$\displaystyle=(-1)^{n-1}~{}\text{\bf Adj}[H(x_{i}-\lambda)]\,.$ (17)

Therefore one only needs to calculate $\text{\bf Det}[H]$ and $\text{\bf
Adj}[H]$, a considerable simplification. We call the method given here, the
“adjugate method,” and it is the most straight forward and systematic for
analytic neutrino oscillation calculations. This happens because the full
Hamiltonian can be split into an easily diagonalized part, $H_{0}$, and an
additional piece, $H_{1}=H-H_{0}$, that is a sparse matrix.

All of these methods, as well as combinations, can be used to find the
required results, however, the last method given above, adjugate method,
requires the least algebraic computation and gives simple expressions for
neutrino oscillations because of the form of the Hamiltonian. This is the
method that will be used in the rest of this paper for the numerator of eq. 5.

A comment about the denominator of adjugate equation, eq. 5, is in order as
well. To evaluate this denominator it appears that you need all the
eigenvalues of H, not just the eigenvalue of the eigenvector one is
calculating. This is illusionary, as one can easily show that

$\displaystyle{\bf\Pi}_{j\neq i}(\lambda_{i}-\lambda_{j})$
$\displaystyle=\frac{\partial}{\partial\lambda}\,{\bf Det}[\lambda
I-H]\,\biggr{|}_{\lambda=\lambda_{i}}=n\lambda^{n-1}_{i}+(n-1)\lambda^{n-2}_{i}d_{1}+\cdots+d_{n-1}\,,$
(18)

see for example ref. Denton et al. (2022). ${\bf Det}[\lambda I-H]$ is just
the characteristic polynomial and the coefficients are given by the $d_{i}$’s
of eq. 8. Thus, if you know the position of a zero of the characteristic
polynomial and the slope of this polynomial at that zero, one has enough
information using eq. 5 to determine that eigenvector i.e. in physicist’s
language the determination of the eigenvector is local in $\lambda$, one does
not need to know all the other eigenvalues. That is, the adjugate equation,
eq. 5, could have been written as

$\displaystyle V_{\alpha i}V^{*}_{\beta i}=\biggr{(}\text{\bf Adj}[\lambda
I-H]_{\alpha\beta}\biggr{/}\frac{\partial}{\partial\lambda}\,{\bf Det}[\lambda
I-H]\biggr{)}\,\biggr{|}_{\lambda=\lambda_{i}}\,,$ (19)

which manifestly only depends on the single eigenvalue of $H$, $\lambda_{i}$,
associated with the eigenvector given by $V_{\alpha i}V^{*}_{\beta i}$. Of
course frequently one may need all eigenvalues, if so, then RHS of eq. 18
maybe the simplest route but it is not necessary for the determination of a
single eigenvector. In addition eq. 18 can also be used as a cross check.

Alternatively, since we require $\sum_{\alpha}|V_{\alpha i}|^{2}=1$, then the
denominator of eq. 5 is clearly equal to the $\text{\bf Tr}[\text{\bf
Adj}[\lambda_{i}I-H]]$, see appendix B, as all of these are equivalent since

${\bf\Pi}_{j\neq
i}(\lambda_{i}-\lambda_{j})=\frac{\partial}{\partial\lambda}\,{\bf
Det}[\lambda I-H]\,\biggr{|}_{\lambda=\lambda_{i}}=\text{\bf Tr}[\text{\bf
Adj}[\lambda_{i}I-H]]=\text{\bf Tr}[{\bf\Pi}_{j\neq i}(H-\lambda_{j}I)].$

###  III.1 Relationship to previous works

The method used in the papers of Kimura et al. (KTY), Kimura et al. (2002a) &
Kimura et al. (2002b), specifically for three neutrinos, is closely related
but not exactly the calculation of $\text{\bf Adj}(\lambda_{i}I-H)$ for a 3x3
matrix. In the KTY papers the adjugate equation, eq. 5, never appears, nor do
the words “adjugate” or “classical adjoint” or ”adjunct” appear, even though
some of the results of their papers are closely related to what we find for
three neutrinos in the next section. In fact, the results given in the KTY
papers are midway888It is easy to show that the numerator of eq. 39 of Kimura
et al Kimura et al. (2002b), the key equation of that paper, can be rewritten
with the eigenvalues $(\lambda_{i},\lambda_{j},\lambda_{k})$ with (i,j,k) all
different, as
$-\lambda_{i}(\lambda_{j}+\lambda_{k})I_{3}+\lambda_{i}H+H^{-1}\text{\bf
Det}(H)=\lambda_{j}\lambda_{k}I_{3}-(\lambda_{j}+\lambda_{k})H+H^{2}=\lambda^{2}_{i}I_{3}+\lambda_{i}(H-\text{\bf
tr(}H)I_{3})+\text{\bf Adj}(H).$ The left expression from KTY, depends
explicitly on all 3 eigenvalues, $(\lambda_{i},\lambda_{j},\lambda_{k})$, the
middle, which comes from ${\bf\Pi}_{j\neq i}(H-\lambda_{j}I)$, Harrison and
Scott (2002), depends explicitly on 2 eigenvalues,
$(\lambda_{j},\lambda_{k})$, whereas the right, which comes from $\text{\bf
Adj}(\lambda_{i}I-H)$, depends explicitly only on the eigenvalue,
$\lambda_{i}$. Here, (2E) is set to 1 for simplicity.  between the
calculations using $\text{\bf Adj}(\lambda_{i}I-H)$ given in the next section
and ${\bf\Pi}_{j\neq i}(H-\lambda_{j}I)$ of ref. Harrison and Scott (2002).
There are multiple paths to extending the KTY method to more than three
neutrinos, whereas $\text{\bf Adj}(\lambda_{i}I-H)$ and ${\bf\Pi}_{j\neq
i}(H-\lambda_{j}I)$ are uniquely determined. The methods based on the Caley-
Hamilton Theorem, ${\bf\Pi}^{n}_{j=1}(H-\lambda_{j}I)=0$, such as papers
Ohlsson and Snellman (2000), Ohlsson and Snellman (1999) & Bustamante (2019),
are more distantly related and are more complicated to implement.

Our method is very simple, general and systematic, and can be applied to more
than 3 neutrinos or to more complicated matter potentials than the method in
KTY. In some sense it can be considered a modification of the above methods
that allows a straight forward generalization to more than 3 neutrinos or more
complicated matter potentials.

##  IV 3 Flavors in Matter: revisited

In this section we have used the adjugate method to calculate of the
eigenvectors, i.e $V_{\alpha i}V^{*}_{\beta i}$, and to reproduces the known
results for 3 flavors. This not only demonstrates the simplicity of this
method but also by reproducing known results confirms its validity for this
important and non-trivial example.

For three flavors in matter the Hamiltonian is simply

$\displaystyle(2E)H_{3x3}$
$\displaystyle=UM^{2}U^{\dagger}+\text{Diag}(a,0,0)\,,$ (20)

where $U$ is the PMNS matrix, $M^{2}$ the mass matrix squared,
$M^{2}=\text{Diag}(m^{2}_{1},m^{2}_{2},m^{2}_{3})$ and
$a=2\sqrt{2}G_{F}N_{e}E$ is the Wolfenstein matter potential, Wolfenstein
(1978), where $E$ is the energy of the neutrino in the rest frame of the
matter. The term proportional to the number density of neutrons, has been
subtract from all diagonal elements of the Hamiltonian, eq. 20, as terms
proportional to the identity matrix in the Hamiltonian do not effect
oscillations. Now

$\displaystyle\text{Det}[(2E)H_{3x3}]$
$\displaystyle=\Pi_{i}\,m^{2}_{i}+a\sum_{i}|U_{ei}|^{2}\,\Pi_{j\neq
i}m^{2}_{j}\,.$ (21)

Useful identities for this calculation can be found in appendix E. Therefore
by replacing $m^{2}_{i}\rightarrow(m^{2}_{i}-\lambda)$ we obtain the
characteristic equation999This is equivalent to
$\lambda^{3}-\lambda^{2}Tr[(2E)H]+\lambda(Tr^{2}[(2E)H]-Tr[(2E)^{2}H^{2}])/2-\text{Det}[(2E)H]=0.$

$\displaystyle\text{\bf Det}[\lambda I-(2E)H_{3x3}]$
$\displaystyle=\lambda^{3}-\lambda^{2}\biggr{(}\sum_{i}m^{2}_{i}+a\biggr{)}+\lambda\biggr{(}\sum_{i>j}m^{2}_{i}m^{2}_{j}+a\sum_{i}m^{2}_{i}(1-|U_{ei}|^{2})\biggr{)}$
$\displaystyle\quad-\biggr{(}\Pi_{i}m^{2}_{i}+a\sum_{i}|U_{ei}|^{2}\,\Pi_{j\neq
i}m^{2}_{j}\biggr{)}=0\,.$ (22)

Note that if one sets $(m^{2}_{1},m^{2}_{2},m^{2}_{3})=(0,\Delta
m^{2}_{21},\Delta m^{2}_{31})$, eq. 22 reduces to the results of Zaglauer-
Schwarzer, Zaglauer and Schwarzer (1988), and there are significant
simplifications as any term with $m^{2}_{1}$ is zero. In this paper we prefer
not to do this so as to maintain the symmetry between the $m^{2}_{i}$. The
exact solutions to this cubic equation are given in appendix F.

To obtain the eigenvectors, we need both the numerator and the denominator of
the adjugate eqn. 5. The denominator can be simply evaluated using $\text{\bf
Det}^{\prime}[\lambda I-(2E)H_{3x3}]$, see eqn. 18, with $\lambda$ being the
eigenvalue of the eigenvector. Note it is linear in the matter potential if
one uses only the eigenvalue of the eigenvector being calculated.

The numerator is easily evaluated from $\text{\bf Adj}[(2E)H_{3x3}$, which is
super simple to calculate and given by

$\displaystyle\text{\bf Adj}[(2E)H_{3x3}]_{\alpha\beta}$
$\displaystyle=\sum_{i}U_{\alpha i}U^{*}_{\beta i}\,\Pi_{j\neq
i}\,m^{2}_{j}\,+a\left\\{\begin{array}[]{ll}0&\alpha=e~{}~{}or~{}~{}\beta=e\\\
\sum_{i}m^{2}_{i}|U_{\tau i}|^{2}&\alpha=\beta=\mu\\\ \sum_{i}m^{2}_{i}|U_{\mu
i}|^{2}&\alpha=\beta=\tau\\\ -\sum_{i}m^{2}_{i}U_{\alpha i}U^{*}_{\beta
i}&\alpha\neq\beta,\alpha\neq e,\beta\neq e,\end{array}\right.\,.$ (27)

Again each sum has three terms (two if $m^{2}_{1}=0$). The matter potential
independent part here is most easily calculated using the identity given in
eqn. 15 whereas the matter dependent part is trivial. Note, $\text{\bf
Adj}[(2E)H_{3x3}]$ is also linear in the matter potential.

The eigenvectors are then calculated as follows:

$\displaystyle V_{\alpha i}V^{*}_{\beta i}$ $\displaystyle=\frac{\text{\bf
Adj}[\lambda_{i}I-(2E)H_{3x3}]_{\alpha\beta}}{\Pi_{l\neq
i}(\lambda_{i}-\lambda_{l})}\equiv\frac{\lambda^{2}_{i}\delta_{\alpha\beta}-\lambda_{i}S_{\alpha\beta}+T_{\alpha\beta}}{\Pi_{l\neq
i}(\lambda_{i}-\lambda_{l})}\,,$ (28)

where $\text{\bf Adj}[\lambda I-(2E)H_{3x3}]$ can be easily obtained from
$\text{\bf Adj}[(2E)H_{3x3}]$ by replacing
$m^{2}_{i}\rightarrow(m^{2}_{i}-\lambda)$. Thus $S_{\alpha\beta}$ and
$T_{\alpha\beta}$ are simply,

$\displaystyle S_{\alpha\beta}$
$\displaystyle=\sum_{i}m^{2}_{i}(\delta_{\alpha\beta}-U_{\alpha i}U^{*}_{\beta
i})+a\delta_{\alpha\beta}(1-\delta_{\alpha e})\quad\text{and}\quad
T_{\alpha\beta}=\text{Adj}[(2E)H_{3x3}]_{\alpha\beta}\,.$ (29)

The diagonal terms ($\beta=\alpha$) are then easily obtained as

$\displaystyle|V_{\alpha i}|^{2}$
$\displaystyle=\frac{\lambda^{2}_{i}-\lambda_{i}S_{\alpha\alpha}+T_{\alpha\alpha}}{\Pi_{l\neq
i}(\lambda_{i}-\lambda_{l})}\quad\text{with}\quad
S_{\alpha\alpha}=\sum_{i}m^{2}_{i}(1-|U_{\alpha i}|^{2})+a(1-\delta_{\alpha
e})$ $\displaystyle\quad\text{and}\quad T_{\alpha\alpha}$
$\displaystyle=\sum_{i}U_{\alpha i}U^{*}_{\beta i}\,\Pi_{j\neq
i}m^{2}_{j}+a\left\\{\begin{array}[]{cl}0&\alpha=e\\\
\sum_{i}m^{2}_{i}|U_{\tau i}|^{2}&\alpha=\mu\\\ \sum_{i}m^{2}_{i}|U_{\mu
i}|^{2}&\alpha=\tau\\\ \end{array}\right.$ (33)

In the adjugate method where we only use the eigenvalue of the eigenvector we
are calculating, then both the numerator and denominator are quadratic
polynomials in the eigenvalue with coefficients that are linear in the matter
potential. Other formulations such as the Lagrange method do not have such a
simple dependence, but of course, they give the same numerical values.

###  IV.1 Matter potential invariant quantities for three flavors

The elements of the Hamiltonian given by eq. 20 are nearly all independent of
the matter potential, except for the $H_{ee}$ component. In this section we
will explain what the physics implications of this independence has on the
oscillation probabilities. First we note that because we can re-phase the
flavor states we must also take this into account when identify these physics
implications.

First, the diagonal elements of the Hamiltonian give the “mass” of the flavor
states:

$\displaystyle m^{2}_{\alpha}=\sum_{i}|U_{\alpha
i}|^{2}m^{2}_{i}+a\delta_{\alpha e}=\sum_{i}|V_{\alpha i}|^{2}\lambda_{i}\,.$
(34)

Note the “mass” of $\nu_{e}$ is the exception here and depends linearly on the
matter potential whereas the ”masses” of $\nu_{\mu}$ and $\nu_{\tau}$ are
independent of the matter potential.

Second the square of the off-diagonal elements of the Hamiltonian
($|H_{\alpha\beta}|^{2}$, $\alpha\neq\beta$) are related to the appearance
oscillation probabilities in the small L/E limit. Therefore in this limit the
coefficients of the $(L/E)^{2}$ terms in the oscillation probabilities, i.e.
CP or T conserving parts, must be invariant for $\alpha\neq\beta$, Harrison
and Scott (2002):

$\displaystyle\sum_{j>k}{\cal R}(V_{\alpha j}V^{*}_{\beta j}V^{*}_{\alpha
k}V_{\beta k})~{}(\Delta\lambda_{jk})^{2}$ $\displaystyle=\sum_{j>k}{\cal
R}(U_{\alpha j}U^{*}_{\beta j}U^{*}_{\alpha k}U_{\beta k})~{}(\Delta
m^{2}_{jk})^{2}\,.$ (35)

This identity identity is easily checked using the expressions given in the
previous subsection and can be used as cross check.

Another identity, which is not independent of the above, is associated with
the disappearance probabilities:

$\displaystyle\sum_{j>k}|V_{\alpha j}|^{2}|V_{\alpha
k}|^{2}(\Delta\lambda_{jk})^{2}$ $\displaystyle=\sum_{j>k}|U_{\alpha
j}|^{2}|U_{\alpha k}|^{2}(\Delta m^{2}_{jk})^{2}$ (36)

This identity guarantees that in $L/E\rightarrow 0$ limit, the disappearance
oscillation probabilities is independent of matter effects i.e. the the
coefficients of the terms of order $(L/E)^{2}$ are equal in matter and vacuum.
This identity follows from the matter potential independence of
$\sum_{\beta\neq\alpha}|H_{\alpha\beta}|^{2}$.

What about the intrinsic CP violating (CPV) or T violating part of the
appearance probabilities? The CP or T violating part of the oscillation
probabilities appears first an $(L/E)^{3}$ term given by

$\displaystyle\Im(V_{\alpha i}V^{*}_{\beta i}V^{*}_{\alpha j}V_{\beta
j})\Delta\lambda_{ij}\Delta\lambda_{ik}\Delta\lambda_{jk}$ (37)

and is also invariant. This is the Naumov-Harrison-Scott (NHS) invariant,
Naumov (1992) & Harrison and Scott (2000) which follows from the invariance of
$\Im(H_{e\mu}H_{\mu\tau}H_{\tau e})$. Note, rephasing of the flavor states
does not change this invariant. $\Re(H_{e\mu}H_{\mu\tau}H_{\tau e})$ gives no
new invariant. The NHS identity can be easily derived from our expressions as
follows: for $\alpha\neq\beta$ and $i\neq j$, we obtain

$\displaystyle V_{\alpha i}V^{*}_{\beta i}V^{*}_{\alpha j}V_{\beta j}$
$\displaystyle=\frac{\lambda_{i}\lambda_{j}|S_{\alpha\beta}|^{2}-\lambda_{i}S_{\alpha\beta}T^{*}_{\alpha\beta}-\lambda_{j}S^{*}_{\alpha\beta}T_{\alpha\beta}+|T_{\alpha\beta}|^{2}}{\Pi_{l\neq
i}(\lambda_{i}-\lambda_{l})~{}\Pi_{k\neq j}(\lambda_{j}-\lambda_{k})}$ (38)

and therefore the Jarlskog invariant, Jarlskog (1985), is given by

$\displaystyle\Im(V_{\alpha i}V^{*}_{\beta i}V^{*}_{\alpha j}V_{\beta j})$
$\displaystyle=\frac{(\lambda_{j}-\lambda_{i})~{}{\Im}(S_{\alpha\beta}T^{*}_{\alpha\beta})}{\Pi_{l\neq
i}(\lambda_{i}-\lambda_{l})~{}\Pi_{k\neq j}(\lambda_{j}-\lambda_{k})}$
$\displaystyle=\frac{\Delta m^{2}_{ij}\Delta m^{2}_{ik}\Delta
m^{2}_{jk}}{\Delta\lambda_{ij}\Delta\lambda_{ik}\Delta\lambda_{jk}}~{}~{}{\Im}(U_{\alpha
i}U^{*}_{\beta i}U^{*}_{\alpha j}U_{\beta j})$ (39)

where (i,j,k) are all different. This important identity as be used in many
places to understand the effects of matter on CP and T violating oscillation
probabilities, see for example Yokomakura et al. (2000), Parke and Weiler
(2001), Nunokawa et al. (2008), Wang and Zhou (2019) and Denton and Parke
(2019).

This next identity follows from the fact that matter potential is invariant
under rotations in the $\mu$ and $\tau$ sectors. It is simple algebra to show
that Kimura-Takamura-Yokomakura (KTY) identity, Kimura et al. (2002b),

$\displaystyle|V_{e1}|^{2}|V_{e2}|^{2}|V_{e3}|^{2}(\Delta\lambda_{21}\Delta\lambda_{31}\Delta\lambda_{32})^{2}$
(40)

is an invariant, i.e. independent of the matter potential. One only has to use
that the $\lambda_{i}$’s satisfy the 3-flavor characteristic equation, eq. 4
to prove this identity. When taking the square root of this identity care must
be taken with the sign (phase) which depends on conventions, see Denton and
Parke (2022).

The Toshev invariant, Toshev (1991), is obtained by dividing the NHS identity
by the square root of the KTY identity, i.e.

$\displaystyle\left(\frac{\Delta\lambda_{ij}\Delta\lambda_{ik}\Delta\lambda_{jk}}{|\Delta\lambda_{ij}\Delta\lambda_{ik}\Delta\lambda_{jk}|}\right)\frac{\Im(V_{\alpha
i}V^{*}_{\beta i}V^{*}_{\alpha j}V_{\beta j})}{|V_{e1}V_{e2}V_{e3}|}\,,$ (41)

where (i,j,k) are all different. Note, the piece of this expression in
$\left(\cdots\right)$ is either $\pm 1$. Some care is needed with the phases
of $V_{ei}$ and signs of $\Delta\lambda_{ij}$ when taking the square root KTY
identity, as noted in reference Denton and Parke (2022). Using the PDG PMNS
matrix with all $\Delta\lambda_{ij}>0$ one obtains the usual Toshev invariant
(the invariance of $\sin 2\theta_{23}\sin\delta$), however departures from the
PDG conventions can introduce an additional sign. Eq. 41 is the PMNS
convention independent Toshev invariant which will be useful if a different
PMNS convention than that of the PDG is used Denton and Pestes (2021).

##  V 4 Flavors in Matter:

For four flavors in matter the Hamiltonian is simply

$\displaystyle(2E)H_{4x4}$
$\displaystyle=UM^{2}U^{\dagger}+\text{Diag}(a,0,0,b)$ (42)

where $U$ is the 4 flavor PMNS matrix, $M^{2}$ the mass matrix squared,
$M^{2}=\text{Diag}(m^{2}_{1},m^{2}_{2},m^{2}_{3},m^{2}_{4})$, “a” is the
$N_{e}$ matter potential with $a=2\sqrt{2}G_{F}N_{e}E_{\nu}$ and “b” is the
$N_{n}$ matter potential with $b=\sqrt{2}G_{F}N_{n}E_{\nu}$ where $N_{e}$ and
$N_{n}$ are the number density for electrons and neutrons respectively 101010
$bI$ has been added to $H$..

The determinant is easily calculated as

$\displaystyle\text{Det}[(2E)\,H_{4x4}]=$
$\displaystyle\Pi_{i}m^{2}_{i}+a\sum_{i}|U_{ei}|^{2}\,\Pi_{j\neq
i}m^{2}_{j}+b\sum_{i}|U_{si}|^{2}\,\Pi_{j\neq i}m^{2}_{j}$
$\displaystyle+ab\sum_{i>j}|U_{ei}U_{sj}-U_{ej}U_{si}|^{2}\,\Pi_{k\neq
i,\,k\neq j}m^{2}_{k}\,.$ (43)

Here we use $|U_{\mu i}U_{\tau j}-U_{\mu j}U_{\tau
i}|^{2}=|U_{ek}U_{sl}-U_{el}U_{sk}|^{2}$ with $(k,l)\neq(i,j)$. Then the
characteristic equation

$\displaystyle\text{Det}[\lambda
I_{4}-(2E)H_{4x4}]=\lambda^{4}-f_{1}\lambda_{3}+f_{2}\lambda^{2}-f_{3}\lambda+f_{4}=0$
(44)

is easily obtained by making the $m^{2}_{i}\rightarrow(m^{2}_{i}-\lambda)$
substitution in 43, which gives

$\displaystyle f_{1}$ $\displaystyle=\sum_{i}m^{2}_{i}+a+b$ (45)
$\displaystyle f_{2}$
$\displaystyle=\sum_{i>j}m^{2}_{i}m^{2}_{j}+a\sum_{i}m^{2}_{i}(1-|U_{ei}|^{2})+b\sum_{i}m^{2}_{i}(1-|U_{si}|^{2})+ab$
(46) $\displaystyle f_{3}$ $\displaystyle=\sum_{i}\,\Pi_{j\neq i}m^{2}_{j}$
$\displaystyle+a\sum_{i>j}m^{2}_{i}m^{2}_{j}(1-|U_{ei}|^{2}-|U_{ej}|^{2})+b\sum_{i>j}m^{2}_{i}m^{2}_{j}(1-|U_{si}|^{2}-|U_{sj}|^{2})$
$\displaystyle+ab\sum_{i}m^{2}_{i}(1-|U_{ei}|^{2}-|U_{si}|^{2})$ (47)
$\displaystyle f_{4}$
$\displaystyle=\text{Det}[(2E)H_{4x4}]=\Pi_{i}m^{2}_{i}+a\sum_{i}|U_{ei}|^{2}\,\Pi_{j\neq
i}m^{2}_{j}+b\sum_{i}|U_{si}|^{2}\,\Pi_{j\neq i}m^{2}_{j}$
$\displaystyle\hskip
85.35826pt+ab\sum_{i>j}|U_{ei}U_{sj}-U_{ej}U_{si}|^{2}\,\Pi_{k\neq i,\,k\neq
j}m^{2}_{k}\,.$ (48)

See appendix E for some useful identities. This is in agreement with ref. Li
et al. (2018),

The solution to the 4 flavor characteristic equation is as follows:

$\displaystyle 0$
$\displaystyle=\lambda^{4}-f_{1}\lambda^{3}+f_{2}\lambda^{2}-f_{3}\lambda+f_{4}$
(49)
$\displaystyle=(\lambda^{2}-.\frac{1}{2}(f_{1}-\sqrt{f_{1}^{2}+4(y-f_{2})}~{})\lambda+\frac{1}{2}(y+\sqrt{y^{2}-4f_{4}}))$
(50)
$\displaystyle\times(\lambda^{2}-\frac{1}{2}(f_{1}+\sqrt{f_{1}^{2}+4(y-f_{2})}~{})\lambda+\frac{1}{2}(y-\sqrt{y^{2}-4f_{4}}))$
(51)

where

$\displaystyle
y^{3}-f_{2}y^{2}+(f_{3}f_{1}-4f_{4})y+f_{4}(4f_{2}-f^{2}_{1})-f^{2}_{3}=0\,.$
(52)

The solution to this cubic equation can be found in the appendix F. Therefore
the solutions are then simply given by the solutions of the quadratics

$\displaystyle\lambda^{2}-.\frac{1}{2}(f_{1}-\sqrt{f_{1}^{2}+4(y-f_{2})}~{})\lambda+\frac{1}{2}(y+\sqrt{y^{2}-4f_{4}})=0$
(53)
$\displaystyle\lambda^{2}-\frac{1}{2}(f_{1}+\sqrt{f_{1}^{2}+4(y-f_{2})}~{})\lambda+\frac{1}{2}(y-\sqrt{y^{2}-4f_{4}})=0\,.$
(54)

The solutions to four flavor cubic, eq. 52, is the challenging part to solving
the four flavor characteristic equation.

To obtain the eigenvectors, we first calculate the $\text{\bf
Adj}[(2E)H_{4x4}]$, which is given by

$\displaystyle\text{\bf Adj}[(2E)H_{4x4}]_{\alpha\beta}=\sum_{i}U_{\alpha
i}U^{*}_{\beta i}~{}\Pi_{j\neq i}m^{2}_{j}$ $\displaystyle\hskip
14.22636pt+a\begin{pmatrix}\begin{matrix}0\end{matrix}&\vline&0~{}~{}~{}~{}~{}~{}0~{}~{}~{}~{}~{}~{}0\\\
\hline\cr\begin{matrix}0\\\ 0\\\
0\end{matrix}&\vline&\begin{matrix}{\large\text{\bf Adj}[{\cal
W}_{\mu,\tau,s}]}\end{matrix}\end{pmatrix}_{\alpha\beta}+b\begin{pmatrix}\begin{matrix}{\large\text{\bf
Adj}[{\cal W}_{e,\mu,\tau}]}\end{matrix}&\vline&\begin{matrix}0\\\ 0\\\
0\end{matrix}\\\
\hline\cr\begin{matrix}0~{}~{}~{}~{}~{}0~{}~{}~{}~{}~{}0\end{matrix}&\vline&0\end{pmatrix}_{\alpha\beta}+ab\begin{pmatrix}0&\vline&0\hskip
28.45274pt0&\vline&0\\\ \hline\cr\begin{matrix}0\\\
0\end{matrix}&\vline&\begin{matrix}{\large\text{\bf Adj}[{\cal
W}_{\mu,\tau}]}\end{matrix}&\vline&\begin{matrix}0\\\ 0\end{matrix}\\\
\hline\cr 0&\vline&0\hskip 28.45274pt0&\vline&0\\\
\end{pmatrix}_{\alpha\beta}$ (55)

where ${\cal W}_{\alpha\beta}$ (${\cal W}_{\alpha\beta\gamma}$) is the matrix
formed by the $\alpha\beta$ ($\alpha\beta\gamma$) rows and columns of $(2E)H$
with a=b=0, i.e. only the vacuum components. Using
$W_{\alpha\beta}\equiv\sum_{i}U_{\alpha i}U^{*}_{\beta i}m^{2}_{i}$,

$\displaystyle{\cal W}_{\mu,\tau}$
$\displaystyle=\begin{pmatrix}W_{\mu\mu}&W_{\mu\tau}\\\
W_{\tau\mu}&W_{\tau\tau}\end{pmatrix}\quad\text{and}\quad{\cal
W}_{\mu,\tau,s}=\begin{pmatrix}W_{\mu\mu}&W_{\mu\tau}&W_{\mu s}\\\
W_{\tau\mu}&W_{\tau\tau}&W_{\tau s}\\\
W_{s\mu}&W_{s\tau}&W_{ss}\end{pmatrix}.$ (56)

Where $\text{\bf Adj}[{\cal W}_{\mu,\tau}]=\text{tr}({\cal
W}_{\mu,\tau})I-{\cal W}_{\mu,\tau}$ and the elements of $\text{\bf Adj}[{\cal
W}_{\mu,\tau,s}]$ are bilinear combinations of the W’s and the identity

$\displaystyle
W_{\alpha\beta}W_{\gamma\delta}-W_{\alpha\delta}W_{\gamma\beta}$
$\displaystyle=\sum_{i>j}m^{2}_{i}m^{2}_{j}(U_{\alpha i}U_{\gamma j}-U_{\alpha
j}U_{\gamma i})(U_{\beta i}U_{\delta j}-U_{\beta j}U_{\delta i})^{*}$ (57)

is useful. Note, if $\alpha=\gamma$ or $\beta=\delta$ this bilinear sum is
zero, so this combination never appears. Again see appendix E.

The eigenvectors are then calculated as follows:

$\displaystyle V_{\alpha i}V^{*}_{\beta i}$ $\displaystyle=\frac{\text{\bf
Adj}[\lambda_{i}I_{4}-(2E)H_{4x4}]_{\alpha\beta}}{\Pi_{j\neq
i}(\lambda_{i}-\lambda_{j})}$ (58)

with the numerator obtained from eq. 55 by making the
$m^{2}_{i}\rightarrow(m^{2}_{i}-\lambda)$ replacement, which also implies a
$W_{\alpha\beta}\rightarrow W_{\alpha\beta}-\lambda\delta_{\alpha\beta}$
replacement. Defining

$\eta_{\alpha\beta}\equiv 1-\delta_{\alpha\beta},$

the opposite of the kronecker delta, we have

$\displaystyle\text{\bf Adj}[\lambda I_{4}-H_{4x4}]_{\alpha\beta}$
$\displaystyle=\lambda^{3}\delta_{\alpha\beta}$
$\displaystyle-\lambda^{2}\biggr{(}\sum_{i}m^{2}_{i}(\delta_{\alpha\beta}-U_{\alpha
i}U^{*}_{\beta
i})+a\,\eta_{e\alpha}\,\delta_{\alpha\beta}+b\,\eta_{s\alpha}\,\delta_{\alpha\beta}\biggr{)}$
$\displaystyle+\lambda\biggr{(}\sum_{i>j}m^{2}_{i}m^{2}_{j}(\delta_{\alpha\beta}-U_{\alpha
i}U^{*}_{\beta i}-U_{\alpha j}U^{*}_{\beta
j})+a\,\eta_{e\alpha}\eta_{e\beta}\,({\cal W}_{\mu,\tau,s}-\text{Tr}[{\cal
W}_{\mu,\tau,s}]I_{3})_{\alpha\beta}$ $\displaystyle\hskip
56.9055pt+b\,\eta_{s\alpha}\eta_{s\beta}\,({\cal
W}_{e,\mu,\tau}-\text{Tr}[{\cal
W}_{e,\mu,\tau}]I_{3})_{\alpha\beta}+ab\,\eta_{e\alpha}\,\eta_{s\alpha}\,\delta_{\alpha\beta}\biggr{)}$
$\displaystyle-\biggr{(}\sum_{i}U_{\alpha i}U^{*}_{\beta i}~{}\Pi_{j\neq
i}m^{2}_{j}+a\,\eta_{e\alpha}\eta_{e\beta}\,\text{\bf Adj}[{\cal
W}_{\mu,\tau,s}]_{\alpha\beta}+b\,\eta_{s\alpha}\eta_{s\beta}\,\text{\bf
Adj}[{\cal W}_{e,\mu,\tau}]_{\alpha\beta}$ $\displaystyle\hskip 56.9055pt-
ab\,\eta_{e\alpha}\eta_{e\beta}\,\eta_{s\alpha}\eta_{s\beta}\,({\cal
W}_{\mu,\tau}-\text{Tr}[{\cal W}_{\mu,\tau}]I_{2})_{\alpha\beta}\biggr{)}\,,$
(59)

where

$\displaystyle(\text{Tr}[{\cal W}_{\mu,\tau,s}]I_{3}-{\cal
W}_{\mu,\tau,s})_{\alpha\beta}$ $\displaystyle=\sum
m^{2}_{i}(\delta_{\alpha\beta}(1-|U_{ei}|^{2})-U_{\alpha i}U^{*}_{\beta
i})\,,$ $\displaystyle(\text{Tr}[{\cal W}_{e,\mu,\tau}]I_{3}-{\cal
W}_{e,\mu,\tau})_{\alpha\beta}$ $\displaystyle=\sum
m^{2}_{i}(\delta_{\alpha\beta}(1-|U_{si}|^{2})-U_{\alpha i}U^{*}_{\beta i}),$
$\displaystyle(\text{Tr}[{\cal W}_{\mu,\tau}]I_{2}-{\cal
W}_{\mu,\tau})_{\alpha\beta}$ $\displaystyle=\sum
m^{2}_{i}(\delta_{\alpha\beta}(1-|U_{ei}|^{2}-|U_{si}|^{2})-U_{\alpha
i}U^{*}_{\beta i}).$

Everything is now available for calculating the eigenvectors, $V_{\alpha
i}V^{*}_{\beta i}$, as well as the oscillation probability using eq. 1.

###  V.1 Matter potential invariant quantities for four flavors and beyond:

As in the the 3x3 case, the CP conserving identities are given by

$\displaystyle\sum_{j>k}{\cal R}(V_{\alpha j}V^{*}_{\beta j}V^{*}_{\alpha
k}V_{\beta k})~{}(\Delta\lambda_{jk})^{2}$ $\displaystyle=\sum_{j>k}{\cal
R}(U_{\alpha j}U^{*}_{\beta j}U^{*}_{\alpha k}U_{\beta k})~{}(\Delta
m^{2}_{jk})^{2}\quad\alpha\neq\beta\,,$ $\displaystyle\sum_{j>k}|V_{\alpha
j}|^{2}|V_{\alpha k}|^{2}(\Delta\lambda_{jk})^{2}$
$\displaystyle=\sum_{j>k}|U_{\alpha j}|^{2}|U_{\alpha k}|^{2}(\Delta
m^{2}_{jk})^{2}$ (60)

These identities are easily derived from the matter potential independence of
$|H_{\alpha\beta}|^{2}$, $\alpha\neq\beta$ in the flavor basis, and are valid
when there is more than one sterile neutrino. In general these identities are
valid whatever the matter potential is only along the diagonal in the flavor
basis. As we will see in the next section, non-diagonal entries in the flavor
basis of matter potential to the Hamiltonian break these identities.

The CP violating identities can be derived using the matter potential
independence of

$\displaystyle\sum_{\gamma\neq\alpha,\beta}\Im(H_{\alpha\beta}H_{\beta\gamma}H_{\gamma\alpha})$
(61)

$\alpha\neq\beta$ and $\alpha,\beta,\gamma=e,\mu,\tau,s$. This is the
generalization of what was done for 3 flavors to 4 (or more flavors). This
leads to a set of matter potential invariants that can be written as

$\displaystyle\sum_{i>j}\Im(V_{\alpha i}V^{*}_{\beta i}V^{*}_{\alpha
j}V_{\beta j})~{}\Delta\lambda_{ij}\Delta\lambda_{ik}\Delta\lambda_{jk}$
$\displaystyle=\sum_{i>j}\Im(U_{\alpha i}U^{*}_{\beta i}U^{*}_{\alpha
j}U_{\beta j})~{}\Delta m^{2}_{ij}\Delta m^{2}_{ik}\Delta m^{2}_{jk}$ (62)

for any fixed k. Note, when i=k or j=k, that term vanishes.

Using the notation111111Some useful properties of the J’s are
$J^{\alpha\beta}_{ij}=-J^{\alpha\beta}_{ji}=-J^{\beta\alpha}_{ij}=J^{\beta\alpha}_{ji}$
and $J^{\alpha\beta}_{ij}=-\sum_{k\neq j}J^{\alpha\beta}_{ik}=\sum_{k\neq
i}J^{\alpha\beta}_{kj}$.  $J^{\alpha\beta}_{ij}(U)\equiv\Im(U_{\alpha
i}U^{*}_{\beta i}U^{*}_{\alpha j}U_{\beta j})$, then for 4 flavors the CP
violating piece of the appearance oscillation probability was given in Sec II
as

$\displaystyle\sum_{i>j}J^{\alpha\beta}_{ij}(V)\sin 2\Delta_{ij}=$
$\displaystyle\quad
4J^{\alpha\beta}_{32}(V)\sin\Delta_{32}\sin\Delta_{21}\sin\Delta_{13}$
$\displaystyle+4J^{\alpha\beta}_{42}(V)\sin\Delta_{42}\sin\Delta_{21}\sin\Delta_{14}$
$\displaystyle+4J^{\alpha\beta}_{43}(V)\sin\Delta_{43}\sin\Delta_{31}\sin\Delta_{14}$
(63)

using k=1. In L/E $\rightarrow 0$ limit, the coefficient of the $(L/E)^{3}$
terms is given by

$\displaystyle
J^{\alpha\beta}_{32}(V)\,\Delta\lambda_{21}\Delta\lambda_{31}\Delta\lambda_{32}+J^{\alpha\beta}_{42}(V)\,\Delta\lambda_{21}\Delta\lambda_{41}\Delta\lambda_{42}+J^{\alpha\beta}_{43}(V)\,\Delta\lambda_{41}\Delta\lambda_{31}\Delta\lambda_{43}$
(64)

and is therefore independent of matter potentials, i.e. is invariant. Such
identities are the 4-flavor equivalents to the NHS ID. There are 6 of them,
but there only 3 independent ones: if one chooses the independent ones to be
$(\alpha,\beta)=(e,\mu),(e,s)~{}\&~{}(\mu,s)$ then the others
$(e,\tau),(\mu,\tau)~{}\&~{}(\tau,s)$ are linear combinations as follows:
since

$\displaystyle J^{e\tau}_{ij}=-J^{e\mu}_{ij}-J^{es}_{ij},\quad
J^{\mu\tau}_{ij}=+J^{e\mu}_{ij}-J^{\mu s}_{ij},\quad J^{\tau
s}_{ij}=-J^{es}_{ij}-J^{\mu s}_{ij}\,,$ (65)

for $i\neq j$. This is confirmed by the fact that there are 3 independent
Dirac CP phases for 4 neutrino flavors. All of these identities have been
confirmed using the expressions from the previous sub-section.

For n-flavor of neutrinos, with n-3 sterile neutrinos, the equivalent NHS
invariant follows from the matter potential independence of

$\sum_{\gamma\neq\alpha,\beta}\Im(H_{\alpha\beta}H_{\beta\gamma}H_{\gamma\alpha}).$

The identities are

$\displaystyle\sum_{i>j}\Delta m^{2}_{ij}\Delta m^{2}_{jk}\Delta
m^{2}_{ki}~{}J^{\alpha\beta}_{ij}(U)$
$\displaystyle=\sum_{i>j}\Delta\lambda_{ij}\Delta\lambda_{jk}\Delta\lambda_{ki}\,J^{\alpha\beta}_{ij}(V)$
(66)

for k=(1, 2, $\cdots$, n) and $\alpha,\beta=(e,\mu,\tau,s_{1},\cdots)$.
Although there are a large number of identities here only (n-1)(n-2)/2 are
independent. A proof of this identity is given in appendix G. Even for one
sterile neutrino eq. 66 is more general than eq. 64, since k can take any
value from 1 to 4.

##  VI 3 flavors with off-diagonal non-standard interactions (NSI)

Here we sketch how the systematic methods presented earlier can be used for
the case when there are non-standard interactions in the off-diagonal
components of the Hamiltonian in the flavor basis121212For non-standard
interactions in the diagonal components the method advocated here is also
simple to apply..

###  VI.1 NSI e$\mu$

The Hamiltonian is

$\displaystyle(2E)H$
$\displaystyle=UM^{2}U^{\dagger}+\left(\begin{array}[]{ccc}a&b&0\\\
b^{*}&0&0\\\ 0&0&0\end{array}\right)\,,$ (70)

where in the usual notation $b\equiv a\,\epsilon_{e\mu}$ in this subsection.
The determinant is given by

$\displaystyle\text{\bf Det}[(2E)H]-\text{\bf Det}[(2E)H_{3x3}]$
$\displaystyle=-b(W_{\mu e}W_{\tau\tau}-W_{\mu\tau}W_{\tau
e})-b^{*}(W_{e\mu}W_{\tau\tau}-W_{\tau\mu}W_{e\tau})-bb^{*}W_{\tau\tau}\,,$
$\displaystyle=\sum_{P}m^{2}_{i}m^{2}_{j}(b\,U^{*}_{ek}U_{\mu
k}+b^{*}\,U_{ek}U^{*}_{\mu k})-bb^{*}\sum_{i}m^{2}_{i}|U_{\tau i}|^{2}$ (71)

where $\text{\bf Det}[(2E)H_{3x3}]$ is given by eq. 21. Then the
characteristic equation can be calculated using the replacement
$m^{2}_{i}\rightarrow(m^{2}_{i}-\lambda)$, giving

$\displaystyle\text{\bf Det}[\lambda I-(2E)H]$ $\displaystyle=\text{\bf
Det}[\lambda I-(2E)H_{3x3}]-\lambda(bW_{\mu
e}+b^{*}W_{e\mu})+bb^{*}W_{\tau\tau}=0\,,$ (72)

where $\text{\bf Det}[\lambda I-(2E)H_{3x3}]$ is given by eq. 22.

The eigenvectors can be determined from $\text{\bf Adj}[(2E)H]$ in a similar
fashion. First,

$\displaystyle\text{\bf Adj}[(2E)H]-\text{\bf Adj}[(2E)H_{3x3}]$
$\displaystyle\quad=-b\begin{pmatrix}0&W_{\tau\tau}&-W_{\mu\tau}\\\ 0&0&0\\\
0&-W_{\tau e}&W_{\mu e}\end{pmatrix}-b^{*}\begin{pmatrix}0&0&0\\\
W_{\tau\tau}&0&-W_{e\tau}\\\
-W_{\tau\mu}&0&W_{e\mu}\end{pmatrix}-bb^{*}\begin{pmatrix}0&0&0\\\ 0&0&0\\\
0&0&1\end{pmatrix}$

where $\text{\bf Adj}[(2E)H_{3x3}]$ is given by eq. 27. Then by using the
$m^{2}_{i}\rightarrow(m^{2}_{i}-\lambda)$ replacement, we have

$\displaystyle\text{ \bf Adj}[\lambda I-(2E)H]$ $\displaystyle.\quad=\text{\bf
Adj}[\lambda I-(2E)H_{3x3}]+\lambda\begin{pmatrix}0&b&0\\\ b^{*}&0&0\\\
0&0&0\end{pmatrix}+(\text{\bf Adj}[(2E)H]-\text{\bf Adj}[(2E)H_{3x3}])$ (73)

where $\text{\bf Adj}[(2E)H_{3x3}]$ and $\text{\bf Adj}[\lambda
I-(2E)H_{3x3}]$ are given by eq. 27 and 28, respectively.

Since, in the flavor basis, $H_{e\tau}$ and $H_{\mu\tau}$ are independent of
the matter potential, therefore we have

$\displaystyle\sum_{j>k}{\cal R}(V_{\alpha j}V^{*}_{\beta j}V^{*}_{\alpha
k}V_{\beta k})~{}(\Delta\lambda_{jk})^{2}$ $\displaystyle=\sum_{j>k}{\cal
R}(U_{\alpha j}U^{*}_{\beta j}U^{*}_{\alpha k}U_{\beta k})~{}(\Delta
m^{2}_{jk})^{2}$ (74)

is invariant for
$(\alpha,\beta)=(e,\tau),(\tau,e),(\mu,\tau)~{}\text{and}(\mu,\tau)$ but not
$(e,\mu)$ or $(\mu,e)$. Adding $(\tau,e)$ and $(\tau,\mu)$ one also obtains

$\displaystyle\sum_{j>k}|V_{\tau j}|^{2}|V_{\tau
k}|^{2}(\Delta\lambda_{jk})^{2}$ $\displaystyle=\sum_{j>k}|U_{\tau
j}|^{2}|U_{\tau k}|^{2}(\Delta m^{2}_{jk})^{2}$ (75)

as an invariant. So in the small L/E, the $(L/E)^{2}$ coefficients of the CP
conserving part of the oscillation probabilities are independent of the matter
potential for $\nu_{e}\leftrightarrow\nu_{\tau}$,
$\nu_{\mu}\leftrightarrow\nu_{\tau}$ and $\nu_{\tau}$ disappearance but not
for $\nu_{e}\leftrightarrow\nu_{\mu}$ or $\nu_{e}$ and $\nu_{\mu}$
disappearance.

For the equivalent of the NHS identity, or CPV invariant, the dependence of
$H_{e\mu}$ on the matter potential means that $\Im(H_{e\mu}H_{\mu\tau}H_{\tau
e})$ also depends on the matter potential and therefore in the small L/E
limit, the coefficient $(L/E)^{3}$ term of the CP violating part of the
appearance oscillation probabilities depends on the matter potential for all
channels, see also Fong (2022). This statement is true if the off-diagonal
matter effect appears in $H_{e\tau}$ or $H_{\mu\tau}$. So any off-diagonal
contribution to the matter potential in the flavor basis renders the NHS
expression matter potential dependent.

Similar results can be obtained for an NSI in the e$\tau$ sector.

###  VI.2 NSI $\mu\tau$

For an NSI component in the $\mu\tau$ sector, the Hamiltonian is given by

$\displaystyle(2E)H$
$\displaystyle=UM^{2}U^{\dagger}+\left(\begin{array}[]{ccc}a&0&0\\\ 0&0&b\\\
0&b^{*}&0\end{array}\right)\,,$ (79)

where in the usual notation $b\equiv a\epsilon_{\mu\tau}$ in this subsection.
The determinant and adjugate of the Hamiltonian are slightly different as
product terms between “a” and “b” appear:

$\displaystyle\text{\bf Det}[(2E)H]$ $\displaystyle=\text{\bf
Det}[(2E)H_{3x3}]$ $\displaystyle-b(W_{ee}W_{\tau\mu}-W_{e\mu}W_{\tau
e})-b^{*}(W_{ee}W_{\mu\tau}-W_{\mu e}W_{e\tau})$
$\displaystyle-a(bW_{\tau\mu}+b^{*}W_{\mu\tau})-bb^{*}(a+W_{ee})$ (80)

where $\text{\bf Det}[(2E)H_{3x3}]$ is given by eq. 21. Similarly

$\displaystyle\text{\bf Adj}[(2E)H]$ $\displaystyle=\text{\bf
Adj}[(2E)H_{3x3}]$ $\displaystyle-b\begin{pmatrix}W_{\tau\mu}&0&-W_{e\mu}\\\
-W_{\tau e}&0&W_{ee}\\\
0&0&0\end{pmatrix}-b^{*}\begin{pmatrix}W_{\mu\tau}&-W_{e\tau}&0\\\ 0&0&0\\\
-W_{\mu e}&W_{ee}&0\end{pmatrix}-a\begin{pmatrix}0&0&0\\\ 0&0&b\\\
0&b^{*}&0\end{pmatrix}-bb^{*}\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&0\end{pmatrix}\,.$ (81)

Both $\text{\bf Det}[\lambda I-(2E)H]$ and $\text{\bf Adj}[\lambda I-(2E)H]$
can be calculated using the $m^{2}_{i}\rightarrow(m^{2}_{i}-\lambda)$
procedure as before.

Using similar arguments as in $(e,\mu)$ case, in the small L/E, the
$(L/E)^{2}$ coefficients of the CP conserving part of the oscillation
probabilities are independent of the matter potential for
$\nu_{e}\leftrightarrow\nu_{\mu}$, $\nu_{e}\leftrightarrow\nu_{\tau}$ and
$\nu_{e}$ disappearance but not for $\nu_{\mu}\leftrightarrow\nu_{\tau}$ or
$\nu_{\mu}$ and $\nu_{\tau}$ disappearance.

If the NSI components appear in more than one off-diagonal component of the
Hamiltonian the calculation of $\text{\bf Det}[\lambda I-(2E)H]$ and
$\text{\bf Adj}[\lambda I-(2E)H]$ follows in a similar fashion but of course
the resulting equations are more complicated. For this case it maybe simpler
instead of using $\text{\bf Adj}[\lambda_{i}I-(2E)H]$ in the numerator of eq.
5 but to use ${\bf\Pi}_{j\neq i}(H-\lambda_{j}I)$ as for three neutrinos only
$H$ and $H^{2}$ appear.

##  VII Conclusion

Neutrino oscillations in matter is important for measuring the remaining
unknown parameters in the neutrino sector as well as for the discovery of
beyond the Standard Model physics in this sector. The LBL experiments in
particular will discover the final piece of the neutrino mass ordering,
$\Delta m^{2}_{32}>\text{or}<0$, determine whether the neutrino mass state
with least $\nu_{e}$ is dominated by $\nu_{\mu}$ or $\nu_{\tau}$, octant of
$\theta_{23}$, and determine the size and sign of CP violation in the neutrino
sector, i.e. determine the Jarlskog invariant. All of the current and future
LBL experiments, T2K, NOvA, HK, DUNE and KNO, have significant matter effects
which must be included in their analysis. Understanding these matter effects
for the Standard Model three flavor paradigm as well as for beyond the
Standard Model scenarios is important for extracting the most physics out of
these expensive experiments.

In this paper we have provided a new simple way to determine the mixing matrix
in constant density matter once the mass eigenvalues have been determined.
This method uses a known result in linear algebra that has not been used in
the neutrino literature before and has general applicability in particle
physics. We call this method the adjugate method as it uses the adjugate of
the Hamiltonian to determine the mixing matrix in matter. This method is
particularly well suited to neutrino oscillations due to the form of the
Hamiltonian but is also very general and can used for any Hamiltonian. Once
the mixing matrix in matter has been determined it is then very simple to
determine the oscillation probabilities for these LBL experiments. We have
also rewritten the CP violating part of the neutrino oscillation probability
in a way that better reflects the physics of this important part of the
oscillation probability, see Sec. II, and is especially useful in our
discussion of four or more flavor neutrino oscillations.

In Sec. IV, we have revisited matter effects for three flavor oscillations
that has previously been addressed by many authors. The derivation here is
particularly simple and straight forward such that we do not need the form of
the PMNS matrix to obtain results. We have also derived all of the useful
matter invariant quantities including the NHS identity eq. 37 and the KTY
identity eq. 40 as well as identities related to the disappearance
probabilities eq. 36 and the CPC part of the appearance probabilities, eq. 35.
The scenario with three active and one sterile is addressed in Sec. V. Again
we derive simple expressions for the mixing matrix in matter and give a number
of invariant quantities especially the generalization of the NHS identity to
four or more neutrinos. This is a sum of three neutrino factors which reduces
to three flavor result as the fourth neutrino is decoupled. In Sec. VI, we
have given a sketch of how to apply this method to Hamiltonians with NSI and
discuss what happens to the invariant quantities with different non-zero off-
diagonal elements in the Hamiltonian in the flavor basis.

The adjugate method for determining the mixing matrix, i.e. the eigenvectors
of the Hamiltonian, is of general applicability but is especially well suited
to determining neutrino oscillation probabilities in constant density matter
due to the form of the Hamiltonian. We have illuminated its use in a variety
of scenarios from the three flavor Standard Model as well as additional
sterile neutrinos and also with non-standard interactions. We have used this
method to illuminate a variety of quantities that are independent of the
matter potential in the various scenarios and discussed the physics of these
invariants. Extension of the adjugate method to non-normal Hamiltonians will
be the subject of a follow up paper.

###### Acknowledgements.

We are very grateful to Peter Denton for discussions, encouragement and
comments on this manuscript. Fermilab is operated by the Fermi Research
Alliance under contract no. DE-AC02-07CH11359 with the U.S. Department of
Energy. This project was initiated at the KITP’s program, “Neutrinos as a
Portal to New Physics and Astrophysics”, and was supported in part by the
National Science Foundation under Grant No. PHY-1748958. This project has
received funding/support from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No
860881-HIDDeN.

##  Appendix A Rewriting the CP violating term in Oscillation Probabilities

.

First we defining $J^{\alpha\beta}_{ij}$ and note some of J’s properties:

$\displaystyle J^{\alpha\beta}_{ij}$ $\displaystyle\equiv{\Im}(V_{\alpha
i}V^{*}_{\beta i}V^{*}_{\alpha j}V_{\beta j})$ $\displaystyle
J^{\alpha\beta}_{ij}$
$\displaystyle=-J^{\alpha\beta}_{ji}=-J^{\beta\alpha}_{ij}=J^{\beta\alpha}_{ji}$
$\displaystyle J^{\alpha\beta}_{kj}$ $\displaystyle=-\sum_{i\neq
k}J^{\alpha\beta}_{ij},\quad J^{\alpha\beta}_{ik}=-\sum_{j\neq
k}J^{\alpha\beta}_{ij},\quad J^{\alpha\beta}_{kk}=0\,.$ (82)

Using the above identities one can show that

$\displaystyle\sum_{i,j}J^{\alpha\beta}_{ij}\sin 2\Delta_{ij}$
$\displaystyle=\sum_{i\neq k,j\neq k}J^{\alpha\beta}_{ij}\sin
2\Delta_{ij}+\sum_{j\neq k}J^{\alpha\beta}_{kj}\sin 2\Delta_{kj}+\sum_{i\neq
k}J^{\alpha\beta}_{ik}\sin 2\Delta_{ik}+J^{\alpha\beta}_{kk}\sin 2\Delta_{kk}$
$\displaystyle\sum_{i,j}J^{\alpha\beta}_{ij}\sin 2\Delta_{ij}$
$\displaystyle=\sum_{i\neq k,j\neq k}J^{\alpha\beta}_{ij}\biggr{(}\sin
2\Delta_{ij}+\sin 2\Delta_{jk}+\sin 2\Delta_{ki}\biggr{)}$
$\displaystyle=4\sum_{i,j}J^{\alpha\beta}_{ij}\biggr{(}\sin\Delta_{ij}\sin\Delta_{jk}\sin\Delta_{ik}\biggr{)}\,.$
(83)

For the last line the trig identity

$\displaystyle\sin 2\Delta_{ij}+\sin 2\Delta_{jk}+\sin
2\Delta_{ki}=-4\sin\Delta_{ij}\sin\Delta_{jk}\sin\Delta_{ki}$ (84)

was used. Therefore, for any fixed k

$\displaystyle\sum_{i>j}J^{\alpha\beta}_{ij}\sin 2\Delta_{ij}$
$\displaystyle=4\sum_{i>j}J^{\alpha\beta}_{ij}\biggr{(}\sin\Delta_{ij}\sin\Delta_{jk}\sin\Delta_{ik}\biggr{)}\,.$
(85)

This is eq. 2.

##  Appendix B Adjugate “Proof”

To demonstrate

$\displaystyle V_{\alpha i}V^{*}_{\beta i}$ $\displaystyle=\frac{\text{\bf
Adj}[\lambda_{i}I-H]_{\alpha\beta}}{\Pi_{j\neq
i}(\lambda_{i}-\lambda_{j})}\,.$ (86)

Consider the case when $H=V\Lambda V^{\dagger}$ with $\Lambda=\text{\bf
Diag}(\lambda_{1},\lambda_{2},\cdots,\lambda_{n})$, then

$\displaystyle\text{\bf Adj}[\lambda_{i}I-H]_{\alpha\beta}$
$\displaystyle=V_{\alpha j}\,\text{\bf
Adj}[\lambda_{i}I-\Lambda]_{jk}\,V_{\beta k}^{*}$ $\displaystyle=\Pi_{j\neq
i}(\lambda_{i}-\lambda_{j})V_{\alpha i}\,V_{\beta i}^{*}\,$ (87)

where we have used $\text{\bf Adj}[UAU^{\dagger}]=\text{\bf
Adj}[U^{\dagger}]\,\text{\bf Adj}[A]\,\text{\bf Adj}[U]=U\,\text{\bf
Adj}[A]\,U^{\dagger}$, when $U$ is unitary matrix. qed.

Alternatively, eq. 86 can be rewritten as follows:

$\displaystyle V_{\alpha i}V^{*}_{\beta i}$
$\displaystyle=\lim_{\lambda\rightarrow\lambda_{i}}\,\left(\frac{\text{\bf
Adj}[\lambda I-H]_{\alpha\beta}}{\Pi_{j\neq
i}(\lambda-\lambda_{j})}\right)=\lim_{\lambda\rightarrow\lambda_{i}}\,\left(\frac{\text{\bf
Det}[\lambda I-H]\,(\lambda I-H)^{-1}_{\alpha\beta}}{\Pi_{j\neq
i}(\lambda-\lambda_{j})}\right)\,.$ (88)

Here, one needs to consider the limit ${\lambda\rightarrow\lambda_{i}}$ as the
matrix $(\lambda_{i}I-H)$ is not invertible as
$\text{\bf Det}[\lambda_{i}I-H]=0$, although the $\text{\bf
Adj}[\lambda_{i}I-H]$ is well defined. Now consider the case when $H=V\Lambda
V^{\dagger}$ $(H^{-1}=V\Lambda^{-1}V^{\dagger})$ with $\Lambda=\text{\bf
Diag}(\lambda_{1},\lambda_{2},\cdots,\lambda_{n})$, so that

$\text{\bf Det}[\lambda
I-H]=\Pi_{j}(\lambda-\lambda_{j})\quad\text{and}\quad(\lambda
I-H)^{-1}_{\alpha\beta}=\sum_{k}V_{\alpha
k}\frac{1}{(\lambda-\lambda_{k})}V^{*}_{\beta k}~{}.$

Then

$\displaystyle\frac{\text{\bf Adj}[\lambda_{i}I-H]_{\alpha\beta}}{\Pi_{j\neq
i}(\lambda_{i}-\lambda_{j})}$
$\displaystyle=\lim_{\lambda\rightarrow\lambda_{i}}\,\left(\frac{\Pi_{j}(\lambda-\lambda_{j})\,\sum_{k}V_{\alpha
k}\frac{1}{(\lambda-\lambda_{k})}V^{*}_{\beta k}}{\Pi_{j\neq
i}(\lambda-\lambda_{j})}\right)$
$\displaystyle=\lim_{\lambda\rightarrow\lambda_{i}}\,\left(\sum_{k}V_{\alpha
k}\frac{(\lambda-\lambda_{i})}{(\lambda-\lambda_{k})}V^{*}_{\beta k}\right)$
$\displaystyle=V_{\alpha i}V^{*}_{\beta i}$ (89)

qed

For any matrix, X, the following identity holds:

$\displaystyle X.\text{\bf Adj}[X]=\text{\bf Adj}[X].X=\text{\bf Det}[X].I\,.$
(90)

Thus by setting $X=\lambda_{i}I-H$, so that RHS is zero, demonstrates that

$\text{\bf Adj}[\lambda_{i}I-H]_{\alpha\beta}\propto V_{\alpha i}V^{*}_{\beta
i},$

since $HV=V\Lambda$ and $V^{\dagger}H=\Lambda V^{\dagger}$. Since we require
$\sum_{\alpha}|V_{\alpha i}|^{2}=1$, the normalization is obtained using

$\text{\bf Tr}[\text{\bf Adj}[\lambda_{i}I-H]]={\bf\Pi}_{j\neq
i}(\lambda_{i}-\lambda_{j})$

giving eq. 86.

##  Appendix C Le Verrier-Faddeev algorithm.

Following ref. Hou (1998), we can write for any $\lambda$ the following:

$\displaystyle\text{Det}(\lambda I-H)$
$\displaystyle=\lambda^{n}+\lambda^{n-1}d_{1}+\lambda^{n-2}d_{2}+\cdots+d_{n}\,,$
(91) $\displaystyle\text{Adj}(\lambda I-H)$
$\displaystyle\equiv\lambda^{n-1}A_{1}+\lambda^{n-2}A_{2}+\cdots+A_{n}\,,$
(92)

where the A’s and d’s are calculated recursively as follows

$\displaystyle A_{1}$ $\displaystyle=I\,,\quad
d_{1}=-\frac{1}{1}\text{Tr}(HA_{1})=-\text{Tr}(H)$ (93) $\displaystyle A_{2}$
$\displaystyle=HA_{1}+d_{1}I=H-\text{Tr}(H)I$ $\displaystyle d_{2}$
$\displaystyle=-\frac{1}{2}\text{Tr}(HA_{2})=\frac{1}{2}(\text{Tr}^{2}(H)-\text{Tr}(H^{2}))\,.$
(94)

If H is 2x2 matrix then $A_{2}=-\text{Adj}(H)$ and $d_{2}=\text{Det}(H)$.
Continuing

$\displaystyle A_{3}$ $\displaystyle=HA_{2}+d_{2}I=H^{2}+d_{1}H+d_{2}I\,,$
$\displaystyle d_{3}$
$\displaystyle=-\frac{1}{3}\text{Tr}(HA_{3})=-\frac{1}{6}(\text{Tr}^{3}(H)-3\text{Tr}(H)\text{Tr}(H^{2})+2\text{Tr}(H^{3})).$
(95)

If H is a 2x2 matrix, $A_{3}=0$ and $d_{3}=0$. For 3x3 matrix then
$A_{3}=\text{Adj}(H)$ and $d_{3}=-\text{Det}(H)$. Next iteration gives

$\displaystyle A_{4}$
$\displaystyle=HA_{3}+d_{3}I=H^{3}+d_{1}H^{2}+d_{2}H+d_{3}I\,,$ $\displaystyle
d_{4}$
$\displaystyle=\frac{1}{24}(\text{Tr}^{4}(H)-6\text{Tr}^{2}(H)\text{Tr}(H^{2})+3\text{Tr}^{2}(H^{2})+8\text{Tr}(H)\text{Tr}(H^{3})-6\text{Tr}(H^{4}))\,.$
(96)

If H is a 3x3 matrix or smaller, $A_{4}=0$ and $d_{4}=0$. For 4x4 matrix then
$A_{4}=-\text{Adj}(H)$ and $d_{4}=\text{Det}(H)$. Then,

$\displaystyle A_{5}$
$\displaystyle=HA_{4}+d_{4}I=H^{4}+d_{1}H^{3}+d_{2}H^{2}+d_{3}H+d_{4}I\,,$
$\displaystyle d_{5}$
$\displaystyle=-\frac{1}{120}[~{}\text{Tr}^{5}(H)-10\text{Tr}^{3}(H)\text{Tr}(H^{2})+15\text{Tr}(H)\text{Tr}^{2}(H^{2})+20\text{Tr}^{2}(H)\text{Tr}(H^{3})$
$\displaystyle\hskip
85.35826pt-20\text{Tr}(H^{2})\text{Tr}(H^{3})-30\text{Tr}(H)\text{Tr}(H^{4})+24\text{Tr}(H^{5})~{}].$
(97)

If H is a 4x4 matrix or smaller, $A_{5}=0$ and $d_{5}=0$. For 5x5 matrix then
$A_{5}=\text{Adj}(H)$ and $d_{5}=-\text{Det}(H)$. Next iteration gives

$\displaystyle A_{6}$
$\displaystyle=HA_{5}+d_{5}I==H^{5}+d_{1}H^{4}+d_{2}H^{3}+d_{3}H^{2}+d_{4}H+d_{5}I$
$\displaystyle d_{6}$
$\displaystyle=\frac{1}{720}[~{}\text{Tr}^{6}(H)-15~{}\text{Tr}^{4}(H)\text{Tr}(H^{2})+45~{}\text{Tr}^{2}(H)\text{Tr}^{2}(H^{2})-15~{}\text{Tr}^{3}(H^{2})$
$\displaystyle\quad\quad+40~{}\text{Tr}^{3}(H)\text{Tr}(H^{3})-120~{}\text{Tr}(H)\text{Tr}(H^{2})\text{Tr}(H^{3})+40~{}\text{Tr}^{2}(H^{3})$
(98)
$\displaystyle\quad\quad-90~{}\text{Tr}^{2}(H)\text{Tr}(H^{4})+90~{}\text{Tr}(H^{2})\text{Tr}(H^{4})+144~{}\text{Tr}(H)\text{Tr}(H^{5})-120~{}\text{Tr}(H^{6})~{}]$

If H is a 5x5 matrix or smaller, $A_{6}=0$ and $d_{6}=0$. For 6x6 matrix then
$A_{6}=-\text{Adj}(H)$ and $d_{6}=\text{Det}(H)$. All of the above satisfy
$\text{\bf Tr}[A_{m+1}]=(n-m)d_{m}$ as required.

For an nxn matrix this continues until

$\displaystyle A_{n}$ $\displaystyle=HA_{(n-1)}+d_{(n-1)}I$
$\displaystyle=\text{Adj}(-H)=(-1)^{n-1}\text{Adj}(H)$ (99) $\displaystyle
d_{n}$ $\displaystyle=-\frac{1}{n}\text{Tr}(HA_{n})=(-1)^{n}\text{Det}(H)\,.$
(100)

The next iteration returns zero for both $A_{n+1}$ and $d_{n+1}$, as follows

$\displaystyle A_{n+1}$
$\displaystyle=HA_{n}+d_{n}I=(-1)^{n-1}\biggr{(}H\text{Adj}(H)-\text{Det}(H)I\biggr{)}=0$
(101) $\displaystyle d_{n+1}$
$\displaystyle=-\frac{1}{n+1}\text{Tr}(HA_{n+1})=0\,$

thus terminating the series.131313Note for any $d_{k}$ if you put all traces
equal to “k”, then $d_{k}=(-1)^{k}$. This is equivalent to replacing H with
$I_{k}$, and is a useful but simple cross check.

This recursion method is simple to program for numerical studies. Also eq. 99
gives an alternative way to calculate the $\text{Adj}(H)$ for an nxn matrix in
terms of powers of H up to (n-1) plus their traces, this is equivalent to the
Caley-Hamilton identity.

It is worth noting that eq. can be written as follows:

$\displaystyle\text{\bf Det}(\lambda I-H)$
$\displaystyle=\lambda^{n}-\lambda^{n-1}\text{\bf
Det}_{1}(H)+\lambda^{n-2}\text{\bf Det}_{2}(H)+\cdots+(-1)^{n}\text{\bf
Det}_{n}(H)\,,$ (102) $\displaystyle\text{\bf Adj}(\lambda I-H)$
$\displaystyle=\lambda^{n-1}I-\lambda^{n-2}\text{\bf
Adj}_{2}(H)+\lambda^{n-3}\text{\bf Adj}_{3}(H)+\cdots+(-1)^{n-1}\text{\bf
Adj}_{n}(H)\,,$ (103)

where $\text{\bf Det}_{m}(H)$ and $\text{\bf Adj}_{m}(H)$ are defined as the
polynomials of $H$, that gives the determinant of H and the adjugate of $H$ is
if $H$ where an m x m matrix, i.e.

$\displaystyle\text{\bf Det}_{1}(H)$ $\displaystyle\equiv\text{\bf
Tr}(H)=\sum_{j}\lambda_{j}\ $ $\displaystyle\text{\bf Det}_{2}(H)$
$\displaystyle\equiv\frac{1}{2}(\text{\bf Tr}^{2}(H)-\text{\bf
Tr}(H^{2}))=\sum_{j1>j2}\lambda_{j1}\lambda_{j2}$ $\displaystyle\text{\bf
Det}_{3}(H)$ $\displaystyle\equiv\frac{1}{24}(\text{\bf Tr}^{3}(H)-3\text{\bf
Tr}(H)\text{\bf Tr}(H^{2})+2\text{\bf
Tr}(H^{3}))=\sum_{j1>j2>j3}\lambda_{j1}\lambda_{j2}\lambda_{j3}$
$\displaystyle\vdots\quad$ $\displaystyle\quad\quad\vdots$
$\displaystyle\text{\bf Det}_{n}(H)$ $\displaystyle\equiv\text{\bf
Det}(H)=\Pi_{j}\lambda_{j}$ (104)

and

$\displaystyle\text{\bf Adj}_{2}(H)$ $\displaystyle\equiv\text{\bf Tr}(H)I-H$
$\displaystyle\text{\bf Adj}_{3}(H)$ $\displaystyle\equiv H^{2}-\text{\bf
Tr}(H)H+\frac{1}{2}(\text{\bf Tr}^{2}(H)-\text{\bf Tr}(H^{2}))I$
$\displaystyle\vdots\quad$ $\displaystyle\quad\quad\vdots$
$\displaystyle\text{\bf Adj}_{n}(H)$ $\displaystyle\equiv\text{\bf Adj}(H)$
(105)

##  Appendix D Proof of $\text{Adj}(\lambda_{i}I-H)={\bf\Pi}_{j\neq
i}(H-\lambda_{j}I)$

Rearrange the LeVerrier-Faddeev expression for $\text{\bf Adj}(\lambda I-H)$
as a power series in the matrices $I,~{}H,~{}H^{2},\cdots,~{}H^{n-1}$:

$\displaystyle\text{\bf Adj}(\lambda I-H)$
$\displaystyle\equiv\lambda^{n-1}A_{1}+\lambda^{n-2}A_{2}+\cdots+A_{n}\,,$
$\displaystyle=\lambda^{n-1}I+\lambda^{n-2}(H+d_{1}I)+\lambda^{n-3}(H^{2}+d_{1}H+d_{2}I)$
$\displaystyle\quad+\cdots+(H^{n-1}+d_{1}H^{n-2}+\cdots+d_{n-1}I)\,.$

Grouping powers of $H$, we have

$\displaystyle\text{\bf Adj}(\lambda I-H)$
$\displaystyle=(\lambda^{n-1}+d_{1}\lambda^{n-2}+d_{2}\lambda^{n-3}+\cdots+d_{n-1})I$
$\displaystyle\quad+(\lambda^{n-2}+d_{1}\lambda^{n-3}+d_{2}\lambda^{n-4}+\cdots+d_{n-2})H$
$\displaystyle\quad\quad\quad\vdots$
$\displaystyle\quad+(\lambda^{2}+d_{1}\lambda+d_{2})H^{n-3}$
$\displaystyle\quad+(\lambda+d_{1})H^{n-2}$ $\displaystyle\quad+H^{n-1}$ (106)

Now for nxn Hamiltonian with eigenvalues $\lambda_{j}$, its easy to show that

$d_{m}=(-1)^{m}\sum_{j1>...>jm}\lambda_{j1}\cdots\lambda_{jm}$

as these are just coefficients of the characteristic equation $\text{\bf
Det}(\lambda I-H)=\Pi_{i}(\lambda-\lambda_{i})=0$.

If $\lambda_{i}$ is one of the eigenvalues of $H$,
$(\lambda_{1},\cdots,\lambda_{n})$, then one can show by induction that

$\displaystyle(\lambda_{i}+d_{1})$ $\displaystyle=-\sum_{j1\neq
i}\lambda_{j1}$ (107) $\displaystyle(\lambda^{2}_{i}+d_{1}\lambda_{i}+d_{2})$
$\displaystyle=\sum_{\begin{array}[]{c}j1>j2\\\\[-2.84526pt] \neq
i\end{array}}\lambda_{j1}\lambda_{j2}$ (110)
$\displaystyle(\lambda^{3}_{i}+d_{1}\lambda^{2}_{i}+d_{2}\lambda_{i}+d_{3})$
$\displaystyle=-.\sum_{\begin{array}[]{c}j1>j2>j3\\\\[-2.84526pt] \neq
i\end{array}}\lambda_{j1}\lambda_{j2}\lambda_{j3}$ (113)
$\displaystyle\quad\vdots$
$\displaystyle(\lambda^{n-1}_{i}+d_{1}\lambda^{n-2}_{i}+\cdots+d_{n-1})$
$\displaystyle=(-1)^{n-1}\sum_{\begin{array}[]{c}j1>...>j(n-1)\\\\[-2.84526pt]
\neq i\end{array}}\lambda_{j1}\cdots\lambda_{j(n-1)}$ (116)
$\displaystyle=(-1)^{n-1}\,\Pi_{j\neq i}\lambda_{j}\,,$ (117)

where $\lambda_{i}$ does not appear on the RHS of any of these identities. The
following identity is use repetitively:

$\displaystyle\sum_{\begin{array}[]{c}j1>...>jm\\\\[-2.84526pt] \neq
i\end{array}}\lambda_{j1}\cdots\lambda_{jm}$
$\displaystyle=\sum_{\begin{array}[]{c}j1>...>jm\\\\[-2.84526pt]
\end{array}}\lambda_{j1}\cdots\lambda_{jm}\hskip
14.22636pt-~{}~{}\lambda_{i}\sum_{\begin{array}[]{c}j1>...>j(m-1)\\\\[-2.84526pt]
\neq i\end{array}}\lambda_{j1}\cdots\lambda_{j(m-1)}\,.$ (123)

Therefore, using the fact that all the $\lambda$’s are eigenvalues of H, we
have that

$\displaystyle\text{\bf Adj}(\lambda_{i}I-H)$
$\displaystyle=H^{n-1}-H^{n-2}\sum_{j\neq
i}\lambda_{j}+H^{n-3}\sum_{\begin{array}[]{c}j1>j2\\\\[-2.84526pt] \neq
i\end{array}}\lambda_{j1}\lambda_{j2}$ (126)
$\displaystyle+\cdots+(-1)^{n-1}\,I_{n}\,\Pi_{j\neq i}\lambda_{j}\,,$
$\displaystyle={\bf\Pi}_{j\neq i}(H-\lambda_{j}I)\,.$ (127)

This expression can also be derived from the inverse of the Vandermonde matrix
141414https://proofwiki.org/wiki/Inverse_of_Vandermonde_Matrix.

Note, the LeVerrier-Faddeev recursion method gives both $\text{\bf
Adj}(\lambda I-H)$ and $\text{\bf Det}(\lambda I-H)$ for arbitrary $\lambda$,
where as $\text{\bf Adj}(\lambda_{i}I-H)={\bf\Pi}_{j\neq i}(H-\lambda_{j}I)$
requires all of $\\{\lambda_{1},\cdots,\lambda_{n}\\}$ to be eigenvalues of H.
Both calculations require (n-1) matrix multiplications, so are of similar
computational cost. However, the LeVerrier-Faddeev recursion method automates
the calculation of $\text{\bf Adj}(\lambda I-H)$ and $\text{\bf Det}(\lambda
I-H)$ and thus require very few lines of code to program.

As a cross check consider the case that $H=V\Lambda V^{\dagger}$, then

$\displaystyle\Pi_{j\neq i}(H-\lambda_{j}I)|_{\alpha\beta}$
$\displaystyle=V~{}\Pi_{j\neq
i}(\Lambda-\lambda_{j}I)~{}V^{\dagger}|_{\alpha\beta}=\sum_{k}V_{\alpha
k}~{}\Pi_{j\neq i}(\lambda_{k}-\lambda_{j})~{}V^{*}_{\beta k}$
$\displaystyle=V_{\alpha i}V^{*}_{\beta i}~{}\Pi_{j\neq
i}(\lambda_{i}-\lambda_{j})\,,$ (128)

as required.

As an observation note that

$\displaystyle\Pi_{j\neq i}(H-\lambda_{j}I)$
$\displaystyle=H^{n-1}-\biggr{(}\sum_{j\neq
i}\lambda_{j}\biggr{)}H^{n-2}+\cdots+(-1)^{n-2}\biggr{(}\lambda_{j1}\lambda_{j2}\cdots\lambda_{j(n-1)}\biggr{)}I.$
(129) $\displaystyle\Pi_{j\neq i}(\lambda_{i}-\lambda_{j})$
$\displaystyle=\lambda_{i}^{n-1}-\biggr{(}\sum_{j\neq
i}\lambda_{j}\biggr{)}\lambda^{n-2}_{i}+\cdots+(-1)^{n-2}\biggr{(}\lambda_{j1}\lambda_{j2}\cdots\lambda_{j(n-1)}\biggr{)}\,,$
(130)

where none of the $j$’s are equal i and the quantites in $(\cdots)$ are the
elementary symmetric polynomials of the (n-1) variables
$(\lambda_{1},\cdots\lambda_{n})$ that do not including $\lambda_{i}$, which
are sometimes written as
$e_{k}(\\{\lambda_{1},\cdots,\lambda_{n}\\}\backslash\lambda_{i})$ for
$k=1,\cdots,n-1$. Both the numerator and denominator of eq. 12 are identical
polynomials of $H$ and $\lambda_{i}$ respectively. As this is required for the
special case $H=\Lambda$, as the numerator and denominator must cancel.

Another informative cross check is given by the Caley-Hamilton identity, which
can be written as $\Pi_{j}(H-\lambda_{j}I)=0$, therefore

$\displaystyle\Pi_{j}(H-\lambda_{j}I)|_{\alpha\gamma}$
$\displaystyle=(H-\lambda_{i}I)_{\alpha\beta}~{}\Pi_{j\neq
i}(H-\lambda_{j}I)|_{\beta\gamma}\propto(H-\lambda_{i}I)_{\alpha\beta}V_{\beta
i}V^{*}_{\gamma i}=0$ $\displaystyle=\Pi_{j\neq
i}(H-\lambda_{j}I)|_{\alpha\beta}(H-\lambda_{i}I)_{\beta\gamma}~{}\propto
V_{\alpha i}V^{*}_{\beta i}(H-\lambda_{i}I)_{\beta\gamma}=0\,,$ (131)

as required.

##  Appendix E Some Useful Identities:

For n-flavors:

$\displaystyle
W_{\alpha\beta}W_{\gamma\delta}-W_{\alpha\delta}W_{\gamma\beta}$
$\displaystyle=\sum_{i>j}m^{2}_{i}m^{2}_{j}(U_{\alpha i}U_{\gamma j}-U_{\alpha
j}U_{\gamma i})(U_{\beta i}U_{\delta j}-U_{\beta j}U_{\delta i})^{*}.$ (132)

For 3-flavors, the unitarity of U gives

$\displaystyle(U_{\alpha i}U_{\gamma j}-U_{\alpha j}U_{\gamma i})$
$\displaystyle=\text{Det}[U]U^{*}_{\sigma k}.$ (133)

where $(\alpha,\gamma,\sigma)$ and (i,j,k) are all different, therefore

$\displaystyle
W_{\alpha\beta}W_{\gamma\delta}-W_{\alpha\delta}W_{\gamma\beta}$
$\displaystyle=\sum_{i>j}m^{2}_{i}m^{2}_{j}U_{\sigma k}U^{*}_{\rho k}$ (134)

where $(\alpha,\gamma,\sigma)$, $(\delta,\beta,\rho)$ and (i,j,k) are all
different.

For 4-flavors, the unitarity of U gives

$\displaystyle(U_{\alpha i}U_{\gamma j}-U_{\alpha j}U_{\gamma i})$
$\displaystyle=\text{Det}[U](U_{\beta k}U_{\delta l}-U_{\beta l}U_{\delta
k})^{*}$ (135)

where $(\alpha,\gamma,\beta,\delta)$ and (i,j,k,l) are all different.
Therefore

$\displaystyle|U_{\alpha i}U_{\gamma j}-U_{\alpha j}U_{\gamma i}|^{2}$
$\displaystyle=|U_{\beta k}U_{\delta l}-U_{\beta l}U_{\delta k}|^{2}$ (136)

since $|\,\text{Det}[U]\,|^{2}=1$.

##  Appendix F Solutions to a cubic equation

The exact solutions to the cubic equation

$\displaystyle\lambda^{3}-A\lambda^{2}+B\lambda-C=0$ (137)

are given by Barger et al. (1980); Zaglauer and Schwarzer (1988)

$\displaystyle\lambda_{1}$
$\displaystyle=\frac{1}{3}A-\frac{1}{3}\sqrt{A^{2}-3B}\left(S+\sqrt{3}\sqrt{1-S^{2}}\right)\,,$
$\displaystyle\lambda_{2}$
$\displaystyle=\frac{1}{3}A-\frac{1}{3}\sqrt{A^{2}-3B}\left(S-\sqrt{3}\sqrt{1-S^{2}}\right)\,,$
(138) $\displaystyle\lambda_{3}$
$\displaystyle=\frac{1}{3}A+\frac{2}{3}\sqrt{A^{2}-3B}\,S\,,$

with

$\displaystyle S$
$\displaystyle\equiv\cos\left\\{\frac{1}{3}\cos^{-1}\left[\frac{2A^{3}-9AB+27C}{2(A^{2}-3B)^{3/2}}\right]\right\\}\,,$
(139)

The terms $A$, $B$, and $C$ are the sum of the eigenvalues, the sum of the
products of the eigenvalues, and the triple product of the eigenvalues
respectively.

As an example of the analytic impenetrability of $S$, eq. 139, set
$A=\sum_{i}m^{2}_{i}$, $B=\sum_{i>j}m^{2}_{i}m^{2}_{j}$ and
$C=\Pi_{i}\,m^{2}_{i}$, then to recover the eigenvalues,
$(m^{2}_{1},m^{2}_{2},m^{2}_{3})$, from eq. 138 is a highly non-trivial
exercise.

##  Appendix G n-flavor NHS ID

Assuming all off-diagonal elements of the Hamiltonian in the flavor basis are
independent of the matter potential, then

$\displaystyle(2E)^{3}~{}\Im\biggr{[}\sum_{\gamma\neq\alpha,\beta}H_{\alpha\beta}H_{\beta\gamma}H_{\gamma\alpha}\biggr{]}$
$\displaystyle=\Im\biggr{[}\sum_{\gamma\neq\alpha,\beta}\biggr{(}\sum_{i}\lambda_{i}V_{\alpha
i}V^{*}_{\beta i}\biggr{)}~{}\biggr{(}\sum_{j}\lambda_{j}V_{\beta
j}V^{*}_{\gamma j}\biggr{)}~{}\biggr{(}\sum_{k}\lambda_{k}V_{\gamma
k}V^{*}_{\alpha k}\biggr{)}\biggr{]}$ (140)

is an invariant. Using

$\sum_{\gamma\neq\alpha,\beta}V^{*}_{\gamma j}V_{\gamma
k}=\delta_{jk}-V^{*}_{\alpha j}V_{\alpha k}-V^{*}_{\beta j}V_{\beta k}.$

Only the $\delta_{jk}$ contributes to the imaginary part as the rest lead to
real terms. With the notation $J^{\alpha\beta}_{ij}(V)\equiv{\Im}(V_{\alpha
i}V^{*}_{\beta i}V^{*}_{\alpha j}V_{\beta j})$, we obtain

$\displaystyle(2E)^{3}~{}\Im(\sum_{\gamma\neq\alpha,\beta}H_{\alpha\beta}H_{\beta\gamma}H_{\gamma\alpha})$
$\displaystyle=\sum_{i,j}\lambda_{i}\lambda^{2}_{j}~{}J^{\alpha\beta}_{ij}(V),$
$\displaystyle=\sum_{i>j}(\lambda_{i}-\lambda_{j})(\lambda_{j}-\lambda_{k})(\lambda_{k}-\lambda_{i})J^{\alpha\beta}_{ij}(V)$
(141)

for some fixed k, usually k=1. The last line is derived using the properties
of $J^{\alpha\beta}_{ij}$ and a similar procedure as in appendix A. No terms
proportional to $J^{\alpha\beta}_{ik}$ or $J^{\alpha\beta}_{kj}$ appear.
Therefore,

$\displaystyle\sum_{i>j}\Delta m^{2}_{ij}\Delta m^{2}_{jk}\Delta
m^{2}_{ki}~{}J^{\alpha\beta}_{ij}(U)$
$\displaystyle=\sum_{i>j}\Delta\lambda_{ij}\Delta\lambda_{jk}\Delta\lambda_{ki}~{}J^{\alpha\beta}_{ij}(V)$
(142)

is the NHS identity for arbitrary number of sterile neutrinos. There are
(n-1)(n-2)/2 terms on both sides of this identity, none of them contain
$J^{\alpha\beta}_{ik}$ or $J^{\alpha\beta}_{kj}$.

## References

  * Wolfenstein (1978) L. Wolfenstein, Phys. Rev. D17, 2369 (1978), [,294(1977)].
  * Mikheyev and Smirnov (1985) S. P. Mikheyev and A. Yu. Smirnov, Sov. J. Nucl. Phys. 42, 913 (1985), [305(1986)].
  * Parke (1986) S. J. Parke, Phys. Rev. Lett. 57, 1275 (1986).
  * Cleveland et al. (1995) B. T. Cleveland, T. Daily, R. Davis, Jr., J. Distel, K. Lande, C. K. Lee, P. Wildenhain, and J. Ullman, Nucl. Phys. B Proc. Suppl. 38, 47 (1995).
  * Ahmad et al. (2002) Q. R. Ahmad et al. (SNO), Phys. Rev. Lett. 89, 011301 (2002), eprint nucl-ex/0204008.
  * Workman et al. (2022) R. L. Workman et al. (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Abe et al. (2011) K. Abe et al. (T2K), Nucl. Instrum. Meth. A 659, 106 (2011), eprint 1106.1238.
  * Ayres et al. (2004) D. S. Ayres et al. (NOvA) (2004), eprint hep-ex/0503053.
  * Abe et al. (2018) K. Abe et al. (Hyper-Kamiokande), PTEP 2018, 063C01 (2018), eprint 1611.06118.
  * Acciarri et al. (2015) R. Acciarri et al. (DUNE) (2015), eprint 1512.06148.
  * Jarlskog (1985) C. Jarlskog, Phys. Rev. Lett. 55, 1039 (1985).
  * Mena and Parke (2004) O. Mena and S. J. Parke, Phys. Rev. D 70, 093011 (2004), eprint hep-ph/0408070.
  * Kelly and Parke (2018) K. J. Kelly and S. J. Parke, Phys. Rev. D 98, 015025 (2018), eprint 1802.06784.
  * King et al. (2020) S. F. King, S. Molina Sedgwick, S. J. Parke, and N. W. Prouse, Phys. Rev. D 101, 076019 (2020), eprint 2001.05505.
  * Nunokawa et al. (2008) H. Nunokawa, S. J. Parke, and J. W. F. Valle, Prog. Part. Nucl. Phys. 60, 338 (2008), eprint 0710.0554.
  * Barger et al. (1980) V. D. Barger, K. Whisnant, S. Pakvasa, and R. J. N. Phillips, Phys. Rev. D 22, 2718 (1980).
  * Zaglauer and Schwarzer (1988) H. W. Zaglauer and K. H. Schwarzer, Z. Phys. C 40, 273 (1988).
  * Harrison and Scott (2000) P. F. Harrison and W. G. Scott, Phys. Lett. B476, 349 (2000), eprint hep-ph/9912435.
  * Ohlsson and Snellman (2000) T. Ohlsson and H. Snellman, J. Math. Phys. 41, 2768 (2000), [Erratum: J.Math.Phys. 42, 2345 (2001)], eprint hep-ph/9910546.
  * Ohlsson and Snellman (1999) T. Ohlsson and H. Snellman, Phys. Rev. D 60, 093007 (1999), eprint hep-ph/9903252.
  * Kimura et al. (2002a) K. Kimura, A. Takamura, and H. Yokomakura, Phys. Lett. B 537, 86 (2002a), eprint hep-ph/0203099.
  * Kimura et al. (2002b) K. Kimura, A. Takamura, and H. Yokomakura, Phys. Rev. D 66, 073005 (2002b), eprint hep-ph/0205295.
  * Barenboim et al. (2019) G. Barenboim, P. B. Denton, S. J. Parke, and C. A. Ternes, Phys. Lett. B 791, 351 (2019), eprint 1902.00517.
  * Denton et al. (2022) P. B. Denton, S. J. Parke, T. Tao, and X. Zhang, Bull. Am. Math. Soc. 59, 31 (2022), eprint 1908.03795.
  * Denton et al. (2020) P. B. Denton, S. J. Parke, and X. Zhang, Phys. Rev. D 101, 093001 (2020), eprint 1907.02534.
  * Hou (1998) S.-H. Hou, SIAM Rev. 40, 706 (1998).
  * Harrison and Scott (2002) P. F. Harrison and W. G. Scott, Phys. Lett. B 535, 229 (2002), eprint hep-ph/0203021.
  * Bustamante (2019) M. Bustamante (2019), eprint 1904.12391.
  * Naumov (1992) V. A. Naumov, Int. J. Mod. Phys. D1, 379 (1992).
  * Yokomakura et al. (2000) H. Yokomakura, K. Kimura, and A. Takamura, Phys. Lett. B 496, 175 (2000), eprint hep-ph/0009141.
  * Parke and Weiler (2001) S. J. Parke and T. J. Weiler, Phys. Lett. B 501, 106 (2001), eprint hep-ph/0011247.
  * Wang and Zhou (2019) X. Wang and S. Zhou, JHEP 05, 035 (2019), eprint 1901.10882.
  * Denton and Parke (2019) P. B. Denton and S. J. Parke, Phys. Rev. D 100, 053004 (2019), eprint 1902.07185.
  * Denton and Parke (2022) P. B. Denton and S. J. Parke, Phys. Rev. D 105, 013002 (2022), eprint 2106.12436.
  * Toshev (1991) S. Toshev, Mod. Phys. Lett. A 6, 455 (1991).
  * Denton and Pestes (2021) P. B. Denton and R. Pestes, JHEP 05, 139 (2021), eprint 2006.09384.
  * Li et al. (2018) W. Li, J. Ling, F. Xu, and B. Yue, JHEP 10, 021 (2018), eprint 1808.03985.
  * Fong (2022) C. S. Fong (2022), eprint 2210.09436.