Looking into analytical approximations for three-flavor neutrino oscillation probabilities in matter

Motivated by tremendous progress in neutrino oscillation experiments, we derive a new set of simple and compact formulas for three-flavor neutrino oscillation probabilities in matter of a constant density. A useful definition of the η-gauge neutrino mass-squared difference Δ∗ ≡ ηΔ31 + (1 − η)Δ32 is introduced, where Δji≡mj2 − mi2 for ji = 21, 31, 32 are the ordinary neutrino mass-squared differences and 0 ≤ η ≤ 1 is a real and positive parameter. Expanding neutrino oscillation probabilities in terms of α ≡ Δ21/Δ∗, we demonstrate that the analytical formulas can be remarkably simplified for η = cos2θ12, with θ12 being the solar mixing angle. As a by-product, the mapping from neutrino oscillation parameters in vacuum to their counterparts in matter is obtained at the order of Oα2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}^2\right) $$\end{document}. Finally, we show that our approximate formulas are not only valid for an arbitrary neutrino energy and any baseline length, but also still maintaining a high level of accuracy.


JHEP12(2016)109
Normal mass ordering (NMO) Inverted mass ordering (IMO) best-fit 3σ range best-fit 3σ range  Table 1. The best-fit values and 3σ ranges of two neutrino mass-squared differences ∆ 21 and ∆ 31 , three mixing angles {θ 12 , θ 13 , θ 23 } and the CP-violating phase δ from a global fit of current experimental data [6].  table 1 for a summary of the latest global-fit results from ref. [6]. Besides precision measurements of the known mixing parameters, the primary goals of ongoing and forthcoming oscillation experiments are to pin down neutrino mass ordering (i.e., the sign of ∆ 31 ), to measure the leptonic CP-violating phase δ, and to determine the octant of θ 23 (i.e., θ 23 < 45 • or θ 23 > 45 • ). In order to study the experimental sensitivities and better understand future experimental results, we should pay particular attention to the Mikheyev-Smirnov-Wolfenstein (MSW) matter effects on the propagation of neutrino beams in a medium [7,8]. Roughly speaking, current and future neutrino oscillation experiments can be categorized into three different types, in which terrestrial matter effects on neutrino oscillations always play an important role.
• Medium-baseline Reactor Neutrino Experiments. The reactor experiments with a baseline length L ≈ 50 km and a neutrino-beam energy E ≈ 4 MeV, such as JUNO [9] and RENO-50 [10], are sensitive to the oscillations driven by both ∆ 21 and ∆ 31 . Hence, they will be able to determine neutrino mass ordering and precisely measure oscillation parameters. It has been found [11] that the Earth matter effects for JUNO are as large as 1%, significantly affecting the determination of sin 2 θ 12 and ∆ 21 , whose precisions are estimated to be 0.54% and 0.24%, respectively [12].
• Long-baseline Accelerator Neutrino Experiments. For the long-basline accelerator experiments T2K [13] (L = 295 km and E ≈ 0.6 GeV), NOνA [14] (L = 810 km and E ≈ 2 GeV) and LBNF-DUNE [15] (L = 1300 km and E ≈ 3 GeV), it is the relative sign between the matter potential for electron neutrinos (or antineutrinos) and ∆ 31 that changes the oscillation probability of ν µ → ν e (or ν µ → ν e ), opening another possibility to pin down neutrino mass ordering. The difference between oscillation probabilities of neutrinos and those of antineutrinos implies leptonic CP violation, which however suffers from a contamination induced by the CP-asymmetric Earth
• Huge Atmospheric Neutrino Experiments. The experiments PINGU [19], ORCA [20], and Hyper-Kamiokande [21] will implement huge ice or water Cherenkov detectors to precisely measure atmospheric neutrinos, for which a wide range of energies (0.1 GeV E 100 GeV) and baseline lengths (10 km L 10 4 km) should be considered. Though neutrinos and antineutrinos cannot be distinguished in these experiments, the MSW resonance in the Earth matter occurs either in neutrino oscillations for NMO or in antineutrino oscillations for IMO. Therefore, matter effects help determine neutrino mass ordering. A 3σ significance can be reached at the ICAL detector of INO, which can also discriminate between ν µ and ν µ events [22].
In principle, for any neutrino energy and baseline length, one can exactly calculate neutrino and antineutrino oscillation probabilities in the Earth matter by numerical methods. However, it is obviously difficult in this way to reveal the underlying physics for neutrino oscillations and to fully understand the numerical results.
For this reason, two theoretical approaches have been suggested to study the Earth matter effects. First, one can establish an exact relation between the effective parameters in matter, i.e., the mixing matrix U (parametrized in terms of three mixing angles θ ij and one CP-violating phase δ) and three neutrino masses m i (or mass-squared differences ∆ ji ≡ m 2 j − m 2 i ), and the intrinsic parameters in vacuum. For example, J ∆ 21 ∆ 31 ∆ 32 = J ∆ 21 ∆ 31 ∆ 32 holds exactly for a constant matter density [23][24][25], where the Jarlskog invariant in vacuum [26,27] is defined by J γ,k αβγ ijk ≡ Im U αi U * αj U * βi U βj with the Greek and Latin letters running over (e, µ, τ ) and (1, 2, 3), respectively, and likewise for J and U . Another example is the Toshev relation sin 2 θ 23 sin δ = sin 2θ 23 sin δ [28] in the standard parametrization. In addition, the notion of unitarity triangles has also been introduced to describe leptonic CP violation [29][30][31][32], and the exact and approximate relations between the unitarity triangles in matter and those in vacuum have been found in refs. [33][34][35][36][37]. Although these exact relations are interesting in themselves, they are in practice not useful to directly explain experimental observations and extract fundamental oscillation parameters.
The second approach is to expand the oscillation probabilities in terms of small perturbation parameters, which can be α ≡ ∆ 21 /∆ * ≈ 0.03 [where ∆ * ≡ η∆ 31 + (1 − η)∆ 32 with 0 ≤ η ≤ 1, cf. eq. (2.6)] and the smallest mixing angle sin θ 13 ≈ 0.147 [38][39][40][41][42][43][44][45]. Another choice is A ≡ A/∆ * , where A is the matter potential from the coherent forward scattering of neutrinos on the background particles and defined as A ≡ 2 √ 2G F N e E, with G F being the Fermi constant, N e the number density of background electrons, and E the neutrino energy. In the case of low neutrino energies or low matter densities, an expansion in A is useful to show the corrections of matter potential to the oscillation probabilities in vacuum.
In the seminal paper by Freund [39], the analytical approximations for three-flavor neutrino oscillation probabilities have been systematically studied, and the formulas are JHEP12(2016)109 valid as long as the oscillation driven by ∆ 21 has not developed and the corresponding MSW resonance is not reached. The latter condition corresponds to A α [39], namely, where ρ is the matter density. For the Earth matter, the electron fraction is Y e ≈ 0.5 and N e = Y e N A [ρ/(1 g cm −3 )] with N A being the Avogadro's number. Although Freund's formulas actually work even for E < 0.34 GeV, it has been shown in ref. [41] and ref. [46] that the series expansion of ≡ (α 2 + A 2 cos 4 θ 13 − 2 Aα cos 2θ 12 cos 2 θ 13 ) 1/2 in terms of α is problematic in the region of low energies or small matter densities, where A → 0. More accurate approximate formulas for low energies E < 1 GeV have been derived in ref. [46] by retaining . However, the analytical results in refs. [41,46] are not applicable for large matter effects and higher neutrino energies. Furthermore, a critical problem for the sin θ 13 expansion is related to the atmospheric resonance A → 1, where the function C ≡ [(1 − A) 2 + 4 A sin 2 θ 13 ] 1/2 cannot be expanded correctly. As we will show later, and C are two key parameters to avoid any difficulties associated with the low-energy solar resonance and the high-energy atmospheric resonance, respectively. In fact, analytical formulas for arbitrary neutrino energies and baseline lengths are derived in refs. [42,43], where the resonances related to ∆ 21 and ∆ 31 have been treated carefully by introducing a few intermediate rotation angles for basis transformations. Thus, the analytical results can be cast into a simple and compact form, in which the eigenvalues of the zeroth-order Hamiltonian and rotation angles, instead of intrinsic mixing parameters, are involved. Since all the existing analytical approximations are not fully satisfactory, we are well motivated to derive a new set of analytical formulas for neutrino oscillation probabilities, which fulfills the following three criteria: 1. They are valid for arbitrary neutrino energies and any baseline length. Such formulas are applicable to atmospheric neutrino experiments.
2. They are expressed in terms of intrinsic oscillation parameters, and in a simple and compact form. Any complicated formulas are not very useful in practice.
3. They give accurate values of oscillation probabilities, under the condition that the first two criteria are met at the same time.
For this purpose, we expand the oscillation probabilities in terms of α, but retain the parameter that corresponds to ( C) in the case of low (high) energies or small (large) matter densities. In addition, an η-gauge neutrino mass-squared difference ∆ * ≡ η∆ 31 + (1 − η)∆ 32 is introduced so as to seek an optimal value of η that greatly simplifies approximate formulas. The remaining part of this paper is organized as follows. In section 2, we briefly review the basic strategy to derive analytical formulas of neutrino oscillation probabilities. We introduce ∆ * ≡ η∆ 31 +(1−η)∆ 32 and demonstrate that it is suggestive of simple analytical formulas for η = cos 2 θ 12 . The oscillation probabilities in the special case of η = cos 2 θ 12 are presented in section 3, where the mapping between effective and intrinsic mixing parameters is also obtained as a by-product. The accuracies of the analytical formulas are examined JHEP12(2016)109 and compared with previous ones. Finally, we summarize our main results in section 4. Some useful formulas are listed in three appendices.

