Eigenvalues: the Rosetta Stone for Neutrino Oscillations in Matter

We present a new method of exactly calculating neutrino oscillation probabilities in matter. We leverage the"eigenvector-eigenvalue identity"to show that, given the eigenvalues, all mixing angles in matter follow surprisingly simply. The CP violating phase in matter can then be determined from the Toshev identity. Then, to avoid the cumbersome expressions for the exact eigenvalues, we have applied previously derived perturbative, approximate eigenvalues to this scheme and discovered them to be even more precise than previously realized. We also find that these eigenvalues converge at a rate of five orders of magnitude per perturbative order which is the square of the previously realized expectation. Finally, we provide an updated speed versus accuracy plot for oscillation probabilities in matter, to include the methods of this paper.


I. INTRODUCTION
While the majority of the parameters in the threeneutrino oscillation picture have been measured, measurements of the remaining parameters will come by leveraging the matter effect in long-baseline experiments such as the currently running T2K and NOvA experiments, the now funded and under construction T2HK and DUNE experiments and the proposed T2HKK and ESSnuSB experiments, [1][2][3][4][5][6]. In this context, only a full three-flavor picture including matter effects is adequate to probe the remaining parameters. Given the time and effort that is going into these experiments, it is paramount that we understand neutrino oscillations in matter as best as we can, both analytically and numerically, so as to maximize the oscillation physics output from these major experiments.
The matter effect is the fact that while neutrino propagation in vacuum occurs in the mass basis, in matter since the electron neutrino experiences an additional potential, they propagate in a new basis. This effect was first identified in 1978 by Wolfenstein [7]. Exact analytic solutions for neutrino oscillation probabilities in constant matter densities are difficult to fully enumerate; a solution using Lagrange's formula appeared in 1980 [8] 1 while the full solution was first written down for three flavors in 1988 by Zaglauer and Schwarzer (ZS) [10]. The exact solution requires solving a cubic equation which, in the general case, has the unsightly and impenetrable cos( 1 3 cos −1 (· · · )) term present in the eigenvalues which are then in nearly every expression involving neutrino oscillations in matter 2 . Given the eigenvalues, there are then several choices of how to map this onto the oscillation probabilities. ZS mapped the eigenvalues onto the effective mixing angles and CP phase in matter; given the phase and the angles, it is then possible to write down the oscillation probabilities in matter using the vacuum expressions and the new phase, angles, and eigenvalues. In 2002, Kimura, Takamura, and Yokomakura (KTY) presented a new mapping from the eigenvalues onto the oscillation probabilities by looking at the products of the lepton mixing matrix that actually appear in the probabilities [12,13], see also [14][15][16][17]. Another formulation of the exact result in the context of the time evolution operator is ref. [18]. Along with these exact expressions, numerous approximate expressions have appeared in the literature in various attempts to avoid the cos( 1 3 cos −1 (· · · )) term, for a recent review see ref. [19].
In this article we use the eigenvector-eigenvalue identity, that has been recently and extensively surveyed in [20,21], to write the exact expressions for the mixing angles in matter in terms of the eigenvalues of the Hamiltonian and its principal minors. The benefit of this approach is two-fold. First, it makes the expressions for the mixing angles in matter, clearer, symmetric, and very simple. Since our approach for the oscillation probabilities in matter is based on the form of the vacuum expressions, the intuition that exists for the vacuum still applies in matter. Second, it allows for a simple replacement of the complicated exact eigenvalues with far simpler approximate eigenvalues in a straightforward fashion. We find that since this approximate approach only relies on approximate expressions for the eigenvalues, it is more accurate than previous methods, including Denton, Minakata and Parke (DMP) [22], with a comparable level of simplicity. We also explore the convergence rate of the eigenvalues in DMP and find that since all odd-order corrections vanish, they converge much faster than expected, at a rate of ∼ 10 −5 .

II. AN EIGENVALUE BASED EXACT SOLUTION
The technique of ZS is to determine expressions for the mixing angles and CP-violating phase in matter ( θ 23 , θ 13 , θ 12 , and δ) as a function of the eigenvalues and other expressions, while KTY derives the general expression for the product of elements of the lepton mixing matrix, U αi U * βj . In this section, we describe a technique of using both approaches.
First, we note that, given the eigenvalues, the mixing angles can be determined from various |U αi | 2 terms. This employs a simpler version of the main result of KTY. Then, to address the CP-violation part of the oscillation probabilities, we use the Toshev identity [23].

