Analytical approximation of the neutrino oscillation matter effects at large θ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 = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ 2\sqrt{2}{G_F}{N_e}E $\end{document}, where Ne 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 θ13, as discovered by the reactor experiments, the running of the mixing angle θ23 and the CP violating phase δ 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 θ13, or in simplicity when higher order terms in θ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][8]. This makes the understanding of matter effects [9][10][11][12] on the neutrino oscillation probabilities an indispensable part of analyzing such experiments. These matter effects can of course be calculated numerically for arbitrary matter profiles, but approximate analytical expressions are useful not only for making initial estimates on the requirements one must place on long-baseline experiments, but in obtaining a deeper understanding of the physics involved.
The exact three-flavor neutrino oscillation probabilities in constant-density matter can be expressed analytically [12][13][14][15][16][17][18][19][20][21][22][23]. This requires the diagonalization of the 3 × 3 effective Hamiltonian in matter whose ee-element in the flavor basis is shifted by a = 2 √ 2G F N e E, where N e is the electron density and E is the neutrino energy. The eigenvalues of the effective Hamiltonian yield the effective neutrino mass-squared differences in matter, 1 while the diagonalization matrix is multiplied with the vacuum neutrino mixing matrix to yield its in-matter counterpart. Many authors adopt the standard vacuum parameterization of the mixing matrix to the matter version, and absorb matter effects into shifts of the mixing angles and CP violating phase, yielding the effective values of these parameters in matter [13,14,16,22]. Thus, the neutrino oscillation probabilities in matter can be obtained from those in vacuum by simply replacing the mass-squared differences, mixing angles, and CP violating phase with their effective values. Unfortunately, the final exact expressions for the neutrino oscillation probabilities obtained this way are too complicated to yield physical insight, especially if re-expressed in terms of the vacuum parameters.
Consequently, various analytical approximations have been devised to better understand the physics potential of various neutrino experiments [12,22,[25][26][27][28][29][30][31][32]. 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. [29]. In some cases the matter-effect parameter a = 2 √ 2G F N e E was also assumed to be small [12,26]. For instance, the formula of Cervera et al. in ref. [27] and that of Ahkmedov et al. in ref. [29] include terms of order O(α 2 ), O(αs 13 ), and O(s 2 13 ). Unfortunately, the accuracies of these formulae suffer when the value of θ 13 is as large as was measured by Daya Bay [33,34] and RENO [35], consistent with both earlier and later determinations by T2K [36], MINOS [37,38], and Double Chooz [39,40]. Given that the current world average of s 13 = sin θ 13 is about 0.15 [41], the terms included are not all of the same order. Asano and Minakata [32] 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. Further improvements in accuracy are possible at the expense of simplicity, as was shown by Freund in ref. [22] where an approximate expression for the oscillation probability P (ν e → ν µ ) including all orders of θ 13 was derived.

JHEP04(2014)047
In previous papers [42,43], we had argued the advantage of not expressing the neutrino oscillation probabilities in matter directly in terms of the vacuum parameters, but to maintain their expressions in terms of the effective parameters in matter which 'ran' with the parameter a = 2 √ 2G F N e E. Further, it was shown that the Jacobi method [44] could be utilized to find approximate expressions for the 'running' parameters in a systematic fashion, leading to fairly simple and compact expressions. In particular, it was shown that the effective values of θ 23 and the CP violating phase δ do not 'run' to the order considered, maintaining their vacuum values at all neutrino energies and baselines. (The non-running of θ 23 and δ has also been discussed in ref. [17].) The a-dependence of the resulting expressions for the oscillation probabilities in terms of the approximate running parameters could be analyzed in a simple manner, facilitating the understanding of matter effects.
The approximation of refs. [42,43] worked extremely 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 for large θ 13 , 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 the approximate 'running' parameters in our approach. Approximate oscillation probabilities are obtained by replacing the mass-squared differences and mixing angles in the vacuum oscillation probabilities with their effective 'running' values. 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. [27], Akhmedov et al. [29], Asano-Minakata [32], and Freund [22] 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.

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, 2 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 2 At baseline length L = 10690 km or longer, the neutrino beam crosses the core-mantle-boundary and experiences a sudden jump in matter density. See ref. [46] for treatments of non-adiabatic transitions.  Figure 1. The dependence of the line-averaged mass density ρ on the baseline length L based on the Preliminary Reference Earth Model [45]. 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. 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 'running' values in matter [9][10][11]: 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 figure 1.

JHEP04(2014)047
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 running mixing angles
While it is possible to write down exact analytical expressions for ∼ U αi and δλ ij , as well as their anti-neutrino counterparts [16], simpler and more transparent approximate expressions are often desirable. One popular approach is to expand the probability formulae in terms of small parameters such as δm 2 21 /|δm 2 31 | and θ 13 . Our approach, however, utilizes the Jacobi method [44]. Instead of obtaining approximations for the probabilities directly, we derived the approximations for the effective mixing parameters. In the following two sections, we list the expressions necessary to calculate the effective running mixing angles and the effective running mass-squared differences for the neutrino and anti-neutrino cases separately. Detailed derivation of our approximation is given in appendix B.

JHEP04(2014)047
where θ 12 and θ 13 are given by tan 2θ 12  while the angle θ 23 and the CP-violating phase δ at kept at their vacuum values [17]. The eigenvalues λ i (i = 1, 2, 3) of H a are also given approximate running expressions: and s 2 12 = sin 2 θ 12 . For the inverted hierarchy case, δm 2 31 < 0, the above expressions simplify to 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 running 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.

Anti-neutrino case
Similarly, in the anti-neutrino case, the mixing matrix can be parameterized by: (2.14) Note that the sign in front of the matter effect parameter a is flipped relative to the neutrino case, so these effective mixing angles will be different. Our approximation is given by Again, θ 23 and δ are unaffected while θ 12 and θ 13 are replaced by their effective running values in matter. The eigenvalues λ i (i = 1, 2, 3) of H a are given approximate running expressions as in the neutrino case. 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 . For the normal hierarchy case, δm 2 31 > 0, the above expressions simplify to Thus, the calculation of matter effects for anti-neutrino beams entails the replacements θ 12 → θ 12 , θ 13 → θ 13 , δm 2 ij → δλ ij = λ i − λ j .

The β-dependence of mixing parameters
We show plots depicting how our various effective parameters run with 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 running parameters as functions of β instead of a. Here β = 0 corresponds to a = |δm 2 31 |, β = −2 to a = δm 2 21 , and so on. The dependence of the effective mixing angles on β are shown in figure 2 and that of the sines of twice these angles in figure 3. The β-dependence of approximate eigenvalues of the effective Hamiltonian are shown in figure 4. 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 mixing angles Figure 2. The dependences of the effective mixing angles on β = − log ε (a/|δm 2 31 |) for the neutrino (a) and antineutrino (b) cases. β = 0 corresponds to a = |δm 2 31 |, and β = −2 to a = δm 2 21 . The βdependences of θ 13 and θ 13 depend on the mass hierarchy: when δm 2 31 > 0 (normal hierarchy, NH) θ 13 increases toward π/2 whereas θ 13 decreases toward zero, while in the δm 2 31 < 0 case (inverted hierarchy, IH), it is the other way around.  . The β-dependences of the sines of twice the effective mixing angles for the neutrino (a) and antineutrino (b) cases. The difference in the behavior of the effective θ 13 mixing angle for normal and inverted hierarchies will allow us to determine which is chosen by nature.

Demonstration of the accuracy of the approximation
In this section, we plot neutrino oscillation probabilities in several scenarios to demonstrate the accuracy of our approximation. As seen in the previous section, 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 anti-neutrino cases can be calculated with the same code by simply flipping the sign of the matter-effect parameter a, cf. eqs. (2.10) and (2.16). The same can be said of λ ± and λ ± defined in eqs. (2.12) and (2.18). In the case of λ ± and λ ± , one also needs to make the swap λ + ↔ λ − but otherwise the code will be essentially the same. For the vacuum values of the mixing angles and mass-squared    We begin by comparing our approximation to eq. (16) of Cervera et al. [27], eq. (3.5) of Akhmedov et al. [29], sum of eqs. (4.2) to (4.4) of Asano and Minakata [32], and eq. (36) of Freund [22]. Note that both Cervera et al.  In figure 5(a), we plot the approximate ν µ → ν e oscillation probabilities calculated using these three approximations against the exact numerical result for the baseline length L = 4000 km. This is the distance used by Asano and Minakata in ref. [32] to demonstrate the strength of their formula. The line-averaged constant Earth matter density 4 for this baseline is 3.81g/cm 3 which has been estimated using the Preliminary Reference Earth Model (PREM) [45]. 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 figure 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. or the Akhmedov et al. formulae. The Freund formula works well in the energy range E 8 GeV, but leads to a kink at E ∼ 8 GeV due to some terms in the expression changing sign at a = |δm 2 13 | cos 2θ 13 . The comparison at a shorter baseline length of L = 810 km, which is the distance from Fermilab to NOνA, is shown in figure 6. There, all five 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 [47], as can be seen in figure 7. There, the Cervera et al. and the Akhmedov et al. formulae greatly overestimate P (ν µ → ν e ), while the Asano-Minakata formula leads to negative probability for E ∼ 4 GeV. The Freund formula is accurate up until E ∼ 7 GeV where a kink occurs at a = |δm 2 13 | cos 2θ 13 . In comparison, our approximation remains accurate for all energies.
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 figures 8 and 9 for the two baselines L = 1300 km and L = 2300 km, respectively. These  3 . This has been estimated using the PREM profile of the Earth [45]. Note that the Asano-Minakata formula gives negative probability for E ∼ 4 GeV.
distances correspond to those between Fermilab and Homestake (1300 km), and CERN and Pyhäsalmi (2300 km) [48]. 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: From figure 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. 5 For the normal hierarchy case, when a ≈ δm 2 31 we have If we expand the running parameters in our eq. (4.4) in powers of the vacuum s13 and α = δm 2 21 /δm 2 31 , the leading order term expressed using the notation of Freund [22] takes the form . This is the same as eq. (38a) of Freund with sin 2 2θ13 replaced by 4s 2 13 , and agrees with the corresponding term of Akhmedov et al. [29]. The enhancement discussed in the main text can be seen to occur atÂ = 1, that is a = δm 2 31 , which is possible only when δm 2 31 > 0. However, the formulae of Cervera et al. [27], Akhmedov et al. [29], Asano and Minakata [32], and Freund [22] compared in the previous section all suffer in accuracy around the resonance regionÂ ≈ 1. This is not the case for our expression, which has a smooth transition across the resonance. The fact that the CP-phase dependent terms are negligible at the relevant energies and baselines is also clear in our approach. Therefore, . (4.7) From figure 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.

The "magic" baseline
The "magic" baseline is the baseline at which the dependence of P (ν e → ν µ ) on the CP violating phase δ vanishes [49]. 6 Looking at eq. (4.2), the only term without δ-dependence 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 figure 4 that δλ 21 ≈ a = 2 √ 2G F N e E, and the above condition reduces to 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 . This is satisfied for L ≈ 7500 km as can be read off of figure 1. Indeed, in figure 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. 6 An illuminating discussion on the physical meaning of the "magic baseline" can be found in ref. [50]. 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 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 figure 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 matter effects. Our approximation was derived utilizing the Jacobi method [44], 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 [33] and RENO [35]. Our formulae are compact and can easily be coded as well as be manipulated by hand. The application of our method to finding the ν e → ν µ , ν τ resonance conversion condition, and that to the determination of the 'magic' baseline [49,50] 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 -16 -

JHEP04(2014)047
neutrino flavor would add extra terms to the effective Hamiltonian, which would require further rotations for diagonalization. This has been discussed previously in ref. [43], and the potential constraints on New Physics from long baseline neutrino oscillations experiments in refs. [51][52][53]. Updates to these works will be presented in future publications.

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

JHEP04(2014)047
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 (A.14) Since ∆ 32 = ∆ 31 − ∆ 21 , (A.15) 8 Note that our notation differs from that of Cervera et al. in ref. [27]. 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. [22] where ∆ = δm 2 31 , and∆ = δm 2 31 L/4E. Huber and Winter in ref. [49] define ∆ = δm 2 31 L/4E, which is also used in ref. [58]. So care is necessary when comparing formulae.

JHEP04(2014)047
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 [59]: The oscillation probabilities for the anti-neutrinos are obtained by replacing U αi with its complex conjugate, which only amounts to flipping the sign of δ in the parametrization of eq. (A.4). It is clear from eq. (A.16) that P (ν α → ν α ) = P (ν α → ν α ), which is to be expected from the CPT theorem. For flavor changing oscillations, only the Jarskog term in eq. (A.17) changes sign.

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 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 [16]. 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 [42,43]  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 [44] 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.27)
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 (cf. eq. (B.17)) we obtain eq. (B.27). The β-dependence of the difference φ − φ is plotted in figure 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.