General formalism
In the framework of three-flavor neutrino oscillations, the effective Hamiltonian responsible for the evolution of neutrino flavor eigenstates in matter is given by (2.1) In the case of a constant matter density, i.e., a constant value of A, we then have two distinct ways to derive the exact oscillation probabilities. First, one can diagonalize the effective Hamiltonian by using a unitary transformation where m i for i = 1, 2, 3 are effective neutrino masses in matter, and U is the effective PMNS matrix, which can also be parametrized in terms of three mixing angles { θ 12 , θ 13 , θ 23 } and one CP-violating phase δ. In terms of these effective parameters, it is straightforward to write down the oscillation probabilities P αβ ≡ P (ν α → ν β ) as follows have been defined in the same manner as for neutrino oscillations in vacuum, and L is the baseline length. The probabilities for antineutrino oscillations ν α → ν β can be obtained by replacing J → − J in eq. (2.3) and A → −A everywhere in the effective parameters. Second, according to the Cayley-Hamilton theorem, the evolution matrix S = e −i H f L of neutrino flavor eigenstates is determined by three eigenvalues of the effective Hamiltonian H f and the matrix elements of H f [47][48][49], namely, where I denotes the 3 × 3 unit matrix and the relevant coefficients are ,

JHEP12(2016)109
with ω i ≡ m 2 i /(2E) being the eigenvalues of H f . The neutrino oscillation probabilities are simply given by P αβ = |S βα | 2 , while the results for antineutrino oscillations can be derived by changing U → U * and A → −A in the effective Hamiltonian.

