Fibonacci Fast Convergence for Neutrino Oscillations in Matter

Understanding neutrino oscillations in matter requires a non-trivial diagonalization of the Hamiltonian. As the exact solution is very complicated, many approximation schemes have been pursued. Here we show that one scheme, systematically applying rotations to change to a better basis, converges exponentially fast wherein the rate of convergence follows the Fibonacci sequence. The results presented here are generally applicable to systems beyond neutrino oscillations as well.


Introduction
Measurements of neutrino oscillations have triggered an immense interest in gaining a better understanding of neutrino oscillations, specifically in the presence of the Wolfenstein matter effect [1]. Due to the complexity of the exact analytic solution for more than two flavors and the presence of the cos( 1 3 cos −1 (· · · )) term, the exact expressions are no more insightful than numerically diagonalizing the Hamiltonian directly. To address this, many approximate formulas have been developed, see ref. [2] for a 2019 review. A recent example is a rotation method known as the Jacobi method [3] which has been applied to neutrino oscillations to calculate the energy eigenvalues and eigenstates to high precision with simple structure and high calculation speed [2,[4][5][6][7][8][9][10][11][12][13][14]. The principle of this method is performing rotations to extinguish the largest off-diagonal elements and resolve the crossings of the diagonal elements of the effective Hamiltonian.
In this paper we further expand upon the properties of the rotation method. We report a phenomenon that in the context of three neutrino oscillations, precision of the approximation will be improved very rapidly with the number of rotations implemented. In details, the order of the error follows the Fibonacci series and grow exponentially. This feature provides a fast way of high precision calculation of neutrino oscillations in matter.
The structure of this paper is listed following. In section 2 we derive our main result. Section 3 gives a numerical test of the theory. Section 4 is the conclusion and summary of this paper's contents. Other materials we believe necessary can be found in the Appendix.

Preliminary rotations
The three neutrino problem in matter requires solving a fully populated complex 3×3 Hermitian matrix for its eigenvectors and eigenvalues. Given those (or just the eigenvalues, see [15][16][17][18][19]) determining the oscillation probabilities is straightforward. Because of the matter effect, the PMNS matrix with the vacuum parameters no longer diagonalizes the Hamiltonian and the eigenvalues are also altered. The Hamiltonian is typically split into a large part and a small part where the small part is proportional to the ratio of ∆m 2 's which is ∼ 3%. Applying perturbation theory at this point suffers from two problems: 1) the zeroth order eigenvalues (diagonal elements) cross at two matter potential values and 2) the perturbative terms (off-diagonal elements) are not small. Recently a rotation method has been applied to overcome the above defects [6,9]. The rotations applied can be used to address both issues: they can eliminate the largest off-diagonal elements in the perturbing Hamiltonian and they cause the level crossings to repel each other.
For an arbitrary n×n Hermitian matrix H, we choose two diagonal elements H pp and H qq and the two corresponding off-diagonal element H pq and H qp . The selected four elements form a 2 × 2 Hermitian submatrix h, It can be diagonalized by a single 2 × 2 complex rotation u = cos α e iβ sin α −e −iβ sin α cos α , where and the new diagonal elements are From Eq. 5 we see that the gap between λ u,v is at least 2|H pq | so any level crossings of the chosen diagonal elements will be resolved. Since the rotation will not increase the scales of any other elements, if the chosen H pq , the pivot, is the largest in absolute value among all the off-diagonal elements of the full matrix, the leading scale of perturbative terms (off-diagonal elements) has been reduced. The above process can be repeated. By implementing a series of rotations one can eliminate all the crossings of the eigenvalues and squeeze the off-diagonal elements as much as desired. This procedure, selecting the largest off-diagonal element, maximizes the precision of the entire matrix. If however, only certain elements of the matrix are necessary for a given calculation, different techniques may be more optimal. In the context of neutrino oscillations, our goal is to provide as unified of a framework as possible to apply equally to all channels.

Fibonacci recursive process
After the preliminary rotations to remove all level crossings (two rotations for the standard neutrino picture), the basis is such that a perturbative expansion of the Hamiltonian is possible [9]. However, in [10] an alternative method to perturbation theory has been discussed, i.e. using more rotations to further reduce the off-diagonal elements to enhance the precision without using perturbation theory. It has been shown that successive rotations will match the precision of perturbation theory. In this subsection we show that implementing additional rotations will be more efficient to achieve very high order precision than perturbative expressions, after a sufficient number of rotations.
For simplicity we focus on a 3 × 3 Hamiltonian, although our results are generally valid. After the preliminary rotations, the Hamiltonian is where H 0 is a diagonal zeroth order Hamiltonian and H 1 is the perturbative part where all diagonal elements vanish. That is, and where 0 < 1 is a small scale, a, b are some positive numbers, and |x|, |y| ∼ O(1). Please note that in this example (H 1 ) 12 is zero since for neutrino oscillations the final rotation to remove the level crossings is the (1-2) rotation. Depending on the initial Hamiltonian before the preliminary rotations, the pair of vanishing off-diagonal elements of H 1 may be different, it will not affect the following deviation.
Next, we assume that b ≥ a > 0, thus a x is the leading order off-diagonal element now so we should implement a rotation in the (1-3) sector. Substitute H 1 into Eq. 4 we get After this rotation the perturbative Hamiltonian in the new basis becomes y cos α 13 y * cos α 13   . For the sake of simplicity we define H 1 to have the order of a-b and H 1 have the order of b-(a + b) where in this definition the first number is smaller (corresponding to the order of the largest off diagonal term). It is easy to see that the rotation angle which can extinguish (H 1 ) 23 must be in order of O( b ) and after that the rotated perturbative Hamiltonian must have order of (a + b)-(a+b+b). This is the famous Fibonacci sequence; that is, the order of the size of the largest off-diagonal element is the sum of the order after each of the previous two rotations. The order of the smallness parameter in the perturbative part of the Hamiltonian will increases exponentially in the number of rotations since the Fibonacci sequence grows exponentially. This means that the diagonal part of the Hamiltonian will converge on the true expression very rapidly. Starting from H 1 ∝ , to achieve precision at the n level, the rotation method takes just O(log n) rotations.
A useful special case of H 1 is a = b = 1. We define H (N ) 1 to be the perturbative Hamiltonian after N rotations (not including the preliminary rotations). Then, we have that the size of the Hamiltonian shrinks exponentially as described by log H Moreover, we notice that all the perturbative Hamiltonian's diagonal elements are zero. Since the first order corrections to the eigenvalues are the diagonal elements of the perturbative Hamiltonain, therefore the order of errors of the eigenvalues will be double of the order of the perturbative Hamiltonian [18].
For a = b = 1, we compare orders of errors of the eigenvalues and eigenvectors given by perturbation expansions and the rotation method in fig. 1. The order of the size of the error in the eigenvalues (eigenvectors) grows like 2F n+1 (F n+1 ) where F n is the Fibonacci sequence defined as F 0 ≡ 0, F 1 ≡ 1, and F n = F n−1 + F n−2 for n > 1. Right: The order of the error of the eigenvectors scales like n + 1 using perturbation theory and F n+1 using rotations.
In Fig. 2 we display the ratio of the perturbative Hamiltonian after N rotations (not counting the prelimilnary rotations) to the zeroth order Hamiltonian. As expected, it is shown that the order spans of each rotation increase with the total number of rotations.

Conclusion
We show the potential application of a rotation method in high precision calculation of the neutrinos oscillations in matter in the context of neutrino oscillations in matter. The fast iteration steps of the method can enhance the zeroth order precision very rapidly. More specifically, a Fibonacci recursive process leads to an exponentially growth of the orders (of some small scale) of the zeroth order eigensystem's errors with number of the rotations. This feature grants an advantage of the rotation method in the range of high precision calculation compared with the perturbation expansion methods. In addition, while the complexity of each additional step in a perturbative expansion grows, the complexity of each additional rotation is constant and simple as shown in Eqs. 2-4.