Analytical Approximation of the Neutrino Oscillation Matter Effects at large $\theta_{13}$

We argue that the neutrino oscillation probabilities in matter are best understood by allowing the mixing angles and mass-squared differences in the standard parametrization to `run' with the matter effect parameter $a=2\sqrt{2}G_F N_e E$, where $N_e$ is the electron density in matter and $E$ is the neutrino energy. We present simple analytical approximations to these `running' parameters. We show that for the moderately large value of $\theta_{13}$, as discovered by the reactor experiments, the running of the mixing angle $\theta_{23}$ and the CP violating phase $\delta$ can be neglected. It simplifies the analysis of the resulting expressions for the oscillation probabilities considerably. Approaches which attempt to directly provide approximate analytical expressions for the oscillation probabilities in matter suffer in accuracy due to their reliance on expansion in $\theta_{13}$, or in simplicity when higher order terms in $\theta_{13}$ are included. We demonstrate the accuracy of our method by comparing it to the exact numerical result, as well as the direct approximations of Cervera et al., Akhmedov et al., Asano and Minakata, and Freund. We also discuss the utility of our approach in figuring out the required baseline lengths and neutrino energies for the oscillation probabilities to exhibit certain desirable features.


Introduction
When performing long-baseline neutrino oscillation experiments on the Earth with accelerator based beams, or when detecting atmospheric neutrinos coming from below, the neutrinos necessarily traverse the Earth's interior [1][2][3][4][5][6][7]. This makes the understanding of matter effects [8][9][10] on the neutrino oscillation probabilities an indispensable part of analyzing such experiments.
Though the exact three-flavor neutrino oscillation probabilities in constant-density matter can be expressed analytically, 1 the exact formula is too complicated to yield physical insight [12][13][14][15]. Consequently, various analytical approximations have been devised to better understand the physics potential of various neutrino experiments [9,[16][17][18][19][20][21]. These approximations relied on expansions in the small parameters α = δm 2 21 /δm 2 31 ≈ 0.03 and/or s 13 = sin θ 13 in one form or another, a systematic treatment of which can be found in Ref. [19]. Of the formulae available, the one derived by Cervera et al. in Ref. [22] has been used frequently due to its simplicity. Unfortunately, the accuracy of the Cervera et al. formula suffers when the value of θ 13 is as large as was measured by Daya Bay [23] and RENO [24], consistent with both earlier and later determinations by T2K [25], MINOS [26,27], and Double Chooz [28,29]. This is due to the fact that the Cervera et al. formula includes terms of order O(α 2 ), O(αs 13 ), and O(s 2 13 ). Given that the current world average of s 13 = sin θ 13 is about 0.15 [30], not all terms of the same order are included. Asano and Minakata [31] have derived the order O(αs 2 13 ) and O(s 4 13 ) corrections to the Cervera et al. formula, but the simplicity of the original expressions is lost.
In previous papers [32,33], we had advocated the utility of the Jacobi method [34] in deriving approximate formulae for neutrino oscillation in matter. The formula presented there maintained the expressions for the oscillation probabilities in vacuum, and only replaced the mixings angles and the mass-squared differences with their effective values in matter. The approximation worked very well except when θ 13 was very small, a possibility that could not be ignored until the Daya Bay/RENO measurements. In this paper, we reintroduce the method with further refinements which improve the accuracy of the approximation, while maintaining its ease of use. This paper is organized as follows. In section 2, we explain our approach to the matter effect problem, and list all the formulae necessary to calculate our approximation. In section 3, we demonstrate the accuracy of our approximation at various baseline lengths, different mass hierarchies, and different values of the CP violating phase δ. Comparisons with the approximations of Cervera et al. [22] and Asano-Minakata [31] are also made. In section 4, we show how simple calculations using our approximation can be used to derive the baselines and energies at which the oscillation probabilities exhibit desirable features. We conclude in Section 5. Detailed derivation of our approximation is given in appendices B and C.  Figure 1. The dependence of the line-averaged mass density ρ on the baseline length L based on the Preliminary Reference Earth Model [35]. The labels on the right edge of the frame indicate the corresponding values of a/E. The green and red dashed lines indicate ρL = 54000 km · g/cm 3 and ρL = 32000 km · g/cm 3 , respectively, which are conditions that will be discussed in section 4.

The Approximation
In the following, we use the conventions and notation reviewed in Appendix A.

Diagonalization of the Effective Hamiltonian
If the matter density along the baseline is constant, the effective Hamiltonian which governs the evolution of neutrino flavor in matter is given by where U is the neutrino mixing matrix in vacuum, and Here, N e is the electron number density, ρ the matter mass density along the baseline, and E the neutrino energy. The above term appearing in the ee-component of H a is due to the interaction of ν e with the electrons in matter via W -exchange, and Eq. (2.2) assumes N e = N p = N n in Earth matter. It also assumes E M W since the W -exchange interaction is approximated by a point-like four-fermion interaction. Z-exchange effects are flavor universal and only contribute a term proportional to the unit matrix to H a , which can be dropped.
If we write the eigenvalues of H a as λ i (i = 1, 2, 3) and the diagonalization matrix as ∼ U , that is then the neutrino oscillation probabilities in matter are obtained by simply taking their expressions in vacuum and replacing the elements of the mixing matrix U and the masssquare differences δm 2 ij with their "effective" values in matter: Note that a is E-dependent, which means that both ∼ U αi and δλ ij are also E-dependent. They also depend on the baseline length L since the average matter density ρ along a baseline varies with L. The L-dependence of the average ρ and the corresponding value of a/E are shown in Fig. 1.
For anti-neutrino beams, the flavor-evolution Hamiltonian in matter is (2.5) In comparison to Eq. (2.1), the CP violating phase δ in U and the matter-effect term a both acquire minus signs. Let us write the eigenvalues of H a as λ i (i = 1, 2, 3) and the diagonalization matrix as U , that is Note that the tilde above U here is flipped to distinguish it from ∼ U in Eq. (2.3). The antineutrino oscillation probabilities in matter are then obtained by making the replacements in the vacuum expressions.

Effective Mixing Angles
As mentioned in the introduction, it is possible to write down exact analytical expressions for ∼ U αi and δλ ij , as well as their anti-neutrino counterparts [12]. However, simpler and more transparent approximate expressions can be obtained using the Jacobi method [34]. In our approach, the matrix ∼ U is approximated by taking the conventional parameterization of the matrix U in vacuum in terms of three mixing angles and a CP-violating phase, namely U = R 23 (θ 23 , 0)R 13 (θ 13 , δ)R 12 (θ 12 , 0) , (2.8) and replacing the values of θ 12 and θ 13 with their effective values in matter: 2 The angle θ 23 and the CP-violating phase δ are kept at their vacuum values. The eigenvalues λ i (i = 1, 2, 3) of H a are also given approximate expressions. Thus, to take matter effects into account when calculating neutrino oscillation probabilities, all that is necessary is to take their expressions in terms of the mixing angles and CP-phase in vacuum as is, and replace the two angles as well as the mass-squared differences with their effective values in matter: θ 12 → θ 12 , θ 13 → θ 13 , δm 2 ij → δλ ij = λ i − λ j . This simplifies the calculation considerably, and allows for a transparent understanding of how matter-effects affect neutrino oscillation by looking at the a-dependence of the effective parameters.
Similarly, in the anti-neutrino case, we approximate Again, θ 23 and δ are unaffected while θ 12 and θ 13 are replaced by their effective values in matter. The eigenvalues λ i (i = 1, 2, 3) of H a are given approximate expressions as in the neutrino case. Thus, the calculation of matter effects for anti-neutrino beams entails the replacements

Definition of β
In the following, we list the expressions necessary to calculate the effective mixing angles and the effective mass-squared differences for the neutrino and anti-neutrino cases separately. Detailed derivation of our approximation is given in Appendix B. We will also show plots depicting how various effective parameters depend on the matter-effect parameter a. Due to the wide separation in scale between δm 2 21 and δm 2 31 , we find it convenient to introduce the parameter β via 3 and plot our effective parameters as functions of β instead of a. β = 0 corresponds to a = |δm 2 31 |, β = −2 to a = δm 2 21 , and so on.

Neutrino Case
For the neutrino case, the effective mixing angles in matter are given by 2 Absorbing matter effects into θ12 and θ13 is similar in spirit to Ref. [18]. 3 We avoid the use of the symbols α or A since they often respectively denote δm 2 21 /δm 2 31 and a/δm 2 31 in the literature.

Anti-Neutrino Case
For the anti-neutrino case, we need to flip the sign in front of the matter effect factor a, and the effective mixing angles are given by tan 2θ 12 = (δm 2 21 /c 2 13 ) sin 2θ 12 (δm 2 21 /c 2 13 ) cos 2θ 12 + a , tan 2θ 13 = (δm 2 31 − δm 2 21 s 2 12 ) sin 2θ 13 (δm 2 31 − δm 2 21 s 2 12 ) cos 2θ 13 + a . (2.16) The β-dependence of these angles are shown in Fig. 2(b), and that of the sines of twice these angles in Fig. 3(b). The three eigenvalues of the effective Hamiltonian are approximated by where the upper(lower) sign is for the normal(inverted) hierarchy, with and c 2 12 = cos 2 θ 12 . The β-dependence of λ + and λ ± for the normal and inverted mass hierarchies are shown in Figs. 4(c) and (d), respectively. That of λ − is shown in Fig. 15(b). For the normal hierarchy case, δm 2 31 > 0, the above expressions simplify to

Comment on Coding
As seen above, our formulae for both the neutrino and anti-neutrino cases are fairly compact and easy to code. In particular, the effective mixing angles for the neutrino and antineutrino cases can be calculated with the same code by simply flipping the sign of the matter-effect parameter a, cf. Eqs. (2.12) and (2.16). The same can be said of λ ± and λ ± defined in Eqs. (2.14) and (2.18). In the case of λ ± and λ ± , one also needs to make the swap λ + ↔ λ − but otherwise the code will be essentially the same. Once the effective mixing angles and the λ's are obtained, one only needs to input them into the code for calculating the vacuum oscillation probabilities, and we are done.

Demonstration of the Accuracy of the Approximation
In this section, we demonstrate the accuracy of our approximation. For the vacuum values of the mixing angles and mass-squared differences, we use the global fit values from Ref. [30] listed in Table 1.   Figure 6. Comparison of the approximation formulae of Cervera et al., Asano-Minakata, and this work at L = 810 km, which is the distance from Fermilab to NOνA. The average matter density for this baseline length is ρ = 2.80 g/cm 3 . The dashed line gives the exact numerical result. [35] We begin by comparing our approximation to those of Cervera et al. [22] and of Asano and Minakata [31]. In Fig. 5(a), we plot the approximate ν µ → ν e oscillation probabilities calculated using the three approximations against the exact numerical result for the baseline length L = 4000 km. This is the distance used by Asano and Minakata in Table 1. Best-fit values of oscillation parameters taken from Ref. [30].  Ref. [31] to demonstrate the strength of their formula. The average Earth matter density for this baseline is 3.81g/cm 3 [35]. We consider the normal hierarchy case, δm 2 31 > 0, with the CP violating phase δ set to zero. The differences between the exact and approximate formulae are plotted in Fig. 5(b). As can be seen, at this baseline, both the Asano-Minakata formula and our approximation work much better than the Cervera et al. formula.
The comparison at a shorter baseline length of L = 810 km, which is the distance from Fermilab to NOνA, is shown in Fig. 6. There, all three approximations work well, with our approximation being the most accurate.
The situation changes at the longer baseline length of L = 8770 km, which is the distance from CERN to Kamioka [36], as can be seen in Fig. 7. There, the Cervera et al. formula greatly overestimates P (ν µ → ν e ), while the Asano-Minakata formula leads to negative probability for E ∼ 4 GeV. In comparison, our approximation remains accurate.
The accuracy of our approximation for both the neutrino and anti-neutrino cases, and both mass hierarchies, for different values of the CP violating phase δ, is demonstrated in Figs. 8 and 9 for the two baselines L = 1300 km and L = 2300 km, respectively. These distances correspond to those between Fermilab and Homestake (1300 km), and CERN and Pyhäsalmi (2300 km) [37]. As is evident, our approximation maintains its accuracy for all energy ranges and mass densities.

Determination of the Mass Hierarchy from ν e Oscillations
Consider the ν e survival probability in matter which is given by where we have assumed that a δm 2 21 so that s 12 ≈ 1 is a good approximation. Similarly, we find:    Fig. 3, it is clear that the factor sin 2 (2θ 13 ) in these expressions behaves quite differently depending on the mass hierarchy. For normal hierarchy sin 2 (2θ 13 ) will peak around a ≈ δm 2 31 but for the inverted hierarchy case it will not. This will become manifest if the factor sin 2 ( ∼ ∆ 32 /2) also peaked at or near the same energy, For the normal hierarchy case, when a ≈ δm 2 31 we have Therefore, . (4.7) From Fig. 1, it is clear that ρL < 54000 (km · g/cm 3 ) as long as the neutrino beam does not enter the core of the Earth, at which point the constant average matter density approximation breaks down. Therefore, in order to take ∼ ∆ 32 /2 as close as possible to π/2 while preventing the beam from entering the Earth's core, we need L ∼ 10000 km.
For instance, if we take L = 10000 km for which ρ = 4.53 g/cm 3 , we have ρL ≈ 45300 km · g/cm 3 . The value of ∼ ∆ 32 /2 at resonance a ≈ δm 2 31 is then π 2 × 45300 54000 = 0.42 π , (4.8) leading to an oscillation peak/dip factor of sin 2 ( ∼ ∆ 32 /2) = 0.94. Using Eq. (2.2), the neutrino beam energy at which a ≈ δm 2 31 is found to be Indeed, in Fig. 10(a) we show the exact ν e survival probabilities at L = 10000 km for both hierarchies, and we can see that the normal hierarchy case dips by over 95% around E ∼ 6.5 GeV. Thus, our rough estimates give a highly reliable result. If we take a somewhat shorter baseline of L = 8770 km, which is the distance between CERN and Kamioka [36], we have ρ = 4.33 g/cm 3 , and ρL ≈ 38000 km · g/cm 3 . The value of ∼ ∆ 32 /2 at resonance a ≈ δm 2 31 is then π 2 × 38000 54000 = 0.35 π , (4.10) leading to an oscillation peak/dip factor of sin 2 ( ∼ ∆ 32 /2) = 0.8, which is still fairly prominent. Using Eq. (2.2), the neutrino beam energy at which a ≈ δm 2 31 is found to be The actual oscillation peak occurs slightly off resonance around E = 6.5 GeV as can already be seen in Fig. 7. Comparison of P (ν e → ν e ) at L = 8770 km with δ = 0 for the normal and inverted hierarchies are shown in Fig. 10(b). P (ν e → ν µ ) is compared in Fig. 11(b). The differences between the normal and inverted hierarchies for both baselines is manifest, indicating that measuring these oscillation probabilities at this baseline would allow us to determine the mass hierarchy easily. (We consider the dependence on the CP violating phase δ in the next section.) Eqs. (4.1), (4.4), and (4.5) also suggest that the measurement may provide a better determination of sin 2 θ 23 .

The "Magic" Baseline
The "magic" baseline is the baseline at which the dependence of P (ν e → ν µ ) on the CP violating phase δ vanishes [38]. 4 Looking at Eq. (4.2), the only term without δ-dependence is the | To make every other term vanish, we must have Therefore, the magic baseline condition is If we are in the energy and mass-density range such that δm 2 21 < a < |δm 2 31 |, we can see from Fig. 4 that δλ 21 ≈ a = 2 √ 2G F N e E, and the above condition reduces to √ 2G F N e L ≈ 2nπ , (4.14) which is the usual magic baseline condition. Using Eq. (2.2), this condition for the n = 1 case becomes that is ρL km · g/cm 3 ≈ 32000 .
(4.16) This is satisfied for L ≈ 7500 km as can be read off of Fig. 1. Indeed, in Fig. 11(a) we plot the bands that P (ν e → ν µ ) at L = 7500 km sweeps for both mass hierarchies when δ is varied throughout its range of [0, 2π]. We can see that the dependence on δ is very weak. However, if we look at Eq. (4.3) carefully, it is clear that all terms that include the CP violating phase δ are multiplied by c 12 which goes to zero when a δm 2 21 . Indeed, this was why δ did not appear in Eq. (4.4). The condition a δm 2 21 demands ρ g/cm 3 which is clearly satisfied around the oscillation peak for the L = 8770 km case just discussed in the previous section. Thus, P (ν e → ν µ ) for this baseline is also only very weakly dependent on δ as shown in Fig. 11(b). We can conclude that, in general, as long as Eq. (4.17) is satisfied, one does not need to be at a specific "magic" baseline to suppress the δ-dependence of P (ν e → ν µ ).

Summary
We have presented a new and simple approximation for calculating the neutrino oscillation probabilities in matter. Our approximation was derived utilizing the Jacobi method [34], and we show in the appendix that at most two rotations are sufficient for approximate diagonalization of the effective Hamiltonian. The two rotation angles are absorbed into effective values of θ 12 and θ 13 . As explained in detail in the appendix, the approximation works when θ 13 = O(ε), where ε = δm 2 21 /|δm 2 31 | = 0.17, a condition which has been shown to be satisfied by Daya Bay [23] and RENO [24]. Our formulae are compact and can easily be coded as well as be manipulated by hand for back-of-the-envelope calculations. The application of our method to finding the ν e → ν µ , ν τ resonance conversion condition, and that to the determination of the 'magic' baseline [38,39] have been demonstrated.
In this paper, only the matter effect due to Standard Model W exchange between the neutrinos and matter was considered. New Physics effects which distinguish between neutrino flavor would add extra terms to the effective Hamiltonian, which would require further rotations for diagonalization. This has been discussed previously in Ref. [33], and the potential constraints on New Physics from long baseline neutrino oscillations experiments in Refs. [40][41][42]. Updates to these works will be presented in future publications [43].

A Conventions, Notation, and Basic Formulae
Here, we collect the basic formulae associated with neutrino oscillation in order to fix our notation and conventions.
The Majorana phases, α 21 and α 31 , only appear in lepton-number violating processes such as neutrinoless double beta decay, and cannot be determined via neutrino→neutrino oscillations. The absolute scale of the neutrino masses also remain undetermined since neutrino oscillation is an interference effect.

A.2 Neutrino Oscillation
If a neutrino of flavor α is created at x = 0 with energy E, then the state of the neutrino at x = 0 is Assing m j E we can approximate so that and we find Therefore, the amplitude of observing the neutrino of flavor β at x = L is given by (dropping the irrelevant overall phase) and Thus, the probability of oscillation from |ν α to |ν β with neutrino energy E and baseline L is given by only two of the three ∆ ij 's in Eq. (A.13) are independent. Eliminating ∆ 32 from Eq. (A.13) for the α = β case yields and for the α = β case we have where J (α,β) is the Jarlskog invariant [49]: Note that our notation differs from that of Cervera et al. in Ref. [22]. There, the symbol ∆ij is defined without the factor of L, that is, ∆ij = δm 2 ij /2E. It also differs from that used by Freund in Ref. [17] where ∆ = δm 2 31 , and∆ = δm 2 31 L/4E. Huber and Winter in Ref. [38] define ∆ = δm 2 31 L/4E, which is also used in Ref. [48]. So care is necessary when comparing formulae.

A.3 Matter Effects
If the matter density along the baseline is constant, matter effects on neutrino oscillations can be taken into account by replacing the PMNS matrix elements and mass-squared differences with their "effective" values in matter: where ∼ U is the unitary matrix that diagonalizes the modified Hamiltonian, The factor a is due to the interaction of the |ν e component of the neutrinos with the electrons in matter via W -exchange: Assuming N e = N p ≈ N n in Earth matter, N e for mass density per unit volume of ρ can be expressed using Avogadro's number N A = 6.02214129 × 10 23 mol −1 as N e = N p ≈ ρ N A /2 = 3.011 × 10 23 /cm 3 × ρ g/cm 3 . (A.25) Thus, putting back powers of c to convert from natural to conventional units, we find For anti-neutrino beams, a is replaced by −a in Eq. (A. 22). Note that a is E-dependent, which means that both ∼ U and ∼ ∆ ij are also E-dependent. It is also assumed that E M W since the W -exchange interaction is approximated by a point-like four-fermion interaction in deriving this expression.

B.1 Setup
As mentioned in the introduction, it is possible to write down exact analytical expressions for ∼ ∆ ij and ∼ U αi [12]. However, simpler and more transparent approximate expressions can be obtained using the Jacobi method as will be shown in the following.
We introduce the matrix Q = diag(1, 1, e iδ ) , (B.1) and start with the partially diagonalized Hamiltonian: The matrix Q serves to rid H a of any reference to the CP violating phase δ. The strategy we used in our previous papers [32,33]  That is, = ε 2 . So care is necessary when comparing formulae.

B.2 Diagonalization of a 2 × 2 Matrix
Recall that for 2 × 2 real symmetric matrices, such as and we obtain where the upper and lower signs are for the cases α < γ and α > γ, respectively. The Jacobi method [34] entails iteratively diagonalizing 2 × 2 submatrices of a larger matrix in the order that requires the largest rotation angle at each step. In the limit of infinite iterations of this procedure, the matrix will converge to a diagonal matrix. In the case of H a given in Eq. (B.2), at most two iterations are sufficient to achieve approximate diagonalization, neglecting off-diagonal elements of order O(ε 2 s 13 |δm 2 31 |), regardless of the size of a. We demonstrate this in this appendix.

B.3.1 Mixing Angles and Mass-squared Differences
Let us first evaluate the sizes of the sines and cosines of the three vacuum mixing angles in comparison to the parameter ε defined in Eq. (B.3). The current best fit values for the mass-squared differences and mixing angles are listed in Table 1
We define where c φ = cos φ , s φ = sin φ , tan 2φ ≡ as 12 sin 2θ 13 δm 2 31 + as 2 The angle φ is in the first quadrant when δm 2 31 > 0, and in the fourth quadrant when δm 2 31 < 0. Then, where the upper(lower) sign corresponds to the normal(inverted) hierarchy case with λ ± ≡ λ + + (δm 2 31 + as 2 13 ) ± λ + − (δm 2 31 + as 2 13 ) 2 + 4a 2 s 2 12 c 2 13 s 2 13 2 . (B.23) The β-dependences of λ ± and φ are shown in Fig. 4 ((a) and (b)), and Fig. 14(a), respectively, for both mass hierarchies. For the normal hierarchy case, δm 2 31 > 0, the values of λ ± away from the level crossing point a ∼ δm 2 31 are approximately for all a. This approximation is obtained by first noting that φ is significantly different from zero only when a δm 2 21 . The expansion of λ + in the denominator of the right-hand-side of Eq. (B.21) in powers of δm 2 21 /a was given in Eq. (B.19). Keeping only the first two terms, and noting also that s 12 ≈ 1 to the same order when a δm 2 21 (c.f. Eq. (??)) we obtain Eq. (B.27). The β-dependence of the difference φ − φ is plotted in Fig. 14(b), and we can see that the disagreement is at most O(ε 4 ). Thus, we replace φ with φ in the following. Now, the effective Hamiltonian after the second rotation was given by Eq. (B.22). Note that all of the non-zero off-diagonal elements include the factor ac 12 , which is never larger than O(ε 2 |δm 2 31 |) regardless of the value of a as discussed above. They also all include a factor of s 13 , which is O(ε) as we have seen in Eq. (B.9). Therefore, all off-diagonal elements of H a are of order O(ε 2 s 13 |δm 2 31 |) = O(ε 3 |δm 2 31 |) or smaller regardless of the size of a. Note that had the value of s 13 been smaller, the sizes of the neglected terms would have been proportionately smaller also. We conclude that, at this point, off-diagonal elements are negligible and a third rotation is not necessary.

B.3.4 Absorption of φ into θ 13
From the above consideration, we conclude that the matrix which diagonalizes H a , Eq. (B.2), is given approximately by V W , and that the effective neutrino mixing matrix becomes (B.28) Using Here, we argue that that is, the 2-3 rotation becomes a 1-3 rotation when commuted through R 12 (θ 12 , 0). This is due to the fact that φ only becomes non-negligible when a δm 2 12 where s 12 ≈ 1 and c 12 ≈ 0, which means and it is straightforward to see that where s φ = sin φ and c φ = cos φ . In the range a δm 2 21 , the angle φ is very small and both R 23 (φ , 0) and R 13 (φ , 0) are approximately unit matrices and Eq. (B.31) is trivially satisfied. Curiously, this approximation breaks down around a ∼ δm 2 31 for the normal hierarchy case when θ 13 is O(ε 2 ) or smaller, as is discussed in appendix C. However, given that the current experimentally preferred value of θ 13 is O(ε), the approximation is valid. Thus, The diagonal phase matrix Q appearing rightmost in the above matrix product can be absorbed into the redefinition of the major phases and can be dropped. Thus, we arrive at our final approximation in which the vacuum mixing angles are replaced by their effective values in matter (B.37) and the eigenvalues of the effective Hamiltonian are given by Note that of the mixing angles, only θ 12 and θ 13 are shifted. θ 23 and δ stay at their vacuum values.

B.4.1 First Rotation
For the anti-neutrino case, the matter effect parameter a acquires a minus sign. Thus, the effective hamiltonian to be diagonalized is The largest off-diagonal element is the 1-2 element. Therefore, our first step is to diagonalize the 1-2 submatrix. Define where c ϕ = cos ϕ , s ϕ = sin ϕ , tan 2ϕ ≡ − ac 2 13 sin 2θ 12 δm 2 21 + ac 2 13 cos 2θ 12 The β-dependences of θ 12 and λ ± are shown in Fig. 15. Note that in contrast to the neutrino case, there is no level crossing. θ 12 decreases monotonically toward zero as a is increased. For a δm 2 21 , s 12 and c 12 behave as and we see that, this time, we have ac 12 ≈ a and as 12 ≈ δm 2 21 s 12 c 12 /c 2 13 = O(ε 2 |δm 2 31 |). These β-dependences of s 12 , c 12 , as 12 , and ac 12 are shown in Fig. 16(a) and (b). In the range a δm 2 21 , λ ± can be expanded as The angle φ is in the fourth quadrant when δm 2 31 > 0, and the first quadrant when δm 2 31 < 0. Using W , we find where the upper(lower) sign corresponds to normal(inverted) mass hierarchy with The β-dependence of λ ± and φ are shown in Fig. 4 ((c) and (d)), and Fig. 17(a), respectively, for both normal and inverted mass hierarchies. For the normal hierarchy case, δm 2 31 > 0, there is no level crossing, and λ ± are well approximated by Level crossing occurs for the inverted hierarchy case, δm 2 31 < 0, in which we have  Fig. 17(b), and it is clear that the difference is negligible. Now, the effective Hamiltonian after the second rotation was given by Eq. (B.51). Note that all of the non-zero off-diagonal elements include the factor as 12 , which is never larger than O(ε 2 |δm 2 31 |) regardless of the value of a as discussed above. They also all include a factor of s 13 , which is O(ε) as we have seen in Eq. (B.9). Therefore, all off-diagonal elements of H a are of order O(ε 2 s 13 |δm 2 31 |) = O(ε 3 |δm 2 31 |) or smaller regardless of the size of a. We conclude that, at this point, off-diagonal elements are negligible and a third rotation is not necessary.
C Commutation of R 13 and R 23 through R 12 In the derivation of our approximation formulae above, Eqs. (B.31) and (B.59) played crucial roles in allowing the second rotation angle to be absorbed into θ 13 . In this appendix, we evaluate the validity of these approximations.

C.1 Neutrino Case
The difference between the two sides of Eq. (B.31) is given by    It is clear that δR will vanish in the two limits a → 0 where s 12 → s 12 , c 12 → c 12 , s φ → 0, and c φ → 1, and a → ∞ where s 12 → 1, c 12 → 0, s φ → c 13 (−s 13 ), and c φ → s 13 (c 13 ) for normal(inverted) hierarchy. The question is whether δR will stay negligible in between as s 12 runs from s 12 to 1, c 12 from c 12 to 0, s φ from 0 to c 13 (normal) or −s 13 (inverted), and c φ from 1 to s 13 (normal) or c 13 (inverted) as shown in Fig. 18. The dependence of the non-zero elements of δR on β = − log ε (a/|δm 2 31 |) is shown in Fig. 19. The bumps at a ∼ δm 2 21 for both hierarchies, and that at a ∼ δm 2 31 for the normal hierarchy, happen due to the θ 12 factor competing with the φ factor as one of them goes through a resonance while the other damps to zero. The heights of the bumps depend on the narrowness of the resonances.
For the case shown in Fig. 19, which was generated with the numbers in Table 1 as input, all elements of δR are O(ε 3 ) or smaller for the entire range of a, with the maximum value of ∼ 0.01 ≈ 2ε 3 occurring in c 12 s φ near a ∼ δm 2 31 in the normal hierarchy case. Since the size of the third rotation angle we neglected in the Jacobi procedure was O(ε 2 s 13 ), Eq. (B.31) is valid to the same order provided s 13 = O(ε).
For smaller values of s 13 , the resonance at a ∼ δm 2 31 would have been narrower, and the peaks in c 12 s φ and c 12 (1 − c φ ) higher. This is illustrated in Fig. 20. In the limit s 13 → +0, s φ and 1 − c φ will become step functions at β ∼ 0, and the maximum height of the peak will be c 12 (a ∼ δm 2 31 ) ≈ s 12 c 12 ε 2 = 0.46 ε 2 = 0.014 , (C.2) as can be discerned from Eq. (B.17). This is the same as the asymptotic value of ac 12 /δm 2 31 discussed earlier. While this value may not seem particularly large, only a factor of 3/2 larger than the peak in Fig. 19(a), it is parametrically O(ε 2 ). On the other hand, the third rotation angle neglected in the Jacobi procedure was O(ε 2 s 13 ). Thus, using Eq. (B.31) would lead to dropping terms that are larger than the ones we keep when s 13 = O(ε 2 ) or smaller. Also, the sudden change in the accuracy of Eq. (B.31) across a ∼ δm 2 31 , as can be seen in Fig. 20, will lead to kinks in the resulting oscillation probabilities.
The dependence of the non-zero elements of δR on β = − log ε (a/|δm 2 31 |) is shown in Fig. 22, which was generated with the numbers in Table 1 as input. We can see that all    elements of δR are O(ε 3 ) or smaller for the entire range of a, with the maximum value of ∼ 0.01 ≈ 2ε 3 occuring in s 12 s φ near a ∼ δm 2 31 in the inverted hierarchy case. Since the size of the third rotation angle we neglected in the Jacobi procedure was O(ε 2 s 13 ), Eq. (B.59) is valid to the same order provided s 13 = O(ε).
For smaller values of s 13 , the resonance at a ∼ δm 2 31 would have been narrower, and the peaks in s 12 s φ and s 12 (1 − c φ ) higher. This is illustrated in Fig. 23. In the limit s 13 → +0, s φ and 1 − c φ will become step functions at β ∼ 0, and the maximum height of the peak will be the same as Eq. (C.2), and the asymptotic value of as 12 /|δm 2 31 |, as can be discerned from Eq. (B.46). This is parametrically O(ε 2 ), while the third rotation angle neglected in the Jacobi procedure was O(ε 2 s 13 ). Thus, using Eq. (B.59) would lead to dropping terms that are larger than the ones we keep when s 13 = O(ε 2 ) or smaller. Also, the sudden change in the accuracy of Eq. (B.59) across a ∼ δm 2 31 , as can be seen in Fig. 23, will lead to kinks in the resulting oscillation probabilities.