η-gauge mass-squared difference
In order to simplify the analytical formulas as much as possible, we tentatively introduce a generic definition of neutrino mass-squared difference where 0 ≤ η ≤ 1 is a real and positive parameter. It is evident that ∆ * reduces to the conventional definitions of atmospheric neutrino mass-squared differnce ∆ 32 for η = 0, ∆ 31 for η = 1, and (∆ 31 + ∆ 32 )/2 for η = 1/2. In the global-fit analysis of neutrino oscillation data, the first two definitions have been used in refs. [4,5] in the IMO and NMO cases, respectively, while the last one has been implemented in ref. [6] for either neutrino mass ordering. Another definition ∆ ee ≡ cos 2 θ 12 ∆ 31 + sin 2 θ 12 ∆ 32 , corresponding to η = cos 2 θ 12 , has been advocated by Parke [50] and demonstrated to be advantageous to reactor antineutrino experiments. Although for quite a different reason, as we will show later, the introduction of ∆ * in eq. (2.6) with η = cos 2 θ 12 turns out to be very useful in simplifying the approximate formulas of oscillation probabilities.
With the help of ∆ * , the effective Hamiltonian H f can be rewritten as where I is the identity matrix of rank three and with α ≡ ∆ 21 /∆ * and A = A/∆ * . Note that the first term on the right-hand side of eq. (2.7) is flavor-independent and thus irrelevant for neutrino oscillations. In the formalism shown in eqs. (2.4) and (2.5), the evolution matrix is now S = e −i∆ * M f L/(2E) and only the eigenvalues of M f need to be calculated. To find out the eigenvalues of M f , it is more convenient to convert into the mass basis in vacuum via M v = U † M f U , where the neutrino mass term, i.e., the first term on the right-hand side of eq. (2.8), becomes diagonal. More explicitly, we have [51] where U ei for i = 1, 2, 3 are three elements in the first row of the PMNS matrix. Therefore, it is expected that a proper choice of η will be helpful in reducing the complexity of three eigenvalues, and thus the final oscillation probabilities. Furthermore, it is interesting to notice that only the matrix elements U ei for i = 1, 2, 3 are involved in M v and |U e3 | 1, so a suitable value of η is anticipated to be mainly associated with U e1 and U e2 , or θ 12 in the standard parametrization.

