Perturbing neutrino oscillations around the solar resonance

Atmospheric neutrinos at low energies, $E \lsim 500$ MeV, is known to be a rich source of information of lepton mixing parameters. We formulate a simple perturbative framework to elucidate the characteristic features of neutrino oscillation at around the solar-scale enhancement due to the matter effect. The clearest message we could extract from our perturbation theory is that CP violation in the appearance oscillation probability is large, a factor of $\sim 10$ times larger than CP violation at around the atmospheric-scale oscillation maximum. Underlying mechanism for it is that one of the suppression factors on the CP phase dependent terms due to smallness of $\Delta m^2_{21} / \Delta m^2_{31}$ are dynamically lifted by the solar-scale enhancement. Our framework has a unique feature as a perturbation theory in which large $\Delta m^2_{31}$ term outside the key 1-2 sector for the solar-scale resonance does not yield sizeable corrections. On the contrary, the larger the $\Delta m^2_{31}$, the smaller the higher order corrections.


Introduction
Physics of neutrinos has been blossomed in the last 20 years after the discovery of neutrino oscillation [1][2][3][4][5], which implies nonzero masses of neutrinos and existence of lepton flavor mixing [6]. It is interesting to observe that pinning down these evidences of neutrino mass required understanding of the matter effect to a different degree of importance. While no or limited importance of the matter effect in KamLAND reactorν e measurement and the atmospheric neutrino observation in Super-Kamiokande, the theory of matter enhanced flavor conversion [7,8] provided the key to understand the results of solar neutrino experiments pioneered by Davis et al., in particular, the charged current/neutral current (CC/NC) ratio measured by SNO.
Another important feature of neutrino oscillation we see is the existence of two distance scales, atm ≡ 4πE ∆m 2 31 ∼ 100 km and solar ≡ 4πE ∆m 2 21 ∼ 3000 km at E = 100 MeV, 1 which so much differ from each other because of the hierarchy of mass squared differences ∆m 2 31 /∆m 2 21 30 in our three-generation world of fundamental leptons. It is interesting to observe that both of them easily fit into our mother planet. This is shown in the left panel of figure 2 in section 6, in which the equi-oscillation probability contours of ν µ → ν e channel is plotted in space of neutrino energy E and baseline L. It is notable that both of the two matter-effect enhanced regions, or "resonances", are clearly visible and equally prominent. Again, it is intriguing to see that the matter density of the earth is designed in such a way that they both live inside the plot as in the left panel of figure 2.
Since the discovery era, all the mixing angles and the two ∆m 2 have been determined by extensive efforts by the experimentalists, see for example, refs. [9][10][11][12][13][14][15][16] for a very limited list. The latest solar neutrino combined results are given in Fig. 13.2 of ref. [16]. Then, what is left inside νSM by now is to measure CP violating phase δ and determination of the neutrino mass ordering. It appears that the strategies for determination of "what is left" are well formulated. The tantalizing hint for δ ∼ 3π 2 from T2K experiment [17] may be confirmed at higher confidence level (CL). The next-generation long-baseline (LBL) accelerator neutrino experiments, T2HK and DUNE [18,19], are designed to perform definitive measurement of δ. Determination of the neutrino mass pattern can be carried out by utilizing the earth matter effect 2 by the same experiments, or for T2HK possibly by an extended version T2HKK [21]. It can also be done by atmospheric neutrino observation by them as well as by IceCube-Gen2/PINGU and/or KM3NeT-ORCA [22,23]. 3 It may be surprising to see that not only the ongoing and the future LBL experiments but also all the remaining ones, except for JUNO [24], rely on measurement at around the atmospheric oscillation maximum, Given the feature that the solar oscillation 1 atm = 103 km and solar = 3.30 × 10 3 km if we use ∆m 2 31 = 2.4 × 10 −3 eV 2 , ∆m 2 21 = 7.5 × 10 −5 eV 2 , E = 100 MeV. 2 A simple understanding of why it can be done in this way is to note the approximate degeneracy in maximum and the solar resonance are as prominent as the atmospheric ones, it must be more natural to pursue the line of solar-scale oscillation physics for the purpose of mixing parameter measurement. Whether right or wrong, it defines our physics motivation leading to the work described in this paper. It is the purpose of this paper to explore physics and possible advantage of measurement at the solar resonance affected region or around the solar oscillation maximum, for short, the region of "solar-scale oscillation". 4 We do it under a hope that the discovery potential in such region will be eventually fully explored experimentally. We must note, however, that this is our first attempt with "pedestrian spirit", in which we try to understand the characteristic features of neutrino oscillation at a qualitative level. Therefore, the readers must be aware that even when we address utility of solar-scale oscillation physics in the context of mixing parameter measurement, our statements will only be at a qualitative level.
We happily note, with full respect to the foregoing works, that such an approach to solar-scale oscillation physics are neither new nor unique. Notably, the extensive foregoing studies by Smirnov and the collaborators on this subject must be mentioned [26][27][28][29]. Yet, we may still miss many other relevant works. The authors of refs. [26][27][28][29] presented detailed studies of atmospheric neutrinos at low energies, and estimated, for example, how large is the effect of CP phase. Despite these existing massive works, we believe that physics of neutrino oscillations at around the solar oscillation maximum, in a broader sense, has a sufficient importance to be explored further. In fact, there are new ongoing researches toward this direction. Detection of low-energy atmospheric neutrinos by the large liquid scintillator detector is under investigation in the JUNO group with a merit of higher light output at low energies [30]. The similar low energy region is studied with a large liquid Ar TPC, and the sensitivity to CP phase is evaluated in the context of DUNE [31].
In this paper, toward understanding physics of the solar-scale oscillation, we formulate a perturbative framework by expanding around the solar-scale enhancement. More specifically, we aim at constructing perturbation theory of neutrino oscillation by targeting the region E a few × 100 MeV and L a few × 1000 km. In a sense it is a continuation of the proposal of low energy ∼100 MeV conventional superbeam to measure δ [32], the cleanest place for CP measurement, but in a new unexplored territory in that work.
What is the most interesting outcome of our perturbative treatment? We will show that CP violation is large at around the solar resonance, a factor of ∼ 10 larger than those expected at the conventional LBL neutrino experiments such as T2HK and DUNE. As is well known, the effect of CP phase δ, not only CP-odd but also CP-even effects, receives the two suppression factors, ≡ ∆m 2 21 /∆m 2 31 and the product of angles c 23 s 23 s 13 c 12 s 12 [33]. The secret of the factor ∼10 enhancement, which gives us a definite advantage in looking for the CP phase effect, is that the former suppression goes away due to the solarscale enhancement. This property makes atmospheric neutrinos at low energies, say E = 100 − 500 MeV, particularly interesting target for hunting lepton CP violation. We believe that our framework also has an interesting novel feature as a perturbation theory of neutrino oscillation. 5 It will be explained in section 3. We construct our perturbative framework, the "solar resonance perturbation theory", in section 4. The formulas for the neutrino oscillation probability is derived and their properties are elucidated in section 5. Numerical accuracy of the probability formula is examined in section 6, which leaves section 7 for the final remarks on our work.
2 Utility and limitation of perturbation theory of neutrino oscillation It appears to the authors that there is no sensible consensus in the community about the values of perturbative treatment of neutrino oscillation. In this section, we try to give a balanced view on utility and limitation of such treatment, even though it is inevitably based on our prejudices. We do it under a hope that forming a consensus on values of perturbative treatment eventually leads to a further development of the field.

