Simple and compact expressions for neutrino oscillation probabilities in matter

We reformulate perturbation theory for neutrino oscillations in matter with an expansion parameter related to the ratio of the solar to the atmospheric Δm2 scales. Unlike previous works, we use a renormalized basis in which certain first-order effects are taken into account in the zeroth-order Hamiltonian. We show that the new framework has an exceptional feature that leads to the neutrino oscillation probability in matter with the same structure as in vacuum to first order in the expansion parameter. It facilitates immediate physical interpretation of the formulas, and makes the expressions for the neutrino oscillation probabilities extremely simple and compact. We find, for example, that the νe disappearance probability at this order is of a simple two-flavor form with an appropriately identified mixing angle and Δm2. More generally, all the oscillation probabilities can be written in the universal form with the channel-discrimination coefficient of 0, ± 1 or simple functions of θ23. Despite their simple forms they include all order effects of θ13 and all order effects of the matter potential, to first order in our expansion parameter.


Introduction
Neutrino oscillation based on the standard three-flavor scheme provides the best possible theoretical framework available to date to describe most of the experimental results obtained in the atmospheric, solar, reactor, and the accelerator neutrino experiments. Although numerically calculated neutrino oscillation probabilities suffice to analyze experimental data, understanding of the framework, in particular the one in matter [1], has not yet reached a sufficient level, in our opinion. Under the assumption of uniform matter density distribution, the exact expressions of the eigenvalues, mixing angles, and the oscillation probabilities in matter have been obtained [2][3][4]. Yet, the results for these quantities are generally too complicated to facilitate understanding of the structure of the three flavor neutrino oscillations in matter primarily due to the complexities of solving the cubic eigenvalue characteristic equation. For a recent comprehensive treatment of neutrino oscillation in the matter, see ref. [5].
Analytic approaches to the neutrino oscillation phenomenon, so far, are mostly based on variety of perturbative frameworks. If the matter effect is small one can treat it as a small perturbation [6]. In the environments in which the matter effect is comparable to the vacuum mixing effect, the only available small expansion parameter known to us is the ratio of the solar-scale ∆m 2 to the atmospheric-scale ∆m 2 ⊕ , ∆m 2 /∆m 2 ⊕ 0.03. sin θ 13 has been often used as an expansion parameter (there are enormous number of references, see e.g., [7]), but it is now known that its value is not so small, sin θ 13 0.15, which is of the order of ∆m 2 /∆m 2 ⊕ . Moreover, expansion around sin θ 13 = 0 misses the physics of the resonance which exists at an energy around E ∼ 10 GeV for earth densities. Therefore, it appears that the suitable perturbative framework is the one with the unique expansion parameter ∆m 2 /∆m 2 ⊕ . This framework was indeed examined in the past, to our knowledge in refs. [7][8][9][10].
In this paper, we present a new framework of perturbative treatment of neutrino oscillation in matter. We follow the reasoning stated above which led to identification of the unique expansion parameter ≈ ∆m 2 /∆m 2 ⊕ . But, unlike the preceding works, we use a "renormalized basis" as the basis of perturbation. That is, we absorb certain terms of order to our "zeroth-order" Hamiltonian around which we perturb. Or, in other word, we take the zeroth-order eigenvalues in matter such that it matches the exact eigenvalues to order , see section 5. We will show that use of the renormalized basis makes the structure of the perturbation theory exceptionally transparent, as we will explain in the next section. It allows us to obtain simple, elegant and compact expressions for the oscillation probabilities, which have a universal form even to first order in our expansion parameter. For example, ν e survival probability takes the form to order as P (ν e → ν e ) = 1 − sin 2 2φ sin 2 (λ + − λ − )L 4E (1.1)

JHEP01(2016)180
where φ is θ 13 in matter, and λ ± denote the eigenvalues of the states which participate the 1-3 level crossing. Despite its extremely simple form, P (ν e → ν e ) in (1.1) takes into account all order effects of θ 13 and the matter potential. Since we will only consider terms up to order , in this paper, our results here are not applicable to the region near the solar MSW resonance [1,11]. Our perturbative framework will be called as the "renormalized helio-perturbation theory" in the rest of this paper. The section plan of this paper is somewhat unusual: in the next section 2 we summarize the key features of the perturbative framework we develop in this paper. Then, in section 3 we describe the principle results of this work including the oscillation probabilities for all channels in matter. This section does not describe the derivations but provides a self contained summary of the results of this paper. Following this section, see section 4, we present a systematic exposition of our perturbative framework and how the results of the section 3 are derived. In the appendices A, B, and C we present, respectively, calculational details of the S matrix, the results of oscillation probabilities to order in the standardized form, and useful relationships to verify the equivalence of various expressions.

Structure of the neutrino oscillation probabilities in vacuum and in uniform matter
We start by giving a precise definition of what we mean by "simplicity and compactness" of the expressions for the neutrino oscillation probabilities. Then, we explain the need for reformulating perturbative treatment, followed by outlining "what is new" in this paper.

Simplicity and compactness
In vacuum and in matter with constant density, it is well known that the neutrino oscillation probabilities for ν β → ν α for three-flavor mixing (i, j = 1, 2, 3) can be written in the following form 1 The flavor mixing matrix elements, V αi , relate the flavor states, ν α , to the eigenstates of the Hamiltonian, ν i , with eigenvalues λ i /2E by ν α = V αi ν i . Equation (2.1) applies both in vacuum and in uniform matter since the matrix V diagonalizes the full Hamiltonian, which includes the Wolfenstein matter potential, a [1]. In vacuum, λ i = m 2 i where m i denotes the mass of the i-th neutrino state and V = U MNS .
Notice that the identical form of the oscillation probabilities in vacuum and in matter, eq. (2.1), apart from replacement of ∆m 2 ij by (λ i − λ j ), implies that physical interpretation of the formulas in matter is as transparent as in vacuum. 1 We have used for the last term in (2.1) ] which follows from unitarity, and the identity (C.2).