η-gauge oscillation probabilities
Now it is time to derive the oscillation probabilities by using eqs. (2.4) and (2.5). First of all, the eigenvalues λ i (for i = 1, 2, 3) of M f or equivalently M v can be obtained by solving the following eigen-equation where the relevant coefficients are The eigenvalues of the effective Hamiltonian have been known for a long time [51][52][53], but it has recently been noticed that the results depend also on neutrino mass ordering [46]. To be explicit, taking λ 1 < λ 2 < λ 3 , we have for the NMO; or for the IMO, where we have defined (2.14) Note that λ i 's are the eigenvalues of M f , and (λ 1 , λ 2 , λ 3 ) correspond to ( m 2 1 , m 2 2 , m 2 3 ) in the NMO case with ∆ * > 0, but to ( m 2 2 , m 2 1 , m 2 3 ) in the IMO case with ∆ * < 0. In the latter case, though λ 3 is the largest eigenvalue, ∆ * λ 3 becomes negative and thus m 2 3 is the smallest one. In addition, it is easy to verify that λ 2 − λ 1 > 0 holds for either neutrino mass ordering.

JHEP12(2016)109
According to eqs. (2.4) and (2.5), it is straightforward to compute the evolution matrix S = e −2iF * M f with F * ≡ ∆ * L/(4E) and thus the oscillation probabilities where the flavor-dependent coefficients ξ αβ i (for i = 1, 2, 3) with α and β running over e, µ and τ can readily be identified from similar equations for M f to those for H f in eqs. (2.4) and (2.5). A further exploration of the right-hand side of eq. (2.15) gives rise to where the identity λ 1 + λ 2 = −(b + λ 3 ) has been implemented. Although we will not show the exact expressions of ξ αβ i 's, some useful properties of them can be implemented to further simplify the oscillation probabilities and the series expansions of ξ αβ i 's with respect to the small parameter α have been collected in appendix A.
In the appearance channel ν α → ν β with α = β, the identity ξ αβ 1 = −ξ αβ 2 holds exactly. Therefore, it is easy to verify that Im[ξ αβ * . Under these conditions, eq. (2.16) will be reduced to where one can observe four different types of oscillation terms. In the disappearance channel ν α → ν α , we have ξ αβ 1 + ξ αβ 2 = 1 and ξ αβ * i = ξ αβ i , and thus arrive at in which only three oscillation terms survive. In order for the oscillation probabilities to be valid for arbitrary neutrino energies and baseline lengths, as we shall see later, it is important to always keep those oscillation terms not expanded at all.

Analytical and numerical results
So far, all the analytical results in the previous section are exact. In this section, we will first expand the eigenvalues λ i 's in terms of α and derive the approximate formulas of neutrino oscillation probabilities in the general η-gauge. Simple and compact formulas in the special case of η = cos 2 θ 12 then emerge in an obvious way. As a by-product, the mapping between effective and fundamental mixing parameters is also obtained. Finally, numerical verifications are carried out to show high precisions of our analytical formulas, in comparison with the exact ones.

Approximate formulas
Let us begin with the series expansion of three eigenvalues. First, to clearly see the relative sizes of λ i 's, we have shown their exact values as functions of the neutrino energy E in figure 1, where the matter density ρ = 2.8 g cm −3 , the electron fraction Y e ≈ 0.5 and η = 1 have been taken for illustration. In addition, the best-fit values of neutrino oscillation parameters from table 1 are adopted. In the left panel, the results for the NMO are given, where one can observe a potential level crossing at E ∼ 0.3 GeV for the solar resonance, and another one around E ∼ 10 GeV for the atmospheric resonance, if a poor approximation to λ i is adopted. In the right panel, since both λ 1 and λ 2 in the IMO case are actually negative, their absolute values are shown together with λ 3 . It is obvious that there is no level crossing in this case between λ 3 and λ 2 , but the level crossing at E ∼ 0.3 GeV for the solar resonance still exists. For antineutrino oscillations in matter, as is well known, the atmospheric resonance will be present in the IMO case, while absent in the NMO case. A correct treatment of these eigenvalues in the regions of resonances is crucial to get well-behaved analytical results.

JHEP12(2016)109
In this work we shall use α as the only expansion parameter and deal carefully with the would-be divergences in the neighborhood of resonances and in the limiting cases (e.g., the vacuum oscillations with A → 0). Looking at the analytical results of λ i 's in eqs. (2.12) and (2.13), one should first expand x 2 − 3y and z in terms of α, and then insert their approximate expressions back into eqs. (2.12) and (2.13). After a straightforward but tedious calculation, we finally get and where c φ ≡ cos φ and s φ ≡ sin φ have been introduced also for φ = 2θ ij . In addition, we have defined a regulator for the atmospheric resonance [39] and C ≡ ( C 2 + Ac 2 Note that C appears in the denominators and will cause divergences in the further expansions in terms of sin 2 θ 13 when A = 1. Therefore, we shall keep the exact form of C in eq. (3.3) in our calculations of the oscillation probabilities.
On the other hand, in the low-energy or vacuum limit with A → 0, we have learned from refs. [41] and [46] that one cannot expand the function mentioned in section 1 in terms of α. A further study shows that this function arises from the difference between two eigenvalues λ 2 and λ 1 , namely, the terms proportional to √ 1 − z 2 in eqs. (2.12) and (2.13). Therefore, we define ≡ λ 2 − λ 1 and expand λ 3 up to the second order of α. Then, λ 1 and JHEP12(2016)109 λ 2 can be obtained from the identity λ 1 + λ 2 = −(b + λ 3 ) and the definition of , namely, where one can clearly observe that the above coefficients will be greatly simplified for η = cos 2 θ 12 . In particular, ρ 2 = 0 implies that the first order correction to λ 3 is vanishing, so the leading-order results are already very precise. Additionally, at the second order, only one term is left in λ 3 . However, this is not the case for λ 1 and λ 2 , as extra contributions come from , which can be determined from In the derivation of eq. (3.6), the identities λ 1 + λ 2 = −(b + λ 3 ) and λ 1 λ 2 = −dλ −1 3 have been used, where both b and d have been given in eq. (2.11). Note that, instead of itself, 2 has been expanded in α in eq. (3.6) where high-order terms of O(α 3 ) have been neglected. As we will show in the next subsection, reduces to in the case of η = cos 2 θ 12 and in the limit of A → 0. Therefore, is the parameter that we should retain in the series expansion. where The expansion of F + is given to O(α), but a few terms proportional to α 2 are kept in the coefficients in front of oscillation terms, as they may become important in some cases. For completeness, we also present the complete expression for P τ µ in eq. (B.1) in appendix B. As is proved in ref. [40], only two oscillation probabilities are independent, say, P µe and P τ µ . The other probabilities can be constructed by making use of unitarity condition α P αβ = β P αβ = 1 and the time-reversal transformation P αβ = P βα (δ → −δ) for a constant matter density. Furthermore, considering that the rotation matrix in the 2-3 sector commutes with the matter potential term in the effective Hamiltonian, we can establish the relations P eτ = P eµ (θ 23 → θ 23 + π/2) and P µµ = P τ τ (θ 23 → θ 23 + π/2).

JHEP12(2016)109
For this reason, only two independent appearance probabilities P µe and P τ µ are shown in this work, while P ee is given as an example in the disappearance channel. Note that the oscillation probabilities are for NMO, and we can get the corresponding results by replacing with − for IMO. Given the above approximate formulas for oscillation probabilities, we also verify that these expansions indeed reduce to those that already exist in the literature. For example, in the low energy range, α, A and are of the same order and can be expanded simultaneously, from which one can arrive at eq. (4.6) in ref. [46]. On the other hand, for the high energy region with A ∼ α, one can safely expand in terms of α and restore the familiar results of Freund [39] and Akhmedov et al. [40].
As we have mentioned, the analytical expressions can be substantially simplified when η = cos 2 θ 12 is adopted, where all the terms proportional to (η − cos 2 θ 12 ) automatically disappear. The resultant simplified formulas will be presented and discussed in the next subsection. Here we show that the choice of η = cos 2 θ 12 is not only advantageous for the analytical simplicity, but also for numerical accuracy. To examine the influence of η on the accuracy of analytical approximations in the oscillation probabilities, we have adopted the best-fit values of neutrino oscillation parameters in table 1. In addition, the matter density ρ ≈ 2.8 g cm −3 for the Earth's crust and Y e ≈ 0.5 are taken for illustration. In order to test the numerical accuracy, we define the absolute error of the analytical approximations of P (ν α → ν β ) as ∆ P αβ for α, β = e, µ, τ , i.e., where ( P αβ ) Exact is calculated by a fully numerical evolution of the neutrino flavor states. Note that an unusual baseline of L = 6500 km is employed in order to make the fine structure of oscillations more prominent. The oscillation probabilities and their absolute errors are given in figure 2, where we can observe that the case of η = cos 2 θ 12 is the most accurate one for almost the entire range of neutrino energies. 2 Comparing with previous analytical approximations of the oscillation probabilities, our results are advantageous in several aspects. First, we have included all the possible leading terms of the whole energy region. Taking the expansion terms α 2 / 2 and α for instance, although α 2 / 2 is a higher-order term than α near the atmospheric resonance, it is significantly enhanced in the low energy range where is small. Thus both are maintained in the expansion. Second, our analytical results keep and C as independent parameters in order to avoid any divergence in the low-energy limit and near the atmospheric resonance, respectively. Third, for the first time, we have presented the analytical results with a generic η value, which is convenient to make a comparison with previous results. We further show that η = cos 2 θ 12 is the best choice in terms of both simplicity and numerical accuracy. 3 2 Note that the spikes along the curves for ∆ P αβ in figure 2 do not mean the best precision but the intersection points of exact and approximate oscillation probabilities, which are caused by the modifications of oscillation frequency and amplitude in the approximate formulas. 3 Although we demonstrate that η = cos 2 θ 12 leads to simpler and more accurate oscillation probabilities, the underlying physical reason is not clear and deserves further studies [54]. We notice that the same mass-squared difference ∆m 2 ee ≡ cos 2 θ 12 ∆ 31 + sin 2 θ 12 ∆ 32 has been shown in ref. [50] to be advantageous for ν e disappearance experiments without matter effects. This observation may provide a clue to better understand the choice of η = cos 2 θ 12 . If η = cos 2 θ 12 is fixed, we can obtain much simpler formulas for relevant oscillation parameters in matter and those for the oscillation probabilities as well. First, let us focus on the two regulators for eliminating possible divergences. As indicated in eq. (3.3), C depends on η implicitly through A and ∆ * , so its expression is not modified. The other one is
In the low-energy limit, A will also be a small parameter, just like α. In this case, it is easy to verify where higher-order terms of O(α 2 A) are omitted. It has been found in ref. [46] that one can keep in the oscillation probabilities, whose low-energy behaviors will then be remarkably improved. In the high-energy limit, it is safe to expand in terms of α, and thus we get Note that in this case is not a small parameter, as the matter effects become important or even dominant, e.g., A 1. For an arbitrary neutrino energy, it is necessary to make use of the full result of in eq. (3.11).
Second, as it is useful to define the oscillation phases F − ≡ ∆ 21 L/(4E) and F + ≡ (∆ 31 + ∆ 32 )L/(4E) in vacuum, we obtain the following analytical approximations of their counterparts in matter with the help of eq. (3.9), namely, reflecting the corrections induced by the Earth matter to neutrino mass-squared differences.
Given the above parameters, the oscillation probabilities for the special case of η = c 2 θ 12 turn out to be which are much simpler than the general formulas in eqs. (3.7) and (3.8). See also the results of P τ µ in eq. (B.2) in appendix B.
To carry out a systematic test of numerical accuracy of analytical approximations, we consider the absolute errors ∆ P αβ defined in eq. (3.10) and the approximate results are now obtained by using the simplified formulas in the case of η = cos 2 θ 12 . The numerical results of ∆ P αβ for a wide range of neutrino energies and baseline lengths have been shown in figure 3, where the sizes of absolute errors are denoted by different colors. Some comments on the numerical calculations are in order: • In figure 3, the matter density of ρ ≈ 2.8 g cm −3 with Y e ≈ 0.5 and the best-fit values of neutrino oscillation parameters from table 1 have been used in numerical calculations. In addition, to avoid fast oscillations at low energies, we have averaged the oscillation probabilities over a Gaussian energy resolution of 1%. The baseline lengths and neutrino energies have been set to be 0.1 km ≤ L ≤ 10 4 km and 1 MeV ≤ E ≤ 100 GeV, respectively. Hence, both current and future oscillation experiments as mentioned in the introduction are essentially covered. As for the atmospheric neutrinos, our assumption of a constant matter density renders it impossible to reveal the structure of parametric resonances [55][56][57][58][59]. However, it suffices to illustrate the numerical difference between our analytical formulas and the exact oscillation probabilities.
• In the lower part of each plot in figure 3, i.e., for L ≤ 1 km, one can observe that the errors are always far below the level of 10 −8 . This can be understood by noticing that the oscillations driven by ∆ 21 have not yet developed for a short baseline. The ∆ 31 -driven oscillations indeed take place for short baseline lengths and low neutrino energies, however, the amplitudes will be suppressed by the smallest mixing angle θ 13 . For higher neutrino energies, we need longer baseline lengths for the ∆ 31 -driven oscillations to develop. The errors in the entire range of baseline lengths and energies are below 10 −3 , demonstrating an excellent agreement between our approximate formulas and the exact ones.  , where the matter density of ρ ≈ 2.8 g cm −3 with the electron fraction Y e ≈ 0.5 and the best-fit values of neutrino oscillation parameters in table 1 have been used. The absolute errors ∆ P αβ (for αβ = ee, µe, τ µ) have been defined in eq. (3.10), and the probabilities are averaged over a Gaussian energy resolution of 1%.

JHEP12(2016)109
• For IMO on the right column of figure 3, one can observe that the discrepancy is at most 10 −5 ∼ 10 −4 , as a consequence of the absence of resonances in this case. For NMO, the region of largest errors always appears around E ≈ 10 GeV and L ≈ 5000 km, where the atmospheric resonance is encountered, while around the region of the solar resonance relatively smaller errors are observed. Such a difference on the size of error at two different resonances may be attributed to the fact that λ 2 and λ 3 are more close to each other at the atmospheric resonance than λ 1 and λ 2 at the solar resonance.
Notice that the approximate formulas of P ee and P µe in eqs. (3.15) and (3.16), together with that of P τ µ in eq. (B.2) in appendix B, are the main results of this work. Given their simplicity and high level of numerical accuracy, one may directly employ them to perform both analytical and numerical studies on neutrino oscillation phenomena in current and upcoming oscillation experiments. We leave such applications for a future work.

Parameter mapping
With those newly obtained approximate formulas for oscillation probabilities, a more accurate mapping of the intrinsic mixing parameters to the effective mixing parameters in matter can actually be established as a by-product. To see this clearly, we first re-express the exact formulas of neutrino oscillation probabilities in eq. (2.3) in terms of effective parameters in matter. Starting with the disappearance channel ν e → ν e , we have where c θ ij ≡ cos θ ij and s θ ij ≡ sin θ ij have been defined as before. Comparing between P ee in eq. (3.17) and P αβ with α = β = e in eq. (2.18), one can immediately realize by identifying the oscillation terms of the same kind. In the derivation of eq. (3.18), we have implemented the following relations 19) which are verified by using λ i = [ m 2 i − (m 2 2 − η∆ 21 )]/∆ * and ∆ ji ≡ m 2 j − m 2 i . Thus far, the results are exact and no approximations have been made. To get more useful results for effective mixing angles θ 12 and θ 13 , we have to expand ξ ee i (for i = 1, 2, 3) with respect to α but keep unchanged as before.