Utility in data analyses
An extremist critical statement which is often made is that there is no need for perturbation theory of neutrino oscillation because the three-flavor neutrino evolution is so simple to solve numerically even with handy personal computers nowadays. In fact, it seems that most of the numerical analyses of the experimental data are done in this way, producing a feeling that numerical computation suffices in the data analyses. 6 Another problem in using the analytic formulas in data analyses would be that they often rely on idealistic approximations such as uniform matter density approximation in the earth. Then, the users cannot be confident about numerical accuracies of the results.

Utility in qualitative understanding of neutrino oscillation
We would like to emphasize our view that utility of perturbative formulation of neutrino oscillation exists in a completely different realm, offering the way of understanding the phenomenon of neutrino oscillation at a qualitative level.
The best example showing this feature is the understanding of parameter degeneracy offered by the approximate formulas of the oscillation probability derived by Arafune et al. [34], and Cervera et al. [35]. 7 The former (the latter) applies to the case of weak (sizeable) matter effects in the LBL accelerator neutrino experiments. Using the formulas it has been recognized that there exists the three sets of parameter degeneracy, the phenomenon 5 Utility of perturbation theory of neutrino oscillation is, however, often an issue of controversy. Since this is a quite involved subject to discuss here in Introduction, we defer it to the next section 2 which is devoted to this specific topics. 6 The situation, however, might change when we have to deal with a huge number of systematic errors and nuisance parameters in precision era of the neutrino experiments, as already hinted, e.g., in Super-Kamiokande atmospheric neutrino analysis. 7 Another examples to which we do not enter into the details is that an explicit computation of probabilities reveals new regularity such as correlation between the mixing parameters. A quite limited list includes the one between θ13 and the NSI parameters [36], and the one between νSM CP phase and the phases which represent the effects of new physics [37].
of multiple solutions for a given set of observed probabilities. They are, so called, the intrinsic θ 13 −δ degeneracy [38], the sign-∆m 2 31 degeneracy (which implies different solutions depending upon the mass orderings) [20], and the θ 23 − δ degeneracy [39]. Even with exact analytic expressions of the oscillation probability derived by Kimura et al. (KTY) [40], and by Zaglauer and Schwarzer (ZS) [41], it would have been difficult to recognize such structure of degeneracy unless these structure-transparent approximate formulas were known.
Notice that the Cervera et al. formula can be derived by perturbative treatment to second order in the small parameters, ≡ ∆m 2 21 /∆m 2 31 and s 13 . The formula is not super-accurate, since the measured value of θ 13 is large, s 2 13 = 0.022 ∼ = 0.031 (see e.g., refs. [42][43][44][45]), and it requires higher order corrections to s 4 13 to match for secondorder accuracy in [33]. Nonetheless, it did not affect utility of the formula to unfold the structure of parameter degeneracy. The experience tells us that numerical accuracy of the perturbative formula is not the whole story.