JHEP01(2016)180
What is less well appreciated is that the expressions of the oscillation probabilities in (2.1) are maximally simple and compact. That is, they contain 5 (including the constant term) functions of L/E which are linearly independent. This property can be easily verified by showing the Wronskian is nonvanishing, which implies that none of the functions can be written in terms of the other functions for all L/E. Hence, in a true mathematical sense, eq. (2.1) and it's equivalents give the simplest and most compact form for the complete set of three flavor oscillation probabilities in vacuum and in uniform matter. 2 A special feature of the oscillation probabilities which is also worth noting is that each term in eq. (2.1) factorizes into the characteristic sin [(λ j − λ i )L/4E] factor and the products of the V matrix elements which control the amplitude of the oscillation. Both the eigenvalues, λ i and the matrix elements of V are independent of the baseline, L, but are functions of the mixing angles, θ's, the ∆m 2 ji , and the product of the energy of neutrino times the matter density via the matter potential. 3 CP and T violation is described by the last term in (2.1), which has the universal, channel-independent form in the three neutrino mixing.

Raison d'être and the requirements for perturbative treatment
With the simplest form of the oscillation probability eq. (2.1) and by knowing both the exact form of the eigenvalues, λ i /2E, and the elements of the V matrix, see [3], one might expect we have all that is needed for theoretical discussions. Unfortunately, the analytic expressions for these eigenvalues, λ i /2E, as well as the V matrix elements are notoriously complex and give no analytic insight into the oscillation physics in uniform matter. This is true even when one of ∆m 2 21 , sin θ 13 , sin θ 12 or a is set equal to zero. In any one of these limits, the characteristic equation for the eigenvalues factorizes. But, the form of the general solution does not simplify trivially to yield the correct eigenvalues, even though it must. The structure of the general solutions of the cubic characteristic equation is the root cause of this rather unusual behaviour. Hence, there is a need for a reformulated perturbative framework so that we can obtain approximate but much simpler expressions for the eigenvalues, λ i , and the mixing matrix elements, V αi , which provide the necessary physics insight.
In this paper, we formulate perturbation theory by which we can calculate the eigenvalues λ i and the elements of V matrix as a simple power series of the small expansion parameter, a renormalized ∆m 2 /∆m 2 ⊕ . At the same time the structure-revealing form of the oscillation probabilities (2.1) is kept intact. While the existing frameworks do not satisfy the latter requirement, the key to the success in our case is due to the correct decomposition of the Hamiltonian into the unperturbed and the perturbed terms. Use of the renormalized zeroth-order basis, allows us to correctly determine the eigenvalues to the appropriate order in our expansion parameter. The resultant renormalized eigenvalues and 2 There are equivalent ways to write this set of oscillation probabilities with 5 independent L/E functions, e.g one could use the trigonometric identity 2 sin 2 x = 1 − cos 2x to replace the sin 2 x. Any other way of writing these oscillation probabilities will have 5 or more linearly independent functions, e.g., using the identity (C.2) increases the number of independent functions by 2.
3 a ∝ ρE with ρ being the matter density. The V matrix elements are determined solely by the elements of the Hamiltonian matrix. It is easy to observe this feature by using the method developed in ref. [4], which is valid for nonuniform matter density as well.

JHEP01(2016)180
the mixing parameters effectively absorb the additional terms that arise in the conventional perturbative frameworks, as extra functions to the minimally required 5 linearly independent functions. This occurs automatically and we believe that such a framework has never before been formulated. 4 Because of the structural simplicity of the L/E dependence, the expressions of the oscillation probabilities calculated by our method are extremely simple and compact, as will be fully demonstrated in the next section.
3 Results of our perturbative expansion for all oscillation probabilities In this section, we describe the main results of this paper without derivations and with only minimal discussion. In later sections we provide the derivation and more in depth discussions. We start with the approximate eigenvalues of the Hamiltonian, the approximate neutrino mixing matrix and then give the oscillation probabilities for all channels to first order in the expansion parameter, , see eq. (3.3) for the precise definition.

Mass squared eigenvalues in matter
In vacuum the three eigenvalues of the full Hamiltonian which governs the neutrino oscillation is given in the form m 2 i /2E, where m i is the mass of i-th mass eigenstate of neutrinos, i = 1, 2, 3. Similarly, in matter we write the three eigenvalues as where the state label runs over i = −, 0, + for the approximate Hamiltonian of three flavor mixing system. To treat the normal and the inverted mass orderings (NO and IO respectively) in a unified way, we define the eigenvalues as follows 5 λ − = 1 2 ∆m 2 ren + a − sign(∆m 2 ren ) (∆m 2 ren − a) 2 + 4s 2 13 a∆m 2 ren + ∆m 2 ren s 2 12 , λ 0 = c 2 12 ∆m 2 ren , (3.1) λ + = 1 2 ∆m 2 ren + a + sign(∆m 2 ren ) (∆m 2 ren − a) 2 + 4s 2 13 a∆m 2 ren + ∆m 2 ren s 2 12 .
In eq. (3.1), the renomalized ∆m 2 ≡ ∆m 2 ren , the expansion parameter , and the Wolfenstein matter potential [1], a, are defined as follows: 6 4 See, however, the clarifying remark at the end of subsection 3.3.6. 5 We note that the eigenvalues in (3.1) above appear in ref. [5]. See section 4.1 for the derivation and a comment on the treatment in [5]. 6 The following notation is used throughout: ∆m 2 ij ≡ m 2 i − m 2 j , sij = sin θij and cij = cos θij where θij are the standard neutrino mixing angles and GF is the Fermi constant, Ne is the number density of electrons, E is the energy of the neutrino, Ye the electron fraction and ρ is the density of matter.