JHEP12(2016)109
Since P ee is independent of θ 23 and the CP-violating phase δ, one should further consider the oscillations in the appearance channels. As an example, we study the oscillation probability in the appearance channel ν µ → ν e , namely, 20) from which additional relations between { θ 23 , J } and the parameters {ξ αβ i , }, similar to those in eq. (3.18) can be found. Using the series expansions of ξ αβ i listed in appendix A and setting η = cos 2 θ 12 , we finally arrive at the mapping for three mixing angles and that for the Jarlskog invariant J 2αJ Note that the leptonic CP violation is now described by the Jarlskog invariant, and the direct relation between δ and the vacuum mixing parameters can be easily obtained using eqs. (3.21) and (3.22). In appendix C we also list the mapping of the three mixing angles and the Jarlskog invariant for a generic value of η. 4E It turns out that the expression in eq. (B.2) can be exactly reproduced, when both of them are matched to the same order of α.
It is worth mentioning that the mapping relations for mixing angles and the Jarlskog invariant have been truncated at the second order of α and serve as excellent approximations to the exact results. For illustration, we have calculated the effective mixing angles {sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 } and the effective Jarlskog invariant J for different neutrino energies. As depicted in figure 4, the exact results are denoted as solid curves (red), while the approximate results based on eqs. (3.21) and (3.22) are represented by dashed curves (blue). One can see that our approximate results are in perfect agreement with the exact ones, and the differences between them are invisible from the plots. For comparison, the numerical results according to the mapping relations found by Freund in ref. [39] are given as dotted curves (green). Significant deviations can be observed in the figures for sin 2 θ 12 and J, which can be explained by the divergence encountered in the low-energy region.
Moreover, given the approximate expressions of effective mixing parameters in eqs. (3.21) and (3.22), one can insert them back into eqs. (3.17), (3.20) and (3.23) and obtain a new set of oscillation probabilities, which we call P ee , P µe and P τ µ , respectively. As the effective mixing parameters are expanded up to O(α 2 ), these new oscillation probabilities will be more accurate in the sense that part of higher-order terms are now included. To illustrate this point, we compute the absolute errors ∆ P αβ according to eq. (3.10) and compare it with ∆ P αβ from figure 2 in the case of η = cos 2 θ 12 . The results are shown in figure 5, where one can find ∆ P αβ (blue solid curves) are almost always one or two orders of magnitude smaller than ∆ P αβ (red dashed curves).

