Why is the neutrino oscillation formula expanded in Δm212/Δm312 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 α = Δm212/Δm312 ≈ 0.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.

The formula was originally derived in [10,11] as a series expansion in α. But the problem is that due to the non-perturbative behavior near the solar resonance, the expansion is expected to be valid only when the neutrino energy is well above the solar resonance, conventional formula will be proven to be accurate below the bound (1.3). 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 [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. As a byproduct of our analysis, a new approximate formula is derived in this paper, with better accuracy. Though the exact form is a little more complicated than (1.1), for practical use in neutrino simulation, it is useful and covers most aspects. This is important considering that simulation of LBL experiments and performing χ 2 -fits require a fast and simple method to compute a large number of oscillation probabilities. Therefore even though the numerical calculation is always viable, there are still many studies on analytic approximation formulae for neutrino oscillation in matter [13,[30][31][32][33][34][35][36][37][38][39][40][41].
This paper is organized as follows. In section 2, we reformulate the original derivation of the formula (1.1) and numerically show the accuracy 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 very accurate. Then in section 3 from the viewpoint of singularities, we show that non-differentiable singularities in many parameters originate from the branch cut and result in the failure of the α-expansion. In section 4 we solve the problem rigorously and then compute the analytical error of (1.1). Based on the calculation in section 4, we also propose an alternative to the conventional formula. Their accuracies are numerically verified, which will be shown in section 5. Finally we conclude in section 6.

The α-expansion and the accidental accuracy
In this section, we first introduce analytic diagonalization of the 3 × 3 effective Hamiltonian, which has early been done by Zaglauer and Schwarzer [42] without any approximation. Then we show the α-expansion of the result from Freund's calculation [11] and compare the approximate result with the exact result (though 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 final result (i.e. the assembled oscillation probability) from these parameters is very accurate.
Neutrino oscillation in matter is subjected to the Schrödinger equation in the 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 in the standard neutrino oscillation framework, Here U and m i 's are neutrino mixing matrix and masses in vacuum respectively. The second term in eq. (2.2) comes from the matter effect. Without the second term (i.e. N e = 0), the JHEP10(2015)090 solution of (2.1) is quite simple since the first term has already been diagonalized. So in vacuum, the transition amplitude of ν α → ν β is just Here S αβ is usually referred to as the S-matrix in neutrino oscillation. For neutrino oscillation in matter, we need to diagonalize (2.2) to obtain the effective mixing matrixŨ and the effective neutrino massesm k , defined as ThenŨ andm k , combined in the way similar to (2.3), gives the S-matrix in matter. The 3 × 3 matrix H can be analytically diagonalized by solving first the eigenvalues and then the corresponding eigenvectors, though the computation is complicated.
The process can be a little simplified if we extract a dimensionless matrix M from where u ≡ (U e1 , U e2 , U e3 ) is the first row of U and is real by proper rephasing U . The cubic equation for the eigenvalues of M d is The eigenvalues of M d (Note that 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ỹ

JHEP10(2015)090
ThenŨ can be computed by solving the corresponding eigenvectors of λ k . The reader may refer to [42] for the full form of eigenvectors. OnceŨ is computed, we can extract effective mixing angles from the standard parametrization ofŨ . All effective parameters (masses and mixing angles) expanded in α are given below [11]: 14) The effective mixing angles in matter are (we use a superscript m to distinguish them from vacuum parameters) So one can replace the corresponding vacuum parameters in (2.21) with the effective parameters in matter given above. This will produce the conventional formula (1.1). Note that the above α-expansion of effective parameters requires not only α 1 but also α/A 1. If we look at the effective mixing angles in (2.17)-(2.20), we find that the α-expansion of sin 2θ m 12 may have a problem at α/A ∼ 1. In (2.18) we see sin 2θ m 12 ∼ α/A which implies the correction from α is amplified by 1/A, so the expansion may be not accurate if A is small. We compare it with the exact value from [42] in figure 1, where the energy range is 0.1 − 1.2GeV and matter density is 1.3g/cm 3 (the case of 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 (the gray region). This is because the unitarity ofŨ 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 that below 0.3 GeV the effective mass square difference ∆m 2 21 = m 2 3 (λ 2 − λ 1 ) from the exact solution (solid curves) can be several times that of the α-expansion (dashed curves), which also implies invalidity of the α-expansion at low energies.  But interestingly, 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 of P is because sin 2θ m 12 does not appear at the leading order (LO) of (1.1), but only at the next-to-leading order (NLO) and the next-to-next-to-leading order (NNLO) which are of order α 1 and α 2 , respectively. This suppresses 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.
It might be expected that the calculation at a higher order of α can explain this by finding some cancellations between errors. However, at a higher order, the accuracy in figures 1, 2 turns out to be improved very little. Actually, as will be discussed in the next section, there is an underlying problem that some variables in the calculation are not JHEP10(2015)090 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.

Non-differentiabilities, singularities and branch cuts in the oscillating system
To reveal the key problem in the expansion, we start from an oversimplified but heuristic problem, series expansion of the following function If α is small, but A is relatively not, then we can expand it in α, Here we see α/A 1 is necessary to make eq. (3.2) accurate. If A is much smaller than α, then we can expand it in A as g(α, A) = α + A 2 /(2α) + · · · . But what if A is close to α? One may think that if A is close to α, then both are small so we can make a double expansion of the function, where c 0 = g(0, 0), c 1 and c 2 are the partial derivatives ∂ α g and ∂ A g at α = A = 0. However, we cannot expand √ α 2 + A 2 in this way because the partial derivatives c 1 and c 2 do not exist (as one can check explicitly). Geometrically this is easy to understand since g = √ α 2 + A 2 is a cone in the α − A − g space. Expansion at the tip of the cone will certainly fail.  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.
Actually the function (3.1) does exist in the eigenvalues of the oscillation system 2 so in the original derivation of (1.1) α/A 1 has to be assumed. If we use formulae derived from such an expansion but simply ignore the condition α/A 1, then we are at the risk of getting total wrong results, such as sin 2θ m 12 ¿1 shown in figure 1. So basically the question is why for the oscillation probablity this problematic expansion works very well. This will be answered next by branch cut singularities.
First we 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 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 different directions at level crossings which implies the behavior near the resonance is quite non-perturbative.
As we suppress A 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 but not differentiable functions of α.
In the left plot of figure 5, we plot λ 1 and λ 2 as functions of α and A. It shows that there is a singularity (here we mean non-differentiable singularity) in the eigenvalues. The singularity is intrinsic and can not be removed by proper ordering of eigenvalues. So this implies that double expansion in α and A does not work.
The intrinsic singularity in figure 5 is the kernel of the non-perturbativity problem in the oscillation system. It actually comes from a branch cut singularity. 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 2 The eigenvalues (and thus the corresponding oscillation parameters) contain more complicated square roots like α 2 + A 2 c 2 13 − 2αAκ where κ 1/3 (see, e.g., calculation in [32]), but the problem caused by α ∼ A in the expansion is essentially the same as the simplified one in (3.1). further in terms of α, A and u according to (2.9). They look like smooth analytic functions everywhere but they are actually not. Note that there is a square root and a cubic root in (2.11). Functions like √ z or z 1/3 have branch cut singularities which make them not analytic 3 on the complex plane.

JHEP10(2015)090
In the right plot in figure 5 we show the two branches of ± √ z where z ≡ x + iy connecting with each other at the branch cut (the line y = 0 for x < 0). At x = 0 which is the end 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 the singularity in the eigenvalues shown in the left plot.

Solution
Identifying that the singularity in the eigenvalues originates from the branch cut singularity makes a 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, α/A will not appear any more because 1/A will be absorbed by some continuous and smooth functions which we will see below.

JHEP10(2015)090
After the singularity problem is solved, we will mathematically show the conventional formula is accurate near the solar resonance.
The solution of the singularity problem can be summarized by the three key steps below: 1. Use the Cayley-Hamilton theorem [43][44][45] to express the S-matrix only in terms of the eigenvalues λ 1,2,3 . 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, though the singularity still exists in λ − since it is a branch cut singularity.
3. It turns out that λ − only appears in cosine function and the following function − which is singularity-free, the singularity is removed.
In short, we will first make the S-matrix only depend on (λ 1 , λ 2 , λ 3 ) and then after the transformation from (λ 1 , λ 2 , λ 3 ) to (λ 2 − , λ + , λ 3 ), the singularity in the S-matrix will be explicitly removed. Next we will show the calculation in detail.
The  where we put a −it to be used later. The coefficient s 0 , s 1 and s 2 can be determined in various methods [43] such as the Lagrange interpolation or the Newton interpolation. They have been computed in [44,45], JHEP10(2015)090 From eq. (2.1) and (2.5), the S-matrix is where 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 an identity matrix, S eµ can be written as two terms which implies we do not need to compute s 0 for the appearance probability.
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, we have λ 1 + λ 2 = −b − λ 3 and λ 1 λ 2 = −dλ −1 3 . Therefore where b, d are apparently free from the singularity [as shown in eq. (2.9)] and λ 3 is also singularity-free (shown in figure 4, see also the proof in the appendix). So λ + and λ 2 − are singularity-free. Note that, however, λ − has a singularity originating from the branch cut singularity, which can be seen from figure 5.

JHEP10(2015)090
Now that we have expressed the S-matrix in terms of (λ 3 , λ + , λ 2 − ) and thus 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 we have not taken any approximation. Then we will use the approximation .

(4.21)
Since the p term and q term have been expressed in terms of singularity-free quantities λ 3 , λ + and λ 2 − , we can use (4.20), (4.21) to compute them. The calculation is straightforward (see the appendix) and the result is The oscillation probability is where P (A) is defined as To derive (1.1), we need to expand the modulus squared of (4.26), where |p term| 2 , |q term| 2 and the cross term are computed in the appendix and the result is (4.30) Now the conventional formula (1.1) can be analytically justified near and below the solar resonance. Here we denote it as P (B) . We can see that the first and last terms of P B just correspond to eqs. (4.28), (4.29) respectively and the middle term in P (B) corresponds to the cross term (4.30). Combine the analytic errors, we have which implies that the conventional formula is accurate up to the O-terms above and the O-terms in (4.25). We see there is no α/A in all these O's, so we draw the conclusion that the bound (1.3) which originally requires α/A 1 can be safely removed. The conventional formula is still accurate without this bound, as long as these O-terms are small.
According to eq. (4.25), the error of P (A) is δP (A) = O(s 4 13 A∆) + O(s 2 13 αA∆ 2 ). Taking T2K as an example, for E = 0.25GeV which is below the conventional domain of validity, we have A 0.02 and ∆ 3.7. This gives s 2 13 αA∆ 2 2 × 10 −4 and s 4 13 A∆ 4 × 10 −5 . The error is very small and the dominant correction would be O(s 2 13 αA∆ 2 ) if we want to improve the accuracy. Actually, if we only concern ourselves with the ∆ 1 region, then O(s 2 13 αA∆ 2 ) is larger than O(s 4 13 A∆) since s 2 13 α s 4 13 . Therefore we expect that typically O(s 2 13 αA∆ 2 ) is the principal source of the error of P (A) . For the same set of parameter values, the four terms in eq. (4.31) have the values, s 2 13 α∆ 3 × 10 −3 , s 13 α 2 ∆ 2 2 × 10 −3 and α 3 A∆ 4 α 4 ∆ 4 1 × 10 −4 . So the dominant errors of P (B) are O(s 2 13 α∆) and O(s 13 α 2 ∆ 2 ). Note that they do not depend on A, which implies that the main source of inaccuracy of P (B) is not due to inadequately accounting for the matter effect contribution, but rather than due to an insufficient expansion of the small phase α∆. Though the phase α∆ is small, terms quadratic in α∆ would not be negligible if we want to improve its accuracy. In conclusion, the dominant errors of P  Figure 6. A comparison of the approximate oscillation formulae with the numerical solution to (2.1) in the T2K case. The left plot shows that all these approximate formulae are very accurate and the errors are negligible for practical use. The right plot shows that the residuals defined as |δP | = |P −P numerical | are consistent with our analytic estimation of the error, which are represented by green and yellow shades for δP (A) and δP (B) respectively.

Numerical verification
As our study of the problem is originally motivated by the fact that T2K covers the solar resonance, we would like to numerically verify our analysis in that case first. The matter density in T2K is ρ = 2.6g/cm 3 [12] so we take the electron density to be N e = 1.3N A /cm 3 under the assumption that for matter Z/A = 1/2 in average. Figure 6 shows that both P (A) and P (B) are accurate enough for practical use while the new formula P (A) has better accuracy than the conventional formula P (B) . We also plot the analytic errors according to (4.25) and (4.31) in the right panel of figure 6, using light green and yellow shades. The actual residuals defined as |δP | = |P − P numerical | where P numerical is the numerical solution are well compatible with analytic estimation, which implies the errors are correctly estimated. Therefore figure 6 verifies both P (A,B) and δP (A,B) in the T2K case.
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 1 and for simplicity we take the same matter density as T2K for all the other experiments, since the neutrino beams in these experiments only go though the earth crust. We see again that in current or future accelerator neutrino experiments, the formulae are accurate enough for practical use and the errors are well described by our analytic estimation.
Although computers are becoming more powerful, it is still desirable to have efficient methods of computation. For example, the χ 2 -fit in a high dimensional parameter space (including both oscillation parameters and experiment parameters) is always extremely time-consuming. When nuisance parameters are being marginalized, the likelihood function has to be invoked an enormous amount of times to complete a sub-process of minimization (only for the frequentist treatment, the Bayesian approach usually needs much more computations). The package GLoBES [48,49]  neutrino oscillation experiments has optimized the diagonalization procedure [50] combining the QL decomposition with additional developed algorithms. This is not necessary for oscillations in constant density matter, where a simple analytical formula performs better (GLoBES allows users to replace the probability engine with a user-defined function). In this case, we recommend the use of P (A) instead. A simple test on Mathematica 8.0 with Intel Core i7 CPU shows that 10 5 evaluations of P (A) and P (B) cost 4 4.7 and 5.6 seconds, which implies P (A) can be computed at a speed not slower than P (B) . Finally there is one issue related to the solar resonance to be discussed. Strictly speaking, the matter effect can be safely regarded as a small perturbative effect only if 2 √ 2G F N e E is much less than ∆m 2 21 (A α), i.e. only the region between the vacuum limit and the solar resonance can be regarded as the truly small-A region, where no physics can be changed greatly by the matter effect. Typically LBL accelerator neutrino experiments are in the region between the solar resonance and the atmospheric resonance which we can refer to as medium-A region. Note that originally A can not be treated perturbatively in the medium-A region [11]. From the small-A region to the medium-A region, the solar mixing will experience a resonance. It is interesting that, according to the formulae we derived , the contribution from the matter effect passes through the resonance gradually without showing any resonances, despite the solar mixing being affected drastically in that region (note that the solar mixing has sizable contributions to these experiments). In other words, the region with a perturbative matter correction can be extended from the small-A region to the medium-A region for current LBL accelerator neutrino experiments.

Conclusion
The conventional formula obtained by an 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 though the effective masses and effective mixing angles computed in the α- 4 For compiled languages which are used in practical simulation such as the C-based GLoBES package, the speed will be about several hundred times faster. But the ratio of the speeds of computing P (A) and p (B) varies little for different machines or different languages.

JHEP10(2015)090
We have shown that the inaccuracies are because the expansion is too close to the branch cut singularity in the eigenvalues. This singularity is inherent in the eigenvalues so it cannot be removed by interchanging eigenvalues. But certain combinations of them such as their 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 proven that the conventional formula is still accurate near the solar resonance.
A new oscillation formula P (A) in (4.26) which might be practically useful is derived when we try to prove the accuracy of the conventional one. Both the conventional and the new formulae are very accurate in various accelerator neutrino experiments for baseline lengths varying from 150km (MOMENT) to 1300km (LBNE), as shown in figures 6 and 7. We have also estimated the analytic errors for these formulae.

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. (4.20). Note that the cubic equation (2.8) has the following identity b + c + d = As 2 13 (α − 1) − 1, (A.1) 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.2) 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,