A. Mixing Angles in Matter
The neutrino oscillation Hamiltonian in matter in the flavor basis is where we have subtracted out an overall 2G F n e E is the Wolfenstein matter potential [7], and the PMNS lepton mixing matrix [24,25] is parameterized, where s ij = sin θ ij , c ij = cos θ ij , and we have shifted the CP-violating phase δ from its usual position on s 13 to s 23 which does not affect any observable. For our numerical studies we use ∆m 2 21 = 7.55 × 10 −5 eV 2 , ∆m 2 31 = 2.5 × 10 −3 eV 2 , s 2 12 = 0.32, s 2 13 = 0.0216, s 2 23 = 0.547, and δ = 1.32π from [26].
Using the eigenvector-eigenvalue identity [20], the squares of the elements of the lepton mixing matrix in matter are simple functions of the eigenvalues of the neutrino oscillation Hamiltonian in matter, λ i /2E for i ∈ {1, 2, 3}, and new submatrix eigenvalues, ξ α /2E and χ α /2E for α ∈ {e, µ, τ }. In general, the square of the elements of the mixing matrix are parameterization independent, where i, j, and k are all different, and the λ i are the exact eigenvalues, see appendix A. This result, eq. 3, can also be directly obtained from KTY as shown in appendix B. This equation is valid for every element of the mixing matrix, even the µ and τ rows, which are relatively complicated in the standard parameterization. Eq. 3 is one of the primary results of our paper. Given the eigenvalues of the Hamiltonian and the eigenvalues of the submatrix Hamiltonian, it is possible to write down all nine elements of the mixing matrix in matter, squared. This result is also quite simple and easy to memorize which is contrasted with the complicated forms from previous solutions [10,13,27].
The submatrix eigenvalues ξ α /2E and χ α /2E are the eigenvalues of the 2 × 2 submatrix of the Hamiltonian, for α, β, and γ all different. Explicit expressions for the Hamiltonian are given in appendix C and the eigenvalues of the submatrices, which require only the solution to a quadratic, are plotted in fig. 1. We note that while solving a quadratic is necessary to evaluate the submatrix eigenvalues, since only the sum and the product of the eigenvalues (that is, the trace and the determinant of the submatrix Hamiltonian) appear in eq. 3 whose numerator can be rewritten as λ 2 i −λ i (ξ α +χ α )+ξ α χ α , the submatrix eigenvalues do not have to be explicitly calculated. The expressions for the sums and products of the eigenvalues are given in appendix C.
Given the standard parameterization of the lepton mixing matrix, this allows us to write all three mixing angles in matter as simple expressions of the eigenvalues, where the hat indicates that it is the mixing angle in matter. While similar versions of eqs. 5 and 6 have previously appeared in the literature [10], eq. 7 is original to this manuscript. The general form of eq. 3 allows us to write any term in the lepton mixing matrix, and thus any mixing angle with considerable ease. In addition, as The two submatrix eigenvalues, ξα and χα, as a function of neutrino energy, are shown in the solid blue curves with α = e, µ, τ in the left, center, and right figures respectively. For comparison, the full matrix eigenvalues λi are shown in dashed red, green, and orange in each panel. When a submatrix eigenvalue (solid) overlaps one of the full matrix eigenvalues (dashed) the corresponding |Uαi| 2 → 0, as seen from the numerator of eq. 3. Note the Cauchy interlacing condition is satisfied, λ1 ≤ ξα ≤ λ2 ≤ χα ≤ λ3, for each α = (e, µ, τ ) and all E, using the convention ξα < χα. See appendix D for further discussion.
we will show in the next section, this also allows us to calculate the CP violating phase quite easily.
In appendix E we show how to use this method in the vacuum (θ 23 , δ)-rotated flavor basis. Further extensions of eq. 3 to an arbitrary number of neutrinos is also discussed in appendix F.

B. CP-Violating Phase in Matter
In order to determine the CP-violating phase in matter, we note that cos δ can be determined, given the other mixing angles in matter, from |U µ1 | 2 (or |U µ2 | 2 , |U τ 1 | 2 , or |U τ 2 | 2 ) from eq. 3. The sign of δ needs to be separately determined. We note that the sign of δ must be the same as the sign of δ. To see this, we employ the Naumov-Harrison-Scott (NHS) identity [28,29], where J = [U e1 U µ2 U * e2 U * µ1 ] = s 23 c 23 s 13 c 2 13 s 12 c 12 sin δ is the Jarlskog invariant [30]. We note that the numerator and denominator in eq. 8 always have the same sign, so sin δ has the same sign as sin δ. That is, the eigenvalues in matter never cross.
In practice, it is simpler to determine sin δ from the Toshev identity [23], and use θ 23 determined in eq. 7. An alternative method for determining the CP violating phase is given in appendix E.