JHEP01(2016)180
This choice of ∆m 2 ren is crucial to the compact formulas for the oscillation probabilities that will be given in this paper. Note also that the sign of ∆m 2 ren signals the mass ordering, both ∆m 2 ren and are positive (negative) for NO (IO). However, for both orderings ∆m 2 ren = ∆m 2 21 > 0, as required by nature. Notice that λ 0 is the same for the both mass orderings, and when we switch from NO to IO we also switch the sign in front of the square root in eq. (3.1). The nicest feature of the sign choice is that the oscillation probability has a unified expression and the solar resonance is in ν − -ν 0 level crossing for the both mass orderings. ∆m 2 ren is equal to the effective atmospheric ∆m 2 measured in a electron (anti-) neutrino disappearance experiment in vacuum, ∆m 2 ee ≡ c 2 12 ∆m 2 31 + s 2 12 ∆m 2 32 [12]. This quantity is identical to ∆m 2 ee recently measured by the reactor θ 13 experiments [13][14][15] up to effects of O(∆m 2 21 /∆m 2 31 ) 2 . Whether the coincidence between ∆m 2 ee and ∆m 2 ren reflects a deep aspect of neutrino oscillation or not will be judged depending upon what happens at second order in . This point as well as the relevance of the other effective ∆m 2 µµ [12], ν µ equivalent of ∆m 2 ee , will be discussed in depth in a forthcoming communication.

The mixing angle θ 13 and mixing matrix in matter
We use the angle φ to represent the mixing angle θ 13 in matter. With the definitions of the eigenvalues (3.1), the following mass-ordering independent expressions for cosine and sine 2φ (see section 4.2) are given by cos 2φ = ∆m 2 ren cos 2θ 13 − a λ + − λ − , sin 2φ = ∆m 2 ren sin 2θ 13 It is easy to show that φ goes from 0 → π/2 as a goes from − ∞ to + ∞ for the NO and as a goes from + ∞ to − ∞ for the IO. In vacuum (a = 0), φ = θ 13 and φ = π/4 at the atmospheric resonance, when a = ∆m 2 ren cos 2θ 13 , for both mass orderings. The mixing matrix in matter, V , relates the flavor eigenstates, ν e , ν µ , ν τ , to the perturbatively defined matter mass eigenstates, ν − , ν 0 , ν + as follows (see section 4.4): where the matrix V is unitary. It is convenient to split V into a zeroth order term, V (0) , and a first order term, V (1) in our expansion, where the zeroth order matrix is given by

JHEP01(2016)180
with δ being the CP violating phase, whereas the first order correction is given by (3.9) A factorized form of the V matrix is given in eq. (4.25).
As an outcome of the consistent perturbative treatment the total V matrix given by (3.7) with (3.8) and (3.9) must be unitary to order . In fact it is, since the following two conditions are satisfied (3.10) Of course, none of what follows is self consistent without unitarity here.
With the matter eigenvalues, λ's , definite by eq. (3.1) and the matter mixing matrix, V , given by eq. (3.7), simple and compact expressions can be easily derived for the oscillation probabilities in matter for all channels, to leading order in , as will be shown in the next section.

Compact formulas for the oscillation probabilities in matter
In this section we start by presenting the shortest path to the oscillation probabilities of the ν e -related channel by using the eigenvalues, λ ±,0 , and mixing matrix, V , given in the previous section to order . Then, we derive the universal expression of the oscillation probabilities which is applicable to all channels.

JHEP01(2016)180
Notice that the formula in eq. (3.11) takes into account the matter effect as well as the effect of s 13 to all orders. Nonetheless, it keeps an exceptional simplicity, an effective twoflavor form in matter which consists of single term with the unique eigenvalue difference λ + −λ − , the feature we believe to be unique in the market. The feature stems from the fact that there is no ν e component at zeroth order in in the "0" state in matter. It is expressed in the zero in the V e0 element of the zeroth-order V matrix as in (3.8), see also section 4.4.
3.3.2 ν e → ν µ and ν e → ν τ appearance channels Now, we discuss the appearance channels ν e → ν µ and ν e → ν τ . We describe here the simplest way to derive the formulas for the oscillation probabilities starting from the V matrix by using unitarity. The oscillation probability P (ν e → ν µ ) can be computed as (3.12). Then, we obtain where the common shorthand notation for the kinematic phase here J r , the reduced Jarlskog factor, is This expression for the ν e → ν µ appearance channel probability is quite compact, despite that it contains all-order contributions of s 13 and a. In particular, it keeps the similar structure as the one derived by the Cervera et al. [7], which retains terms of order 2 but is expanded by s 13 only up to second order. This method of computing P (ν e → ν µ ) in the above offers the shortest path to the expression of the oscillation probability which is manifestly free from the apparent singularity as λ − → λ 0 because 1/(λ − − λ 0 ) always appears adjacent to sin [(λ − − λ 0 )L/4E]. We will refer this method as the "shortcut method" in the rest of this paper.

JHEP01(2016)180
The expression P (ν e → ν µ ), eq. (3.14), is not quite of the form of eq. (2.1) because of . This can be easily remedied, by using the following identity which leads to exactly the form of eq. (2.1). However, if this is used then some of the terms are singular when λ − = λ 0 , yet the total expression is equivalent to eq. (3.14) and is finite. This is the reason why we prefer the form of eq. (3.14). The oscillation probability for ν e → ν τ channel can be obtained by all of the following three methods: (1) the similar calculation by the shortcut method, (2) using the unitarity relation P (ν e → ν τ ) = 1 − P (ν e → ν e ) − P (ν e → ν µ ), (3) using the relation P (ν e → ν τ ) = P (ν e → ν µ : c 23 → −s 23 , s 23 → c 23 ) [10], whose derivation is sketched in section 4.1. We note that J r changes sign by the transformation. In general the oscillation probability in arbitrary channel can be obtained by (1) the shortcut method, or by (2) rewriting the expressions in appendix B using the formulas in appendix C.

