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 $\eta$-gauge neutrino mass-squared difference $\Delta^{}_* \equiv \eta \Delta^{}_{31} + (1-\eta) \Delta^{}_{32}$ is introduced, where $\Delta^{}_{ji} \equiv m^2_j - m^2_i$ for $ji = 21, 31, 32$ are the ordinary neutrino mass-squared differences and $0 \leq \eta \leq 1$ is a real and positive parameter. Expanding neutrino oscillation probabilities in terms of $\alpha \equiv \Delta^{}_{21}/\Delta^{}_*$, we demonstrate that the analytical formulas can be remarkably simplified for $\eta = \cos^2 \theta^{}_{12}$, with $\theta_{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 ${\cal O}(\alpha^2)$. 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.


Introduction
Thanks to a number of elegant neutrino oscillation experiments in the past few decades, it has been well established that neutrinos are massive and lepton flavors are mixed [1]. In the framework of three-flavor neutrino oscillations, the lepton flavor mixing is described by a 3 × 3 unitary matrix U, i.e., the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [2,3], which is conventionally parametrized in terms of three mixing angles {θ 12 , θ 13 , θ 23 where s θ ij ≡ sin θ ij and c θ ij ≡ cos θ ij for ij = 12, 13, 23 have been defined. 1 The global-fit analysis of solar, atmospheric, reactor and accelerator neutrino oscillation experiments [4,5,6] yields three mixing angles θ 12 ≈ 33 • , θ 13 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 1 If massive neutrinos are Majorana particles, two extra CP-violating phases are needed to parameterize the PMNS matrix, but they are irrelevant for neutrino oscillations.  and ∆ 31 , three mixing angles {θ 12 , θ 13 , θ 23 } and the CP-violating phase δ from a global fit of current experimental data [6]. 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 matter. In addition, the neutrino super-beam experiments ESSνSB (L ≈ 500 km and 0.2 GeV E 0.6 GeV) and MOMENT (L ≈ 150 km and 0.15 GeV E 0.20 GeV) have also been proposed to measure the CP-violating phase with relatively low energy neutrinos and short baseline lengths [16,17,18] .
• 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. 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 CPviolating 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.
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 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 ≡ 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] 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 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 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 where I denotes the 3 × 3 unit matrix and the relevant coefficients are , , 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 Ref. [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. (8) 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. (9) is flavor-independent and thus irrelevant for neutrino oscillations. In the formalism shown in Eqs. (6) and (7), 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 where the neutrino mass term, i.e., the first term on the right-hand side of Eq. (10), becomes diagonal. More explicitly, we have [51] M 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. (6) and (7). 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 [52,53,51], 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 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.
According to Eqs. (6) and (7), 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. (6) and (7). A further exploration of the right-hand side of Eq. (17) 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. (18) 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

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 Fig. 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.
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. (14) and (15), one should first expand x 2 − 3y and z in terms of α, and then insert their approximate expressions back into Eqs. (14) and (15). 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 θ 13 ) 1/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. (23) 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. (14) and (15).
Therefore, we define ǫ ≡ λ 2 − λ 1 and expand λ 3 up to the second order of α. Then, λ 1 and λ 2 can be obtained from the identity λ 1 + λ 2 = −(b + λ 3 ) and the definition of ǫ, namely, where the higher-order terms of O(α 3 ) have been omitted. Note that Eq. (24) is valid for both NMO and IMO. The corresponding coefficients ρ i (for i = 1, 2, 3) in Eq. (24) can be directly computed by making use of Eqs. (13), (21) and (22). More explicitly, we have 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. (26), 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. (13). Note that, instead of ǫ itself, ǫ 2 has been expanded in α in Eq. (26) 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.
Having obtained λ i 's, we can calculate ξ αβ i 's according to Eqs. (6) and (7). Their analytical expressions have been collected in Appendix A. After inserting ξ αβ i 's and λ i 's into Eqs. (19) and (20), we finally obtain the approximate formulas of oscillation probabilities 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. (56) 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). 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 Fig. 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

Special Case 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. (23), C depends on 2 Note that the spikes along the curves for ∆ P αβ in Fig. 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 masssquared 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 .  Table 1 have been adopted and a baseline of L = 6500 km is employed. The left panel is for NMO and the right panel is for IMO.
η implicitly through A and ∆ * , so its expression is not modified. The other one is where Eq. (26) with η = cos 2 θ 12 has been used.
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. (31).
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. (29), 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 and P µe ≃ s 2 2θ 13 s 2 which are much simpler than the general formulas in Eqs. (27) and (28). See also the results of P τ µ in Eq. (57) in Appendix B.
To carry out a systematic test of numerical accuracy of analytical approximations, we consider the absolute errors ∆ P αβ defined in Eq. (30) 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 Fig. 3, where the sizes of absolute errors are denoted by different colors. Some comments on the numerical calculations are in order: • In Fig. 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 Fig. 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.
• For IMO on the right column of Fig. 3, one can observe that the discrepancy is at most 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. (35) and (36), together with that of P τ µ in Eq. (57) 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.  Table 1 have been used. The absolute errors ∆ P αβ (for αβ = ee, µe, τ µ) have been defined in Eq. (30), and the probabilities are averaged over a Gaussian energy resolution of 1%.

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. (5) 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. (37) and P αβ with α = β = e in Eq. (20), one can immediately realize by identifying the oscillation terms of the same kind. In the derivation of Eq. (38), we have implemented the following relations 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.
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, from which additional relations between { θ 23 , J } and the parameters {ξ αβ i , ǫ}, similar to those in Eq. (38) 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 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. (41) and (42). In Appendix C we also list the mapping of the three mixing angles and the Jarlskog invariant for a generic value of η.
As a cross check, we further use the relations derived in Eqs. (41) and (42) It turns out that the expression in Eq. (57) 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   Table 1 have been used.
Moreover, given the approximate expressions of effective mixing parameters in Eqs. (41) and (42), one can insert them back into Eqs. (37), (40) and (43) and obtain a new set Figure 5: Numerical comparison between ∆ P ′ αβ and ∆ P αβ , where η = cos 2 θ 12 is fixed and the other input parameters are the same as in Fig. 2.
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  and cannot be expanded in terms of α in the low energy range with A α. Thus, we keep ǫ intact in the calculations.
• Our calculations employ a generic η-gauge neutrino mass-squared difference ∆ * and derive the η-gauge oscillation probabilities as shown in Eqs. (27), (28) and (56)  Moreover, as demonstrated in Fig. 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. (35), (36) and (57) 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 Fig. 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. (41) and (42) 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. (37), (40) and (43) for P ′ ee , P ′ µe and P ′ τ µ , respectively. The accuracy of the effective mixing parameters is proved in Fig. 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 Fig. 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 aforementioned criteria and can be readily applied to future oscillation experiments. We leave such applications for a separate work in the near future.