C. Oscillation Probabilities in Matter
Finally, these can all be combined into any oscillation probability. For the primary appearance channel at NOVA, T2K, DUNE, T2HK(K), ESSnuSB [1][2][3][4][5][6], or any other long-baseline neutrino experiment, the oscillation probability can be written in the following compact form, where the upper (lower) sign is for neutrinos (antineutrinos) and The above determination of the mixing angles and CP violating phase in matter allow for the simple determination of the oscillation probability in matter for the ν µ → ν e appearance channel, or any other channel via the vacuum oscillation probabilities. Therefore, any physics intuition already obtained for vacuum oscillation probabilities is easily transferred to oscillation probabilities in matter.
In addition to all the appearance channels this approach also works in a straightforward fashion for the disappearance channels as well. For disappearance the oscillation probabilities in matter can be written as, Thus the coefficients, | U αi | 2 | U αj | 2 , can be read off as simple functions of the eigenvalues and the submatrix eigenvalues, eq. 3, without any need to even convert to the mixing angles in matter.

III. APPROXIMATE EIGENVALUES
While the form of the mixing angles in matter presented above is exact, it still relies on the complicated expression of the eigenvalues. It has been previously shown, however, that the eigenvalues can be extremely well approximated via a mechanism of changing bases as demonstrated by Denton, Minakata, and Parke (DMP) [22], see also refs. [31][32][33]. While expressions for the differences of eigenvalues in DMP are quite compact [34], the expressions in eqs. 5-7 require the individual eigenvalues so we list those here as well. Beyond the zeroth order expressions, it is possible to derive higher order terms through perturbation theory [22] or through further rotations [33]. This approach leads to a smallness parameter that is no larger that c 12 s 12 ∆m 2 21 ∆m 2 31 ∼ 1.5% and is zero in vacuum confirming that the exact solution is restored at zeroth order in vacuum, see eq. 23 below.

A. Zeroth Order Eigenvalues
The zeroth-order eigenvalues are extremely precise with a fractional error in the difference of the eigenvalues of < 10 −5 at DUNE. We define eigenvalues of two intermediate steps. First, the eigenvalues of the un-rotated Hamiltonian, after a constant U 23 (θ 23 , δ) rotation 3 , λ a = a + s 2 13 ∆m 2 ee + s 2 12 ∆m 2 21 , λ c = c 2 13 ∆m 2 ee + s 2 12 ∆m 2 21 .
Next, after an O 13 rotation, we have and Finally, after an O 12 rotation, the eigenvalues through zeroth order are and Here x represents an approximate expressions for the quantity x in matter. At this order, θ 23 and δ, are unchanged from their vacuum values. Eqs. 18 to 22 define the zeroth order approximation. We note that φ and ψ are an excellent approximation for θ 13 and θ 12 respectively, [22], see ref. [33] for the explicit higher order correction terms. These are effective two-flavor approximations to θ 12 and θ 13 while eqs. 5 and 6 are the full three-flavor exact expressions. We further discuss the similarity in these expressions in sec. F.