The general form of the ν α → ν β oscillation probabilities
The expressions of the oscillation probabilities in the ν µ − ν τ sector can be derived by one of the methods mentioned in the previous subsection. Then, one observes a remarkable feature that all the oscillation probabilities P (ν α → ν β ) (including the ν e sector) can be written in a universal form: where the eight coefficients A αβ ij , B αβ ij , C αβ and S αβ are given in table 1. Notice that they are 0, ± 1, or the simple functions of θ 23 .
In table 1, we observe that three relations hold between the coefficients These identities hold due to the invariance of the oscillation probabilities under the following transformation φ → π/2 + φ and λ + ↔ λ − .
Order 0 : Order cos δ: Order sin δ: This invariance must be respected by the oscillation probabilities because the two cases in eq. (3.19) are both equally valid ways of diagonalizing the zeroth-order Hamiltonian. Look at eq. (3.5) to confirm that it is invariant under (3.19). Then, the former two identities in (3.18) trivially follow, but the last one requires use of the kinematic relationship 7 We note that θ 23 dependence in the oscillation probabilities, apart from that in J r , is confined into the five independent coefficients. All the other parameters, ∆m 2 's, θ 12 , θ 13 , δ, E and a, are contained in the remaining part of the probabilities, which takes the universal form for all the oscillation channels.
The antineutrino oscillation probabilities P (ν α →ν β ) can be easily obtained from the neutrino oscillation probabilities. One can show that by taking complex conjugate of the evolution equation for anti-neutrinos of energy E that it is equivalent to the evolution equation for neutrino with energy equal to −E. Therefore, (3.20) Here, we give some critical observations regarding the formulas (3.17) and the associated table 1: • The oscillation probability P (ν e → ν e ), to first order in , is of a simple two-flavor form with the matter-modified mixing angle θ 13 and the mass squared difference given by ∆m 2 ren . As far as we know, having single two-flavor form as the expression of P (ν e → ν e ) in matter at this order is the unique case among all the perturbative frameworks so far studied. 7 Notice that the relation B αβ +− = C αβ needs to be satisfied only to order 0 because these terms are already suppressed by .

JHEP01(2016)180
• Formulating the perturbation theory of neutrino oscillation with automatic implementation of the oscillation probability formulas (2.1), which is to be carried out in section 4, is the principal result of this paper. However, the formula (3.17) goes one step further beyond eq. (2.1) by displaying the channel-independent universal function of L/E, up to the coefficients A ± , B ± , etc. which depend only on θ 23 .
• With the universal form of the probabilities, unitarity can be trivially checked just by adding up the appropriate columns of table 1 with proper care of using the correct signs in the last row for the intrinsic CP-violating terms. 8 One can easily confirm it by adding, for example, the ν e → ν e , ν e → ν µ , and ν e → ν τ columns to give zero in each row.
• If one uses the identity (3.16) to transform the sin term to the right-hand side of (3.16), then the transformational invariance given by eq. (3.19) will be manifest and eq. (3.17) will become the form of eq. (2.1). However, in the limit λ − → λ 0 some terms will diverge even though the total expression is perfectly finite.
• In the universal form, the amplitude of a given oscillation frequency, to leading (zeroth) order in , is controlled by φ. The amplitudes of the frequencies For the channels involving ν e , only the (λ + − λ − ) frequency appears whereas in the ν µ − ν τ sector all frequencies appear and are of equal magnitude near the θ 13 resonance, when φ ≈ π/4, but every where else either the frequency (λ + − λ 0 ) ≈ ∆m 2 ren or (λ − − λ 0 ) ≈ ∆m 2 ren dominates.
Finally, we give more clarifying remarks on generic features of the oscillation probability formulas presented in this section 3.
• The term proportional to sin δ in P (ν α → ν β ), in eq. (3.17), is equal to up to an overall sign. This is because the quantity in {· · · } is just the Jarlskog factor in matter due to the Naumov-Harrison-Scott identity, [17,18].
• Due to the above theorem sin δ terms are always accompanied by the Jarlskog factor J r = c 12 s 12 c 23 s 23 c 2 13 s 13 . Though less well known, it can be proven that cos δ terms must be proportional to the reduced Jarlskog factor J r /c 2 13 [19]. It can also be shown that both of the δ-dependent terms must come with the suppression factor of [17][18][19]. It is nice to observe that these general properties are realized in our formulas.
• Association of the factor and the all angle factors to the δ-dependent terms indicates that they are genuine three flavor effects. In our general formula (3.17), all the terms JHEP01(2016)180 with explicit factor are of this δ-dependent terms. Notice that other corrections of order , though they do not show up explicitly, are implicitly contained in ∆m 2 ren , λ ±,0 and φ which do not contain any δ or θ 23 dependence.
• In vacuum, a = 0, the above oscillation probabilities reproduce the standard results, to first order in . The form is somewhat unusual but we have checked that the expressions are identical. Also, we have checked that in the limit L/E → 0, where matter effects are negligible, the coefficients of the (L/E) n , for n ≤ 3, terms in eq. (3.17), are identical to those in vacuum for all oscillation channels.