Utility in having framework of systematic inclusion of sub-leading effects
Another point for utility of perturbative formulation of the standard three-flavor neutrino oscillation is that it offers a systematic way of including the sub-leading effects beyond the νSM.
One might argue that since such effect is small in size one can perturb with the small sub-leading effect itself. However, it is not practical even for constant matter density, given the complexity of the KTY or ZS solutions as the leading term. Therefore, one need a suitable perturbative framework of three neutrino evolution to include beyond the νSM subleading effect. Examples of the effects under discussion include non-standard interactions (NSI) [46,47], and non-unitarily due to possible new physics at high-or low-energy scales [48][49][50][51]. Perturbative treatments with such effect were tried, e.g., in [36] and [37] for NSI and non-unitarily, respectively. We note that having formalism for including such subleading effects systematically may become more important as accuracy of measurement further improves.

Limitation of perturbative formula
Yet, we must keep in mind that perturbative expression of the oscillation probability has many limitations. One of them, which is most important in our opinion is that the kinematical region, e.g., in energy E and baseline L space, of validity of the formula is quite restricted depending upon how the perturbation theory is formulated.
An example of this feature is that perturbation theory to treat neutrino oscillation in region around the "atmospheric-scale resonance" fails to make proper treatment of the eigenvalue level crossing in the "solar-scale resonance" region [52], and vice versa. 8 However, we must emphasize that it is perfectly normal in perturbation theory, nothing wrong with it by itself. Nevertheless, it must be kept in mind that we must know the region of validity of each perturbative framework. 8 They serve for understanding the features of neutrino oscillations in LBL accelerator or atmospheric neutrino experiments with detector at around the atmospheric oscillation maximum ∆m 2 31 L/4E = π/2 (e.g., [52]), and the solar oscillation maximum ∆m 2 21 L/4E = π/2 (this paper), respectively.
It may be appropriate to mention here that there is an approach to remedy the problem of perturbative treatment described above, as initiated by the authors of ref. [53], who performed a 1-2 space rotation first and then a 1-3 space rotation to approximately diagonalize the Hamiltonian. It was followed by Denton et al. who did the two rotations in the opposite order [54]. Here, we do not discuss this approach further, regarding it a method for approximate diagonalization of Hamiltonian rather than perturbation theory, because the leading order term itself has a sufficient accuracy for most purposes.

Summary
To summarize our discussion in this section, with full of the prejudices, we quote the following three statements: • A perturbative formula of neutrino oscillation probability may be useful for data analyses, if it is reasonably simple and accurate, and can handle a realistic situation such as adiabatically varying matter density.
• A perturbative formula can reveal qualitative features of neutrino oscillation as illustrated by an example of how the phenomenon of parameter degeneracy was uncovered.
• Perturbation theory has an inherent limitation that the applicability region of the formula is restricted to a particular region depending on how it is formulated. It is crucially important to know the region of validity.