Summary
In this work we have taken a deep look into analytical approximations for three-flavor neutrino oscillation probabilities in matter of a constant density and presented a new set of simple and compact formulas. A useful definition of the η-gauge neutrino mass-squared difference ∆ * ≡ η∆ 31   expansions of α (i.e., α ≡ ∆ 21 /∆ * ). The approximate oscillation probabilities are valid for arbitrary neutrino energies and any baseline length. Among different choices of η, it turns out that the case of η = cos 2 θ 12 is the best one in terms of both simplicity and numerical accuracy. These formulas are particularly useful for the future long baseline accelerator JHEP12(2016)109 P ee , P µe and P τ µ respectively. Given the expressions of ρ i (i = 1, 2, 3) in eq. (3.5) and 2 in eq. (3.6), the analytical results of oscillation probabilities are greatly simplified for η = cos 2 θ 12 , where all the terms proportional to (η − cos 2 θ 12 ) automatically disappear. Moreover, as demonstrated in figure 2 for different values of η, the choice of η = cos 2 θ 12 is the most accurate one for almost the entire range of neutrino energies.
• Fixing the gauge at η = cos 2 θ 12 , the oscillation probabilities are presented in eqs. (3.15), (3.16) and (B.2) for P ee , P µe and P τ µ respectively, constituting the main results of this work. Regarding the accuracy of these analytical approximations, a careful study is performed in figure 3 for the neutrino energies from 10 −3 GeV to 10 2 GeV and the baseline length range 10 −1 km ≤ L ≤ 10 4 km. One can observe that in the NMO case the errors in the entire range of baseline lengths and neutrino energies are below 10 −3 , while for IMO below 10 −4 . The largest errors appear in NMO around E ≈ 10 GeV and L ≈ 5000 km, where the atmospheric resonance is encountered and the small energy splitting between λ 2 and λ 3 slows down the convergence of the series expansions.
• As a by-product a more accurate mapping of the intrinsic mixing parameters to the effective mixing parameters in matter is established in eqs. (3.21) and (3.22) for three mixing angles and the Jarlskog invariant, respectively. With the effective mixing parameters, one can obtain a new set of oscillation probabilities in eqs. (3.17), (3.20) and (3.23) for P ee , P µe and P τ µ , respectively. The accuracy of the effective mixing parameters is proved in figure 4 for the whole energy range including the regions of the solar and atmospheric resonances. For the new set of oscillation probabilities, one can find from figure 5 that the accuracy of P αβ will be one or two orders of magnitude better than P αβ because some higher-order terms are also properly included.
• Finally, in the low energy range, α, A and are of the same order and can be expanded simultaneously, from which one can arrive at eq. (4.6) in ref. [46]. On the other hand, for the high energy region with A ∼ α, one can safely expand in terms of α and restore the familiar results of Freund [39] and Akhmedov et al. [40] For future long-baseline accelerator and atmospheric neutrino experiments, with the goals of determining the neutrino mass ordering and measuring the leptonic CP violating phase, a set of compact and simple analytical approximations of oscillation probabilities in matter is very helpful. These analytical oscillation probabilities should be directly connected with the fundamental oscillation parameters and be valid for arbitrary neutrino energies and any baseline length. In this sense, our analytical approximations in this work meet all the afore-mentioned criteria and can be readily applied to future oscillation experiments. We leave such applications for a separate work in the near future.
A Expressions for the ξ αβ i terms In this appendix, we present expressions for the ξ αβ i terms, which are coefficients in front of various oscillation terms in eq. (2.15). For this purpose, we employ the Cayley-Hamilton theorem to express the evolution matrix S = e −2iF * M f into a form similar to eq. (2.4) with

JHEP12(2016)109
H f replaced by M f . Correspondingly, ω i in eq. (2.5) are now the eigenvalues of M f , i.e., λ i , and e −iω i L read as e −2iF * λ i . Then, with the explicit form of M f given in eq. (2.8), we are able to obtain the expressions of various ξ αβ i for various oscillation channels P αβ , according to the definitions of ξ αβ i in eq. (2.15). As shown in eq. (2.16), the final oscillation probabilities P αβ only depend on certain combinations of ξ αβ i , we therefore just show the analytical expansions for those relevant ones. For P ee , we have ξ ee 1 + ξ ee 2 = 1, ξ ee i = ξ ee * i and ξ ee For P µe , we have ξ µe 1 = −ξ µe 2 and (A.12)