Range of applicability of the formulas
To discuss the range of applicability of our expressions, it is useful to first consider the vacuum expressions to first order in the expansion parameter . For all channels, the expansion of the vacuum oscillation probabilities to first order in does not include terms proportional to sin 2 ∆ 21 which starts at second order in , where ∆ ji ≡ ∆m 2 ji L/4E. When, ∆ 31 = π/2, that is at the first atmospheric oscillation maximum, 2 ∆ 2 31 ≈ 0.002 which is small for the channels ν e → ν x where x = e, µ, or τ since 1 − P (ν e → ν e ), P (ν e → ν µ ) and P (ν e → ν τ ) are all of order sin 2 2θ 13 ≈ 0.1. However, at the second atmospheric oscillation maximum, ∆ 2 31 = 3π/2 and 2 ∆ 2 31 ≈ 0.02, which is significant compared to the sin 2 2θ 13 term. So in vacuum our first order expansion is only a good approximation for ∆ 2 31 π or L/E 1000 km/GeV for these ν e channels. For the other channels, ν µ → ν µ , ν τ → ν τ and ν µ → ν τ , our first order expansion is a good approximation to somewhat beyond L/E = 1000 km/GeV because the leading terms are not suppressed by the smallness of sin 2 2θ 13 .
Then, what about the validity in matter? In section 5 we will argue that our perturbative description is valid outside the solar resonance. Notice that the region without validity (no guarantee for approximation being good) is rather wide and includes the vacuum because the solar resonance width |∆a| = √ 3(sin 2θ 12 / cos 2 θ 13 )∆m 2 21 is larger than the solar resonance position a = (cos 2θ 12 / cos 2 θ 13 )∆m 2 21 . We expect then that our helioperturbation theory works for the matter potential a larger than a few tenth of |∆m 2 ren |. To give the reader a sense of the precision of our approximation we have plotted in figure 1, the contours of equal probability for the exact and the approximate solutions for the channels ν e → ν µ , ν e → ν e and ν µ → ν µ . The right (left) half plane of each panel of figure 1 corresponds to the neutrino (anti-neutrino) channel. As expected, for large values of the matter potential, |a| > 1 3 |∆m 2 ren | we find we have no restrictions on L/E, to have a good approximation to the exact numerical solutions. Whereas for small values of the matter potential, |a| < 1 3 |∆m 2 ren | we still need the restriction L/E 1000 km/GeV. 9 JHEP01(2016)180 Figure 1. The iso-probability contours for the exact (solid blue) and approximate (dashed red) oscillation probabilities for upper left, ν e → ν e , upper right, ν e → ν µ and lower, ν µ → ν µ . The upper (lower) half plane is for normal ordering (inverted ordering), whereas positive (negative) L/E is for neutrinos (antineutrinos). For treatment of antineutrinos, see section 3.3.3. The order of the contours given in the title is determined from the line L/E=0. The discontinued as one crosses Y e ρ|E| = 0 is because we are switching mass orderings at this point. In most of parameter space the approximate and exact contours sit on top of one another so the lines appear to alternate blue-red dashed. Note that, for L/E >1000 km/GeV and |Y e ρE| < 5 g.cm −3 .GeV, the difference between the exact and approximate contours becomes noticeable at least for ν e → ν e and ν e → ν µ .

JHEP01(2016)180
We note that most of the settings for the ongoing and the proposed experiments, except possibly for the one which utilizes the second oscillation maximum, fall into the region L/E 1000 km/GeV. To improve the accuracy to larger values of L/E, especially for values of |a| < 1 3 |∆m 2 ren |, second order perturbation theory in is needed, which will be the subject of a future publication.

An example of the power of our oscillation probabilities
In this section we present an example of the power of our compact expression of the oscillation probabilities in understanding the physics of oscillations to first order in our expansion parameter . For simplicity we consider the ν e disappearance channel. In figure 2 the ν e disappearance probability is shown, as a function of neutrino energy E, for baselines, L, of 3000 km (upper panel) and 5000 km (lower panel). As in figure 1, the solid blue curves are drawn by using the exact expression of P (ν e → ν e ), and the dashed red curves by our compact formula in eq. (3.11). The black dotted curves are for the vacuum case.
We first notice that our approximate formula agrees well with the exact expression in particular at higher energy. Secondly, because of the long baselines, the matter effect is sizeable, producing not only a large shift in the position of the first oscillation minimum but also an significant increase in the depth of the minimum. These two effects are correlated by the energy dependent quantity, (λ + − λ − ) which can be seen graphically in figure 3: the depth of the oscillation minimum changes from whereas the energy at which the the first oscillation minimum occurs changes from Because of the simplicity of our expression for P (ν e → ν e ), eq. (3.11), these shifts are accurate to first order in the expansion parameter . This simple understanding of the features of P (ν e → ν e ) is new to this paper.

Comparison with the existing perturbative frameworks
As we emphasized in section 2, our machinery has an advantage over the existing perturbative frameworks by having the minimum number of terms composed of sin [(λ j − λ i )L/4E] (i, j = 1, 2, 3) in the oscillation probabilities. This contrasts with the features of the existing perturbative frameworks in which much larger number of terms than those minimally necessary as in (2.1) are produced. They include, typically, the terms with either extra L/E dependences or different frequencies in the sine functions, or often both, which easily obscures the physical interpretation.
To give a feeling to the readers on how simple and compact our formulas are, we compare our expressions to the ones in the existing literatures to the same order in expansion. For definiteness, we pick the ones in ref. [10] to make the comparison, the most recent one among the reference list given in section 1. Our expression of P (ν e → ν e ) in (3.11) which has only a single term (ignoring unity) may be compared with eqs. (4.6) and (4.7) which consist of total 3 terms. With regard to P (ν e → ν µ ), if we count numbers of terms with different L/E dependence we have one δindependent and 2 δ-dependent terms in (3.14), total 3. Whereas eqs. (4.8) and (4.9) in [10] have total 7 terms, 3 δ-independent and 4 δ-dependent ones. Our expression of P (ν µ → ν τ ) in (3.17) has 3 δ-independent and 2 δ-dependent terms, adding up to total 5. On the other hand, eqs. (4.10) and (4.11) in [10] contain total 8 δ-independent and 10 δ-dependent terms (2 sin δ and 8 cos δ terms), which adds up to 18. So not only do our expressions for the oscillation probabilities have less than half the number of L/E structures, but the form of the L/E dependence is made manifest and identical to the vacuum form.

JHEP01(2016)180
We note that when our formulas are expanded by ∆m 2 21 /∆m 2 31 , they agree, of course, with the existing ones calculated previously. Therefore, our formalism may be regarded as a systematic way of organizing the equivalent expressions to order ∆m 2 21 /∆m 2 31 into neater, structure-revealing ones. The benefit for having such simple expressions is that the physical interpretation of the formulas is transparent, as emphasized in section 2.

