Compact Perturbative Expressions for Oscillations with Sterile Neutrinos in Matter

We extend a simple and compact method for calculating the three flavor neutrino oscillation probabilities in uniform matter density to schemes with sterile neutrinos, with favorable features inherited. The only constraint of the extended method is that the scale of the matter potential is not significantly larger than the atmospheric $\Delta m^2$, which is satisfied by all the running and proposed accelerator oscillation experiments. Degeneracy of the zeroth order eigensystem around solar and atmospheric resonances are resolved. Corrections to the zeroth order results are restricted to no larger than the ratio of the solar to the atmospheric $\Delta m^2$. The zeroth order expressions are exact in vacuum because all the higher order corrections vanish when the matter potential is set zero. Also because all the corrections are continuous functions of matter potential, the zeroth order precision is much better than $\Delta m^2_\odot/\Delta m^2_\text{atm}$ for weak matter effect. Numerical tests are presented to verify the theoretical predictions of the exceptional features. Moreover, possible applications of the method in experiments to check the existence of sterile neutrinos are discussed.


I. INTRODUCTION
Since the discovery of neutrino oscillations, [1], which determined that neutrinos are massive particles, many studies of neutrino scenarios beyond the three-flavor Standard Model have been performed.
One promising solution to the origin of the neutrino masses is a theoretical scheme with additional sterile neutrinos. In such a scheme neutrino oscillations will be modified because of the additional mixing with sterile neutrinos. In matter, calculations of neutrino propagation will be significantly more complicated since the sterile neutrinos also change the Wolfenstein matter effect term [2] in the Hamiltonian. There have been some analytical derivations of the matter effect in a 3+1 scenario, i.e. one sterile neutrino [3] in addition to the 3 active ones. However, the exact analytical solutions are impossible for more than one sterile neutrino because a quintic or even higher order equation will be encountered. Consequently, alternative perturbation approaches should be considered.
A satisfying perturbative framework, regardless of the existence of sterile neutrinos, is expected to possess the following properties: the expansion parameter is small; crossings of zeroth order eigenvalues are avoided anywhere; the approximated values go to the exact ones in vacuum. Recently, a compact perturbative framework achieving all the objectives above was developed by Denton, Minakata and Parke (DMP) to calculate the propagation of neutrinos in matter under the assumption of the standard three-flavor scheme [4][5][6]. * parke@fnal.gov; 0000-0003-2028-6782 † xining@uchicago.edu; 0000-0001-8959-8405 The main focus in this paper is to extend the principle and method of the DMP framework to schemes with sterile neutrinos when the scale of matter potential a is smaller or comparable to ∆m 2 atm , which is the case of all running and proposed accelerator neutrino oscillation experiments. The expansion parameter [5,7,8], which will be retained by the extension, is ≡ ∆m 2 21 /∆m 2 ee 0.03, ∆m 2 ee ≡ cos 2 θ 12 ∆m 2 31 + sin 2 θ 12 ∆m 2 32 .
The perturbative Hamiltonian will have no diagonal elements, and all its off-diagonal elements are proportional to and vanish in vacuum. Crossings of the zeroth order active eigenvalues will be resolved by a series of real or complex rotations; whereas crossings of the large sterile eigenvalues will not be considered since this will happen only if the matter effect is extremely large.
The structure of this paper is listed as follows: in Section II, we derive details of the rotations. This gives the zeroth order PMNS matrix and eigenvalues. The perturbative Hamiltonian is also determined by the rotations. In Section III, we discuss the higher order corrections by perturbative expansions after the rotations. A numerical test will also be presented to verify the predicted precision. In Section IV we use these perturbative expressions to calculate the oscillation probabilities of different channels and baselines. Moreover, potential applications of the method are discussed. We compare our method to some former works in Section V. Section VI is the conclusion. All other remarks and supplementary materials that are useful can be found in the Appendices.

