Why Is The Neutrino Oscillation Formula Expanded In $\Delta m_{21}^{2}/\Delta m_{31}^{2}$ Still Accurate Near The Solar Resonance In Matter?

The conventional approximate formula for neutrino oscillation in matter which is obtained from the expansion in terms of the ratio of mass square differences $\alpha=\Delta m_{21}^{2}/\Delta m_{31}^{2}\approx0.03$, first proposed by Cervera, et al and Freund, turns out to be an accurate formula for accelerator neutrino experiments. Originally it required the neutrino energy to be well above the solar resonance to validate the expansion but it is found to be still very accurate when the formula is extrapolated to the resonance, which is practically important for the T2K experiment. This paper shows that the accuracy is guaranteed by cancellations of branch cut singularities and also, for the first time, analytically computes the actual error of the formula. The actual error implies that the original requirement can be safely removed in current experiments.


Introduction
The hierarchy of neutrino mass square differences ∆m 2 21 ∆m 2 31 simplifies 3-neutrino oscillation to 2-neutrino-dominated oscillation in many practical cases [1]. For example in short-baseline(SBL) reactor neutrino experiments such as Daya Bay [2] and RENO [3], we usually simply assume the oscillation generated by ∆m 2 21 is negligible. For current longbaseline(LBL) accelerator neutrino experiments such as T2K [4,5], MINOS [6] and NOvA [7,8], the oscillation is still ∆m 2 31 -dominated since the ratio L/E of the baseline length and neutrino energy is of order O(1/∆m 2 31 ). If the approximation is not good enough, we can expand the neutrino oscillation probability in terms of the mass hierarchy parameter α defined as α ≡ ∆m 2 21 /∆m 2 31 ≈ 0.03. (1.1) For example, ν µ → ν e oscillation can be expanded in α as P νµ→νe ∼ = P 0 + P cos δ + P sin δ + P 3 (1.2) where P 0 , P cos δ + P sin δ and P 3 are proportional to α 0 , α 1 and α 2 respectively. In vacuum with small sin θ 13 approximation they are where ∆ ≡ ∆m 2 31 L/(4E) and J CP = s δ c 2 13 s 13 c 12 s 12 c 23 s 23 is the Jarlskog invariant. All other parameters are defined in the same convention of PDG [9], including the short denotations such as c 13 ≡ cos θ 13 and s δ ≡ sin δ. The leading order(LO) P 0 enables ν e appearance experiments at T2K and MINOS to measure θ 13 (θ 23 in this term can be measured by ν µ disappearance [4] from the same neutrino beams.). P cos δ and P sin δ are CP even and odd terms including the Dirac CP phase δ. Due to sin θ 13 suppression, actually the leading order(LO) term is only about 3 times the NLO term P cos δ + P sin δ when δ = ±π/2 . Therefore measurement of θ 13 at T2K and MINOS largely depends on the CP phase δ and this in turn implies these experiments have sensitivities to δ if θ 13 is fixed from SBL reactor neutrino experiments such as Daya Bay and RENO. Actually this is how T2K excluded δ between 0.19π and 0.80π at 90% C.L.(i.e. δ = π/2 is disfavored) in their recent paper [5].
Although the α−expansion in eq.(1.3) is useful in discussing physical implications such as the correlation of θ 13 and δ in T2K, it does not include the matter effect [10][11][12] which will contribute roughly 10% correction in T2K and 20% in MINOS. So one may expect a similar α-expansion including the matter effect. However, when the matter effect comes in, the α−expansion will be confronted with a non-trivial problem, i.e. the non-perturbativity near the solar resonance 1 . Anyway, traditionally we do have a similar formula in matter which can be found in PDG [9] and was first obtained in [13,14],  (1.5) and N e is the electron number density in matter, about 1.4cm −3 N A in the earth crust. But, the problem is, due to the non-perturbative behavior near the solar resonance, the expansion including matter effect is expected valid only well above the solar resonance, i.e. the validity of the formula requires a lower bound on the neutrino energy, E 0.34GeV ∆m 2 21 7.6 × 10 −5 eV 2 (1. 6) This was emphasized in ref. [14] ( see also comments on eqs. (14.45)- (14.49) in PDG [9]), because the approximation α/A 1 was used when the formula was derived. We will reformulate the traditional derivation of the formula in section 2 to show the problem more explicitly but here we take the solar mixing angle θ 12 as a good example to show the problem. The effective sin 2θ 12 in matter, denoted as sin 2θ m 12 , expanded by α to the first order, is about [14] sin 2θ m The solar resonance is at A = α cos 2θ 12 ≈ 0.4α so near the solar resonance sin 2θ m 12 ∼ α/A is quite likely to be larger than 1. As will be shown in section 2, sin 2θ m 12 > 1 does appear when the energy is lower than 0.34GeV, which makes the calculation totally invalid. Originally, sin 2θ m 12 in the calculation was expected not only less than 1 but also small, i.e. sin 2θ m 12 1, otherwise the unitarity of the effective mixing matrix would be badly violated. So we should not expect the traditional oscillation formula to be accurate below the bound (1.6) or even near the bound. The formula should only be used by experiments with energies well above the bound.
For MINOS experiment [6], most data were taken at 1-5GeV and the spectrum peaked at 3 GeV, so there is no serious problem. But for T2K, the energy range is 0.1-1.2 GeV and the spectrum peaks at 0.6 GeV [5]. So at the spectrum peak of T2K, eq.(1.4) might be very inaccurate. Furthermore in a part of the current measured range 0.1-0.34 GeV, the calculation might be totally wrong. Recently there was a proposal of a new accelerator neutrino experiment called MOMENT [15] with an even lower neutrino energy beam (200MeV) and a shorter baseline(150 km) which makes the problem more serious if the formula is applied to MOMENT.
However, T2K actually has used this formula in their recent paper [16] because eq.(1.4) presents unexpectedly good accuracy near the solar resonance(for a systematic analysis of accuracies, see [1]). And for the MOMENT experiment, one can check numerically that this formula shows a better accuracy than that in T2K.
Therefore some questions are raised: 1. Why is the traditional approximate formula in matter so accurate near the solar resonance?
2. Is there any underlying mechanism such as some cancellations at a higher order, to guarantee the accuracy of the final result which is assembled from those very inaccurate pieces?
3. If (1.6) is not a real bound of its validity, then to what extent is the formula accurate or valid?
In this paper, these questions will be answered. We will analytically prove that the formula is still accurate near and below the solar resonance as long as O(s 2 13 α)+O(s 13 α 2 )+O(α 3 A)+ O(α 4 ) + O(s 4 13 A) is small, which answers both the first and the third questions. It also implies that the lower bound (1.6) actually can be safely removed and thus theoretically justifies the use of this formula in T2K. As for the second question, a naive higher order calculation still has to face the same problem that a lot of α/A will appear and the unitarity is badly violated. Actually we will show that due to a singularity 2 in the eigenvalues, even infinite order calculation may fail.
The way to rigorously solve this problem involves a deeper understanding of the relation between branch cuts and the well known level crossing [17,18] in the matter effect. The former is less noticed but essentially related to level crossings. Note that the three eigenvalues of the oscillating system come from the same cubic equation but they are different. The differences originate from the different branches in the square roots and cubic roots in the general solutions of a cubic equation. At level crossings two of the eigenvalues are very close to each other which makes the problem quite non-perturbative and this just corresponds to the starting point of the branch cuts, which are called branch cut singularities. The branch cut singularities are the origins of all non-perturbativities in the oscillating system. We will remove the singularity corresponding to the solar resonance in our analytic calculation by transformation of the eigenvalues to some singularity-free variables and compute the S-matrix using the Cayley-Hamilton theorem. In this way the traditional formula will be proved accurate below the bound (1.6). The relation between the branch cut singularities and level crossings will be discussed in detail and thus improve our understanding of the matter effect in neutrino oscillation [19][20][21][22][23][24][25][26][27][28][29][30][31][32].
As a benefit of the analytic proof, we will propose some alternative formulae to the traditional one, with simple forms but better accuracies. This makes a practical sense considering that technically, simulation on the LBL accelerator neutrino oscillations and the χ 2 -fit on data desires a fast and simple calculation of the oscillation probability, due to the highly repetitive evaluation. 3 Besides, a simple formula is also useful in looking for the physical implications, for example, the correlation of θ 13 and δ as we just discussed and the parameter degeneracy problem [36][37][38]. Therefore even though the numerical calculation is always viable, there are still many studies on simple and accurate formulae for neutrino oscillation in matter [1,22,[39][40][41][42][43][44][45][46][47][48].
This paper is organized as follow. In section 2, we reformulate the traditional derivation of the formula and numerically show the accuracies of the α−expansion in the case of T2K. We will see that the α−expansion for some effective parameters is actually invalid below 0.34GeV in T2K while the final result of the probability is unexpectedly accurate. Then in section 3 from the viewpoint of singularities, we show there are non-differentiable singularities in many parameters which originate from the branch cut singularity and result in the failure of the α−expansion. In section 4 we solve the problem rigorously and 2 In this paper, the term singularity does not necessarily mean a point where a function explodes such as 1/x at x = 0, but also somewhere a function is discontinuous or continuous but not differentiable. Most singularities in this paper are branch cut singularities or non-differentiable singularities. For more information, see, for e.g. http://en.wikipedia.org/wiki/Singularity_(mathematics). 3 Actually the GLoBES package (designed for simulation of neutrino oscillation experiments) [33,34] numerically rediagonalizes the effective Hamiltonian using some algorithms including the famous QL decomposition and further developed algorithms [35] to improve the speed of computing the probabilities. But in T2K and many current accelerator neutrino experiments, neutrinos only go though the crust of the earth so the matter density is assumed constant. Then a simple analytic formula with a good accuracy is more practical.
then compute the analytical error of the traditional formula. Based on the calculation in section 4, we also propose some alternative formulae to the traditional formulae. All the approximate formulae and their accuracies are summarized in section 5.
2 The α-expansion and its problem near the solar resonance In this section, we first introduce the analytic calculation of effective neutrino masses and mixing angles from the rediagonalization of the 3 × 3 effective Hamiltonian, which has early been done by Zaglauer and Schwarzer [49] without any approximation. Then we show the α-expansion result of these parameters from Freund's calculation [14] and compare the approximate result with the exact result (complicated but numerically programmable) to see how much it deviates from the exact result. We will show that the α-expansion result of effective neutrino parameters are quite inaccurate and even invalid near the solar resonance but the combined result (i.e. the finally assembled oscillation probability) from these parameters is very accurate. The 3-neutrino oscillation in matter can be described by the Schrödinger equation in the 3-dimensional flavor space, where |ν(L) denotes the flavor state of the evolving neutrino at a distance of L from the source and H is the Hamiltonian represented by the 3 × 3 matrix Here the first part of (2.2) is already diagonalized by the PMNS matrix U and m i 's are neutrino masses. The second part comes from the matter effect [10][11][12]. Without the second part (i.e. N e = 0), since the first part has already been diagonalized, the solution of (2.1) is quite simple. The transition amplitude in vacuum of ν α → ν β is just where ∆ ≡ ∆m 2 31 L/(4E) and α ≡ ∆m 2 21 /∆m 2 31 and S αβ is usually called the S-matrix in neutrino oscillation. Considering the matter effect, we have to rediagonalize (2.2) to obtain the effective mixing matrixŨ and the effective neutrino massesm k , defined as ThenŨ andm k produce the S-matrix in matter, in a way similar to (2.3).
The 3×3 matrix H can be analytically diagonalized by first solving the eigenvalues (just solving a cubic equation) and then solving the corresponding eigenvectors, despite that the computation is complicated. This has early been done in [49], without any approximation.

Diagonalization of H can be a little simplified if we extract a dimensionless matrix M from
Then M d has a simpler form for diagonalization, where u ≡ (U e1 , U e2 , U e3 ) is a row vector of U and is real by proper rephasing U . The cubic equation to solve the eigenvalues of M d is The eigenvalues of M d ( also of M since M has the same eigenvalues as M d ) solved from eq.(2.8) are where k = 0, 1, 2 and Then the effective neutrino masses defined in eq.(2.4) are given bỹ which can be expressed explicitly in terms of α and A according to eqs.(2.9), (2.10) and (2.11). ThenŨ can be computed by solving the corresponding eigenvectors of λ k . But fully writing down the complicated result makes little sense here. The reader may refer to [49] for it. However all the effective parameters (masses and mixing angles) expanded in terms of α have relatively simpler forms. The α-expansion of the eigenvalues in eq.(2.10) are [14]: The effective mixing angles in matter are(we use a superscript m to distinguish them from vacuum parameters) Note that the above α−expansion of effective parameters requires not only α 1 but also α/A 1. Using the above expansion and an additional approximation θ 13 1, one can replace the corresponding vacuum parameters in (1.3) with these effective parameters in matter, to obtain the traditional formula [14].
The formula (1.4) was originally expected to be used well above the solar resonance which is the case of the neutrino factory with a baseline of thousands of kilometers and a neutrino beam higher than several GeV. In such a case, the α-expansion is of course valid. However it has been found to be accidentally accurate near and even below the solar resonance(α/A 1) . Furthermore in the limit A → 0, the formula in matter (1.4) reduces to the vacuum limit [13] which makes people guess that it may be accurate near the solar resonance. Despite the fact that the formula is numerically accurate and has a correct vacuum limit, a lower bound (1.6) is still set both in Freund's paper [14] and PDG [9], maybe due to the lack of theoretical justifications of the validity at this region rigorously.
If we have a close look at each parameter, we first find that the α−expansion of sin 2θ m 12 may have problems at α/A ∼ 1. In (2.18) we see sin 2θ m 12 ∼ α/A which implies the O(α) correction here is amplified by 1/A, so the expansion is not accurate. We compare it with the exact value from [49] in figure 1, where the energy range is 0.1 − 1.2GeV and matter density is 1.3g/cm 3 , all same as the T2K experiment.
From figure 1 we see the expansion formula fits the exact solution well only at E > 0.5 GeV but deviates from it quickly when E < 0.5GeV. More seriously, when the energy goes below 0.35 GeV then sin 2θ m 12 will be larger than 1 which is an unphysical value(the grey region). This is because the unitarity of the diagonalization matrix is badly violated.
The other parameters suffering from the same problem are λ 1 and λ 2 . We plot them with the exact solutions in figure 2. We see the effective mass square difference ∆m 2 21 = m 2 3 (λ 2 − λ 1 ) from the exact solution (solid curves) can be several times the difference from the α-expansion(dashed curves) below 0.3 GeV, which implies the invalidity of the α-expansion of the effective masses at low energies.  But, to our surprise, the formula of oscillation probability assembled from these inaccurate and even invalid pieces is very accurate, as shown in figure 3 where we use the same energy range and matter density as figure 1 and figure 2.
One argument might be that, the accuracy is because sin 2θ m 12 does not appear at the LO (P 0 ) but only at the NLO(P cδ + P sδ ) and NNLO(P 3 ), which suppress the effect of the inaccuracy from sin 2θ m 12 shown in figure 1. But in figure 3 we see that at the second and third peaks(from right to left), the NLO and NNLO are comparable to the LO so the accuracy can not be explained by the NLO suppression. For example, at the third peak, sin 2θ m 12 computed from the α−expansion is more than two times the true value in figure 1, so one would expect that P cδ + P sδ ∝ sin 2θ m 12 deviates from the true value by a factor of 2 and P 3 ∝ sin 2 2θ m 12 by a factor of 4, so P will deviates from the true value by a difference of 2(P cδ + P sδ ) + 4P 3 − (P cδ + P sδ + P 3 ) ∼ 0.12, much larger than the actual error∼ 0.01.  Figure 3. This figure shows that the traditional formula (P 0 + P cδ + P sδ + P 3 , blue solid curve) is very accurate when used in T2K, in contrast with figure 1 and figure 2 where those effective parameters used to compute the probability are very inaccurate. The three components are also plotted (dashed curves) to show that every term in the traditional formula (1.4) is indispensable for it to be accurate in T2K.
It might be expected that the calculation at a higher order of α can explain this by finding some cancellations between errors. However, one will find that, at a higher order, there are still some α/A terms so it improves very little of the accuracy. Actually, as will be discussed in the next section, there is an underlying problem called non-differentiability that some variables in the calculation are not differentiable at α = A = 0. For these variables, the expansion series including α/A can not even converge if α/A 1. That is why a higher order calculation can not solve the accuracy problem.
3 Non-differentiabilities, singularities and branch cuts in the oscillating system The expansion with respect to α has a more serious problem which has not been noticed before. That is, the eigenvalues are not differentiable at some points which are called singularities and will destroy the mathematical prerequisite for series expansion. First let us take a look at the functions λ k (α) which are defined by the exact solution of eigenvalues (2.10) and vary with α, as shown in figure 4. Note that we consider λ's as functions of α rather than the more common picture where they vary with A (or E), since they are expanded with respect to α.
The left plot in figure 4 shows that the eigenvalues can be very close to another at level crossings(corresponding to resonances), but they never really go across another. They turn to a different direction at a level crossing which implies the behavior near a resonance is quite non-perturbative.
As we suppress A more and more close to zero, the curvatures at those turning points in figure 4 become larger and larger. Finally the curvatures go to infinite when A → 0, which makes the curves turn suddenly at some points, then some singularities emerge. The right plot in figure 4 shows the A → 0 limit. In this limit, the eigenvalues are continuous  Figure 4. The three eigenvalues given by eq.(2.10) as a function of α. In the left plot, the matter effect parameter is A = 0.1, corresponding to E = 1.2GeV. In the right plot we take the limit A → 0.
but not differentiable functions of α. Note that the Taylor expansion of a function f (x) with respect to x requires that f (x) should be smooth, i.e. infinitely differentiable.
Here is a natural question that might be asked. Why not redefine λ 1 and λ 2 so that they are interchanged when they pass the level crossing? It seems the problem can be solved in this way since the eigenvalues will not change their directions drastically when passing the level crossing, if the interchange is made.
But this is not feasible, explained as follow, 1. This can be done only when A is exactly zero. Once when neutrinos get into matter, there are always gaps between those curves in figure 4, no matter how close they are. Therefore the curves do not really cross if A = 0 and there is no actual crossing point.
2. If we impose a point to interchange the two eigenvalues, then take a look at figure 2 which covers the main energy region we are interested in. It clearly shows that a naive interchange of the eigenvalues at any point only makes it worse.
Strictly speaking, the α-expansion in matter will not suffer from the non-differentiablility since in matter A always deviates from 0 and the eigenvalues are smooth functions of α as long as A = 0. But in practical calculation, if A is small, then the expansion is close to non-differentiable which makes the expansion series converge very slowly.
To make the problem more clear, we raise another question. As both α and A in T2K are very small (A ∼ α), why not doubly expand the formula in both parameters? The fact is, even though there are some double expansions of neutrino oscillation in matter [1] such as the α − θ 13 double expansion, there is no such approximate formula based on α − A double expansion. The oscillation probability P can be regarded as a function of α and A. The function P (α, A) is actually a continuous (this is obvious by intuition) and differentiable(this will be shown later) function of α and A. Thus in principle it can be expanded as P = P (0, 0) + c 1 α + c 2 A + O(α, A) where c 1 and c 2 are two constant, equal to∂P/∂α and ∂P/∂A at (0, 0). But in actual calculation the double expansion is not directly feasible because the eigenvalues can not be doubly expanded in α and A, due to the non-differentiability. That is why there is no such approximate formula so far. (But in section 4 this can be done.) In figure 5, we plot that the eigenvalues λ 1 and λ 2 as functions of α and A. It shows there is a singularity (here we mean non-differentiable singularity) and the singularity is intrinsic and can not be removed by proper ordering of eigenvalues. So this also answers both the two questions, i.e. why not interchange λ 1 and λ 2 after level crossing and why not doubly expand the P (α, A).
The intrinsic singularity in figure 5 is the kernel of the non-perturbativity problem in the oscillating system. Where is it from? From eqs. (2.10) and (2.11) we see that λ 1,2,3 can be analytically expressed in terms of b, c and d and then further in terms of α, A and u according to (2.9). There is no place for a singularity to appear except for the square root and cubic root in (2.11). We remind the reader that for complex functions, √ z or z 1/3 have singularities since they are not analytic functions 4 on the whole complex plane.
Typically, the functions √ z and z 1/3 are defined by taking the principal branch(for e.g √ i = √ e iπ/2 = e iπ/4 , not √ e −i3π/2 = e −3iπ/4 ) and it connects with the other branches at the branch cuts. A full picture of square root containing all possible branches is shown in figure 6. The two branches( √ z and − √ z) connect with each other at the branch cut(the line y = √ z = 0 with x < 0) and at the starting point of the branch cut there is a singularity. As one can check from (2.10) and (2.11) , the branch cut singularity just corresponds to

Solution
Identifying that the singularity in the eigenvalues originates from the branch cut singularity makes the crucial step to solve the problem, since all branch cut singularities can be easily removed if the multi-branches collapse to one. For example, the branch cut singularity in √ z will disappear when it is squared, i.e. ( √ z) 2 makes the two branches collapse and is singularity-free. Once the singularity is removed, the α/A appears as a nuisance in the traditional calculation will not appear here any more because 1/A will be absorbed by a continuous and smooth function which we will see soon. After the singularity problem is solved, we prove the accuracy of the the traditional formula at the end of this section.
The solution of the singularity problem can be summarized by the three key steps below: 1. Using the Cayley-Hamilton theorem [50][51][52] to express the S-matrix only in terms of the eigenvalues λ 1,2,3 in section 2. The eigenvectors are not needed.
2. The singularity still exists in the eigenvalues but can be partially removed by the transformation where it will be shown that λ + and λ 2 − are singularity-free but the singularity still exists in λ − since it is a branch cut singularity.
3. It turns out λ − only appears in cosine function and the following function where x is proportional to λ − . Note that f (x) can be well defined at x = 0. Since f (x) = 1 − x 2 6 + x 4 120 + ... and cos x = 1 − x 2 2 + x 4 24 + ... are actually functions of x 2 , so the singularity in λ − is removed due to that λ 2 − does not have the singularity.
The Cayley-Hamilton theorem is a famous theorem in linear algebra which states that if p(λ) is the characteristic polynomial of a matrix A (for example the left-handed side of eq.(2.8)), then "substituting" the matrix A for λ in this polynomial results in the zero matrix, i.e. p(A) = 0. Take the example of eq.(2.8), this means M 3 + bM 2 + cM + d = 0 or where we put a −it here to be used later. The coefficient s 0 , s 1 and s 2 can be determined in various methods [50] such as the Lagrange interpolation or the Newton interpolation. They have been computed in [51,52], where From eq.(2.1) and (2.5), the S-matrix is So we can identify t = 2∆ to use eq.(4.4) directly. Now the transition amplitude of ν µ → ν e is S eµ = s 0 I 12 +s 1 M 12 +s 2 (M 2 ) 12 but because I is the identity matrix, S eµ can be written as two terms which implies we do not need to compute s 0 for the appearance probability Although the denominator ∆ λ in s 1 and s 2 is possible to be zero, but we will show next that s 1 and s 2 are not divergent at ∆ λ = 0 and smooth (differentiable) everywhere. For example the singularity from λ 1 − λ 2 = 0 can be removed by the transformation λ ± = 1 2 (λ 1 ± λ 2 ). This singularity is the only one that confronts us in the energy range of current experiments.
According to Vieta's formulas for a cubic equation, where b, d are obviously free from the singularity as shown by eq.(2.9) and λ 3 is also singularity-free according to figure 4 (see also the proof in the appendix). So λ + and λ 2 − are singularity-free. Note that, however, λ − has a singularity here(continuous but nondifferentiable), originated from the branch cut singularity. After the transformation λ ± = 1 2 (λ 1 ± λ 2 ), s 1 and s 2 are (4.14) We see they depend only on λ 3 , λ + , λ − and d which are all continuous and smooth functions except for λ − . But since λ − only appears in cos(λ − t) and f (λ − t) which are actually functions of λ 2 − ( note that cos(x) = 1 − x 2 2 + x 4 24 + ... and f (x) = 1 − x 2 6 + x 4 120 + ...), we come to the conclusion that s 1 and s 2 are continuous and smooth functions of α and A. Now that we have removed the singularity, expansion in α will not suffer from any problems. The α/A appears in section 2 will not appear any more if we use eqs.(4.14-4.15) to compute the probability.
Define We call the two terms in eq.(4.17) as p term and q term respectively. From eqs.(4.14-4.15) we have So far all the calculations have been done without any approximation. We simplify the p term and q term in the appendix with some terms dropped into O(s 3 13 A), O(s 13 αA) and O(αs 2 13 A). The result is The neutrino oscillation probability in matter for ν µ → ν e is just where P (A) is defined as which implies that the traditional formula is accurate up to O(s 2 13 α)+O(s 13 α 2 )+O(α 3 A)+ O(α 4 ) and the O-terms in (4.23). We see there is no α/A in all these O's and the middle procedure, so we draw the conclusion that the bound (1.6) which requires α/A 1 can be safely removed. The traditional formula is still accurate without this bound, as long as  Table 1. A summary of all approximate formulae for neutrino oscillation in matter derived in this paper and their analytic errors. For numerical errors, see figure 7 and 8. Note that P (C) is the traditional formula proposed by Cervera, et al [13] and Freund [14], not a new formula from this paper.

Other approximate oscillation formulae and their accuracies
When taking A → 0 in eq.(4.24) we get the vacuum limit of P (A) Since all the O-terms in eq.(4.23) vanish in vacuum, eq.(5.2) should be equal to the exact vacuum formula without any approximation. One of our concerns is that how large is the correction that the matter effect contributes to the neutrino oscillation probability. In the appendix we calculate the difference of the probabilities in matter and in vacuum, which is where P mat+ is the matter enhancement on the probability, then we have another approximate formula, named as P (B) , P (B) = P vac + P mat+ (5.5) which has the form of the vacuum probability plus the correction from the matter effect. Similar to the principle of obtaining P (C) from P (A) , we have the α-expansion of P (B) , named as P (D) ,  Figure 7. Compare the approximate oscillation formulae with the numerical result in the T2K case, where P (C) is the traditional formula proposed by Cervera, et al [13] and Freund [14]. The upper plot shows that all these approximate formulae are very accurate and the errors are negligible for practical use. The lower plots shows the errors defined as |δP | = |P − P numerical |. , as calculated in the appendix. So far we have derived four approximate formulae, named as P (A) , P (B) , P (C) and P (D) . We make a summary of them in table.1 as well as their errors, denoted as δP .
All these formulae have different advantages. P (A) is suitable for fast and accurate numerical computation and P (B) shows explicitly the correction from the matter effect. P (C) and P (D) are useful in understanding the correlation of neutrino parameters because it shows clearly the leading order which contains θ 13 and θ 23 and the NLO with CP even and CP odd contributions.
As our study of the problem is originally motivated by the use of the traditional formula in T2K, we also concern about the accuracies of the other formulae in T2K. We plot all these formulae in figure 7. The matter density in T2K is ρ = 2.6g/cm 3 [16] so we have the electron density N e = 1.3N A /cm 3 if assuming Z/A (Z: atomic number, A: mass number) of the matter is 1/2 in average. Figure 7 shows that all of them are accurate enough for practical use while the formula P (A) which is the first formula we obtained in the above calculation has the best accuracy. The accuracy of P (B) is next to P (A) while the traditional one P (C) proposed by Refs. [13,14] and P (D) have much larger errors. The reason may be that both P (A) and P (B) have correct vacuum limits without approximation while the vacuum limits of P (C) and P (D) are approximate to the exact vacuum oscillation.
Besides T2K, we also show the accuracies of these formulae in other accelerator neutrino experiments. The information of the baselines and neutrino energies are listed in table 2 and for simplicity we take the same density as T2K for all the other experiments, since the neutrino beams in these experiments only go though the earth crust.

Conclusion
The traditional formula obtained by expansion in the mass hierarchy parameter α = ∆m 2 21 /∆m 2 31 ≈ 0.03 turns out to be very accurate near the solar resonance, as shown in figure 3 and figure 7 while the effective masses and effective mixing angles computed in the α−expansion are very inaccurate or even invalid at this region, as shown in figure 1 and figure 2. So it is interesting that the inaccuracies in the middle procedure cancel each other out at the end.
We show that the inaccuracies are because the expansion is too close to the nondifferentiable singularity in the eigenvalues which originates from the branch cut. The non-differentiable singularity is intrinsic in the eigenvalues so it cannot be removed by interchanging eigenvalues. But some combinations of them such as the sum of the eigenvalues λ 1 + λ 2 do not have the singularity and the oscillation probability only depends on these singularity-free combinations. By computing the probability in this way, we have analytically proved that the traditional formula is still accurate near the solar resonance as long as O(s 2 13 is small. We also propose some other approximate formulae which have simple forms, as summarized in table 1. All these formulae are practically useful as they have both simple forms and good accuracies in various accelerator neutrino experiments where the baseline length varies from 150km (MOMENT) to 1300km (LBNE), as shown in figures 7 and 8.

A Some details of analytic calculations
A.1 Simplify the p term and q term Here we show how to simplify the p term and q term step by step. All the approximations in the calculation should be analytically treated so we use O() instead of ≈. The first result we will derive is eq.(A.1) can be easily understood from figure 4 where we see the largest eigenvalue is almost a constant near the solar resonance. Note that the cubic equation (2.8) has the following identity b + c + d = As 2 which provides a fast way to compute λ 3 as follow. We assume λ 3 = 1 + x with x 1 and replace the λ in the cubic equation (2.8) with 1 + x. Then the leading order vanishes and the next-to-leading order(NLO) gives 3x + x(2b + c) + O(x 2 ) = As 2 13 (1 − α) (A.3) which implies x = O(s 2 13 A) while the explicit form of x is not important here. Therefore from eq.(4.13) we have so the q term can be greatly simplified. whereĀ ≡ (A + α) 2 − 4Aαc 2 12 c 2 13 andĀ 2 = 4λ 2 − .