Formulating the renormalized helio-perturbation theory
In this section, we formulate the helio-to-terrestrial ratio perturbation theory, for short the helio-perturbation theory, which has the unique expansion parameter two ways of deriving the oscillation probabilities, the S matrix method and the wave function method. Here we sketch both of them, leaving technical or computational parts into appendices A and B. The meaning of the agreement between results obtained by both the S matrix and the wave function methods will be discussed at the end of this section. The S matrix describes neutrino flavor changes ν β → ν α after traversing a distance L, ν α (L) = S αβ ν β (0), (4.2) and the oscillation probability is given by When the neutrino evolution is governed by the Schrödinger equation, i d dx ν = Hν, S matrix is given as where T symbol indicates the "time ordering" (in fact "space ordering" here). In the standard three-flavor neutrinos, Hamiltonian is given by where the symbols are defined in an earlier footnote. For the case of constant matter density, the right-hand side of (4.4) may be written as e −iHx . We recapitulate here the earlier footnote: in (4.5) ∆m 2 ji ≡ m 2 j −m 2 i where m i denotes the mass of i-th mass eigenstate neutrinos. Position dependent function a(x) ≡ 2 √ 2G F N e (x)E is a coefficient for measuring the matter effect on neutrinos propagating in medium of electron number density N e (x) [1] where G F is the Fermi constant and E is the neutrino energy.
The neutrino flavor mixing matrix U is usually taken to be the standard form U PDG given by Particle Data Group. We, however, prefer to work in a slightly different basis, for this paper, in which the flavor mixing matrix has a form (with the obvious notations s ij ≡ sin θ ij etc. and δ being the CP violating phase) under the understanding that the left phase matrix in the first line in eq. (4.6) is to be absorbed into the ν τ neutrino wave functions. Being connected by the phase matrix rotation, it is obvious that our U and U PDG give rise to the same neutrino oscillation probabilities.

Choosing the basis for the renormalized helio-perturbation theory
It is convenient to work with the tilde basis defined asν α = (U † 23 ) αβ ν β , in which the Hamiltonian is related to the flavor basis one as where U 23 is defined in eq. (4.6). The S matrix in the flavor basis is related to the S matrix in the tilde basisS as Notice that the matter term in the Hamiltonian (4.5) is invariant under the U 23 rotation. Hence, the dynamics of neutrino propagation in matter is governed by only the two mixing angles θ 12 and θ 13 [22], which is also independent of δ thanks to our convention of U in (4.6). From the first equation in (4.8), the relationships P (ν e → ν τ ) = P (ν e → ν µ : c 23 → −s 23 , s 23 → c 23 ) etc. simply follow as used in section 3.3.
The simplest formulation of helio-perturbative treatment in the tilde-basis includes decomposition ofH into the zeroth and the first order terms in the expansion parameter where ∆m 2 ren ≡ ∆m 2 31 −s 2 12 ∆m 2 21 and ≡ ∆m 2 21 /∆m 2 ren , as defined in (3.2) and (3.3).H(x) with (4.10) and (4.11) is identical with the tilde-Hamiltonian in (4.9). Note the simplicity of the perturbing Hamiltonian,H 1 , especially the zeros on the diagonal elements. The diagonalization ofH 0 leads to the eigenvalues given in section 3.1, eq. (3.1), and because of the zeros along the diagonal ofH 1 , there are no corrections to these eigenvalues at first order in a perturbative expansion.
The authors of [5] treat the order effect in the Hamiltonian as a renormalization of the matter potential, whereas we regard it as a renormalization of ∆m 2 31 .

JHEP01(2016)180
Our treatment can be easily generalized to cases with matter density variation as far as the adiabatic approximation holds. See e.g., [10] for such a treatment. However, we prefer to remain, for simplicity and compactness of the expressions, into the formulas with constant matter density approximation in the rest of this paper.

Hat basis
We transform the HamiltonianĤ =Ĥ 0 +Ĥ 1 , from the"tilde" basis to the "hat" basis, usingĤ where the unperturbed HamiltonianĤ 0 is diagonal, We take the following form of unitary matrix U φ to diagonalizeH 0 : (4.14) The expressions of the zeroth order eigenvalues λ − , λ 0 , and λ + , are given in eq. (3.1). Similarly, cosine and sine φ are given in eq. (3.5). Also the perturbing Hamiltonian,Ĥ 1 , retains it's simple form thanks to that the U φ rotation keeps "zeros" inH 1 unchanged both on the diagonal and in the top-right and bottom-left corners, In fact,Ĥ 1 is identical toH 1 with θ 13 replaced by (θ 13 − φ).

S matrix and the oscillation probability
The S matrix in the flavor basis is related to the S matrix in the tilde and the hat bases as where we have used explicitly the fact that the matter density is constant:

JHEP01(2016)180
Then, Ω(L) obeys the evolution equation Then, Ω(x) can be computed perturbatively as Collecting the formulas the S matrix can be written as Thus, we are left with perturbative computation of Ω(L) with use ofȞ 1 in (4.20) to calculate the S matrix. With the S matrix in hand it is straightforward to compute the oscillation probabilities by using (4.3). We leave these tasks to appendices A and B.

Mass eigenstate in matter: V matrix method
In this section we calculate the V -matrix directly using our perturbation theory. If we switch off the perturbationĤ 1 , the mass eigenstates in matter, to lowest order, are given by the hat-basis wave functionν i , which are the eigenstates ofĤ 0 in (4.12), and sincê H 0 is diagonal, we haveν Thus, the V matrix is given to zeroth order by V (0) = U 23 U φ whose explicit form is given in section 5.3, and also in eq. (3.8).
In order to obtain the mass eigenstates in matter to first order in , ν i =ν (0) i +ν (1) i , let us compute the first order correction to the hat basis wave functions. Using the familiar perturbative formula for the perturbed wave functionŝ withĤ 1 in (4.12), and the λ i 's are given by the eigenvalues ofĤ 0 , see (3.1). Then the mass eigenstate in matter ν i can be written to first order in as: Using (4.23), this equation is of the form ν i = V † ν α which can be inverted to easily obtained the V -matrix given in eq. (3.7),