II. ROTATIONS TO DERIVE ZEROTH ORDER APPROXIMATIONS AND PERTURBATIVE HAMILTONIAN
The principle of the method in [4][5][6] is that by implementing a series of rotations of the Hamiltonian, one can disentangle the crossings of the diagonal elements and diminish the off-diagonal elements to arbitrary scales. In particular: 1. In the given Hamiltonian in flavor basis, find the sector with leading order (largest absolute value) off-diagonal element, then perform a rotation to diagonalize this sector.
2. Use the rotated Hamiltonian to replace the initial one and repeat the process until all the off-diagonal elements are smaller than the expected scale and the diagonal element crossings are eliminated.
In principle, the above process is not designated to any specific dynamical system and is also applicable to the schemes with sterile neutrinos. However, this scheme must be implemented with considerable care otherwise the resulting analytic expressions becoming extremely long and complicated. First, one has to carefully choose the extension to the PMNS matrix to include sterile neutrinos as the standard choice here is far from optimal. Second, one has to decide whether or not one deals with all level crossings of the diagonal elements of the Hamiltonian or restrict the range of applicability of the result. We address these issues in depth in the following sub-sections.

A. PMNS matrix in vacuum
If we assume a 3 + N scheme, i.e. there are N sterile neutrinos in the scheme, the Hamiltonian in the flavor basis will be where a and b are Wolfenstein's matter potentials [2]: For earth matter the neutron number density N n is approximately equal to the electron number density N e so that b ≈ a/2. The PMNS matrix U PMNS in vacuum, which relates the flavor basis and the mass basis, is the product of a series of (complex) rotations [9,10]. In the Standard Model, the convention is chosen to be U SM PMNS ≡ U 23 U 13 U 12 . In the 3 + N scheme there will be extra rotations mixing with sterile neutrinos. It is natural to require that the convention is equivalent to that of the 3νSM case if all the extra rotations are trivial. Therefore, we will keep the relative positions of the three rotation matrices in the active sector when defining the PMNS matrix with sterile neutrinos.
Also it is observed that both the second and the third row vanish in the matter potential term in Eq. 2, thus we will keep U 23 as the first rotation in the PMNS matrix so the R.H.S of Eq. 2 will be independent of the 2-3 mixing parameters if we perform the U 23 rotation. The last step to determine the convention of the PMNS matrix is finding places after the U 23 for the rotations mixing with the sterile neutrinos. By trying different choices to simplify the calculation processes, we adopt the following convention of the PMNS matrix: (4) where U sterile is the product of all the rotations mixing with sterile neutrinos. This choice leads to significant reductions in the complexity of the calculations and the resulting expressions. Physics, of course, is independent of this choice.
In the following sections, we will use the 3+1 scheme as an example to develop the expressions for the schemes with sterile neutrinos. In particular, we choose 1 Current global fits [11][12][13] suggest |U i4 | ∼ 0.1, so in this paper we assume that U sterile 1 + O( √ ), which means that s i4 ∼ O( √ ) for i = 1, 2, 3. The small parameter is defined in Eq. 1.
The convention in Eq. 4 is different from the usual one used by many papers in which U sterile comes before (i.e. on the left side of) all the three rotations in the active sector (see e.g., [14]). We will derive the relations of the mixing angles and phases connecting both conventions in Appendix A.

B. U23 and U sterile rotations
We first define a rotated basis |ν by |ν f is the flavor basis. After the rotations, the Hamiltonian becomes In the above equation M 2 (b) ≡ ∆m 2 41 + b c 2 14 c 2 24 c 2 34 ,H is a 3 × 3 submatrix in the active sector and inH M all the elements not in the 4th column or row vanish.
Based on the scales we can distribute the elements of H into two parts, i.e.H =H 0 +H 1 .
The leading order term is where and the diagonal elements, which can be approximations to the eigenvalues are In the first order termH 1 , all the diagonal elements vanish, and the off-diagonal elements are (H 1 ) 12 = 2E c 12 s 12 c 13 ∆m 2 ee + b k 12 c 24 c 2 34 e −iδ34 , None-zero elements ofH M are listed below (the Hamiltonian is a Hermitian matrix) Since s i4 ∼ O( √ ), it is easy to see thatH M ∼ O( √ ). AlthoughH M is not as small as O( ), it will be a part of the perturbative Hamiltonian. However, this does not mean that the first order corrections must be as large  Table I. as O( √ ). The mass of the heavy sterile neutrino will be an alternative parameter which controls scales of the correction terms. More specifically, in a perturbative expression, all non-zero elements ofH M are divided by M 2 . For large M 2 the quotient gives a small term in the perturbation expansion. Another condition that is necessary forH M being a perturbative Hamiltonian is that it consists of terms proportional to a and b, which means that it vanishes in vacuum. This is crucial because we require the perturbative expressions to be exact in vacuum.