JHEP04(2014)047
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, where we have defined θ 13 ≡ θ 13 + φ . 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 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. 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 asymptotic values are thus λ + → δm 2 21 , λ − → −ac 2 13 c 2 12 , in the limit a → 0, and λ + → c 2 12 δm 2 21 , λ − → −ac 2 13 , in the limit a → ∞.

B.4.2 Second rotation
After the first rotation, the effective hamiltonian was given by eq. (B.42). When a < δm 2 21 , both non-zero off-diagonal elements are of order O(εa) < O(ε 3 |δm 2 31 |). In contrast to the neutrino case, as a increases beyond δm 2 21 , the angle θ 12 approaches 0, and it is the 1-3 element that becomes the larger of the two. Therefore, a 1-3 rotation is needed next.
We define where c φ = cos φ , s φ = cos φ , tan 2φ ≡ − ac 12 sin 2θ 13 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 The β-dependence of λ ± and φ are shown in figure 4 ((c) and (d)), and figure 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  The difference between φ and φ is shown in figure 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.

JHEP04(2014)047
Here, we argue that that is, the 1-3 rotation passes through R 12 (θ 12 , 0). This is due to the fact that φ only becomes non-negligible when a δm 2 12 where s 12 ≈ 0 and c 12 ≈ 1, which means thus any matrix will commute with R 12 (θ 12 , 0). 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.59) is trivially satisfied. The accuracy of this approximation is discussed in appendix C. Therefore, The 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 θ 12 → θ 12 = θ 12 + ϕ , θ 13 → θ 13 = θ 13 + φ , θ 23 → θ 23 , δ → δ , (B.64) 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.  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 figure 18. The dependence of the non-zero elements of δR on β = − log ε (a/|δm 2 31 |) is shown in figure 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 figure 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 figure 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 figure 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 figure 20, will lead to kinks in the resulting oscillation probabilities.

C.2 Anti-neutrino case
The difference between the two sides of eq. (B.59) 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 → 0, c 12 → 1, s φ → −s 13 (c 13 ), and c φ → c 13 (s 13 ) for normal(inverted) hierarchy. The question is whether δR will stay negligible in between as s 12 runs from s 12 to 0, c 12 from c 12 to 1, s φ from 0 to −s 13 (normal) or c 13 (inverted), and c φ from 1 to c 13 (normal) or s 13 (inverted) as shown in figure 21.
The dependence of the non-zero elements of δR on β = − log ε (a/|δm 2 31 |) is shown in figure 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 figure 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 figure 23, will lead to kinks in the resulting oscillation probabilities.   Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.