JHEP01(2016)180
This can be used to directly compute the oscillation probabilities as was performed in section (3.3), or by inserting the V matrix elements into the expressions of the oscillation probabilities in (4.6). This V matrix method was used to calculate the matter effect correction in the oscillation probabilities [20]. For a different approach toward simpler expressions of the V matrix elements see [21].
In sections 4.3 and 4.4, we have sketched the two different methods for calculating the oscillation probabilities, the S-matrix method (section 4.3) and the V -matrix method (section 4.4). The results obtained by the both methods agree with each other, and the expressions of the oscillation probabilities are of the form given in eq. (2.1).
Let us make a comment on what happens if we use a different decomposition of the Hamiltonian into unperturbed and perturbed parts as in (4.9). Then, the perturbed Hamiltonian has nonzero diagonal terms, and this produces first order corrections to the eigenvalues. If we use these expressions for the eigenvalues and expand by the small expansion parameter we obtain terms that do not exist in (2.1), such as the ones proportional to sin [(λ j − λ i )L/2E]. Of course, the S-matrix method will give the same expressions. Therefore, absence of the explicit first-order correction to the eigenvalues is essential to keep perturbative expressions of the oscillation probabilities of the form (2.1) to order . In our case it is guaranteed by vanishing diagonal elements of the perturbed Hamiltonian (4.11).

More about the renormalized helio-perturbation theory
In this section, we critically examine the framework of the renormalized helio-perturbation theory. Despite a drawback of the current framework (which is to be described below) we argue that our perturbation theory works apart from the vicinity of the solar resonance crossing.

Exact versus zeroth-order eigenvalues in matter
The three eigenvalues of the Hamiltonian are written as λ i 2E , where i runs over 1, 2, 3 for the exact eigenvalues, and i = −, 0, + for the zeroth-order eigenvalues in our perturbative framework. In figure 3 the λ i are plotted as a function of a, the Wolfenstein matter potential, for both the exact and the zeroth-order ones given by eq. (3.1). It is clear in figure 3 that our zeroth-order eigenvalues fail to treat the solar-∆m 2 scale level crossing correctly. As one goes through the solar resonance in the exact solution the two eigenvalues involved, the red and green, repel one another, whereas in our perturbative solution the two corresponding eigenvalues cross with each other.
We point out here that the drawback is the common feature among the similar perturbative treatments of neutrino oscillation available in the market. In this section we argue, for the first time, that despite the drawback, our perturbative framework successfully treat the flavor content of these two states in region reasonably far from the solar resonance.
We first note that, despite the feature, the atmospheric and solar resonances occur at the correct values of the matter potential with our zeroth-order eigenvalues. The atmospheric resonance occurs when λ + − λ − is a minimum, which is at a = ∆m 2 ren cos 2θ 13 JHEP01(2016)180 . The exact eigenvalues are depicted with colored lines, green for ν 1 , red for ν 2 , and flight blue for ν 3 . The eigenvalues calculated by using our renormalized helio-perturbation theory are drawn by the black dashed lines whose state labels are marked on the figures. The approximate eigenvalues for ν 0 and ν − cross at the solar resonance, whereas the exact eigenvalues for ν 1 and ν 2 repelled at solar resonance. To make these features visible we used a solar ∆m 2 21 three times as large as the measured value. and the minimum difference between λ + and λ − is ∆m 2 ren sin 2θ 13 as expected. The solar resonance occurs when λ − − λ 0 is a minimum; this occurs when a = ∆m 2 ren cos 2θ 12 / cos 2 θ 13 and the minimum difference is zero! This is not the value determined by the full Hamiltonian which is ≈ ∆m 2 ren sin 2θ 12 . Therefore, while our perturbative scheme treats the atmospheric resonance correctly to order , it misses the effects of the solar resonance.

Eigenvalues in vacuum and in the asymptotic regions a → ±∞
In this and the next subsections, we will give the arguments to indicate that our renormalized helio-perturbation theory works apart from the vicinity of the solar resonance despite the issue mentioned above.
We first show that the zeroth-order eigenvalues given in (3.1) agrees with the exact ones to order in the asymptotic regions a → ±∞. We use the characteristic equation of the full Hamiltonian (4.5) to derive (i = 1, 2, 3) whose order terms are different from the exact values, λ 1 = 0, λ 2 = ∆m 2 21 , and λ 3 = ∆m 2 31 . To understand the meaning of the failure at order we need to discuss what happens to the flavor content of the matter mass eigenstates as the matter potential changes from below to above the solar resonance. This will be done in the next subsection.
Similarly, the asymptotic behaviour of the angle φ, i.e., θ 13 in matter can be easily worked out. With the definition of φ in (3.5), it is easy to show that φ takes on the following values as a is varied from − ∞ to + ∞. In the case of normal mass ordering,
It reflects a natural view that physics at a → +∞ for the normal mass ordering corresponds to the one at a → −∞ for the inverted mass ordering at least to leading order in .

Neutrino mixing matrix in matter and the ν flavor content at a → ±∞
Since the two levels cross at the solar resonance in our perturbative treatment, one may expect that our treatment fails completely beyond the solar resonance, i.e., in the region with matter density higher than the resonance. However, we will show in this subsection that the flavor contents of the three eigenstates are correctly reproduced at least in the asymptotic region. That is, the zeroth-order V matrix describes correctly the asymptotic behaviour of the exact eigenstates in matter. Suppose we denote the exact flavor mixing matrix in matter as the (almost standard) U matrix defined in (4.6). Let us discuss the case of NO first. At a → −∞ U is nothing but V (0) in (3.8). Notice that s 23 in matter is frozen to its vacuum value for < 0.1, and s 12 0 at a → −∞ [3]. At a → +∞, s 12 1 and c 12 0 the matter U matrix is identical to (3.8) if we interchange ν 1 and ν 2 apart from re-phasing factor −1 for the new second mass eigenstate.
To make the meaning of this feature clearer we write down here the flavor content of the states ν i (i = 1, 2, 3) at a → ±∞. Noticing that ν i = (V † ) iα ν α , the flavor composition at a → −∞ is given at zeroth order by Whereas the composition at a → +∞ is given by The flavor compositions given in (5.8) and (5.9) imply that the flavor content of the lower two mass eigenstate in matter is correctly described in our perturbative framework despite the failure of describing the solar level crossing. Note, the combinations of ν µ and ν τ in the (· · · ) and {· · · } in (5.8) and (5.9) are orthogonal with each other.
In the case of IO, essentially the same discussion goes through. The asymptotic behavior of θ 12 at a → ±∞ is the same as that in NO. The asymptotic behavior of θ 13 at a → ±∞ for NO is mapped into the one at a → ∓∞ for IO. But, this is already taken care of by the definition of φ given in (3.5). Therefore, the same U matrix in matter as the one JHEP01(2016)180 in NO are obtained at a → ±∞. Then, the same flavor compositions as in (5.8) and (5.9) follow for IO. Again it is consistent with those we expect from the level crossing diagram.
We add that our ν + state always corresponds to ν 3 state. At a → +∞ in NO and at a → −∞ in IO the electron neutrino component is all in this ν + = ν 3 state. At a → −∞ in NO and at a → +∞ in IO the electron neutrino component is all in ν − (ν 1 for NO, ν 2 for IO) state. In vacuum for both NO and IO, V (0) e+ = s φ = s 13 , which is the correct value for the ν 3 = ν + state.
To summarize: despite that the eigenvalues calculated by our helio-perturbation theory do not show the correct behavior at around the solar level crossing, the flavor composition of the states are correctly represented by our zeroth-order states. It is the very reason why our perturbative formulas work at much higher and lower values of Y e ρE compared to the one at the solar resonance. This point escaped detection in the previous literature to our knowledge. We have seen that they work practically whole regions except for the vicinity of the solar resonance.