C. U13 rotation
Now the dominating off-diagonal term (except the ones inH M ) comes from the (1-3) sector ofH 0 . Because of the complex phase δ 34 , the rotation will not be real. Let us assume that the rotation is U 13 (θ 13 , α 13 ), whereθ 13 is the real rotation angle and α 13 2 is the complex phase. After this rotation the neutrino basis becomes where Since the 4th index is not engaged in the rotation, we can just focus on the first three indices and define a 3 × 3 submatrix U 13 to be the active sectors of U 13 , i.e.
with λ ± and λ 0 to be determined. SimultaneouslyH 1 becomesĤ It can be shown that The real rotation angle and the complex phase can be determined by The elements ofĤ 1 are Here we are not using the usual phase symbol δ since α 13 is not an effective physical phase in matter. In Appendix B it can be eliminated by implementing a pure phase transformation of the neutrino basis.
+b k 12 c 24 c 2 34s13 e i α13 + k 23 c 34c13 e −iδ34 e iδ24 , The Hamiltonian in the sterile sector becomeŝ At the end of this subsection we define a real parameter and a phase α Obviously ∼ and (Ĥ 1 ) 23 = e iα ∆m 2 ee /2E. It is not hard to see that in the Standard Model = | sin (θ 13 − θ 13 ) s 12 c 12 |, which reconciles with the one defined in [5]. The two new defined parameters will frequently emerge in the following sections. Since in vacuum, a, b = 0,θ 13 = θ 13 and α 13 = 0, must be zero then, as shown in Fig. 1. This guarantees that the perturbative expressions will be exact in vacuum.  Table I.

D. U12 rotation
As pointed out in [5], to resolve the λ 1 and λ 0 crossing at the solar resonance, one more rotation that diagonalizes the (1-2) sector is necessary. Again, since (Ĥ 1 ) 12 is complex, the rotation cannot be real in general. We assume that the rotation in (1-2) sector is U 12 (θ 12 , α 12 ), and after this rotation, the neutrino basis becomes where Similar to the case of the (1-3) rotation, we can again define a 3 × 3 submatrix U 12 by Now we require the U 12 (θ 12 , α 12 ) to diagonalize the (1-2) sector ofĤ. After the rotation the sub-Hamiltonian iš whereȞ 0 andȞ 1 are in zeroth and first order respectively, i.e.
The diagonal elements ofȞ 0 are The real rotation angle and the complex phase can be determined by Values of sin 2θ 13 and sin 2θ 12 are plotted in Fig. 2. After this (1-2) rotation, crossings of the first two diagonal elements λ 1,2 have been resolved, as shown in the top panels of Fig. 3. They will be the zeroth order eigenvalues in the following perturbation expansions in the next section. The difference between 3+1 and 3νSM is small in both panels of Fig. 2 and the bottom panels of Fig. 3 but not insignificant.
The Hamiltonian in the sterile sector now iš FromH M toȞ M , we implemented two rotations in (1)(2)(3) and (1-2) sectors. Because the active and sterile sectors were not mixed by the two rotations, the elements are still combinations of the terms proportional to s i4 ∼ O( √ ). Elements ofȞ M can be found in Appendix C.