B. Second Order Eigenvalues
After performing the rotations that lead to the eigenvalues in eqs. 20 and 21, the smallness parameter is ≡ sin(φ − θ 13 )s 12 c 12 ∆m 2 21 /∆m 2 ee < 1.5% (23) and is zero in vacuum since φ = θ 13 in vacuum. Because of the nature of the DMP approximation, the zeroth order eigenvalues in eqs. 20 and 21 already contain the first order in corrections. That is, the first order corrections are just the diagonal elements in the perturbing Hamiltonian which are all zero by construction. The second order corrections are simply, λ λ 2 .
It is useful to note that λ 1 and λ 2 are related by the 1-2 interchange symmetry [22]. The 1-2 interchange symmetry says that all oscillation observables in matter are independent of the following transformations, and c ψ s ψ → −c ψ s ψ . (27) It is clear that λ is invariant under this interchange, and λ (2) 2 follows directly from λ (2) 1 and the interchange. As with φ in eq. 19, ψ is an excellent approximation for θ 12 . The fractional precision of the eigenvalues at zeroth and second order are shown in fig. 2. Since φ and ψ are good approximations for θ 13 and θ 12 respectively, and since θ 23 and δ don't vary very much in matter, one could imagine using the vacuum probabilities with the approximate eigenvalues and replacing only θ 13 and θ 12 with φ and ψ respectively. This is exactly the DMP approach at zeroth order. Thus one way to quantify the improvement of this approach over DMP is to compare the precision with which we can approximate θ 13 and θ 12 with either φ and ψ which result from a two-flavor rotation (see eqs. 19 and 22) or with eqs. 3, 6, and 5. We have numerically verified that the full three-flavor approach to calculating the mixing angles improves the precision on the mixing angles in matter (and thus the oscillation probabilities) compared with the two-flavor approach that leads to φ and ψ.
Next, we note that for similar reasons that the first order corrections vanish, λ (1) i = 0, all the odd corrections vanish within the DMP framework 4 . That is, λ (k) i = 0 for all i ∈ {1, 2, 3} and any k odd, see appendix G. While it would appear that, given a perturbing Hamiltonian ∝ that the precision would converge as , this shows that, in fact, the precision converges considerably faster at 2 . This result had not been previously identified in the literature.
We now compare the precision of sine of the mixing angles and CP violating phase in matter using the approximate eigenvalues through zeroth order and second order to the exact expressions in fig. 3. Using the zeroth order eigenvalues to evaluate the angles and the phase is quite precise even at zeroth order, at the 1% level or much better. Adding in the second order corrections dramatically increases the precision by about four orders of magnitude for neutrinos and six orders of magnitude for anti-neutrinos, consistent with the fact that is ∼ 10 −2 in the limit as E → ∞ and ∼ 10 −3 in the limit as E → −∞. We also see that we recover the exact answers in vacuum, a trait that many approximation schemes do not share [19].
Next, we show the precision of the appearance oscillation probability for DUNE in fig. 4 5 . The scaling law of the precision remains the same as previously shown and we have verified that it continues at the same rate to even higher orders. In fact, as we continue to higher orders we find that all the odd corrections to the eigenvalues vanish, see appendix G.
Finally, in an effort to roughly quantify the "simplicity" of our results, we computed the speed with which we can calculate one oscillation probability as shown in fig. 5. For comparison we have included many other approximate and exact expressions as previously explored in [19]. For the sake of openness, the nu-pert-compare code used for each of these is publicly available https:// github.com/PeterDenton/Nu-Pert-Compare [35]. We note that while our new results using DMP eigenvalues are not as fast as others, adding in higher order corrections is extremely simple, as indicated in eqs. 24-26 which give rise to an impressive six orders of magnitude improvement in precision for almost no additional complexity. All points use δ = −0.4π except for OS and Exp where δ = 0, for a detailed discussion see ref. [19].
Computational speed is a useful metric not only for simplicity but also for long-baseline experiments which must compute oscillation probabilities many times when marginalizing over a large number of systematics and standard oscillation parameters. In addition, performing the Feldman-Cousins method of parameter estimation is known to be extremely computationally expensive [36].

IV. CONCLUSIONS
In this article we have used the eigenvector-eigenvalue identity to develop a new way to write neutrino oscillation probabilities in matter, both exactly and with simpler approximate expressions. The primary new result involves determining the mixing angles in matter which has the benefit in that intuition gained about vacuum oscillations still applies to oscillations in matter. The CP violating phase in matter is then determined in a straightforward fashion from θ 23 and the Toshev identity [23]. Given the mixing angles and the CP violating phase in matter, writing down the oscillation probability in matter follows directly from the simple vacuum expression.
This technique benefits from the simplicity of using the expression for the oscillation probability in vacuum, and the clear, compact, expressions for the mixing angles in matter, and can be applied to any oscillation channel. There is also the fact that by explicitly writing the mixing angles as simple functions of the eigenvalues, they can be replaced with simple approximate expressions, such as those derived by Denton, Minakata, and Parke (DMP) [22]. This new technique presented here is more precise than that in DMP, order by order, since this result is effectively complete to all orders in the eigenvectors and only requires correction to the eigenvalues. The primary new results of this article are, • Eq. 3, reproduced here, which presents a simple, clear and easy to remember way to determine the norm of the elements of the mixing matrix and hence the mixing angles in matter given the eigenvalues. Then the oscillation probabilities can be calculated in a straightforward fashion using the CP-violating phase in matter from the Toshev identity.
• The form of eq. 3 allows for the direct substitution of approximate eigenvalues, such as those from DMP. As shown here for the first time, the DMP eigenvalues converge extremely quickly, ∼ 10 −5 per step since all the odd order corrections to the eigen- We have plotted the fractional precision at the first oscillation maximum for DUNE at δ = −0.4π versus the time to compute one oscillation probability on a single core. Our results are labeled DPZ and are in orange. ZS [10] and Diag are two exact solutions where Diag represents an off-the-shelf linear algebra diagonalization package. Two other exact solutions from [18] are labeled OS (using the Cayley-Hamilton method) and Exp (exponentiating the Hamiltonian) which do not account for CP violation (δ = 0). We only plot expressions that reach at least 1% precision at the first oscillation maximum. The remaining expressions are: MP [32], AM [37], MF [38], AKT [31], and DMP [22]. For a detailed discussion see ref. [19].
• The form of eq. 3 is trivially generalizable to any number of neutrinos.
Given the formalism presented here, we have a clear and simple mechanism for calculating the oscillation probabilities in matter either exactly or approximately. This method hinges on the eigenvalues, therefore the eigenvalues are the Rosetta Stone for neutrino oscillations in matter.