Concluding remarks
We have developed a new perturbative framework of neutrino oscillations, which allows us to derive compact expressions of the oscillation probabilities in matter, which we examined in this paper to order ≡ ∆m 2 21 /∆m 2 ren ∆m 2 21 /∆m 2 31 . Although vacuum like in their forms, our compact results keep all-order effects of s 13 and the matter potential a, while using as the unique expansion parameter.
The characteristic features of our perturbative framework (dubbed as the renormalized helio-perturbation theory), in contrast to the ones previously studied by various authors, are as follows: • We use the renormalized basis (4.10) defined with the atmospheric mass-squared difference corrected by an order quantity, ∆m 2 ren ≡ ∆m 2 31 − s 2 12 ∆m 2 21 . The decomposition of the Hamiltonian into the unperturbed (4.10) and perturbed terms (4.11), is done in such a way that there is no diagonal entries in the perturbed one. It makes the eigenvalues of the zeroth order Hamiltonian the correct ones to order . Usage of this decomposition is crucial to obtain the form of the oscillation probabilities (2.1) akin to the form in vacuum. This facilitates an immediate physical interpretation that the frequencies of the oscillations are determined by (λ + − λ − ), (λ + − λ 0 ), and (λ − − λ 0 ), where the λ's are the eigenvalues of the zeroth order Hamiltonian, see eq. (3.1). They determine the oscillation pattern, the functional form of L/E dependence, for all the oscillation probabilities as in eq. (2.1) to this order. • From the universal form of the oscillation probabilities, eq. (3.17), we readily observe the followings: (1) Each one of the zeroth-order terms in the oscillation probabilities, in all channels, takes the same form as the corresponding two-flavor oscillation probability in vacuum but with use of the eigenvalues and θ 13 in matter. In the special case of ν e → ν e , only the term with the particular frequency ∝ λ + − λ − is non-zero, so the L/E dependence is of the two flavor form. The dominance of this frequency also occurs in the oscillation probability for ν e → ν µ and ν e → ν τ .
(2) All the first-order correction terms with the explicit factor (≈ ∆m 2 /∆m 2 ⊕ ) in the oscillation probabilities are proportional to either cos δ or sin δ. They are also suppressed by the angle factor ∝ c 12 s 12 c 23 s 23 s 13 . Hence, they are the genuine three-flavor effects. The explicit expressions of the δ-dependent terms in our formulas are in agreement with the general theorems.
Further comments are made in section 3.3.3.
Despite the success of our current framework in almost all regions including the one around the first oscillation maximum, it has a clear drawback. It fails to accommodate the physics at around the solar-scale resonance. This is related to the unphysical feature that the two eigenvalues (λ − and λ 0 ) cross with each other at the solar resonance, the universal fault in all the perturbative treatment which include the expansion parameter in the market. However, we have offered, for the first time to our knowledge, an explanation why our framework works beyond the resonance despite the failure in treating the solar level crossing. We hope that we can return to the problem as a whole in the future.
As we noticed at the end of section 3.1 our renormalized ∆m 2 ren is identical, to order , to the effective ∆m 2 measurable in an (anti-) ν e disappearance experiment in vacuum. It is a tantalizing question whether it is just a coincidence, or is an indication of something deep.
It is quite possible that the features of the current framework mentioned above naturally generalizes to higher orders. That is, one can demand that the oscillation probabilities of the general form in (2.1) calculated by the V -matrix method, with the carefully chosen zeroth-order eigenvalues, agree with the ones from the S-matrix method and are correct to certain order in . Our result in this paper may be regarded as an existence proof of this concept to order . Since we know that this is true in the exact form of the oscillation probabilities (assuming adiabaticity) [3], it is likely to be correct in each order in perturbation theory. It may or may not require higher order renormalization in ∆m Then, the oscillation probabilities are simply given by P (ν β → ν α ; L) = |S αβ | 2 , eq. (4.3).

B Expressions of neutrino oscillation probabilities
With the expressions of S matrix elements obtained in appendix A, and using eq. (4.3) it is straightforward to calculate the neutrino oscillation probabilities. Similarly, one can insert the V matrix elements given in (3.7) into (2.1) to obtain the equivalent results. In this appendix we only give the resulting expressions, but of all the oscillation channels under the viewpoint that unitarity is not to be imposed but must be proven to show the consistency of the calculation. To compare these with the results of section (3.3), use of the identities given in appendix C may be of some help.

JHEP01(2016)180
The only comment worth to give here is about the simple method for transformation c 23 → −s 23 and s 23 → c 23 to obtain P (ν e → ν τ ) from P (ν e → ν µ ), or P (ν τ → ν τ ) from P (ν µ → ν µ ), which is utilized in section 3. Though we work with the rephased flavor mixing matrix defined in (4.6), the transformations produces S eτ from S eµ , and S τ τ from S µµ , up to an overall phase, see eq. (A.6).