E. Crossings of M 2
In principle, there are still some possible crossings of the diagonal elements, namely the crossings to the fourth diagonal element. Since both the (1-3) and the (1)(2) rotations are in the active space (first three rows and columns), the fourth element is still since ∆m 2 41 is much larger than the active eigenvalues in vacuum. Thus, the crossings to M 2 can only happen with very high neutrino energy, as shown in the top panels of Fig. 3. From the figure we can see that if Y e ρ = 1.4 g · cm −3 , for the earth's crust, the neutrino energy must be O(1) TeV. Considering the energy scales of the current and future accelerator based oscillation experiments, we are therefore not considering the energy region of these additional crossings, so they will not effect our result. For much higher energy experiments these additional level crossings would have to be dealt with using matter additional rotations.

F. Summary of the rotations
NowȞ 0 's diagonal elements, λ 1,2,3 , do not cross (crossings to M 2 will not happen in the energy region of interest). All the off-diagonal elements in the active sectors are of scale . We will distribute all the diagonal elements to the zeroth order Hamiltonian and all the offdiagonal elements to the perturbative Hamiltonian, i.e.
Since all possible degeneracies have been removed in the energy scale which we are interested in, we are free to implement a perturbation expansion to achieve even better accuracy. The process of reducing errors by performing rotations and perturbative expansions is summarized in Fig. 4.
For the scenario with more than one sterile neutrino, although it is more complicated, the rotation method developed here is still applicable. Using the same convention of Eq. 4 to define the PMNS matrix and implement the rotations in the sequence of U 23 → U sterile → U 13 → U 12 as above.

III. PERTURBATIVE EXPRESSIONS
Since all the crossings of the zeroth order eigenvalues have been resolved (except for the crossings with M 2 , which are not in the energy region of interest) by the rotations and all the off-diagonal elements are small, we can now calculate the higher order corrections to the eigenvalues and eigenvectors by perturbation methods .
We define V to be the exact PMNS matrix in matter. It can be related to the zeroth order U m PMNS by where W n is nth order correction. The exact eigenvalues are where λ 1,2,3 are defined in Eq. 30 and λ 4 = M 2 , λ (n) i is the nth order correction.
First order corrections to the eigenvalues are First order corrections to the eigenstates are determined by W i defined in Eq. 36, which are The detailed first and second order formulas of the perturbation expansions can be found in Appendix D. In general, with crossings of the zeroth order eigenvalues ruled out, perturbative expansions can go to arbitrary precision. However, numerical tests will suggest that it is sufficient to terminate the approach at second order.

A. Numerical precision test
We now test the accuracy of our perturbative expressions. We choose the ν µ → ν e channel and 1300km baseline of DUNE to do the numerical test. The density of the earth crust is chosen to be Y e ρ = 1.4 g·cm −3 , b = a/2 and all the mixing parameters are listed in Table I. The exact oscillation probabilities can be figured out by [3] or given by a computer algebra system 3 . The results are shown in Fig. 5. The error in the zeroth order expression is expected to be no more than ∼ 10 −2 , which is confirmed by the red curve in the plot; the green curve depicts the error of the first order perturbative expansion, which is under 2 ∼ 10 −4 ; to second order, the error further declines to 3 ∼ 10 −6 , which also coincides with the prediction. In Fig. 5 the expectation values are obtained by averaging over the fast oscillation terms, i.e. the terms with angular velocities proportional to (λ 4 − λ i ). More specifically, 3 Only considering the 3+1 scheme, an analytical solution is still possible since one just need solve a quartic equation; but it is not the case for schemes with more sterile neutrinos  Table I.
Based on the numerical results, we confirm that at least the second order perturabtive expansion is significantly more accurate than any experimental results [15][16][17][18][19].