ACKNOWLEDGMENTS
We want to thank Terence Tao for many useful and interesting discussions on the eigenvector-eigenvalue identity. PBD [8,10,39], the exact eigenvalues in matter are, 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, while S contains the cos( 1 3 cos −1 (· · · )) terms, where a ≡ 2E √ 2G F n e is the matter potential, E is the neutrino energy, G F is Fermi's constant, and n e is the electron number density.
As an example of the analytic impenetrability of S, setting a = 0 and recovering the vacuum values for the eigenvalues, (0, ∆m 2 21 , ∆m 2 31 ), is a highly non-trivial exercise.

Appendix B: Derivation From KTY
Since only the product of elements of the PMNS matrix are necessary to write down the oscillation probabilities, we start with the definition of the product of two elements of the lepton mixing matrix in matter from eq. 39 in KTY, ref. [13], where x is the quantity x evaluated in matter and the λ i 's are the exact eigenvalues in matter 6 , see appendix A. We note that a similar approach was used in [27]. The matrix p is just the Hamiltonian in matter, p αβ = (2E)H αβ = i λ i U αi U * βi , see appendix C. The other term, q, is given by We evaluate eq. B1 in the case of α = β.
where α, β, and γ are all different, as are i, j, and k. We can then write the numerator as (λ i − ξ α )(λ i − χ α ) where ξ α and χ α satisfy That is, ξ α and χ α are the eigenvalues of H α , the 2 × 2 submatrix of the Hamiltonian, We also note that eq. 3 leads to the following identity also presented in ref. [13], which is the NHS identity [28,29] divided by the Toshev identity [23]. That is, this quantity is independent of the matter potential. For the full PMNS matrix in matter, we combine U 23 (θ 23 , δ) with V 23 (α) into U 23 ( θ 23 , δ), as follows The equation for sin δ is exactly the Toshev identity [23]. The full PMNS matrix in matter is then by U P M N S = U 23 ( θ 23 , δ) U 13 ( θ 13 ) U 12 ( θ 12 ) .
where V ij = i|V |j . The rightmost term is clearly zero since V ii = 0. In the first summation if j = k, V jk = 0 for the same reason. If not then the numerator contains the product of all three off-diagonal terms. Since one of these terms is zero in DMP (V 12 ) we have that λ In fact, all λ (m) i = 0 for m odd. A brief proof is that when m is odd λ (m) i is a summation of terms proportional either to a diagonal element of V or the product V 12 V 23 V 31 .
The above conclusion is a special case of a more general statement. Let's consider an n × n perturbing Hamiltonian V n×n for which there is an r such that 0 ≤ r ≤ n wherein i.e. V n×n has the block form V n×n = 0 r×r V r×(n−r) V (n−r)×r 0 (n−r)×(n−r) .
If the above condition is satisfied, the mth order eigenvalue corrections λ (m) i = 0 for m odd. DMP is the case of n = 3, r = 2.
We list the even order corrections through 6 th order. We note that we only need to write down the λ For three neutrinos, the Jacobi method ensures that the conditions that the odd corrections to the eigenvalues can always be met by rotating one off-diagonal element of the perturbing Hamiltonian to zero. The necessary conditions can also be met for four neutrinos by rotating two off-diagonal elements in disconnected sectors (say, U 12 and U 34 ). For general matrices this condition cannot be met for more than four neutrinos.