What is new in the solar resonance perturbation theory?
First of all, we want to place our perturbative framework, "solar resonance perturbation theory", into a context in our discussion in the previous section. It is a perturbation theory formulated to treat the solar resonance region in the scope of atmospheric neutrino observation at low energies, typically E < ∼ 500 MeV, or possibly future accelerator LBL experiments in the similar energy region. It was partly motivated as an complement to "atmospheric resonance perturbation theory" [52], which perturbs around atmosphericscale enhancement. See also refs. [35,[55][56][57] for earlier references. Now, what is new in our perturbation theory as a framework itself? First, we must note that this is a somewhat counterintuitive way of formulating perturbation theory. We aim at perturbing the neutrino evolution by taking the solar 1-2 sector as the non-perturbative part of the Hamiltonian. A potential issue here is that ∆m 2 31 is larger than ∆m 2 21 by a factor of ∼ 30, which means that there is much bigger effect outside the 1-2 zeroth order part of the Hamiltonian. It necessitates a formulation of perturbation theory in such a way that the much larger ∆m 31 term does not disturb the zeroth-order effect. In fact, we will do a much better job, which is done in the following two ways: In leading order in expansion the ∆m 2 31 term decouples, and in the first-order corrections, the ∆m 2 31 shows up only in energy denominators. Consequently, we have the effective expansion parameter A exp = c 13 s 13 a/∆m 2 31 10 −3 for E = 100 MeV [see (4.28)], where a denotes the matter potential. That is, the larger the ∆m 31 , the smaller the higher order corrections.
Then, what is new in the outcome of our perturbative formulation in physics context? The most important feature we can clearly observe by analytic means is that CP violation effect is much larger than the one we expect to see in the conventional accelerator LBL experiments, such as T2HK and DUNE. It is larger by a factor of ∼10 at the level of the appearance oscillation probability, the fact best illustrated by using the bi-probability plot [20] as shown in figure 1. The mechanism behind such large amplification is that one of the suppression factors ∝ ∆m 2 21 /∆m 2 31 on δ dependent terms is dynamically lifted, as will be shown in section 5.3.
In consistent with our own suggestion we examined the region of applicability of our solar resonance perturbation theory. We identified, as a necessary condition (see section 6), the region of validity of our formula as L > ∼ atm = 4πE/∆m 2 31 . The smallness of the first order correction terms that comes from the large hierarchy ∆m 2 31 ∆m 2 21 can be confirmed as a small corrections to the zeroth-order oscillation probability, as will be seen in figure 3. To our surprise we have found that the oscillation probability formula works in a reasonable accuracy for L = 300 km and E < ∼ 1 GeV, suggesting its applicability to the T2K and T2HK experiments.

Formulating perturbation theory of neutrino oscillation around the solar resonance
The standard three-flavor neutrino evolution in matter can be described by the Schrödinger equation in the flavor basis, i d dx ν = Hν, with Hamiltonian where E is neutrino energy and ∆m 2 ji ≡ m 2 j − m 2 i . In (4.1), U denotes the standard 3 × 3 MNS lepton flavor mixing matrix which relates the flavor neutrino states to the vacuum mass eigenstates as ν α = U αi ν i , where α runs over e, µ, τ , and the mass eigenstate indix i runs over 1, 2, and 3. With the obvious notations s ij ≡ sin θ ij etc. and δ for lepton version of Kobayashi-Maskawa CP violating phase, U is given by where we have defined U matrix in a convention used in [52], which is slightly different (but physically equivalent) from U PDG of Particle Data Group [42].
The functions a(x) in (4.1) denote the Wolfenstein's matter potential [8] due to charged current (CC) reactions Here, G F is the Fermi constant, N e is the electron number density in matter. ρ and Y e denote, respectively, the matter density and number of electron per nucleon in matter. For simplicity and clarity we will work with the uniform matter density approximation in this paper. But, it is not difficult to extend our treatment to varying matter density case if adiabaticity holds.

Relevant kinematical region
We are interested in exploring the region where neutrino energy E = (1 − 5) × 100 MeV and baseline L = (1 − 10) × 1000 km. We note that the region is around the vacuum solar oscillation maximum, (4.4) But, since the matter potential is comparable to the vacuum effect (represented by ∆m 2 21 ) in this region, the atmospheric neutrinos at low energies, ∼ a few ×100 MeV, are fully affected by the earth matter effect. We will see later that the role of effective expansion parameter is played by the quantity A exp = c 13 s 13 a/∆m 2 31 10 −3 for E = 100 MeV [see (4.28)].

Intermediate basis, or the check basis
To formulate the solar resonance perturbation theory we transform from the flavor basis to an intermediate basis, or the check basiš with Hamiltoniaň where U 23 rotation on the second term is performed with no effect. Here, we have introduced simplified notations (i, j = 1, 2, 3) Using the parametrization of U matrix in (4.2)Ȟ can be written aš

Formulating perturbation theory with hat basis
We use the "renormalized basis" such that the zeroth-order and the perturbed Hamiltonian takes the formȞ =Ȟ 0 +Ȟ 1 : We then transform to the "hat basis" where U ϕ is parametrized as U ϕ is determined such thatĤ is diagonal. The condition reads cos 2ϕ = cos 2θ 12 − c 2 13 r a cos 2θ 12 − c 2 13 r a 2 + sin 2 2θ 12 , sin 2ϕ = sin 2θ 12 cos 2θ 12 − c 2 13 r a 2 + sin 2 2θ 12 , (4.14) where The eigenvalues of the zeroth order HamiltonianȞ 0 in (4.10) is given by Notice that one can show that Then, the Hamiltonian in the hat basis is given byĤ =Ĥ 0 +Ĥ 1 wherê where we should note thatĤ 1 =Ȟ 1 .
If we treat anti-neutrino channel, the transformations δ → −δ and a → −a suffice. We denote the angle ϕ and the eigenvalues h i transformed by the above transformations asφ andh i . Notice that our framework is constructed such as to allow both the normal (∆m 2 31 > 0) and the inverted (∆m 2 31 < 0) mass orderings. Moreover, h 1 , h 2 and ϕ do not depend on the mass orderings because ∆m 2 21 > 0, whereas h 3 does. Therefore, unlike the case of helio-to-terrestrial ratio ( ≡ ∆m 2 21 /∆m 2 31 ) perturbation theory, no essential classification of the mass ordering is necessary in our present formalism.

Calculation ofŜ matrix
We formulate perturbation theory in the "hat" basis. TheŜ matrix is given bŷ where T symbol indicates the "space ordering", but it simplifies for uniform matter density as shown in the second equality in eq. (4.19).
To calculateŜ(x) we define Ω(x) as (4.20) Using i d dxŜ =Ĥ(x)Ŝ, Ω(x) obeys the evolution equation where Then, Ω(x) can be computed perturbatively as and theŜ matrix is given byŜ Noticing that H 1 is given by Ω(x) can be calculated by using (4.23) with (4.26). By using (4.24),Ŝ matrix is given to first order in H 1 bŷ (4.27) Notice that the large element outside "solar" 1 − 2 sector ofȞ 0 appears solely in the energy denominator. 9 Therefore, as we mentioned earlier, the effective expansion parameter becomes which is very small. It gives us a feeling that the leading term in the expansion must give a fair approximation, which will be confirmed in section 6. Notice also thatŜ respects T invariance,Ŝ ij =Ŝ ji , as it should.

Calculation of S matrix and the oscillation probability
Since the relationship between the flavor basis and the hat basis is given by the S matrix in flavor basis can be readily calculated as In appendix A, we briefly describe the computation and give the resulting expressions of S matrix elements to first order in expansion. Then, the oscillation probability P (ν β → ν α ) at baseline x is given by Give the expressions of S matrix elements, it is straightforward to compute the oscillation probability to first order in expansion. All the expressions of the oscillation probability except for the ones which require ν τ beam are given either in the next section 5, or in appendix B. 10

The oscillation probability and its symmetry
In this section we present the explicit form of the zeroth-order oscillation probability in the ν µ → ν e and ν µ → ν µ channels. To simplify the expression we define the reduced Jarlskog factor in matter J mr ≡ c 23 s 23 c 2 13 s 13 c ϕ s ϕ = J r cos 2θ 12 − c 2 13 r a 2 + sin 2 2θ 12 which is proportional to the reduced Jarlskog factor in vacuum, J r ≡ c 23 s 23 c 2 13 s 13 c 12 s 12 . We have used eq. (4.14) in the second equality in (5.1). Notice that not only sin δ but also cos δ terms in the oscillation probability in the ν e -related sector (ν µ − ν τ sector) must be proportional to J r (J r /c 2 13 ), as proved in [33], and our expressions are consistent with this result. We use the notation x as baseline distance in this section, as well as in appendices A and B.
What is a characteristic difference between the zeroth-order oscillation probability P (ν µ → ν e ) (0) and the one around the atmospheric oscillation maximum? For the latter, if necessary, we refer the expression in ref. [52] for definiteness. There are at least two important differences: • The main term of the oscillation probability is proportional to c 2 23 in our solarresonance perturbation theory [see the first term in (5.2)], while it is proportional to s 2 23 in atmospheric-resonance perturbation theory. This property will play crucial role in resolving the θ 23 − δ degeneracy. A brief discussion of the θ 23 − δ degeneracy with our formula is given in appendix C.
• In our solar-resonance perturbation theory the relevant dynamical variables at the solar-scale enhancement, the eigenvalues h 1 , h 2 , and ϕ, depend on ∆m 2 21 but not on ∆m 2 31 or ∆m 2 32 . It means that there is no essential difference in important features between the normal and the inverted mass orderings. It will simplify the treatment of parameter degeneracy because the sign-∆m 2 31 ambiguity essentially decouples.

Symmetry of the oscillation probability
One recognizes that a symmetry exists in eqs. (5.2) and (5.3). That is, P (ν µ → ν e ) (0) and P (ν µ → ν µ ) (0) are both invariant under the transformation In fact, it is easy to confirm that the symmetry exists in all the oscillation channels. The invariance can be understood as the one under the transformation s ϕ → +c ϕ , c ϕ → −s ϕ , or at the ϕ level under which produces the eigenvalue exchange h 1 ↔ h 2 owing to the relations sin(ϕ − θ 12 ) → cos(ϕ − θ 12 ) and cos(ϕ − θ 12 ) → − sin(ϕ − θ 12 ) under the transformation (5.5). It is interesting to note that the symmetry involves ϕ, the matter-affected mixing angle θ 12 which descries the 1-2 level crossing in matter. Since the transformation (5.5) describes a reparametrization of ϕ, nature of the symmetry can be understood as a "dynamical symmery", not a symmetry of Hamiltonian. Nonetheless, recognizing it is useful for a consistency check of the calculation.
In fact, one can show that a similar symmetry exists in the expression of the oscillation probability in "atmospheric resonance" perturbation theory formulated in [52]. See appendix D for details.

CP violation around the solar resonance
The most important message in this paper is that CP violation effect is large in the solar oscillation enhanced region, say E a few × 100 MeV and L a few × 1000 km, for which our perturbative framework serves. Roughly speaking, it is ∼10 times larger than the CP odd term in the oscillation probability estimated for the ongoing and the upcoming LBL experiments [10,12,18,19].
To show the point, let us first look at the CP-odd sin δ term in P (ν µ → ν e ) in vacuum. It takes the form as 1, the size of sin δ term is given by 8J r sin 2 ∆ 31 x 2 4J r , where we take average over the fast atmospheric scale oscillations. Therefore, the magnitude of CP odd term at around the solar oscillation maximum is 15 times larger than that at the atmospheric oscillation maximum. It is the latter quantity that is planned to be measured in the future LBL experiments. Therefore, in a nutshell, the reason why the effect of CP phase δ is large at around the solar-scale enhanced region is that the factor sin ∆ 21 x 2 , which would become a suppression factor ∆ 21 x 2 ∆m 2 21 /∆m 2 31 = at around the atmospheric oscillation maximum ∆ 31 x 2 1, no more acts as a suppression factor, sin ∆ 21 x 2 1. That is, the suppression of the CP-odd term due to small is dynamically lifted by the solar-scale enhancement.
In practice, such neutrino beam passes through inside the earth, for which the matter effect must be taken into account. Here, we use the zeroth-order oscillation probability P (ν µ → ν e ) (0) in (5.2) to give a rough estimation of size of the δ dependent term. To illuminate the point, we take average over the atmospheric-scale short wavelength oscillations. 11 Then, we obtain for P (ν µ → ν e ) (0) P (ν µ → ν e ) (0) = 1 2 s 2 23 sin 2 2θ 13 + c 2 13 sin 2 2ϕ c 2 23 − s 2 23 s 2 13 sin 2 (h 2 − h 1 )x 2 This is essentially the two-flavor oscillation probability in matter near the resonance. We roughly estimate the value of the sin δ term around the solar resonance. For simplicity, let us tune the energy and baseline such that it is on resonance, cos 2θ 12 −c 2 13 r a = 11 We use 0, or ϕ = π 4 , and sin(h 2 − h 1 )x is maximal, (h 2 − h 1 )x = π/2. We use the best fit values of the mixing parameters given in [43] for the normal mass ordering, which leads to sin 2θ 12 = 0.925 J r ≡ c 23 s 23 c 2 13 s 13 c 12 s 12 = 0.0334, and J mr = J r / sin 2θ 12 = 0.0361. Then, the magnitude of the sin δ term at the solar resonance is given by 2J mr = 2J r / sin 2θ 12 = 0.0722 (5.9) while the cos δ term vanishes. Apart from cos 2ϕ which vanishes at this particular point, the solar resonance point ϕ = π 4 , the magnitude of the cos δ term is also 2J mr , the same as the sin δ term given in (5.9). Therefore, apart from the particular point, generically the CP-even cos δ term has the same order of magnitude as the CP-odd term.
We want to compare the value of CP odd term in (5.9) to the magnitudes of the sin δ term to be seen in the future LBL experiments such as T2K/T2HK, NOνA, and DUNE [10,12,18,19], whose set up are close to the atmospheric oscillation maximum. We use, for a rough comparison, the value calculated above assuming vacuum oscillation, 8 J r 0.0089, which gives a fair, albeit crude, estimation. 12 Therefore, the value of CP odd term at the solar resonance in (5.9) is 8 times larger than the CP odd term in vacuum estimated at the atmospheric oscillation maximum. One must keep in mind that this result is obtained by using (5.8), which is obtained by averaging over the short wavelength atmospheric-scale oscillations. Since its magnitude easily reaches to a few tenth of the averaged probability it is difficult to obtain the unique value for the amplification factor. Therefore, it is fair to say that it is ∼ 10. 13 The contrast of the sizes of CP phase dependence between the solar-scale oscillation enhanced region and the conventional LBL setting, NOνA as just an example, can be best represented by the bi-probability plot in P (ν µ → ν e ) − P (ν µ →ν e ) space [20], as seen in figure 1. In this example the ellipses in the solar-scale enhanced region are about 10 times larger than NOνA's. Notice that the features of bi-probability ellipse in the solar-scale enhanced region, in particular their shapes, are extremely sensitive to the values of the mixing parameters, in particular ∆m 2 31 [60], due to coexistence of the atmospheric-and the solar-scale oscillations. However, we do not enter into the complexity in this paper, restricting ourselves to a simple but conservative estimation of the amplification factor of CP phase effect.
The similar discussion on the size of the CP-odd term can be easily extended to that of the ν e → ν µ channel, because they have the same magnitude in vacuum and in matter due to T-invariance. Similarly, the same argument must go through for the ν e → ν τ channel, because the CP-odd term has the same size as that of the ν e → ν µ channel in vacuum and in matter due to unitarity. 12 The matter effect does not affect so much the sizes of δ dependent terms. The simplest way to see this is to look into the neutrino-antineutrino bi-probability plot [20]. The matter potential shifts the CP ellipse to right-down (left-up) direction in the normal (inverted) mass ordering, but keeping its size almost as it is for relatively short baseline L ∼ 1000 km. 13 Another example of amplification of the δ dependent terms is the one at the second oscillation maximum [58,59]. One may regard this amplification of a factor of ∼3 is due to "imperfect growth" of the solar-scale oscillation.  Figure 1: The bi-probability diagrams in P (ν µ → ν e )−P (ν µ →ν e ) space [20] are presented for the solar-scale enhanced region (large ellipses); E = 0.1 GeV and L = 2000 km and for NOνA setting (small ellipses); E = 2 GeV and L = 810 km. The blue and red lines are for the normal and inverted mass orderings, respectively. The values of the mixing parameters are taken from ref. [43], i.e., the best fit values for the normal mass ordering.
6 How accurate are the oscillation probability formulas?
In this section we investigate numerically the oscillation probability formulas derived by the solar resonance perturbation theory to know how accurate they are. Though the numerical accuracy is not our primary concern, as we discussed in section 2, it is better to check how good is our expectation to the accuracy we spelled out in section 4. We examine first overall (dis-) agreement between our leading order formula and the exact expression of P (ν µ → ν e ), which is presented in figure 2. Then, in figure 3, we display accuracies of the formulas in more details by comparing our leading order and the leading plus firstorder formulas to the exact energy dependences of P (ν µ → ν e ) and P (ν µ → ν µ ) at several baseline distances.
In the left panel of figure 2, plotted are the equi-probability contours of P (ν µ → ν e ) in energy E and distance L plane, displayed by color-graded regions as well as the contours superimposed on it. It presents a global view of P (ν µ → ν e ) which serves for the pedagogical discussions as given in section 1. It illuminates, in particular, the two enhanced regions of the solar-and the atmospheric-scale MSW enhancement. Whereas in the right panel of figure 2, presented is the ratio P (ν µ → ν e ) (0) − P (ν µ → ν e ) exact /P (ν µ → ν e ) exact to show how good is the leading order formula. Or, more precisely, it is meant to be a measure for revealing where is the applicability region of our framework. One can see that the region of agreement, the left-sided yellow colored region, essentially coincides with the region of solar-scale enhancement, the target region of our perturbative treatment, as it should be.
Though not shown, the agreement is much better in the disappearance channels, Figure 2: In the left panel, the oscillation probability P (ν µ → ν e ) is shown as the equiprobability contours superimposed on a color grade plot in neutrino energy E and baseline L plane. In the right panel, the ratio of deviation of the zeroth-order formula to the exact probability, P (ν µ → ν e ) (0) /P (ν µ → ν e ) exact − 1, is plotted in the same style. In both panels the matter density is taken to be a constant, ρ = 4.0 g/cm 3 . The values of the mixing parameters are taken from the best fit values of ref. [43], including δ = 215 degree.
P (ν e → ν e ) and P (ν µ → ν µ ), < ∼ 1% deviation from the exact results. For P (ν e → ν τ ) the agreement is slightly better than P (ν µ → ν e ). While the region of validity of our framework essentially traces the resonance region, there is a visible departure, or widening, of the (yellow-white) validity region toward shorter baseline to L ∼ 100 km at low energies E < ∼ 400 MeV. It has an interesting consequence that our formula may be applicable for shorter baseline, e.g., the T2K/T2HK baseline, L = 300 km. Then, we should try to understand theoretically what would be the condition for validity of our formula.
Our approach toward this understanding is still "empirical" one, which starts from the right panel of figure 2. By approaching from the higher-energy right blue region, we encounter a line which can roughly be expressed as for the region of validity. That is, at distance scale much larger than atm , the neutrino evolution is insensitive to short-wavelength atmospheric oscillations, and hence our perturbative formula applies. On the contrary, in the alternative region L E < ∼ 4π ∆m 2

31
, the energy scale is so high for neutrinos such that it can probe the ∆m 2 31 effect, which renders approximation of our perturbative theory inaccurate. For this reason our leading order perturbative formula ceases to be a good approximation in region L < ∼ atm .
ν α ) (1) (α = e, µ) are given by the green dotted and the red dashed lines, respectively. Whereas, the results of exact calculation are given by the blue solid lines. In each cluster in figure 3, the four panels are for distances L = 300 km, 1000 km, 2000 km, and 10000 km. Both in the appearance and disappearance channels, P (ν µ → ν e ) and P (ν µ → ν µ ), the red dashed lines (the zeroth plus first order formulas) perfectly overlap with the blue solid lines, the exact results, to 1% level or below in the region displayed in figure 3. In the appearance channel the zeroth-order formula shows some deviation in region E > ∼ 300 MeV. But, in the disappearance channel even the zeroth-order formula overlaps well with the exact result. The careful readers might have detected a tiny deviation of the first order formulas from the exact result at the peaks and at the bottoms of dips at relatively higher energies, E > ∼ 500 MeV. Therefore, our perturbative framework works as expected in the target region of our formulation, E a few × 100 MeV and L a few × 1000 km. We also found that our first order formula P (ν µ → ν e ) (0) + P (ν µ → ν e ) (1) works very well at low energies for the T2K/T2HK baseline, L = 300 km.