IV. OSCILLATION PROBABILITIES AND DETECTING STERILE NEUTRINOS
In this section we will discuss a possible application of the perturbative expressions above for detecting sterile neutrinos. The principle of the approach is that one can calculate the theoretical predictions of the oscillation probabilities in different schemes and compare them with the experimental results. Usually for a given baseline and neutrino energy the predictions from different schemes are close, therefore it is essential to figure out sufficiently accurate expressions for the oscillation probabilities. A similar discussion can be found in [28].
In a scheme with N sterile neutrinos, the neutrino oscillation probabilities for ν α → ν β (α, β ∈ {e, µ, τ }) are where λ (ex) i are exact eigenvalues. We can chose the zeroth order results as an approximation, i.e. we adopt where U m PMNS is defined in Eq. 35 and where λ 1,2,3 are defined in Eq. 30 and λ 4 = M 2 (b). For the mass of the sterile neutrino, since it is significantly larger than the active ones, the oscillations related to it will be averaged out. Former and running experimental facilities have provided parameter fitting results of neutrino oscillations for different schemes. With these parameters, for future baselines, one can predict the probabilities in different schemes and this is a potential approach to determine the existence of sterile neutrinos [28]. We present the probabilities given by the 3+1 scheme and the differences of the probabilities | P 3+1 − P 3νSM |, in different channels, in Fig. 6, 7 and 8. The probabilities in the Standard Model are given by [29,30]; the 3+1 scheme is calculated by the 0th order rotation method developed in this paper. All the parameters are given in Table I.  The exact probability (expectation value) in the 3+1 scheme, which is plotted by the gray solid (black solid) curve, can be calculated by [3]. As a contrast, the dashed black line is showing the probabilities in the Standard Model, with Yeρ = 1.4 g · cm −3 .
In the figures we can identify several regions in which the differences are significantly larger than errors of the perturbation expansions. For example, in the ν µ → ν e channel, around the band of L/E 700 (km/GeV), | P 3+1 − P 3νSM | may be larger than 0.02, the differences will be even larger than 0.05 if L/E 1500 km/GeV and the baseline is longer than 500 km. In this channel baselines of T2K/HyperK, NOVA and DUNE (esti- predicted by the 3+1 scheme; differences of the probabilities predicted by the standard three-flavor scheme and the 3+1 scheme are presented in the right plot. P3+1 in both figures are computed by the 0th order rotation method developed in this paper. Parameters used are given in Table I. See Fig. 6 for neutrino flux energies of the listed facilities. mated) are marked [20][21][22]. For the channel of ν µ → ν µ , shifts from the 3νSM will be more than 0.05 with L/E 1000(km/GeV) and the baseline is longer than 1000 km. For the ν µ → ν τ channel, the scale of the greatest difference is larger than 0.16 if L/E 500(km/GeV) or 1500(km/GeV). Future experiments may measure the oscillation probabilities with baselines and neutrino energies in the region of interest predicted above and compare the results with the numerical outcomes.

