Neutrino Oscillations in Matter using the Adjugate of the Hamiltonian

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.


I. INTRODUCTION
The Wolfenstein matter effect [1] has been and continues to be very important for exploring neutrino mixing and masses.It is the solution to the solar neutrino anomaly [2,3], discovered by Davis et.al. [4] and in particular the SNO experiment [5] determined |U e2 | 2 (sin 2 θ 12 ≈ 0.3) and that the sign 1 of ∆m 2  21 was positive , i.e. that the state with most ν e , ν 1 , is lighter than state with the next most ν e , ν 2 .The current long baseline (LBL) experiments, T2K [7] and NOvA [8], as well as the future experiments, Hyper-Kamiokande (HK) [9], Korean Neutrino observatory (KNO) [9] and DUNE [10], will determine the parameters |U µ3 | 2 (≈ sin 2 θ 23 ) and ∆m 2  32 including its sign.If |U µ3 | 2 > 1/2 (|U µ3 | 2 < 1/2) then ν µ (ν τ ) dominates the neutrino mass eigenstate with the least ν e , ν 3 .The sign of ∆m 2  32 determines the atmospheric mass ordering, whether ν 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 [11] (or the phase δ 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, [12].Although the density of the earth's crust varies somewhat over the baseline of all the above LBL experiments it has been shown in [13,14] 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 [15].
The neutrino oscillation probabilities for three neutrino flavors in constant medium 1 Our labelling of the mass eigenstates is determined by 2 and the form of the Pontecorvo, Maki, Nakagawa, Sakata (PMNS) matrix, U, as in the PDG [6].
have been exactly calculated 2 by many authors: Barger et al [16], Zaglauer & Schwarzer [17], Harrison & Scott [18], Ohlsson & Snellman [19,20], Kimura et al [21,22] as well as many approximation schemes reviewed in [23].All of the exact methods first determine the squared masses of the neutrinos in the medium by solving the characteristic equation, Det(λI − H) = 0, where λ'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, 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. 2Use of efficient numerical packages, while excellent for generating figures and doing analysis, they do not provide detailed physics insight or understanding.
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, λ i , of the Hamiltonian in the flavor bases and V αi V * βi of the normalized eigenvectors, V αi , no matter how complicated the Hamiltonian.The λ i 's and the V αi 's are the medium equivalents of m 2 i and elements of the PMNS matrix U αi in vacuum.In general the Hamiltonian is not the same for neutrinos and anti-neutrinos, so the λ i and V αi V * βi need to be calculated for both.With these qualities the neutrino oscillation probabilities are given by with ∆ ij ≡ (λ i − λ 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 α = β.The usual way of writing this term, as in the PDG, 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 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 λ i and V αi V * β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 λ i is well known, the closely associated method for calculating V αi V * β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 λ i 's are known, there is a simple matrix expression using the adjugate of a matrix that relates the λ i 's and the Hamiltonian in flavor basis to the V αi V * βi 's.

III. DIAGONALIZING A HERMITIAN HAMILTONIAN
Let H be an n × n Hermitian matrix that we want to diagonalize as follows: where Λ = Diag(λ 1 , λ 2 , • • • , λ n ) and V αi is a unitary matrix such that i-th column of V is a normalized eigenvector of H with eigenvalue λ i , since HV = V Λ.
It is well known that the eigenvalues are obtained be solving the characteristic equation For n ≤ 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 th eigenvector are easily calculated using where Adj is the adjugate of the matrix 3 , see [24], i.e. the transpose of the co-factor matrix 3 See 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 and references therein 4 .For neutrino oscillations one needs exactly the quantities V αi V * β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 αi V * βi given by eq. 5 is independent of this arbitrary choice but can be used to calculate the magnitude, |V αi |, and relative phase of each component of an eigenvector, since this phase is given by V αi V * βi /(|V αi ||V βi |).A physicist's style sketch of why eq.5 is correct is given in appendix B.
With λ i 's and the V αi V * βi 's constructed this way, eq. 3 is satisfied and 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: [25] & [24], where h α is the principal minor of H with the α-th row and column deleted.λ l (h α ) are the eigenvalues of h α , i.e. the solutions to Det[λ l I − h α ] = 0. Constructed this way i.e. they are normalized appropriately.Eq. 7 is mathematically equivalent to eq. 5 when α = β, 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 [24], 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 αi V * βi , for neutrino oscillation calculations. 4Note, 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. 5While the first of these is easily checked from eq. 5, the latter is most easily checked in numerical Therefore, for the evaluation of the eigenvalues we need, Det(λI − H), and for the eigenvectors, Adj(λI − H) is also needed, see eq. 4 & eq.5.We first note that both are polynomials in λ of degree n and (n-1) respectively, i.e.
where the d i 's are scalars and the A i 's are n × n Hermitian matrices.
There are at least five of different algebraic methods that can be used to calculate we review them here: 1. Using the Leibniz formulae 6 to calculate the determinant and elements of the adjugate of (λI −H).However this standard method does not, in general, give the simplest results without significant additional manipulations.
2. Le Verrier-Faddeev algorithm: a more general method is the very simple and recursive Le Verrier-Faddeev algorithm 7 , see also [26], given by The next iteration gives for any λ.
3. Lagrange method: the results of the Le Verrier-Faddeev algorithm can be rearranged, see appendix D, as This is the Lagrange method that has been used in many papers: [27], [38].The LHS of this equation only depends on the eigenvalue λ i whereas the RHS depends on all the other eigenvalues, λ j 's with j = 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.
etc.This implies that the denominator of eq. 5 can also be written as 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 λ k for both the determinant and the adjugate of (λ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 ≤6, we will now argue that using the adjugate of (λI − H) is simpler, due to the form of the Hamiltonian, as we will see in the next item.
4. Adjugate method: for neutrino oscillation calculations one can use the special form of the Hamiltonian in the flavor basis as follows: where H 0 is the vacuum Hamiltonian in the flavor basis, so that, where The adjugate of H 0 is simply given by which follows from Therefore, replacing and Adj[H] one obtains the required results, since Therefore one only needs to calculate Det[H] and 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 see for example ref. [24].Det[λ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 λ, one does not need to know all the other eigenvalues.
That is, the adjugate equation, eq. 5, could have been written as which manifestly only depends on the single eigenvalue of H, λ i , associated with the eigenvector given by V αi V * β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 α |V αi | 2 = 1, then the denominator of eq. 5 is clearly equal to the Tr[Adj[λ i I − H]], see appendix B, as all of these are equivalent since

A. Relationship to previous works
The method used in the papers of Kimura et al. (KTY), [21] & [22], specifically for three neutrinos, is closely related but not exactly the calculation of Adj(λ 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 midway 8 between the calculations using Adj(λ i I − H) given in the next section and Π j =i (H − λ j I) of ref. [27].There are multiple paths to extending the KTY method to more than three neutrinos, whereas Adj(λ i I − H) and Π j =i (H − λ j I) are uniquely determined.The methods based on the Caley-Hamilton Theorem, Π n j=1 (H − λ j I) = 0, such as papers [19], [20] & [28], 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. 8It is easy to show that the numerator of eq.39 of Kimura et al [22], the key equation of that paper, can be rewritten with the eigenvalues (λ i , λ j , λ k ) with (i,j,k) all different, as The left expression from KTY, depends explicitly on all 3 eigenvalues, (λ i , λ j , λ k ), the middle, which comes from Π j =i (H − λ j I), [27], depends explicitly on 2 eigenvalues, (λ j , λ k ), whereas the right, which comes from Adj(λ i I − H), depends explicitly only on the eigenvalue, λ i .Here, (2E) is set to 1 for simplicity.

IV. 3 FLAVORS IN MATTER: REVISITED
In this section we have used the adjugate method to calculate of the eigenvectors, i.e V αi V * β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 where U is the PMNS matrix, M 2 the mass matrix squared, 3 ) and a = 2 √ 2G F N e E is the Wolfenstein matter potential, [1], 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 Useful identities for this calculation can be found in appendix E. Therefore by replacing Note that if one sets (m 2 1 , m 2 2 , m 2 3 ) = (0, ∆m 2 21 , ∆m 2 31 ), eq.22 reduces to the results of Zaglauer-Schwarzer, [17], 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 see eqn. 18, with λ 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 Adj[(2E)H 3x3 , which is super simple to calculate and given by 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, Adj[(2E)H 3x3 ] is also linear in the matter potential.
The eigenvectors are then calculated as follows: where Adj[λI − (2E)H 3x3 ] can be easily obtained from Adj[(2E)H 3x3 ] by replacing m 2 i → (m 2 i − λ).Thus S αβ and T αβ are simply, The diagonal terms (β = α) are then easily obtained as and 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.

A. 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: Note the "mass" of ν e is the exception here and depends linearly on the matter potential whereas the "masses" of ν µ and ν τ are independent of the matter potential.
Second the square of the off-diagonal elements of the Hamiltonian (|H αβ | 2 , α = β) 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 α = β, [27]: 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: This identity guarantees that in L/E → 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 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 and is also invariant.This is the Naumov-Harrison-Scott (NHS) invariant, [29] & [18] which follows from the invariance of ℑ(H eµ H µτ H τ e ).Note, rephasing of the flavor states does not change this invariant.ℜ(H eµ H µτ H τ e ) gives no new invariant.The NHS identity can be easily derived from our expressions as follows: for α = β and i = j, we obtain and therefore the Jarlskog invariant, [11], is given by 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 [30], [31], [15], [32] and [33].
This next identity follows from the fact that matter potential is invariant under rotations in the µ and τ sectors.It is simple algebra to show that Kimura-Takamura-Yokomakura (KTY) identity, [22], is an invariant, i.e. independent of the matter potential.One only has to use that the λ 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 [34].
The Toshev invariant, [35], is obtained by dividing the NHS identity by the square root of the KTY identity, i.e.
where (i,j,k) are all different.Note, the piece of this expression in (• • • ) is either ±1.Some care is needed with the phases of V ei and signs of ∆λ ij when taking the square root KTY identity, as noted in reference [34].Using the PDG PMNS matrix with all ∆λ ij > 0 one obtains the usual Toshev invariant (the invariance of sin 2θ 23 sin δ), however departures from the PDG conventions can introduce an additional sign.Eq. 34 is the PMNS convention independent Toshev invariant which will be useful if a different PMNS convention than that of the PDG is used [36].

V. 4 FLAVORS IN MATTER:
For four flavors in matter the Hamiltonian is simply where U is the 4 flavor PMNS matrix, M 2 the mass matrix squared, The determinant is easily calculated as Here we use is easily obtained by making the m 2 i → (m 2 i − λ) substitution in 36, which gives See appendix E for some useful identities.This is in agreement with ref. [37], The solution to the 4 flavor characteristic equation is as follows: where 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 The solutions to four flavor cubic, eq.45, is the challenging part to solving the four flavor characteristic equation.
To obtain the eigenvectors, we first calculate the Adj[(2E)H 4x4 ], which is given by where W αβ (W αβγ ) is the matrix formed by the αβ (αβγ) rows and columns of (2E)H with a=b=0, i.e.only the vacuum components.Using Where Adj[W µ,τ ] = tr(W µ,τ )I − W µ,τ and the elements of Adj[W µ,τ,s ] are bilinear combinations of the W's and the identity is useful.Note, if α = γ or β = δ this bilinear sum is zero, so this combination never appears.Again see appendix E.
The eigenvectors are then calculated as follows: with the numerator obtained from eq.48 by making the m 2 i → (m 2 i − λ) replacement, which also implies a W αβ → W αβ − λδ αβ replacement.Defining the opposite of the kronecker delta, we have where Everything is now available for calculating the eigenvectors, V αi V * βi , as well as the oscillation probability using eq. 1.
using k=1.In L/E → 0 limit, the coefficient of the (L/E) 3 terms is given by 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 (α, β) = (e, µ), (e, s) & (µ, s) then the others (e, τ ), (µ, τ ) & (τ, s) are linear combinations as follows: since for i = 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 The identities are for k=(1, 2, • • • , n) and α, β = (e, µ, τ, s 1 , • • • ).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.59 is more general than eq.57, 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 basis 12 .

A. NSI eµ
The Hamiltonian is where in the usual notation b ≡ a ǫ eµ in this subsection.The determinant is given by where Det[(2E)H 3x3 ] is given by eq. 21.Then the characteristic equation can be calculated where Det[λI − (2E)H 3x3 ] is given by eq.22.
For the equivalent of the NHS identity, or CPV invariant, the dependence of H eµ on the matter potential means that ℑ(H eµ H µτ H τ 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 [38].This statement is true if the off-diagonal matter effect appears in H eτ or H µτ .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τ sector.

B. NSI µτ
For an NSI component in the µτ sector, the Hamiltonian is given by where in the usual notation b ≡ aǫ µτ in this subsection.The determinant and adjugate of the Hamiltonian are slightly different as product terms between "a" and "b" appear: where Det[(2E)H 3x3 ] is given by eq. 21. Similarly Both Det[λI − (2E)H] and Adj[λI − (2E)H] can be calculated using the m 2 i → (m 2 i − λ) procedure as before.
Using similar arguments as in (e, µ) 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 ν e ↔ ν µ , ν e ↔ ν τ and ν e disappearance but not for ν µ ↔ ν τ or ν µ and ν τ disappearance.
If the NSI components appear in more than one off-diagonal component of the Hamiltonian the calculation of Det[λI − (2E)H] and Adj[λ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 Adj[λ i I − (2E)H] in the numerator of eq. 5 but to use Π j =i (H − λ 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, ∆m 2  32 > or < 0, determine whether the neutrino mass state with least ν e is dominated by ν µ or ν τ , octant of θ 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.30 and the KTY identity eq.33 as well as identities related to the disappearance probabilities eq.29 and the CPC part of the appearance probabilities, eq.28.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.
where we have used Adj Alternatively, eq.B1 can be rewritten as follows: Here, one needs to consider the limit λ → λ i as the matrix (λ i I − H) is not invertible as Det[λ i I − H] = 0, although the Adj[λ i I − H] is well defined.Now consider the case when For any matrix, X, the following identity holds: Thus by setting X = λ i I − H, so that RHS is zero, demonstrates that giving eq.B1.
For an nxn matrix this continues until The next iteration returns zero for both A n+1 and d n+1 , as follows Tr(HA n+1 ) = 0 thus terminating the series. 13his recursion method is simple to program for numerical studies.Also eq.C9 gives an alternative way to calculate the 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: where Det m (H) and 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.
If λ i is one of the eigenvalues of H, (λ 1 , • • • , λ n ), then one can show by induction that This expression can also be derived from the inverse of the Vandermonde matrix 14 .
Note, the LeVerrier-Faddeev recursion method gives both Adj(λI − H) and Det(λI − H) for arbitrary λ, where as Adj(λ i I − H) = Π j =i (H − λ j I) requires all of {λ 1 , • • • , λ 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 Adj(λI − H) and Det(λI − H) and thus require very few lines of code to program.
As a cross check consider the case that H = V ΛV † , then as required.
As an observation note that the numerator and denominator of eq. 12 are identical polynomials of H and λ i respectively.
As this is required for the special case H = Λ, as the numerator and denominator must cancel.
Another informative cross check is given by the Caley-Hamilton identity, which can be written as Π j (H − λ j I) = 0, therefore 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.F3, set A = i m 2 i , B = i>j m 2 i m 2 j and C = Π i m 2 i , then to recover the eigenvalues, (m 2 1 , m 2 2 , m 2 3 ), from eq.F2 is a highly non-trivial exercise.
is an invariant.Using Only the δ jk contributes to the imaginary part as the rest lead to real terms.With the notation J αβ ij (V ) ≡ ℑ(V αi V * βi V * αj V βj ), we obtain for some fixed k, usually k=1.The last line is derived using the properties of J αβ ij and a similar procedure as in appendix A. No terms proportional to J αβ ik or J αβ kj appear.Therefore, 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 αβ ik or J αβ kj .

4 )
, "a" is the N e matter potential with a = 2 √ 2G F N e E ν and "b" is the N n matter potential with b = √ 2G F N n E ν where N e and N n are the number density for electrons and neutrons respectively 10 .