Concluding remarks
In this paper, we have constructed a perturbative framework, the "solar resonance perturbation theory", whose leading order describes the solar-scale MSW enhancement in a good approximation. We aim at illuminating physics of atmospheric or accelerator neutrinos for the regions E a few × 100 MeV and L a few × 1000 km. Despite the keen interests expressed by people and good amount of works done [26][27][28][29], it appears to us that experimental exploration of the physics in this region has not been performed in a sufficient depth. We hope that our discussion at a qualitative level may be useful to trigger interests of wider class of people in the community.
The clearest message which is born out from our treatment is that the effect of CP violating phase is large, larger by a factor of ∼10 than that expected in the conventional LBL neutrino experiments such as T2HK and DUNE, which utilize the atmospheric oscillation maximum. The CP phase effect is large at around the solar resonance or oscillation maximum because there is no suppression factor ∆ 21 x 2 ∆m 2 21 /∆m 2 31 = which does exist at around the atmospheric oscillation maximum. Instead, sin ∆ 21 x 2 1 at around the solar oscillation maximum.
Probably, the most practical way of exploring the region of solar resonance / oscillation maximum is to observe atmospheric neutrinos at low energies, ideally E < ∼ 500 MeV. This possibility is becoming more and more feasible because of construction of massive detectors such as Hyper-K [18], DUNE [19], JUNO [24] (see [30]), and INO [62]. It is worth to watch to what extent the energy threshold of gigantic size neutrino detectors, IceCube-Gen2/PINGU and KM3NeT-ORCA [22,23], can be lowered. Yet, the major issues here seems to be how to improve accuracies of the atmospheric neutrino flux prediction, and reduce uncertainties in cross section measurement [63]. In this context, it is very interesting to know precisely the performance of the large-volume liquid scintillator detector JUNO for the detection of low energy atmospheric neutrinos [30].
Our framework and the resulting oscillation formulas based on constant matter density approximation might have more direct relevance if accelerator LBL neutrino experiment is a viable possibility to explore this region. It is not only because the uniform density approximation may apply but also because it enables one to fully enjoy maximal CP violating effect by tuning the beam energy and baseline to the solar-scale enhanced region. Probably, a possible setup closest to the reality may be to use JPARC neutrino beam pointed to the second Hyper-K detector in Korea, the T2HKK project [21], which has the baseline L 1100 km. We note that more than 1/3 of the ν e appearance signals are in low energy bins, E < ∼ 600 MeV. See Fig. 13 in ref. [21]. Therefore, it might be worthwhile to look at this possibility in detail. Interestingly, our first-order perturbative formula for the oscillation probability works very well for T2K/T2HK and T2HKK setups in the relevant region of energy and baseline, as we saw in section 6.
We are aware that our methodological view for determination of the remaining unknowns in the lepton mixing and the neutrino mass pattern described in section 1 is too biased. For example, JUNO [24] tries to determine the mass ordering without recourse to the matter effect. Before reactor neutrino measurement of θ 13 became popular, the majority of people thought about using accelerator neutrino beam to determine θ 13 , and essentially nobody talked about precision measurement of θ 13 by reactor neutrinos [64]. This example illuminates the danger we might have when we rely too much on the common knowledge, or a prejudice. In this sense it is worthwhile, we believe, to develop alternative methods for determining the unknowns in lepton flavor mixing. It also should have merit in the context of paradigm test (see e.g., [48][49][50][51]), given that the effect of unitarity violation is large in the target region of this paper [51]. The direction we described in this paper is just one of the alternative possibilities.