V. COMPARE TO EXISTING APPROXIMATION FORMULAS
Approximation methods to calculate neutrino oscillations in matter in 3+1 scheme have been studied by many researchers, for example see [28,[31][32][33][34]. All these works chose to ignore the Hamiltonian's fourth row/column (except for the fourth diagonal element) in zeroth order approximations; thus the problem was reduced to 3×3 case. However, solving a three dimensional eigensystem is still not simple (many approximation methods developed for the 3νSM scheme). Fong, Minakata and Numokawa (FMN) adopted the exact 3-dim solutions [28] which was complicated, see [29,35]. Klop and Palazzo (KP) used one more approximation method for the 3 × 3 submatrix [31] which introduced extra errors. Fig. 9 compares fractional errors of the rotations method developed in this paper (PZ) with 0th order FMN and KP, assuming baselines of T2K/Hyper-K (T2K/HK), NOVA, HyperK-Korea (T2HKK) and DUNE. Compared with KP, just to 0th order PZ is significantly more precise for almost all baselines and energy ranges. Based on Fig. 9, FMN's precision is similar to PZ's for most baselines and energy regions, however, we can still identify PZ's advantage for T2K's baseline or low energy ( 1GeV) regions. We compare computing time of the different methods in Table II, since any specific computing time heavily depends on performance of the computer, we provide a list of relative computing time, i.e 0th order PZ's computing time is set to one unit time. The speed of a numerical method (using Eigen 3.3.7, http://eigen.tuxfamily.org) is also included in the comparison. A similar comparison of the methods for the 3νSM scheme can be found in [36]. Table II shows that compared with the rotation method developed in this paper (PZ), only the KP method is faster, however, its advantage in speed will be offset by the relatively poor precision. The FMN method is comparable in time consumed to the exact analytical solution. For experimentalists, the speed of a numerical method for evaluating the oscillation probability is relevant because it effects the time and computing resources consumed by large multi-dimensional parameter scans.
Besides simplicity and better precision, the rotation method of this paper also gives explicit expressions of zeroth order eigenvalues and mixing angles and CP phases with high precision which are not covered by any former references.

VI. CONCLUSIONS
A compact and simple technique for calculating neutrino oscillation probabilities in matter for schemes with sterile neutrino has been developed from the extension of an analogous method for the 3νSM model [5]. The extended method is appropriate to conditions in which the Wolfenstein matter potentials defined in Eq. 3 are not significantly larger than ∆m 2 atm , meaning that it may be applied to all the current and proposed accelerator neutrino oscillation experiments. The zeroth order eigensystem of the Hamiltonian in the active space (i.e. the three dimensions included in 3νSM) derived by the method is non-degenerate. Meanwhile, numerical study shows that crossings of the zeroth order eigenvalues involving the  [28]; the red curve (KP) is from Klop and Palazzo [31]. Parameters used for this sample calculation are listed in Table I. Relative speed of the methods can be found in Table II [31]; the analytical solution is given by [3]; the numerical method is referred to Eigen 3.3.7, http://eigen.tuxfamily.org.
sterile one only happen with large matter potential (high neutrino energy for the earth's crust), which is out of the paper's scope of discussion. An additional crucial advantage of the method developed in this paper is that errors of the zeroth order results are small when the matter potentials are ≤ ∆m 2 ee and vanish in vacuum because the matter potential terms are factors of all the perturbative terms.
We implement a series of complex or real rotations to kill the leading order off-diagonal elements and resolve crossings of the diagonal elements of the Hamiltonian. The rotation angles and phases are the zeroth order mixing parameters of the effective PMNS matrix in matter. In the rotated Hamiltonian, the diagonal elements are the zeroth order eigenvalues, whereas the off-diagonal elements are the perturbing Hamiltonian. Based on this arrangement, perturbation expansions are performed after the rotations to achieve better accuracy. When the matter effect is comparable to the vacuum mixing effect, i.e. matter potentials defined in Eq. 3 are com-parable to ∆m 2 31 , the expansion parameter is no larger than ∆m 2 21 /∆m 2 31 0.03; when the matter effect itself is weak, the perturbative Hamiltonian will be higher order because it consists of terms proportional to a or b .
Finally numerical tests show that the rotation method developed in this paper balances precision with computing speed nicely; comparisons with a numerical method and some previous approximate methods have been applied to demonstrate the rotation method's advantages as presented in this paper. For the first order perturbation expansion of this paper, absolute errors of the oscillation probabilities are shown to be no more than 10 −4 . This precision is sufficient to distinguish the schemes with sterile neutrinos from the 3νSM model, which makes the method developed by this paper suitable to explore the existence of sterile neutrinos.
convention of the complex phases we need to implement a phase transformation. Firstly, we multiply the 1st row by e −i α12 and the 1st column by e i α12 ; then the 3rd row is multiplied by e i(α13−α12) and the 3rd column is multiplied by e i(−α13+α12) . Finally all the complex phases are absorbed into U 23 , U 24 and U 34 . The zeroth order phases areδ 12 = 0, δ 13 = 0, δ 23 = δ 23 − α 13 + α 12 , δ 24 = δ 24 + α 12 , Since all diagonal elements have been absorbed into the zeroth order Hamiltonian, by Eq. 38, the first order corrections to the eigenvalues, which are the diagonal elements of the perturbative Hamiltonian, are zero, i.e. λ (1) i = 2E (Ȟ 1 ) ii = 0. (D1) As for the eigenvectors, to first order, are derived from Eq. 39 as follows:

Second order corrections
The second order corrections to eigenvalues are or explicitly as where ∆λ ij ≡ λ i − λ j . The second corrections to eigenvectors can be calculated by the corrections to the PMNS matrix: (D5) We list the elements of W 2 below: (W 2 ) 11 = −( ∆m 2 ee ) 2s