A Flavor basis S matrix
For convenience of nomenclature we define another intermediate basis, the "tilde" basis. The definition ofS matrix and its relation to S matrix are given bỹ Using (4.27), the explicit expressions ofS matrix elements are given bỹ Then, finally, U 23 rotation can be performed to convert the above expression ofS αβ to S matrix elements:

B The expressions of neutrino oscillation probability
The oscillation probability can be calculated as We calculate P (ν β → ν α ; x) to first order in the effective expansion parameter, and the results are presented in the following way: The explicit expressions of the leading-order term, P (ν β → ν α ) (0) , and the first-order corrections P (ν β → ν α ) (1) are given below.
If precision of the measurement of P (ν µ → ν e ) at the solar resonance region and in the usual LBL setting becomes good enough, whose latter is described in ref. [61], then the most accurate way of resolving the θ 23 − δ degeneracy would be to combine these two sets of measurements. That is because of c 2 23 and s 2 23 dependence of the main oscillation terms of the perturbative formulas around the solar-and atmospheric-resonance regions, respectively. Of course, the quantitative analysis is called for to justify this statement.

D Symmetry of the oscillation probability in "Simple and Compact" formulas
We refer appendix B in ref. [52] for explicit expressions of the oscillation probabilities in the "atmospheric resonance" perturbation theory. It is then easy to recognize that there is a symmetry under the transformations involving angle φ as (see [52] for notations) which may be summarized as The angle φ is the matter-dressed θ 13 angle, which describes the 1-3 level crossing in matter. The nature of the symmetry is the same as ϕ → ϕ + π 2 symmetry in the "solar resonance" perturbation theory, as described in section 5.2.