Neutrino Propagation in Matter

We describe the effects of neutrino propagation in the matter of the Earth relevant for experiments with atmospheric and accelerator neutrinos and aimed at the determination of the neutrino mass hierarchy and CP-violation. These include (i) the resonance enhancement of neutrino oscillations in matter with constant or nearly constant density, (ii) adiabatic conversion in matter with slowly changing density, (iii) parametric enhancement of oscillations in a multi-layer medium, (iv) oscillations in thin layers of matter. We present the results of semi-analytic descriptions of flavor transitions for the cases of small density perturbations, in the limit of large densities and for small density widths. Neutrino oscillograms of the Earth and their structure after determination of the 1-3 mixing are described. A possibility to identify the neutrino mass hierarchy with the atmospheric neutrinos and multi-megaton scale detectors having low energy thresholds is explored. The potential of future accelerator experiments to establish the hierarchy is outlined.


I. INTRODUCTION
Neutrinos are eternal travelers: once produced (especially at low energies) they have little chance to interact and be absorbed. Properties of neutrino fluxes: flavor compositions, lepton charge asymmetries, energy spectra of encode information. Detection the neutrinos brings unique knowledge about their sources, properties of medium, space-time they propagated as well as on neutrinos themselves.
Neutrino propagation in matter is vast area of research which covers variety of different aspects: from conceptual ones to applications. This includes propagation in matter (media) with (i) different properties (unpolarized, polarized, moving, turbulent, fluctuating, with neutrino components, etc), (ii) different density profiles, and (iii) in different energy regions. The applications cover neutrino propagation in matter of the Earth and the Sun, supernova and relativistic jets as well as neutrinos in the Early Universe.
The impact of matter on neutrino oscillations was first studied by Wolfenstein in 1978 [1]. He marked that matter suppresses oscillations of the solar neutrinos propagating in the Sun and supernova neutrinos inside a star. He considered a hypothetical experiments with neutrinos propagating through 1000 km of rock, something that today is no longer only a thought but actual experimental reality. Later Barger et al [2] have observed that matter can also enhance oscillations at certain energies. The work of Wolfenstein was expanded upon in papers by Mikheev and Smirnov [3][4][5], in particular, in the context of the solar neutrino problem. Essentially two new effects have been proposed: the resonant enhancement of neutrino oscillations in matter with constant and nearly constant density and the adiabatic flavor conversion in matter with slowly changing density. It was marked that the first effect can be realized for neutrinos crossing the matter of the Earth. The second one can take place in propagation of solar neutrinos from the dense solar core via the resonance region inside the Sun to the surface with negligible density. This adiabatic flavor transformation, called later the MSW effect, was proposed as a solution of the solar neutrino problem.
Since the appearance of these seminal papers, neutrino flavor evolution in background matter were studied extensively including the treatment of propagation in media which are not consisting simply of matter at rest, but also backgrounds that take on a more general form. For instance, in a thermal field theory approach [6], effects of finite temperature and density can be taken readily into account. If neutrinos are dense enough, new type of effects can arise due to the neutrino background itself, causing a collective behavior in the flavor evolution. This type of effect could have a significant impact on neutrinos in the early Universe and in central parts of collapsing stars.
There has been a great progress in treatments of neutrino conversion in matter, both from an analytical and a pure computational points of view. From the analytical side, the description of three-flavor neutrino oscillations in matter is given by a plethora of formulas containing information that may be hard to get a proper grasp of without introducing approximations. Luckily, given the parameter values inferred from experiments, various perturbation theories and series expansions in small parameters can be developed. In this review we will explain the basic physical effects arXiv:1306.2903v1 [hep-ph] 12 Jun 2013 important for the current and next generation neutrino oscillation experiments and provide the relevant formalism. We present an updated picture of oscillations and conversion given the current knowledge on the neutrino oscillation parameters.
In this paper we focus mainly on aspects related to future experiments with atmospheric and accelerator neutrinos. The main goals of these experiments are to (i) establish the neutrino mass hierarchy, (ii) discover CP-violation in the lepton sector and determination of the CP-violating phase, (iii) precisely measure the neutrino parameters, in particular, the deviation of 2-3 mixing from maximal, and (iv) search for sterile neutrinos and new neutrino interactions.
Accelerator and atmospheric neutrinos propagate in the matter of the Earth. Therefore we mainly concentrate on effects of neutrino propagation in the Earth, i.e., in usual electrically neutral and non-relativistic matter. We update existing results on effects of neutrino propagation in view of the recent determination of the 1-3 mixing.
The review is organized as follows: In Sec. II we consider properties of neutrinos in matter, in particular, mixing in matter and effective masses (eigenvalues of the Hamiltonian); we derive equations which describe the propagation. Sec. III is devoted to various effects relevant for neutrino propagating in the Earth. We consider the properties of the oscillation/conversion probabilities in different channels. In Sec. IV we explore the effects of the neutrino mass hierarchy and CP-violating phase on the atmospheric neutrino fluxes and neutrino beams from accelerators. Conclusions and outlook are presented in Sec. V.

II. NEUTRINO PROPERTIES IN MATTER
We will consider the system of 3 flavor neutrinos, ν T f ≡ (ν e , ν µ , ν τ ), mixed in vacuum: Here U P M N S is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix [7][8][9] and ν T m ≡ (ν 1 , ν 2 , ν 3 ) is the vector of mass eigenstates with masses m i (i = 1, 2, 3). We will use the standard parameterization of the PMNS matrix, which is the most suitable for describing usual matter effects. In Eq. (2) U ij (θ ij ) are the matrices of rotations in the ij−planes with angles θ ij and I δ ≡ diag(1, 1, e δ ).
In vacuum the flavor evolution of these neutrinos is described by the the Schrödinger-like equation where M is the neutrino mass matrix in the flavor basis and E is the neutrino energy. Eq. (3) is essentially a generalization of the equation E ≈ p + m 2 /2E for a single ultra relativistic particle. According to Eq. (3), the Hamiltonian in vacuum can be written as where M 2 diag ≡ M † M = diag m 2 1 , m 2 2 , m 2 3 and we take the masses m i to be real 1 .

A. Refraction and matter potentials
The effective potential for a neutrino in medium V f can be computed as a forward scattering matrix element V f = Ψ|H int |Ψ . Here Ψ is the wave function of the system of neutrino and medium, and H int is the Hamiltonian of interactions.
At low energies, the Hamiltonian H int is the effective four fermion Hamiltonian due to exchange of the W and Z bosons: γ µ (1 − γ 5 )ν {ēγ µ (g V + g A γ 5 )e +pγ µ (g p V + g p A γ 5 )p +nγ µ (g n V + g n A γ 5 )n} , where g V and g A are the vector and axial vector coupling constants. In the Standard Model the matrix of the potentials in the flavor basis, is diagonal: V f = diag(V e , V µ , V τ , 0...). For medium the matrix elements of vectorial components of vector current are proportional to velocity of particles of medium. The matrix elements of the axial vector current are proportional to spin vector. Therefore for non-relativistic and unpolarized medium (as well as for an isotropic distribution of ultra relativistic electrons) only the γ 0 component of the vector current gives a non-zero result, which is proportional to the number density of the corresponding particles. Furthermore, due to conservation of the vector current (CVC), the couplings g p V and g n V can be computed using the neutral current couplings of quarks. Thus, taking into account that, in the Standard Model, the neutral current couplings of electrons and protons are equal and of opposite sign, the NC contributions from electrons and protons cancel in electrically neutral medium. As a result, the potential for neutrino flavor ν a is where n e and n n are the densities of electrons and neutrons, respectively. Only the difference of potentials has a physical meaning. Contribution of the neutral current scattering to V is the same for all active neutrinos. Since V a (a = µ, τ or a combination thereof) is due to the neutral current scattering, in a normal medium composed of protons neutrons (nuclei) and electrons, V µ − V τ = 0. Furthermore, the difference of the potentials for ν e and ν a is due to the charged current scattering of ν e on electrons (ν e e → ν e e) [1]: The difference of potentials leads to the appearance of an additional phase difference in the neutrino system: φ matter ≡ (V e − V a )t ≈ V x. This determines the refraction length, the distance over which an additional "matter" phase equals 2π,: Numerically, l 0 = 1.6 · 10 9 cm 1 g/cm 3 n e m N , (9) where m N is the nucleon mass. The corresponding column density d ≡ l 0 n e = √ 2π/G F is given by the Fermi coupling constant only.
For antineutrinos the potential has an opposite sign. Being zero in the lowest order the difference of potentials in the ν µ − ν τ system appears at a level of 10 −5 V due to the radiative corrections [10]. Thus in the flavor basis in the lowest order in EW interactions the effect of medium on neutrinos is described byV = diag(V e , 0, 0) with V e given in Eq. (7).
The potential has been computed for neutrinos in different type of media, such as polarized or heavily degenerate electrons, in [11][12][13].
B. Evolution equation, effective Hamiltonian, and mixing in matter

Wolfenstein equation
In the flavor basis, the Hamiltonian in matter can be obtained by adding the interaction term to the vacuum Hamiltonian in vacuum [1, 3-5, 14, 15]: In Eq. (10) we have omitted irrelevant parts of the Hamiltonian proportional to the unit matrix. The Hamiltonian for antineutrinos can be obtained by the substitution There are different derivations of the neutrino evolution equation in matter, in particular, strict derivations starting from the Dirac equation or derivation in the context of quantum field theory (see [16] and references therein).
Although the Hamiltonian H f describes evolution in time, with the connection x = vt ≈ x = ct, Eq. (12) can be rewritten as idν f /dx = (H 0 +V )ν f with V = V (x), so it can be used as an evolution equation in space.
Due to the strong hierarchy of ∆m 2 and the smallness of 1-3 mixing, the results can be qualitatively understood and in many cases quantitatively described by reducing 3ν-evolution to 2ν-evolution. The reason is that the third neutrino effectively decouples and its effect can be considered as a perturbation. Of course, there are genuine 3ν− phenomena such as CP-violation, but even in this case the dynamics of evolution can be reduced effectively to the dynamics of evolution of 2ν−systems. The evolution equation for two flavor states, ν T f = (ν e , ν a ), in matter is − cos 2θ sin 2θ sin 2θ cos 2θ + where the Hamiltonian is written in symmetric form.

C. Mixing and eigenstates in matter
The mixing in matter is defined with respect to ν im -the eigenstates of the Hamiltonian in matter H f . As usually, the eigenstates are obtained from the equation where H im are the eigenvalues of H f . If the density, and therefore H f , are constant, ν im correspond to the eigenstates of propagation. Since H f = H 0 , the states ν im differ from the mass states, ν i . For low density, n → 0, the vacuum eigenstates are recovered: ν im → ν i . If the density, and thus H f , changes during neutrino propagation, ν im and H im should be considered as the eigenstates and eigenvalues of the instantaneous Hamiltonian: . The mixing in matter is a generalization of the mixing in vacuum (1). Recall that the mixing matrix in vacuum connects the flavor neutrinos, ν f , and the massive neutrinos, ν mass . The latter are the eigenstates of Hamiltonian in vacuum: ν H = ν mass . Therefore, the mixing matrix in matter is defined as the matrix which relates the flavor states with the eigenstates of the Hamiltonian in matter ν T H = (ν 1m , ν 2m , ν 3m ): From Eq. (13) we find that Furthermore, the Hamiltonian can be represented in the flavor basis as Inserting this expression as well as the relation ν jm = U m * αj ν α , which follows from Eq. (14) Equation (14) can be inverted to ν H = U m † ν f , or in components ν im = U m * αi ν α , α = e, µ, τ . According to this, the elements of mixing matrix determine the flavor content of the mass eigenstates so that |U m αi | 2 gives the probability to find ν α in a given eigenstate ν im . Correspondingly, the elements of the PMNS matrix determine the flavor composition of the mass eigenstates in vacuum.

D. Mixing in the two neutrino case
In the 2ν-case, there is single mixing angle in matter θ m and the relations between the eigenstates in matter and the flavor states reads ν e = cos θ m ν 1m + sin θ m ν 2m , ν a = cos θ m ν 2m − sin θ m ν 1m .

5
The angle θ m is obtained by diagonalization of the Hamiltonian (12) (see previous section): where R is the resonance factor. In the limit V → 0, the factor R → 1 and the vacuum mixing is recovered. The difference of eigenvalues H im equals This difference is also called the level splitting. or oscillation frequency, which determines the oscillation length: The matter potential and ∆m 2 always enter the mixing angle and other dimensionless quantities in the combination where l 0 is the refraction length. This is the origin of the "scaling" behavior of various characteristics of the flavor conversion probabilities. In terms of the mixing angle in matter the Hamiltonian can be rewritten in the following symmetric form

Resonance and level crossing
According to Eq. (19) the effective mixing parameter in matter, sin 2 2θ m , depends on the electron density and neutrino energy through the ratio (21) of the oscillation and refraction lengths, x = l ν /l 0 ∝ EV . The dependence sin 2 2θ m (V E) for two different values of the vacuum mixing angle, corresponding to angles from the full three flavor framework, is shown in Fig. 1. The dependence of sin 2 2θ m on E has a resonant character [3]. At l ν = l 0 cos 2θ (23) the mixing becomes maximal: sin 2 2θ m = 1 (R = sin 2 2θ). The equality in (23) is called the resonance condition and it can be rewritten as 2EV = ∆m 2 cos 2θ. For small vacuum mixing the condition reads: Oscillation length ≈ Refraction length. The physical meaning of the resonance is that the eigenfrequency, which characterizes a system of mixed neutrinos, ω = 2π/l ν = ∆m 2 /2E, coincides with the eigenfrequency of the medium, 2π/l 0 = 1/V . The resonance condition (23) determines the resonance density The width of resonance on the half of height (in the density scale) is given by 2∆n R e = 2n R e tan 2θ. Similarly, for fixed n e one can introduce the resonance energy and the width of resonance in the energy scale. The width can be rewritten as ∆n R e = n 0 sin 2θ, where n 0 ≡ ∆m 2 /2 √ 2EG F . When the vacuum mixing approaches maximal value, θ → π 4 the resonance shifts to zero density: n R e → 0, the width of resonance increases converging to fixed value: ∆n R e → n 0 . In a medium with varying density, the layer in which the density changes in the interval n R e ± ∆n R e is called the resonance layer. In this layer the angle θ m varies in the interval from π/8 to 3π/8.
For V V R , the mixing angle is close to the vacuum angle: θ m ≈ θ, while for V V R the angle becomes θ m ≈ π/2 and the mixing is strongly suppressed. In the resonance region, the level splitting is minimal [17,18], therefore the oscillation length, as the function of density, is maximal.

E. Mixing of 3 neutrinos in matter
To a large extent, knowledge of the eigenstates (mixing parameters) and eigenvalues of the instantaneous Hamiltonian in matter allows the determination of flavor evolution in most of the realistic situations (oscillations in matter of constant density, adiabatic conversion, strong breaking of adiabaticity). The exact expressions for the eigenstates FIG. 1: Resonance in neutrino mixing. The dependence of sin 2 2θmij on the product V E for vacuum mixing: sin 2 2θ12 = 0.851, ∆m 2 21 = 7.59 · 10 −5 eV 2 (red) and sin 2 θ13 = 0.0241, ∆m 2 31 = 2.47 · 10 −3 eV 2 (green). The left semi-plane corresponds to antineutrinos. The behavior of θ23 with vacuum value sin 2 2θ23 = 0.953 is included for completeness. The dashed lines are the predictions from a strict two-flavor approximation while the solid thin lines are the results of numerical diagonalization of the full three-flavor system. The upper panels show the case of the normal mass hierarchy and the lower panels -the inverted hierarchy.
and eigenvalues [19,20] are rather complicated and difficult to analyze. Therefore approximate expressions for the mixing angles and eigenvalues are usually used. They can be obtained performing an approximate diagonalization of H f which relies on the strong hierarchy of the mass squared differences: Without changing physics, the factor I −δ in the mixing matrix can be eliminated by permuting it with U 12 and redefining the state ν 3 . Therefore, in what follows, we use U P M N S = U 23 I δ U 13 U 12 . Here we will here describe the case of normal mass hierarchy: ∆m 2 31 > 0, ∆m 2 32 > 0. Subtracting from the Hamiltonian the matrix proportional to the unit matrix m 2 1 /2EI, we obtain

Propagation basis
The propagation basis,ν = (ν e ,ν 2 ,ν 3 ) T , which is most suitable for consideration of the neutrino oscillations in matter is defined through the relation Since the potential matrix is invariant under 2-3 rotations the matrix of the potentials is unchanged and the Hamiltonian the propagation basis becomesH It does not depend on the 2-3 mixing or CP-violation phase, and so the dynamics of the flavor evolution does not depend on δ and θ 23 . These parameters appear in the final amplitudes when projecting the flavor states onto propagation basis states and back (27) Here all the off-diagonal elements contain small parameters r ∆ and/or s 13 . Notice that, for the measured oscillation parameters, s 2 13 ∼ r ∆ .

Mixing angles in matter
The Hamiltonian in Eq. (29) can be diagonalized performing several consecutive rotations which correspond to developing the perturbation theory in r ∆ . After a 1-3 rotatioñ over the angle θ m 13 determined by tan 2θ m 13 = sin 2θ 13 cos 2θ 13 the 1-3 element of (29) vanishes. The expression (31) differs from that for 2ν mixing in matter by a factor (1 − s 2 12 r ∆ ), which increases the potential and deviates from 1 by ξ ≡ s 2 12 r ∆ ≈ 10 −2 .
After this rotation the Hamiltonian in the ν basis (30) becomes where and x ≡ 2EV /∆m 2 31 . For ξ = 0, these elements are reduced to the standard 2ν expressions. In the limit of zero density, x → 0, h 11 = ξ = s 2 12 r ∆ and consequently the 11 element of the Hamiltonian equals H 11 = s 2 12 ∆m 2 12 /2E.
In the lowest r ∆ approximation one can neglect the non-zero 2-3 element in Eq. (32). The state ν 3 then decouples and the problem is reduced to a two neutrino problem for (ν 1 , ν 2 ). The eigenvalue of this decoupled state equals The diagonalization of the remaining 1-2 sub-matrix is given by rotation where θ m 12 is determined by Here h 11 and θ m 13 are defined in Eqs. (33) and (31), respectively. The eigenvalues equal According to this diagonalization procedure in the lowest order in r ∆ the mixing matrix in matter is given by where mixing angles θ m 12 and θ m 13 are determined in Eqs. (36) and (31), respectively. The 2-3 angle and the CP-violation phase are not modified by matter in this approximation. The eigenvalues H 1m and H 2m are given in Eq. (37) and H 3m is determined by Eq. (34).
The 2-3 element of matrix (32) vanishes after additional 2-3 rotation by an angle θ 23 ∼ r ∆ : which produces corrections of the next order in r ∆ . With an additional 2-3 rotation the mixing matrix becomes where i.e., the combination sin δ sin 2θ 23 is invariant under inclusion of matter effects. Furthermore, θ m 23 ≈ θ 23 and δ m ≈ δ up to corrections of the order O(r ∆ ). The results described here allow to understand behavior of the mixing parameters sin 2 2θ mij in the EV region of the 1-3 resonance and above it (see Fig. 1).
In Fig. 2 we present dependence of the flavor content of the neutrino eigenstates on the potential. The energy level scheme, the dependence of the eigenvalues H im on matter density, is shown in Fig. 3. The energy levels in matter do not depend on δ or θ 23 , but they do depend on the 1-3 and 1-2 mixing.
In the case of normal mass hierarchy, there are two resonances (level crossings). whose location is defined as the density (energy) at which the mixing in a given channel becomes maximal.
1. The H-resonance, in the ν e − ν τ channel, is associated to the 1-3 mixing and large mass splitting. According to Eq. (31) θ m 13 = π/4 at 2. The L-resonance at low densities is associated to the small mass splitting and 1-2 mixing It appears in the ν e − ν µ channel, where ν e and ν e differ by small (at low densities) rotation given by an angle ∼ θ 13 (see eq. (31)). According to Eq. (36) the position of the L-resonance, θ m 12 = π/4 is given by c 2 12 r ∆ = h 11 , where h 11 is defined in Eq. (33). This leads to For antineutrinos (V E < 0 in Fig. 3), the oscillation parameters in matter can be obtained from the neutrino parameters taking V → −V and δ → −δ. The mixing pattern and level scheme for neutrinos and antineutrinos are different both due to the possible fundamental violation of CP-invariance and the sign of matter effect. Matter The flavor contents of the eigenstates of the Hamiltonian in matter as functions of EV . The vertical width of the band is taken to be 1, then the vertical sizes of the colored parts give |Uei| 2 (red), |Uµi| 2 (green), |Uτi| 2 (blue). The right and left panels correspond to neutrinos and anti-neutrinos, respectively. We take the best fit values of [21] with δ = 0. Variations of δ change the relative νµ− and ντ − content. The dashed red line shows a shift of border between νµ− and ντ − flavors for δ = π. The upper (lower) panel corresponds to normal (inverted) mass ordering.
violates CP-invariance and the origin of this violation stems from the fact that usual matter is CP-asymmetric: in particular, there are electrons in the medium but no positrons.
In the case of normal mass hierarchy there is no antineutrino resonances (level crossings), and with the increase of density (energy) the eigenvalues have the following asymptotic limits:

III. EFFECTS OF NEUTRINO PROPAGATION IN DIFFERENT MEDIA
A. The evolution matrix The evolution matrix, S(t, t 0 ), is defined as the matrix which gives the wave function of the neutrino system ν(t) at an arbitrary moment t once it is known in the initial moment t 0 : Inserting this expression in the evolution equation (12), we find that S(t, t 0 ) satisfies the same evolution equation as ν(t): The elements S(t, t 0 ) αβ of this matrix are the amplitudes of ν β → ν α transitions: S(t, t 0 ) αβ ≡ A(ν β → ν α ). The transition probability equals P αβ = |S(t, t 0 ) αβ | 2 . The unitarity of the evolution matrix, S † S = I, leads to the following relations between the amplitudes (matrix elements) The first and the second equations express the fact that the total probability of transition of ν α to everything is one, and the same holds for ν β . The third and fourth equations are satisfied if With these relations the evolution matrix can be parametrized as The Hamiltonian for a 2ν system is T-symmetric in vacuum as well as in medium with constant density. In medium with varying density the T-symmetry is realized if the potential is symmetric. Under T-transformations S βα → S αβ , and the diagonal elements S αα do not change. Therefore according to (48) the T-invariance implies that S βα = −S * βα , or Re S βα = 0, i.e., the off-diagonal elements of the S matrix are pure imaginary.

B. Neutrino oscillations in matter with constant density
In a medium with constant density and therefore constant potential the mixing is constant: θ m (E, n) = const. Consequently, the flavor composition of the eigenstates do not change and the eigenvalues H im of the full Hamiltonian are constant. The two neutrino evolution equation in matter of constant density can be written in the matter eigenstate basis as where H diag ≡ diag(H 1m , H 2m ). This system of equations splits and the integration is trivial, ν im (t) = e −iHimt ν im (0). The corresponding S-matrix is diagonal:S where φ m ≡ 1 2 ω m x is the half-oscillation phase in matter and a matrix proportional to the unit matrix has been subtracted from the Hamiltonian.
The S matrix in the flavor basis (ν e , ν a ) is therefore Then, for the transition probability, we can immediately deduce where φ m = πx/l m with being the oscillation length in matter. The dependence of l m on the neutrino energy is shown in Fig. 4. For small energies, V E ∆m 2 , the length l m l ν . It then increases with energy and for small θ reaches the maximum l max m = l 0 / sin 2θ at E max = E R / cos 2 2θ, i.e., above the resonance energy. For E → ∞ the oscillation length converges to the refraction length l m → l 0 . A useful representation of the S matrix for a layer with constant density follows from Eq. (52): where σ is a vector containing the Pauli matrices and n ≡ (sin 2θ m , 0, − cos 2θ m ). The dynamics of neutrino flavor evolution in uniform matter are the same as in vacuum, i.e., it has a character of oscillations. However, the oscillation parameters (length and depth) differ from those in vacuum. They are now determined by the mixing and effective energy splitting in matter: sin 2 2θ → sin 2 2θ m , l ν → l m .

C. Neutrino polarization vectors and graphic representation
It is illuminating to consider dynamics of transitions in different media using graphic representation [22][23][24]. Consider the two flavor neutrino state, ψ T = (ψ e , ψ a ). The corresponding Hamiltonian can be written as where σ = (σ 1 , σ 2 , σ 3 ), H is the Hamiltonian vector H ≡ (2π/l m ) · (sin 2θ m , 0, cos 2θ m ) and l m = 2π/∆H m is the oscillation length. The evolution equation then becomes Let us define the polarization vector P In terms of the wave functions, the components of P equal The z-component can be rewritten as P z = |ψ e | 2 − 1/2, therefore P e ≡ |ψ e | 2 = P z + 1/2 and from unitarity P a ≡ |ψ a | 2 = 1/2 − P z . Hence, P z determines the probabilities to find the neutrino of in a given flavor state. The flavor evolution of the neutrino state corresponds to a motion of the polarization vector in the flavor space. The evolution equation for P can be obtained by differentiating Eq. (58) with respect to time and insertingψ andψ † from evolution equation (57). As a result, one finds that If H is identified with the strength of a magnetic field, the equation of motion (60) coincides with the equation of motion for the spin of electron in the magnetic field. According to this equation P precesses around H.
With an increase of the oscillation phase φ (see Fig. 5) the vector P moves on the surface of the cone having axis H. The cone angle θ a , the angle between P and H depends both on the mixing angle and on the initial state, and in general, changes in process of evolution, e.g., if the neutrino evolves through several layers of different density. If the initial state is ν e , the angle equals θ a = 2θ m in the initial moment.
The components of the polarization vector P are nothing but the elements of the density matrix ρ = σ · P. The evolution equation for ρ can be obtained from (60) The diagonal elements of the density matrix give the probabilities to find the neutrino in the corresponding flavor state.

D. Resonance enhancement of oscillations
Suppose a source produces flux of neutrinos in the flavor state ν µ with continuous energy spectrum. This flux then traverse a layer of length L with constant density n e . At the end of this layer a detector measures the ν e component of the flux, so that oscillation effect is given by the transition probability P µe . In Fig. 6 we show dependence of this probability on energy for thin and thick layers. The oscillatory curves are inscribed in to the resonance envelope sin 2 2θ m . The period of the oscillatory curve decreases with the length L. At the resonance energy, oscillations proceed with maximal depths. Oscillations are enhanced up to P > 1/2 in the resonance range (E R ±∆E R ) where ∆E R = tan 2θE R (see Sec. II D 1). This effect was called the resonance enhancement of oscillations.

E. Three neutrino oscillations in matter with constant density
The oscillation probabilities in matter with constant density have the same form as oscillation probabilities in vacuum and the generalization of Eq. (51) is straightforward. In the basis of the eigenstates of the Hamiltonian the evolution matrix equalsS and for the elements of the S matrix in the flavor basis we obtain . Removing e −2iφ2m and using the unitarity of the mixing matrix in matter we have In particular, for the amplitudes in matter involving only ν e and ν µ , we obtain [[do we use this? add more?]] F. Propagation in a medium with varying density and the MSW effect

Equation for the instantaneous eigenvalues and the adiabaticity condition
In non-uniform media, the density changes along neutrino trajectory: n e = n e (t). Correspondingly, the Hamiltonian of system depends on time, H = H(t), and therefore the mixing angle changes during neutrino propagation: θ m = θ m (n e (t)). Furthermore, the eigenstates of the instantaneous Hamiltonian, ν 1m and ν 2m , are no longer the "eigenstates" of propagation. Indeed, inserting ν f = U (θ m )ν m in the equation for the flavor states [c.f., Eq. (3)] we obtain the evolution equation for eigenstates ν im whereθ m ≡ dθ m /dt. The Hamiltonian for ν im (68) is non-diagonal, and consequently, the transitions ν 1m ↔ ν 2m occur. The rate of these transitions is given by the speed with which the mixing angle changes with time. According to Eq. (68) [3,25], |θ m | determines the energy of transition ν 1m ↔ ν 2m and |H 2m − H 1m | gives the energy gap between the levels. The off-diagonal elements of the evolution equation Eq. (68) can be neglected ifθ m is much smaller than other energy scales in the system. The difference of the diagonal elements of the Hamiltonian is, in fact, the only other energy quantity and therefore the criterion for smallness ofθ m iṡ This inequality implies a slow enough change of density and is called the adiabaticity condition. Defining the adiabaticity parameter as [22,25] as the adiabaticity condition can be written as γ 1. For small mixing angle, the adiabaticity condition is most crucial in the resonance layer where the level splitting is small and the mixing angle changes rapidly. In the resonance point, it takes the physically transparent form [3]: where l R m ≡ l ν /sin 2θ is the oscillation length in resonance, and ∆r R ≡ ne dne/dr R tan 2θ is the spatial width of the resonance layer. According to this condition at least one oscillation length should be obtained within the resonance layer.
In the case of large vacuum mixing, the point of maximal adiabaticity violation [26,27] is shifted to density, n e (av), larger than the resonance density: n e (av) → n B > n R . Here n B = ∆m 2 /2 √ 2G F E is the density at the border of resonance layer for maximal mixing. Outside the resonance and in the non-resonant channel, the adiabaticity condition has been considered in [28,29].

G. Adiabatic conversion and the MSW effect
If the adiabaticity condition if fulfilled andθ m can be neglected, the Hamiltonian for the eigenstates becomes diagonal. Consequently, the equations for the instantaneous eigenstates ν im split as in the constant density case. The instantaneous eigenvalues evolve independently, but the flavor content of the eigenstates changes according to the change of mixing in matter. This is the essence of the adiabatic approximation: We neglectθ m in evolution equation but do not neglect the dependence of θ m on density. The solution can be obtained immediately as in symmetric form. The only difference from the constant density case is that the eigenvalues now depend on time and therefore integration appears in the phase factors. The evolution matrix in the flavor basis can be obtained by projecting back from the eigenstate basis to the flavor basis with the mixing matrices corresponding to initial and final densities: If the initial and final densities coincide, as in the case of neutrinos crossing the Earth, we obtain the same formulas as in constant density case: with the mixing angle taken at the borders (initial or final state). In particular, the survival probability equals Averaging over the phase, which means that the contributions from ν 1 and ν 2 add incoherently, gives The mixing in the neutrino production point θ 0 m is determined by density in this point, n 0 e , and the resonance density. Consequently, the picture of the conversion depends on how far from the resonance layer (in the density scale) a neutrino is produced. Strong transitions occur if the initial and final mixings differ substantially, which is realized when initial density is much above the resonance density and final one is below the resonance density and therefore neutrinos cross the resonance layer.
According to Eq. (73) the oscillation depth equals D = | sin 2θ m sin 2θ 0 m |. Both the averaged probability (75) and the depth (73) are determined by the initial and final densities and do not depend on the density distribution along the neutrino trajectory. Essentially they are determined by the ratios y ≡ n/n R in the initial and final moments. This is a manifestation of the universality of the adiabatic approximation result.
In contrast, the phase do depend on the density distribution and the period of oscillations (the latter is given is by the oscillation length in matter). So, it is the phase that encodes an information about the density distribution.
The probability depends on t via the phase φ m (t) and also via the mixing angle θ m (t). Two degrees of freedom are operative and P dependence on time is an interplay of two effects: oscillations, associated to the phase φ m (t), and the adiabatic conversion related to change of θ m . Depending on initial condition n 0 e , the relative importance of the two effects is different. If neutrinos are produced far above the resonance, n 0 e n R e , the initial mixing is strongly suppressed, θ 0 m ≈ π/2. Consequently, the neutrino state, e.g. ν e , consists mainly of one eigenstate, ν 2m , and furthermore, one flavor ν e , dominates in ν 2m . Since the admixture of the second eigenstate is very small, oscillations (interference effects) are strongly suppressed. Thus, here the non-oscillatory flavor transition occurs when the flavor of whole state (which nearly coincides with ν 2m ) follows the density change. At zero density ν 2m = ν 2 , and therefore the probability to find the electron neutrino (survival probability) equals [3] The final probability, P = sin 2 θ, is the feature of the non-oscillatory transition (as pure adiabatic conversion). Deviation from this value indicates the presence of oscillations, see Eq. (73). If neutrinos are produced not too far from resonance, e.g. at n 0 e > n R e , the initial mixing is not suppressed. Although ν 2m is the main component of the neutrino state, the second eigenstate, ν 1m , has appreciable admixture; the flavor mixing in the neutrino eigenstates is significant, and the interference effect is not suppressed. Here we deal with the interplay of the adiabatic conversion and oscillations.
Production in the resonance is a special case: If θ 0 m = 45 • , the averaged probability equalsP = 1/2 independently of the final mixing. This feature is important for determining the oscillation parameters. Strong transitions (P > 1/2) occur when neutrinos cross resonance layer. These features are realized for solar neutrinos when propagating from their production region inside the Sun to the surface of the Sun. The adiabatic propagation occurs also in a single layer of the Earth (e.g. in the mantle).

H. Adiabaticity violation
For most of applications the adiabaticity is either well satisfied (neutrinos in the Sun or supernovae), or maximally broken due to sharp (instantaneous) density change (neutrinos in the Earth, neutrinos crossing the shock wave fronts in supernova). In the former case the evolution is described by the adiabatic formulas. In the latter case description is also simple, one just needs to match the flavor conditions at the borders between layers: find the flavor state before the density jump and then use it as an initial state for the evolution after the jump. The intermediate case of the adiabaticity breaking can be realized for neutrinos in the mantle of the Earth, for high energy neutrinos propagating in the Sun (neutrinos from annihilation of hypothetical WIMPs) or for sterile neutrinos with very small mixing.
If the density changes rapidly,θ m is not negligible in (68) and the adiabaticity condition (70) is not satisfied. The transitions ν 1m ↔ ν 2m become noticeable and therefore the admixtures of the eigenstates in a given propagating state change. The S matrix in the flavor basis is given by whereS is the evolution matrix in the basis of instantaneous eigenstates. Then the ν e − ν e transition probability P ee ≡ |S f (x, 0) ee | 2 equals where P 21 ≡ |S 21 | 2 is the probability of ν 2m → ν 1m transitions and P int is an interference term which depends on the oscillation phase. The averaged probability (P int = 0) equals [30] If the initial density is much larger than the resonance density, then θ m (0) ≈ π/2 and cos 2θ m (0) = −1. In this case the averaged probability can be rewritten as Violation of adiabaticity weakens transitions if cos 2θ m (t) > 0, thus leading to an increase of the survival probability.
In the adiabatic case S 11 = e iφm , S 21 = 0, and therefore S 2 11 + S * 2 11 = 2 cos 2φ m (x), so that Eq. (77) is reduced to (73). In the graphic representation (Fig. 5), the neutrino vector moves on the surface of the cone (phase change) and the axis of the cone rotates according to the density change. The cone angle θ a changes as a result of violation of the adiabaticity).
There are different approaches to compute the flop probability P 21 . In the adiabatic regime the probability of transition between the eigenstates is exponentially suppressed P 12 ∼ exp (−π/2γ) with γ given in Eq. (70) [30,31]. One can consider such a transition as penetration through a barrier of height H 2m − H 1m by a system with the kinetic energy dθ m /dt. This leads to the Landau-Zener probability: where h ≡ n(dn/dr) −1 [32]. In the case of weak adiabaticity violation, one can develop an adiabatic perturbation theory which gives the results as a series expansion in the adiabaticity parameter [33].
I. Theory of small matter effects

Minimal width condition
If the vacuum mixing angle is small, there exists a lower limit on amount of matter needed to induce significant flavor change due to matter effect. The amount of matter is characterized by the column density of electrons along the neutrino trajectory: We can define d 1/2 as the column density for which the oscillation transition probability surpasses 1/2 for the first time in the course of propagation. Then it is possible to show that [34] for all density profiles. Furthermore, the minimum, d min , is realized for oscillations in a medium of constant density equal to the resonance density. The relation (83) is known as the minimal width condition. This condition originates from an interplay between matter effects and vacuum mixing: The acquired matter phase, √ 2G F d, must be large. At the same time, since matter effects by themselves are flavor conserving, also vacuum mixing is required in order to induce flavor conversion. The smaller vacuum mixing, the large width is required.

Vacuum mimicking
Vacuum mimicking [35], which states that regardless of the matter density, the initial flavor evolution of neutrino state is similar to that of vacuum oscillations. Consequently for small baselines, L, it is not possible to see matter effect and any such effect appears in higher order of L. Indeed, consider the evolution matrix where T denotes time ordering of the exponential. For small values of L, it can be expanded as If initial neutrino state has definite flavor, the amplitude of flavor transition is given by the off-diagonal element of H(x) which does not depend on matter potential. The matter contribution to H(x) is diagonal. Therefore the flavor transitions depends on the matter density only at higher order in L. This result holds true as long as L l m or when the phase of oscillation is small [36].
This can be seen explicitly in the case of medium with constant density where expanding the oscillatory factor for small oscillation phase we have the transition probability Note that vacuum mimicking only occurs if the initial neutrino state is a flavor eigenstate [36]. If the initial neutrino is in a flavor-mixed state, e.g. in a mass eigenstate, then matter will affect this state already at lowest order in L. This situation is realized in several settings involving astrophysical neutrinos propagating through the Earth, e.g., solar and supernova neutrinos, where the neutrinos arrive at the Earth as mass eigenstates. The mimicking is not valid if there are non-standard flavor changing interactions, so that matter effect appears in the off-diagonal elements of the Hamiltonian.

Effects of small layers of matter
If the minimal width condition is not satisfied, that is d = nx G −1 F , the matter effect on result of evolution is small. This inequality can be written as V x 1 which means that the oscillation phase is small. In this case the matter effect can be considered as small perturbation of the vacuum oscillation result even if the MSW resonance condition is satisfied.
The reasons for the smallness of the matter effect are different depending on the energy interval. Consider a layer of constant density with the length x. There are three possibilities (i) E E R , (E R is the resonance density) -nearly vacuum oscillations in low density medium take place. Matter effect gives small corrections to the oscillation depth and length which are characterized by 2V E (ii) E ∼ E R -modification of oscillation parameters is strong, however l R ν ∼ l ν / sin 2θ ∼ 2π/(V sin 2θ). Consequently, x/l R ν = xV sin 2θ/2π 1. Oscillations are undeveloped due to smallness of phase. (iii) E E R -matter suppresses oscillation depth by a factor E R /E 1. Since the oscillation length equals l m ≈ 2π/V , one obtains x/l m = xV /2π 1. Hence in this case the distance is very small and oscillation effect in the layer has double suppression. J. Propagation in multilayer medium

Parametric effects in the neutrino oscillations
The strong transitions discussed in the previous sections require the existence of large effective mixing, either in the entire medium (constant density) or at least in a layer (adiabatic conversion). There is a way to get strong transition without large vacuum or matter mixings. This can be realized with periodically or quasi periodically changing density [24,37] when the conditions of parametric resonance are satisfied. Although the flavor conversion in a layer which corresponds to one period is small, strong transitions can build up over several periods. For large mixing even a small number of periods is enough to obtain strong flavor transitions.
The usual condition of parametric resonance is that the period of density change T n is an integer times the effective oscillation length l m [38]: or l T /l m = k. Such an enhancement has been considered first for modulation of the profile by sine function [39]. This may have some applications for intense neutrino fluxes when neutrino-neutrino interactions become important. The solvable case, which has simple physical interpretation, is provided by the castle wall profile, for which the period l T is divided into two parts l 1 and l 2 (l 1 + l 2 = l T ) with densities n 1 and n 2 , respectively (n 1 = n 2 and, in general, l 1 = l 2 ). Thus, the medium consists of alternating layers with two different densities [37,[40][41][42][43][44][45][46].
For the "castle wall" profile, the simplest realization of the parametric resonance condition is reduced to equality of the oscillation phases acquired by neutrinos over the two parts of the periods [41]: The enhancement of transition depends on the number of periods and on the amplitude of perturbation, which determines the swing angle (the difference of the mixing angles in the two layers, ∆θ ≡ 2θ 1m − 2θ 2m ). For small ∆θ a large transition probability can be achieved after many periods. For large "swing" angle, even a small number of periods is sufficient.

Parametric enhancement, general consideration.
In general the condition (88) is not necessary for the enhancement or even for maximal enhancement. First, consider the oscillation effect over one period. The corresponding evolution matrix is given by the product where S k (k = 1,2) is the evolution in layer k given by Eq. (55). For brevity we will write it as S k = c k I − is k (σ · n k ), k = 1, 2, where c k ≡ cos φ k , s k ≡ sin φ k and φ k is the half-phase acquired in layer k: Here θ mk is the mixing angle in layer k. Insertion of S k from (55) into (89) gives [37] where Y ≡ c 1 c 2 − s 1 s 2 (n 1 · n 2 ), X = s 1 c 2 n 1 + s 2 c 1 n 2 − s 1 s 2 [n 1 × n 2 ].
Explicitly: (n 1 · n 2 ) = cos(2θ m1 − 2θ m2 ) and [n 1 × n 2 ] = sin(2θ m1 − 2θ m2 )e y . Using unitarity of S T , which gives X 2 + Y 2 = 1, one can parametrize X and Y with a new phase Φ as Y ≡ cos Φ and X ≡ sin Φ. Then the evolution matrix S T can be written in the form S T = cos Φ − i sin Φ(σ ·X) = e −i(σ·X)Φ , whereX ≡ X/X. Consequently, the evolution matrix after n periods equals It is simply accounted for by an increase of the phase: Φ → nΦ. This is the consequence of the fact that the evolution matrices over all periods are equal and therefore commute. If the evolution ends at some instant t which does not coincide with the end of a full period, i.e., t = nT + t , then S(t) = S(t )S n . The transition probability computed with Eq. (92) is It has the form of the usual oscillation probability with phase nΦ and depth (X 2 1 +X 2 2 )/X 2 . The oscillations described by Eq. (93) are called the parametric oscillations. Under condition − X 3 = s 1 c 2 cos 2θ m1 + s 2 c 1 cos 2θ m2 = 0, which is called the parametric resonance condition, the depth of oscillations (93) becomes 1 and the transition probability is maximal when nΦ = π/2 + πk, where k is an integer. There are different realizations of the condition (94) which imply certain correlations among the mixing angles and phases. The simplest one, c 1 = c 2 = 0, coincides with Eq. (88).

Parametric enhancement in three layers
For small number of layers an enhancement of flavor transition can occur due to certain relations between the phases and mixing angles in different layers. This in turn impose certain conditions on the parameters of the layers: their densities and widths. The conditions are the similar to the parametric resonance condition and this enhancement is called the parametric enhancement of flavor transitions. These conditions can be satisfied for certain energies and baselines for neutrinos propagating in the Earth.
Consider conditions for maximal enhancement of oscillations for different number of layers. It is possible to show [47] that they are generalizations of the conditions in one layer which require that (i) the depths of oscillations is 1 we call it the amplitude condition and the oscillation phase is φ = π/2 + πk -the phase condition.
Consider first the case of one layer with (in general) varying density (it can correspond to the mantle crossing trajectories in the Earth). The resonance condition for constant density case, cos 2θ m = 0, can be written according to Eqs. (22) and (49) 22 . Let us find the conditions for extrema for density profiles consisting of two layers. We have S (2) = S 2 S 1 , where S 12 = α 2 β 1 + β 2 α * 1 , and α i , β i for each layer have been defined in Eq. (49). The sum of the two complex numbers in the transition amplitude S 12 has the largest possible result if they have the same phase: arg(α 2 β 1 ) = arg(β 2 α * 1 ), which can also be rewritten as arg(α 1 α 2 β 1 ) = arg(β 2 ) .
This condition is called the collinearity condition [47]. It is an extremum condition for the two-layer transition probability under the constraint of fixed transition probabilities in the individual layers. In other words, if the absolute values |β i | of the transition amplitudes are fixed while their arguments are allowed to vary, then the transition probability reaches an extremum when these arguments satisfy Eq. (95). The conditions for maximal transition probability for three layers can be found in the following way. The 1-2 elements of the evolution matrix S (3) equals In the case of neutrino oscillations in the Earth, the third layer is just the second mantle layer, and its density profile is the reverse of that of the first layer. The evolution matrix for the third layer is therefore the transpose of that for the first one [48], i.e., α 3 = α 1 , β 3 = −β * 1 , and the expression for S 12 can be written as S Note that β 2 is pure imaginary because the core density profile is symmetric. Therefore the amplitude S 12 in Eq. (97) is also pure imaginary, as it must be because the overall density profile of the Earth is symmetric as well. If the collinearity condition for two layers (95) is satisfied, then not only the full amplitude S 12 , but also each of the four terms on the right hand side of Eq. (97) is pure imaginary. If the collinearity condition is satisfied for two layers, then it is automatically satisfied for three layers. This is a consequence of the facts that the density profile of the third layer is the reverse of that of the first layer and that the second layer has a symmetric profile. The conditions described here allow to reproduce very precisely all main structures of the oscillograms of the Earth (see sect. IV A).

K. Oscillations of high energy neutrinos
At high energies or in high density medium when V > ∆m 2 /2E, we can use ∆/V ≡ ∆m 2 /4EV as a small parameter and develop a perturbation theory using its smallness. However, in most situations of interest, the neutrino path length in matter L is so large that ∆ · L > ∼ 1. Therefore the vacuum part of the Hamiltonian cannot be considered as a small perturbation in itself and the effect of ∆ on the neutrino energy level splitting should be taken into account. For this reason we split the Hamiltonian as H =H 0 + H I with where ω m is the oscillation frequency (20) and ≡ (2∆ cos 2θ − V + ω m )/2∆ sin 2θ ≈ ∆ V sin 2θ 1. The ratio of the second and the first terms in the Hamiltonian (98) is given by the mixing angle in matter θ m : 2∆ sin 2θ/ω m = sin 2θ m . Therefore for sin 2θ m 1 the term H I can be considered as a perturbation. Furthermore, ∼ sin 2θ m , so the diagonal terms in H I can be neglected in the lowest approximation.
The solution for S matrix can be found in the form S = S 0 · S I , where S 0 is the solution of the evolution equation with H replaced by H 0 [see Eq. (71)]. The matrix S I then satisfies the equation whereH I ≡ S −1 0 H I S 0 is the perturbation Hamiltonian in the "interaction" representation. Eq. (99) can be solved by iterations: S I = I + S (1) I + ..., which leads to the standard perturbation series for the S matrix. For neutrino propagation between x = 0 and x = L we have, to the lowest non-trivial order, The ν e ↔ ν a transition probability P 2 = [S(L)] ae is given by For density profiles that are symmetric with respect to the center of the neutrino trajectory, V (x) = V (L − x), Eq. (101) gives where z = x − L/2 is the distance from the midpoint of the trajectory and φ(z) is the phase acquired between this midpoint and the point z. The transition probability P 2 decreases with the increase of neutrino energy essentially as E −2 . The accuracy of Eq. (101) also improves with energy as E −2 . Inside the Earth, the accuracy of the analytic formula is extremely good already for E > ∼ 8 GeV. When neutrinos do not cross the Earth's core (cos Θ > −0.837), and so experience a slowly changing potential V (x), the accuracy of the approximation (101) is very good even in the MSW resonance region E ∼(5 -8) GeV.
The above formalism applies in the low energy case as well, with only minor modifications: the sign of H 0 in Eq. (98) has to be flipped, and correspondingly one has to replace ω m → −ω m in the definition of . The expressions for the transition probability in Eqs. (101) and (102) remain unchanged.

L. Effects of small density perturbations
Let us consider perturbation around smooth profile for which exact solution is known. The simplest possibility that has implications for the Earth matter profile is the constant density with additional perturbation: V (x) =V +∆V (x). Correspondingly, the Hamiltonian of the system can be written as the sum of two terms: whereH ≡ω − cos 2θ sin 2θ sin 2θ cos 2θ , ∆H ≡ ∆V (x) 2 Here,θ = θ m (V ) is the mixing angle in matter andω = ω m (V ) is half of the energy splitting (half-frequency) in matter, both with the average potentialV . We will denote byS(x) the evolution matrix of the system for the constant density case H(x) =H. The expression forS(x) is given in Eq. (52) with θ m =θ and φ m ( The solution of the evolution equation with Hamiltonian (103) [47] is of the form where K 1 (x) satisfies |K 1 (x) ab | 1. Inserting Eq. (105) into the evolution equation, one finds the following equation for K 1 (x) to the first order in ∆H(x) and K 1 (x): − cos 2θ sin 2θ sin 2θ cos 2θ + sin 2θ cos 2φ G(θ) + sin 2θ sin 2φ σ 2 .
For practical purposes it is useful to have an expression for S which is exactly unitary regardless of the size of the perturbation. For this we rewrite Eq. (108) as follows: where sin ξ = ∆J √ (∆J) 2 +(∆I) 2 and = sin 2θ · (∆J) 2 + (∆I) 2 . Thus, S =S + ε S and we replace it by S = cos εS + sin ε S .
Here both S andS are unitary matrices, and due to their specific form the combination on the right-hand-side of Eq. (110) is exactly unitary. For a symmetric density profile with respect to the midpoint of the trajectory, the term ∆J is absent. From Eqs. (52), (108) and (110) we immediately get the transition probability P = cos ε sin 2θ sin φ + sin ε cos 2θ 2 ≈ sin 2 2θ sin φ + ∆I cos 2θ 2 , where ε ≡ sin 2θ ∆I and φ ≡ φ(L) =ωL. Here the first term in the square brackets describes oscillations in constant density matter with average potentialV 1 .

M. Oscillation probabilities and their properties
It is convenient to consider the neutrino flavor evolution in the propagation basisν = (ν e ,ν 2 ,ν 3 ) T , defined in Eq. (27). In this basis propagation is not affected by the 2-3 mixing and CP-violation. The dependence on these parameters appears when one projects the initial flavor state on the propagation basis and the final state back onto the original flavor basis. The propagation basis states are related to the mass states as Since the transformations which connectν and ν f , do not depend on matter potential and therefore distance, the statesν satisfy the the evolution equation i dν dt =Hν, with the HamiltonianH defined in Eq. (28).

S-matrix and oscillation amplitudes
A number of properties of the oscillation probabilities can be obtained from general consideration of matrix of the oscillation amplitudes. We introduce the evolution matrix (the matrix of amplitudes) in the propagation basis as Then according to Eq. (27) the S matrix in the flavor basis equals In this part, we use the notation A ij for the amplitudes in the propagation basis and S ij for the amplitudes in the flavor basis. In terms of the propagation-basis amplitudes (113) the S matrix in the flavor basis can be written as where K µµ ≡ s 23 c 23 (e −iδ A23 + e iδ A32) K µτ ≡ c 2 23 e −iδ A23 − s 2 23 e iδ A32 K τ µ = K µτ (δ → −δ,2 ↔3) The scheme of transitions is shown in Fig. 7. There is certain hierarchy of the amplitudes which can be obtained immediately from the form of the Hamiltonian in the propagation basis (29): i.e., A23 and A32 are the smallest amplitudes. In the propagation basis there is no fundamental CP-or T-violation. Therefore for a symmetric density profile with respect to the middle point of trajectory (as in the case of the Earth) the neutrino evolution is T-invariant which yields Consequently, for K αβ we obtain These terms proportional to small amplitudes A23 and A32 are of the order O(s 2 13 ). For a symmetric density profile, from Eqs. (115), (118) and (119) one finds for the probabilities P αβ ≡ |S βα | 2 : P µµ = |c 2 23 A22 + s 2 23 A33 + 2 s 23 c 23 cos δA23| 2 , P µτ = |s 23 c 23 (A33 − A22) + (cos 2θ 23 cos δ + i sin δ)A23| 2 . For antineutrinos the amplitudes can be obtained from the results presented above substituting Notice that the amplitudes of transitions (121) and (122), that involve ν e , are given by linear combinations of two propagation-basis amplitudes. The other flavor amplitudes depend on three propagation-basis amplitudes.

Factorization approximation and amplitudes for constant density
As follows immediately from the form of the HamiltonianH in Eq. (29), in the limits ∆m 2 21 → 0 or/and s 12 → 0 the stateν 2 decouples from the rest of the system, and consequently, the amplitude A e2 vanishes. In this limit, A e3 (as well as A33 and S ee ) is reduced to a 2ν amplitude which depends on the parameters ∆m 2 31 and θ 13 : A A (∆m 2 31 , θ 13 ) ≡ A e3 (∆m 2 21 = 0). The corresponding probability equals P A ≡ |A A | 2 . In the limit s 13 → 0 the stateν 3 decouples while the amplitude A e3 vanishes and the amplitude A e2 reduces to a 2ν amplitude depending on the parameters of the 1-2 sector, ∆m 2 21 and θ 12 . Denoting this amplitude by A S we have A S (∆m 2 21 , θ 12 ) ≡ A e2 (θ 13 = 0). We will use the notation P S ≡ |A S | 2 . This consideration implies that to the leading non-trivial order in the small parameters s 13 and r ∆ the amplitudes A e2 and A2 e are reduced to two neutrino probabilities and depend only on the "solar" parameters, whereas the amplitudes A e3 and A3 e -only on the "atmospheric" parameters: The approximate equalities in Eq. (126) are called the factorization approximation. Due to the level crossing phenomenon the factorization approximation (126) is not valid in the energy range of the 1-3 resonance where the 1-3 mixing in matter is enhanced. In the case of a matter with an arbitrary density profile, one can show, using simple power counting arguments, that the corrections to the factorization approximation for the amplitude A e2 are of order s 2 13 , whereas the corrections to the "atmospheric" amplitude A e3 are of order r ∆ [50], in agreement with our consideration for constant density. The amplitude A e3 does not in general have a 2-flavor form, once the corrections to the factorization approximation are taken into account.
Using the expressions for U m ei and U m µi in terms of the mixing angles in the standard parametrization, we can rewrite Eq. (65) as Here φ m 31 = φ m 32 + φ m 21 . Since to a good approximation θ m 23 ≈ θ 23 and δ m ≈ δ (see Sec. II E) [20,51] where A cst 23 ≡ i e iφ m 21 sin θ m 13 sin 2θ m 12 sin φ m 21 .
Notice that A cst 22 has exactly the form of the corresponding 2ν amplitude driven by the solar parameters. The amplitude A cst 33 also coincides to a very good approximation with the corresponding 2ν amplitude driven by the atmospheric parameters. In the approximation θ m 23 ≈ θ 23 and δ m ≈ δ the amplitudes (130), (131) and (132) can be identified with the corresponding amplitudes in the propagation basis.
3. Properties of the flavor oscillation probabilities 1). ν e − ν e channel. The total probability of the ν e disappearance equals 1 − P ee = P eµ + P eτ = P e2 + P e3 . (133) The probability P ee does not depend on the CP-violating phase and the 2-3 mixing in the standard parametrization. The interference of the solar and atmospheric modes in P ee originates mainly from P e3 ≡ |A e3 | 2 . The survival probability then equals P ee = 1 − P eµ − P eτ = 1 − P A − P S . At high energies, where the effects of the 1-2 mixing and mass splitting in P are suppressed, the probability is P ee ≈ 1 − P eτ ≈ 1 − P A .
Since the amplitude A e2 is suppressed at high energies due to the smallness of the 1-2 mixing in matter, in the lowest approximation we have The maximal value of the probability equals P µe s 2 23 . According to Eqs. (121) and (122) the oscillation probabilities P τ e and P eτ can be obtained from the corresponding probabilities P µe and P eµ through the substitution s 23 → c 23 , c 23 → −s 23 [52]. The interference term has the opposite signs for channels including ν τ as compared with those with ν µ , which can be obtained from the unitarity condition P ee + P µe + P τ e = 1 and the fact that P ee does not depend on δ.
4). ν µ − ν τ channel. For symmetric matter density profiles the probability of ν µ → ν τ oscillations is given in Eq. (124). It can be rewritten as The amplitude depends on δ through the terms proportional to cos δ and sin δ, and therefore P µτ contains both CPand T-even and odd terms. One can show that the δ-dependent interference terms, which are proportional to sin δ and cos δ, satisfy the relation P δ µτ = −P δ µe − P δ µµ . In the limit ∆m 2 21 → 0 we obtain

A. Propagation of neutrinos through the Earth
Flavor neutrino evolution in the Earth is essentially oscillations in a multi-layer medium with slowly changing density in the individual layers and sharp density change on the borders of layers. For energies E > 0.1 GeV, possible short-scale inhomogeneities of the matter distribution can be neglected and the density profile experienced by neutrinos is symmetric with respect to the midpoint of the trajectory: Here L = 2R ⊕ | cos θ z | is the length of the trajectory inside the Earth, R ⊕ = 6371 km is the Earth radius and θ z is the zenith angle related to the nadir angle as Θ ν = π − θ z . For 0 ≤ Θ ν ≤ 33.1 • neutrinos cross both the mantle and the core of the Earth, whereas for larger values of the nadir angle they only cross the Earth's mantle. The column density of the Earth along the diameter equals d Earth = n(x)dx, which is bigger than the minimal width; the size of the Earth is comparable with the neutrino refraction length. For the 1-2 channel, the adiabaticity is well satisfied for all energies. We can therefore use the adiabatic approximation. The results of the evolution are determined by the mixing at the surface of the Earth and by the adiabatic phase. In the 1-3 channel the adiabaticity is broken at the resonance. Thus, the constant density approximation with the average density works well in this regime. For energies below the resonance the matter effect becomes small and the constant density approximation and the adiabatic approximation give very similar results.
For the core crossing trajectories, the profile consists of three layers in the first approximation: (i) mantle (with increasing density); (ii) core (with a symmetric profile) and (iii) second mantle layer (with decreasing density). This second mantle layer is T-inverted with respect to the first. In this approximation the profile can be considered as three layers of constant effective densities. As such, it looks like a part (1.5 period) of the castle wall profile. Consequently, the parametric enhancement of oscillations, and in particular, the parametric resonance can be realized.

Neutrino oscillograms of the Earth
A comprehensive description of effects of neutrino passage through the Earth can be obtained in terms of neutrino oscillograms. The oscillograms are defined as lines of equal probabilities (or certain combinations of probabilities) in the E ν − cos θ z plane. In Fig. 8, we show the oscillograms for the oscillation probabilities P eµ and P µµ , as well as the corresponding probabilities for antineutrinos [44,47,[53][54][55][56].
The structure of the oscillograms is well defined and unique, and reflects the structure of the Earth as well as the properties of the neutrinos themselves. In a sense, the oscillograms are the neutrino images of the Earth. In contrast to usual light, there are several different images in different flavors as well as in neutrinos and antineutrinos.
The positions of all main structures of the oscillograms are determined by different realizations of the amplitude condition and the phase condition. These are generalizations of the condition for maximal flavor transitions in the case of vacuum oscillations or oscillation in uniform matter. Recall that, in the latter case, P = 1 requires • sin 2 2θ m = 1; the amplitude condition, which is nothing but the MSW resonance condition, and • φ = π/2 + πk; the phase condition.
At E > 1 GeV the main structures of oscillograms are generated by the 1-3 mixing. They include: 1. The MSW resonance pattern (resonance enhancement of the oscillations) for trajectories crossing only the mantle, with the main peak at E ν ∼ (5 − 7) GeV. The position of the maximum is given by the MSW resonance condition: whereV 1 (Θ ν ) is the average value of the potential along the trajectory characterized by Θ ν . The phase condition becomes 2φ(E ν , Θ ν ) = 2ω(V , E ν )L(Θ ν ) = π and the intersection of the resonance and the phase condition lines gives the absolute maximum of P A . Combining these conditions gives the coordinates of the peak: cos Θ ν = 0.77 and E R = 6 GeV.
2. Three parametric ridges in the domain of core-crossing trajectories | cos θ z | > 0.87 and E ν > 3 GeV. The parametric ridges differ by the oscillation phase acquired in the core, φ 2 : -Ridge A lies between the core resonance (at Θ ν ∼ 0 • ) and the mantle resonance regions, E ν ≈ 3 − 6 GeV. The phase in the core is φ 2 < ∼ π. This ridge merges with the MSW resonance peak in the mantle. -Ridge B is situated at E ν ≥ 5 GeV. For the lowest energies in the ridge and Θ ν ∼ 0, the half-phase in the core equals φ 2 ∼ (1.2 − 1.3)π. -Ridge C is located at E ν > 11 GeV in the matter dominated region, where the mixing, and consequently, oscillation depth are suppressed.
4. The regular oscillatory pattern at low energies with "valleys" of zero probability and ridges in the mantle domain and a more complicated pattern with local maxima and saddle points in the core domain.
In Fig. 9, we show graphic representations of oscillations which correspond to salient features of the oscillograms. For energies E ν < 1 GeV the main structures are induced by the 1-2 mixing with small corrections due to 1-3 vacuum oscillations. Neglecting effect of θ 13 we have 1 − P ee = |A e2 | 2 ≡ P S . The probabilities of the modes including ν e are expressed in terms of a unique probability P S . The 1-2 pattern differs from the pattern for the 1-3 mixing due to the large value of the 1-2 mixing. The oscillation length at the resonance is smaller than that for the 1-3 mixing, l R m = l ν / sin 2θ 12 ∼ l ν . The resonance energy is shifted to smaller values both due to ∆m 2 21 ∆m 2 31 and because of the factor cos 2θ 12 ≈ 0.4: E R 12 = ∆m 2

21
2V cos 2θ 12 . HereV is the average value of the potential. The adiabaticity is better satisfied than for the 1-3 mixing case and therefore the oscillation probability in the mantle is determined by the potential near the surface of the EarthV averaged over a distance of the order of the first oscillation length. The oscillation length in matter l m monotonically increases with energy, approaching the refraction length l 0 ≡ 2π/V (c.f., Fig. 4). The jump of the mixing angle at the mantle-core boundary is small. Thus, the sudden distortion of the oscillation patterns at Θ ν = 33 • is not as significant as it is for the 1-3 mixing, in particular below the 1-2 resonance energy. These features allow to understand the structure of the oscillograms. In the mantle domain (Θ ν > 33 • ) the oscillation pattern for neutrinos is determined by the resonance enhancement of oscillations. There are three MSW resonance peaks above 0.1 GeV, which differ from each other by value of the total oscillation phase. The outer peak (Θ ν ≈ 82 • ) corresponds to φ ≈ π/2, the middle (Θ ν = 60 • ) to φ ≈ 3π/2, and the inner (Θ ν ≈ 40 • ) to φ = 5π/2. Recall that such a large phase can be acquired due to the smaller resonance oscillation length in comparison to that of the 1-3 mixing case, for which only one peak with φ = π/2 can be realized. The resonance energy is given by Eq. (43), and for the surface potential we find The ratio of the 1-2 and 1-3 resonance energies equals E R 12 /E R 13 ≈ 1 50 . The estimate (143) is valid for the two outer peaks. For the peak at Θ ν = 40 • ,V is larger and, accordingly, the resonance energy is slightly smaller. The width of the 1-2 resonance is large and therefore the regions of sizable oscillation probability are more extended in the E ν direction as compared to the oscillations governed by the 1-3 mixing and splitting.
The resonance energy in the core is E R 12 ≈ 0.04 GeV. Therefore for E ν > (0.10 − 0.15) GeV the 1-2 mixing in the core is substantially suppressed by matter. The peak at E ν 0.2 GeV and Θ ν 25 • is due to the parametric enhancement of the oscillations. It corresponds to the realization of the parametric resonance condition when the oscillation half-phases equal approximately φ mantle ≈ π/2 and φ core ≈ 3π/2 (note that the total phase is ≈ 5π/2 and this parametric ridge is attached to the 5π/2 MSW peak in the mantle domain).

2.
Oscillograms for the inverted mass hierarchy The main change compared to the normal hierarchy is that the 1-3 resonance structure now appears in the antineutrino channel. The level crossing scheme is modified in comparison to NH. In the neutrino channel there is only the 1-2 resonance.
In the approximation of ∆m 2 21 = 0, the neutrino oscillograms for the inverted hierarchy coincide with the antineutrino oscillograms for the normal hierarchy and vice-versa, provided that ∆m 2 31 is taken to be the same in both cases [57]: The inclusion of the 1-2 mixing and mass splitting breaks this symmetry.
B. CP-violation effects

Interference and CP-violation
The survival probability P ee does not depend on the CP-violating phase δ neither for oscillations in vacuum nor in matter [58,59]. This is the consequence of the facts that δ can be removed by transforming to the propagation basis and that ν e is not affected by this transformation. For oscillations in vacuum, or in matter with symmetric density profiles, the other two survival probabilities, P µµ and P τ τ , are T-even quantities dependent on δ only through terms proportional to cos δ and cos 2δ [60]. In contrast to this, for oscillations in a matter with non symmetric density profiles, these probabilities also acquire terms proportional to sin δ and sin 2δ.
Introducing the phase φ ≡ arg(A * e2 A e3 ) and omitting small terms proportional to |A23| 2 = O(s 6 13 ) we obtain where D 23 ≡ 1 2 sin 4θ 23 cos δ Re A * 23 (A33 − A22) is proportional to the small deviation of the 2-3 mixing from maximal one. Notice that D 23 enters P δ µµ and P δ µτ with opposite signs while P δ µe does not depend on D 23 at all. D 23 is CP-even. The sum of these interference terms is zero.
For the other channels, P δ αβ = P −δ βα . For antineutrinos, according to (125), the probabilities have the same form as the corresponding probabilities derived above with a changed sign of δ and the amplitudes computed with the opposite sign of the potential. Thus, the δ dependent parts in all the channels are expressed in terms of two combinations of the propagation basis amplitudes, |A e2 A e3 | and D 23 .

Magic lines and CP-domains
To better assess the effect of δ, one can consider the difference of the oscillation probabilities for two different values of the CP-phase ∆P CP αβ (δ) ≡ P αβ (δ)−P αβ (δ 0 ). In practice, this quantifies how well the phase δ fits with some assumed true value δ 0 . The structure of the oscillograms for ∆P CP αβ (δ) can be understood in terms of the grids of magic lines and interference phase lines along which ∆P CP αβ (δ) ≈ 0. For the ν µ → ν e oscillation probability, the equality is exact and the condition ∆P CP µe = 0 is equivalent to This equality is satisfied if at least one of the following three conditions is fulfilled The last condition implies φ(E ν , Θ ν ) = (δ + δ 0 )/2 + πl . Under the conditions (149), the equality (148) is satisfied identically for all values of δ. In these cases the transition probability does not depend on the CP-phase. Since the amplitudes A e2 and A e3 are complex quantities, these conditions can be satisfied in isolated points of the (Θ ν , E ν ) plane only. In contrast to this, in the factorization approximation A e2 = A S and A e3 = A A both the conditions are fulfilled along certain lines in the oscillograms. This occurs because the amplitudes A S and A A take a 2-flavor form. In the basis of neutrino states where the corresponding 2 × 2 Hamiltonians are traceless, both A A and A S are pure imaginary because of the symmetry of the Earth's density profile [48].
Let us consider the equalities A S = 0 and A A = 0 using the constant density approximation: 1. The condition A S (E ν , Θ ν ) = 0 is satisfied when sin φ S (E ν , Θ ν ) = 0, which leads to At energies E ν > ∼ 0.5 GeV which are much higher than the 1-2 mixing MSW resonance in the mantle and in the core of the Earth one has ω m 21 ≈ V and the condition (151) becomes This expression is energy independent and determines the baselines for which the "solar" contribution to the probability vanishes [62]. In the plane (Θ ν , E ν ) it represents nearly vertical lines at fixed Θ ν . There are three solar magic lines which correspond to n = 1 (in the mantle domain) Θ ν ≈ 54 • and n = 2, 3 (in the core domain) [62] Θ ν ≈ 30 • and 12 • . The existence of a baseline (L ≈ 7600 km) for which the probability of ν e ↔ ν µ oscillations in the Earth is approximately independent of the "solar" parameters (∆m 2 21 , θ 12 ) and of the CP-phase δ was first pointed out in [63] and later discussed in e.g., [62,[64][65][66][67][68][69]. This baseline was dubbed "magic" in [64].
2. The atmospheric magic lines are determined by the condition A A (E ν , Θ ν ) = 0 [62]. Along these lines, the "atmospheric" contribution to the amplitudes of ν µ ↔ ν e and ν τ ↔ ν e transitions vanishes and the probabilities of oscillations involving ν e orν e do not depend on CP-phase. In the constant density approximation, the condition A A = 0 is satisfied when sin φ A = 0 (φ A = πk, k = 1, 2, . . . ) or explicitly For energies which are not too close to the 1-3 MSW resonance, it reduces to which corresponds to the bent curves in the (Θ ν , E ν ) plane. For very large energies, where ∆m 2 31 /2E V , the atmospheric lines approach the same vertical lines as the solar magic lines (152).
3. The condition (150) determines the interference phase lines in the (Θ ν , E ν ) plane. In the constant density approximation φ ≈ −φ m 31 . Consequently in the energy range between the two resonances we have i.e., in the first approximation φ does not depend on the matter density. From Eq. (150) we then obtain Thus, in the factorization approximation, the conditions (149) and (150) define three sets of lines (grid of magic lines) in the oscillograms (see Fig. IV B 2), which play crucial roles in understanding the CP violation effects. Along the lines, the probabilities P µe , P eµ P τ e and P eτ do not depend on the CP-phase in the first order approximation. The other probabilities depend on the phase weakly.
From Fig. IV B 2, we can see that the magic lines described above do not coincide exactly with the lines of ∆P CP µe = 0 which bound the CP-domains. Furthermore, interconnections of the latter occur. This is due to the breakdown of the factorization approximation.

C. Determination of hierarchy with accelerator experiments
An accelerator neutrino experiment has a fixed baseline which corresponds to a vertical line with the length determined by the available energy spectrum. In the oscillogram of Fig. 11 we have included such lines for a handful of accelerator experiments. Furthermore, this energy spectrum is usually peaked at certain energy (or narrow energy range) resulting in the experiment being most sensitive to the oscillation probability at that specific energy. An accelerator neutrino experiment would typically run for several years in neutrinos or antineutrinos before switching polarity and therefore getting information both on P αβ andP αβ . The goal of such a search is to observe in which channel the oscillation probability is suppressed and in which it is enhanced. If a neutrino experiment could run at energy similar to the resonant one and at a baseline of several thousand kilometers, then this determination would be quite simple. However, as can be seen from the oscillogram, accelerator neutrino experiments are confined to relatively shallow trajectories with rather poor oscillatory pattern, and this severely limits their capabilities leading to various degeneracies. In particular, lack of knowledge of the mass hierarchy is part of the famous eightfold degeneracy, which arises as follows. Assume we have access to the values of oscillation probabilities P µe andP µe at a given baseline L and energy E only. Then there exists three types of ambiguities that give rise to the same values of the probabilities in different parts of the parameter space (mixing angles, CP phase, signs of mass differences).
1. Sign (hierarchy) degeneracy: This is the degeneracy due to the unknown neutrino mass hierarchy. Changing the mass hierarchy, it is often possible to find a point in parameter space that predicts the same oscillation probabilities.
Since each of these degeneracies is twofold, an overall degeneracy is eightfold: 2 3 = 8. The first two of these degeneracies can be illustrated in a bi-probability plot of Fig. 12. As follows from this figure, even if both the probabilities (for a given neutrino energy) are known with infinite accuracy, we can not identify the hierarchy within the pink region.
For known mass hierarchy (e.g. normal one) a given value of θ 13 fixes ellipse in the plot along which the CP phase varies. Increasing θ 13 moves the ellipse up and to the right in the plot. Therefore for every point on an ellipse, there will be another ellipse corresponding another value θ prime 13 , which crosses this point and therefore θ prime 13 reproduces the same oscillation probabilities. For example, in the left intersection of the black and white ellipse (Fig.. 12) both combinations of θ 13 and δ correspond to those precise oscillation probabilities and there are also values of θ 13 and δ that will reproduce them in the inverted hierarchy. For the right intersection, the intrinsic degeneracy is still present, while the sign degeneracy is resolved. It should be remembered that this type of figure is just an illustration. In real experiment the neutrino energy spans over wide range, the oscillation probabilities would not be exactly known and strictly this type of consideration becomes invalid.
In order to see how these degeneracies manifest themselves in an experimental setup, we show the oscillation probability P µe as a function of the baseline length in Fig. 13. While the 295 km baseline is too short for matter effects to be very significant, as the baseline increases matter effects start being more and more important. In particular, when the oscillation phase maximum occurs at an energy similar to that of the matter resonance, as is the case of 7500 km baseline, we can see the enhancement of the transition probability in the neutrino channel for the normal hierarchy and the suppression in the inverted. In a simple two-flavor scenario, the amplitude of P µe at the resonance is one by definition in the normal mass hierarchy case. At the same time, the oscillation amplitude in the inverted hierarchy at the same energy is given by sin 2 2θ = sin 2 2θ 1 + 3 cos 2 2θ where the last equality holds for small θ. On the other hand, if the neutrino energy is far below the resonance in order to accumulate a significant oscillation phase, such as in the left and middle panels, then the oscillation amplitude will  be effectively given by The reason that the 810 km baseline separates the hierarchies better than the 295 km one is based mainly on the fact that the oscillation maximum can be reached for higher energies due to the longer baseline, and thus, the relative difference between probabilities for the two hierarchies increases. Also note that the oscillation probabilities for the 7500 km baseline is not very dependent on the CP-violating phase δ. This is due to the so-called magic baseline effect, which has been discussed before. In order to successfully determine the neutrino mass hierarchy in a single accelerator experiment, two conditions are of major importance: 1) The baseline must be long enough to allow for a significant value of V E in order to separate the neutrino and anti-neutrino oscillation probabilities. To separate the mass hierarchy determination from the effects of the CP-phase, this separation must be large enough to avoid overlap of the probabilities within the experimental uncertainties. 2) The statistics must be high enough and the systematics low enough in order to make the split statistically significant. The literature contains several proposals for long baseline experiments with baselines of several thousands of kilometers in order to satisfy these conditions. However, as we will discuss later, the large value of θ 13 also provides us with an opportunity to pin down the value of δ. Such measurements require the presence of interference terms which will be small at the very long baselines, and instead medium long baselines around 1000 km, such as the 810 km baseline shown in Fig. 13, may be preferable due to the significant δ dependence of probabilities.

CP-violation effects and the mass hierarchy
The figure Fig. 13 shows a significant dependence of the probabilities on the CP-violating phase δ, especially at small baselines. We are mainly interested in the oscillation probability at the first or second oscillation maximum, where an experiment would typically be placed. In these baselines L the ν µ − ν e oscillation probability (the "golden channel") can be expanded in the small quantity ∆m 2 21 L/2E which gives [70] P eµ ≈ s 2 23 P 2f + c 13 sin 2θ 13 sin 2θ 12 sin 2θ 23 where P 2f is the two-flavor oscillation probability discussed earlier. In Eq. 159 we have neglected terms of the second (and higher) order in ∆m 2 21 L/2E 2 . as well as the matter effect on ∆m 2 31 . It is the the second term that is responsible for creating the band of different oscillation probabilities displayed in Fig. 13, and hence, for creating the sign degeneracy in accelerator neutrino experiments. The appearance of the sin(V L/2) term is an inheritance from the magic baseline oscillations and will vanish the δ-dependent term when V L = 2π. Furthermore, we can observe that this term contains all of the mixing angles in the same way as the Jarlskog invariant, which is expected due to the CP-dependence of the term.
D. Determination of hierarchy with atmospheric neutrinos

Neutrino fluxes
The original flux of atmospheric neutrinos contains incoherent components of ν e , ν µ and the corresponding antineutrinos, while the original ν τ flux is negligible. We introduce Φ 0 e and Φ 0 µ , the electron and muon neutrino fluxes, as well asΦ 0 e andΦ 0 µ , the electron and muon antineutrino fluxes, at the detector in the absence of oscillations. The flavor ratios There is a mild neutrino-antineutrino asymmetry: the neutrino fluxΦ 0 µ /Φ 0 µ ≈ 0.8 − 0.9. All the fluxes (at E > 1 GeV) decrease rapidly with energy Φ 0 α ∝ E −k , k = k(E) = 3 − 5, and an azimuthal dependence shows up at low energies.
The flux of neutrinos of flavor ν α at a detector, with oscillations taken into account, is given by It is assumed here that experiments do not distinguish the neutrino and antineutrino events and corresponding signals are summed up. The fine-binned distribution of events (168) is shown in Fig. 14. For illustration we use the effective volume of PINGU with 22 additional strings [73]. which increases from ∼ 2 Mt at E ν = 2 GeV to 20 Mt at E ν = 20 GeV. The pattern of the event number distribution follows the oscillatory picture due to the ν µ − ν µ mode of oscillations with a certain distortion in the resonance region. The maxima and minima are approximately along the lines of equal oscillation phases E ν ∼ φ 32 ∆m 2 32 | cos θ z |R ⊕ (where R ⊕ is the Earth radius), with distortion in the resonance region E ν = (4 − 10) GeV. In the high density bins, the number of events reaches 200 and the total number of events is about 10 5 .
The expression for the density of events (169) can be written as where Similarly one can determine the number of events for inverted mass hierarchy. Let us introduce the N-I hierarchy asymmetry for the ij-bin in the (E ν − cos θ z ) plane as The moduli of the asymmetry (171) are the measures of statistical significance of the difference of the number of events for the normal and inverted mass hierarchies: S ij = |A ij |.
The strongest effect of hierarchy change is in the strips along the constant phase lines in the energy interval E ν = (4 − 12) GeV, where these lines are distorted by the matter effect. Here the asymmetry changes sign with the zenith angle and number of intervals with the same sign asymmetry increases with decrease of energy. The ν τ → τ → µ events can be considered as background events and treated within ∼ 5% systematic errors.

Measurements
According to Fig. 14, the hierarchy asymmetry of the ν µ events has opposite signs in different parts of the oscillogram. Thus, the integration over E ν and cos θ z substantially reduces the sensitivity to the hierarchy. Due to this, a relatively good reconstruction of the neutrino energy and direction are required to identify the hierarchy. The uncertainties of the reconstruction of energy σ E and angle σ θ should be comparable to or smaller than the sizes of the domains with the same sign of the asymmetry. The oscillograms for the reconstructed neutrino energy E r ν and angle θ r z can be obtained by smearing of the E ν − cos θ z oscillograms with the reconstruction functions of the width σ E and angle σ θ .
Small uncertainties σ E and σ θ require rather precise measurements of the energy E µ and direction θ µ of the muon, as well as energy of the accompanying hadron cascade E h . Then the neutrino energy equals E r ν = E µ + E h . The reconstruction of the neutrino direction is more involved. In the first approximation, one can use θ ν ≈ θ µ with a spread which decreases with energy: σ θ ∼ A m p /E ν (A = O(1)). Knowledge of the hadron cascade energy allows to reduce this uncertainty. Further improvements could be possible if some information about geometry of the cascade is available. A possibility to separate (at least partially) the neutrino and antineutrino samples would significantly improve sensitivity to the mass hierarchy, as well as to CP-violation.
All this imposes conditions on the detector characteristics. According to Fig. 14, the most sensitive region to the hierarchy is around the resonance and above: E = (5 − 15) GeV. The number of events in Super-Kamiokande is too small, but (upgraded) ice and underwater detectors of the multi-megaton (∼ 10 Mt) scale could collect around the order of 10 5 ν µ events a year in this range so that a high statistics study becomes possible.
A small enough spacing between the PMTs (∼ 10 − 20 m between strings and 3 -5 m in the vertical direction) will allow the reduction of the threshold down to a few GeV and perform reasonably good measurements of the muon and hadron cascade characteristics. Very high statistics will also allow the resolve the problem of parameter degeneracy: effects qualitatively similar to the mass hierarchy effect can be obtained by small (within 1σ interval) variations of ∆m 2 32 and θ 23 . The effect of an unknown CP-phase is small. High statistics would allow to resolve the degeneracy problem by selecting specific regions in the E ν − cos θ z for the analysis, where effects of ∆m 2 32 are suppressed in comparison to the hierarchy effects or averaged out as a result of specific integration. High statistics also allows to perform an analysis of the data using ∆m 2 32 and θ 23 as fit parameters. This will open a possibility to determine the mass hierarchy and measure these parameters simultaneously.
Note that other experimental techniques using atmospheric neutrinos may also prove valuable for determination of the mass hierarchy. In particular, experiments that can separate neutrinos from anti neutrinos on an event basis need a significantly lower number of events to obtain the same sensitivity. Thus, such detectors can be smaller in size as compared to the neutrino telescopes. In this context, a magnetized iron calorimeter, such as the India based Neutrino Observatory [74], could also provide an important contribution to the determination of mass hierarchy. The capabilities of detectors using charge identification were studied in [75].

Interplay between accelerator and atmospheric neutrinos
The atmospheric neutrino data can also be used to compliment the data from accelerator neutrino experiments in order to extract the most information possible. As was demonstrated in [75], the atmospheric neutrino determination of the neutrino mass hierarchy can be significantly affected by the addition of external priors and, in particular, may lead to different sensitivity to the neutrino mass hierarchy for in the cases of true normal or inverted hierarchy. However, once external input on the neutrino oscillation parameters is included by considering also other experiments, the room to mimic the true oscillation pattern in the wrong hierarchy becomes much more restricted and the sensitivity to the hierarchy increases. Adding the accelerator experiments' own sensitivity to the mass hierarchy, a measurement may be possible even for the current generation of accelerator experiments by the addition of detector capable of lepton charge identification. This has been discussed in [76] and the prospects of using a magnetized iron calorimeter detector to augment the current generation of accelerator experiments are a 2-4σ determination of the mass hierarchy within 10 years of data taking, depending on the true value of the oscillation parameters and the characteristics of the detector.

V. DISCUSSION AND CONCLUSIONS
In this paper, we have described the effects of neutrino propagation in matter relevant for experiments with atmospheric and accelerator neutrinos and aimed at the determination of the neutrino mass hierarchy and CP-violation. Thus, to a large extent, we have focused on neutrino propagation in the Earth matter.
1. At relatively low energies, the dominant effect of neutrino interactions with matter is the elastic forward scattering, which is described by an effective potential. Neutrino evolution in matter is then described by a Schrödinger-like equation including this effective potential. The potential differences for neutrinos of different types influence the flavor evolution of the system of mixed neutrinos.
In the majority of realistic situations, neutrinos propagate in normal (unpolarized non-relativistic) matter with nearly constant or slowly changing density. 2. Matter modifies the neutrino flavor mixing and changes the eigenvalues of the Hamiltonian of propagation. This is equivalent to a modification of the dispersion relations of neutrinos. The influence of matter on mixing of neutrinos has a resonance character. At energies or densities for which the eigenfrequency of the neutrino system with mixing ω ij = ∆ 2 ij /2E equals approximately the eigenfrequency of the medium 2π/l 0 , the mixing in matter becomes maximal. Large mixing shifts the position of the resonance to lower values of the potential. At usual densities, there are two resonances related to the two mass squared differences ∆m 2 21 and ∆m 2 31 between the neutrino mass eigenstates. The resonances are realized in oscillation channels involving electron neutrinos.
3. In many practical situations, knowledge of neutrino mixing in matter and the eigenstates of the Hamiltonian in matter allows to immediately find the results of the neutrino flavor evolution. This includes neutrino oscillations in matter with constant density and also adiabatic conversion of neutrinos, where the averaged oscillation results can be written down immediately. In the non-averaged case, the problem is reduced to finding the oscillation phase (integrating the energy splittings over distance). In this sense the Nature has implemented the most (computationally) simple setups. The very convenient presentation of mixing in matter can be obtained as series expansion in the ratio of the two mass squared differences, r ∆ , (perturbative diagonalization of the effective Hamiltonian), which allows to understand a number of subtle results.
The simplest and physically transparent description of dynamics of neutrino flavor evolution can be obtained in the propagation basis (in the case of the standard parameterization). In this basis, the CP-violating phase and 2-3 mixing do not influence the evolution and the amplitudes of transitions do not depend on δ or θ 23 . The dependence on these parameters appear as a result of projecting the states of the propagation basis back to the flavor states at production and detection.
In many practical cases the 3ν evolution can be reduced to evolution of two neutrino systems with certain corrections. 4. There are two practically important cases: (i) neutrino propagation in matter with constant or nearly constant density and (ii) neutrino propagation in matter with slowly (adiabatically) changing density. 5. In the case of constant density, flavor evolution has a character of oscillations with parameters determined by mixing and mass splitting in matter. The oscillations are an effect of a phase difference increase in the course of neutrino propagation. The resonance enhancement of oscillations is realized in an energy region around E R .
If the density is approximately constant, then the results can be obtained by using perturbation theory in the deviation of the density distribution from a constant one. The accuracy improves if the density profile is symmetric with respect to the middle point of the neutrino trajectory, as is realized for neutrinos crossing the Earth.
A simple and rather precise semi-analytical description of neutrino oscillations in matter with varying density can be obtained in the limits of small density, V < ∆m 2 ij /2E, and high density V ∆m 2 ij /2E. The latter gives a very accurate description of neutrino flavor evolution in the Earth at E > (8 − 10) GeV. 6. In a medium with slowly changing density, adiabatic conversion takes place. This effect is related to the change of mixing in matter due to density change. Adiabaticity implies that there are no transitions among the eigenstates of the instantaneous Hamiltonian during propagation.
The strongest flavor transformation is realized when the initial density is much larger, and the final one is much lower than the resonance density. In this case, the initial state (and due to adiabaticity, the state at any other moment of evolution) practically coincides with one of the eigenstates. Therefore, oscillation effects are absent and non-oscillatory flavor conversion takes place. This is realized for supernova neutrinos and approximately -for high energy solar neutrinos. In general, if the initial mixing is not strongly suppressed, an interplay of adiabatic conversion and oscillations occurs.
Adiabatic transformations are also realized for neutrinos with energy ≤ 1 GeV propagating in the mantle of the Earth. In particular, this means that the oscillation depth at the detector is determined by mixing at the surface of the Earth and not by the mixing at average density.
Until now, the mater effects have been observed in solar neutrinos and, indirectly, in atmospheric neutrinos and there is good chance that they will be observed by new generation of the accelerator and atmospheric neutrino experiments. 7. Strong flavor transition can be realized without enhancement of mixing. This occurs in matter with periodic or quasi-periodic density change when the parametric resonance condition is fulfilled. For small mixing strong transition requires a large number of periods.
A similar enhancement can take place in matter with several layers of different densities. Here the enhancement occurs when a certain correlation between the oscillation phases in each layer and amplitudes of oscillations determined by mixing is present. The case of a medium with 3 layers (1.5 periods) is of practical interest for neutrinos crossing both the mantle and the core of the Earth.
For a multilayer medium two conditions must be satisfied to have strong transitions: the amplitude (collinearity) and the phase conditions. 8. For neutrinos crossing a small amount of matter, such as accelerator experiments with baselines up to (1−2)·10 3 km, the column density of matter is small and, according to the minimal width condition, the matter effect on oscillations is small regardless of energy, vacuum mass splitting, and neutrino mixing. Furthermore, if the oscillation phase is small, then mimicking of vacuum oscillations occurs. 9. A comprehensive description of the neutrino flavor transitions in the Earth is given in terms of neutrino oscillograms of the Earth. After the recent determination of the 1-3 mixing, the structure of oscillograms is well fixed. The salient features of oscillograms at high energies (due to 1-3 mixing) are the MSW resonance peak in the mantle domain, three parametric ridges and the MSW peak in the core domain. At low energies (due to 1-2 mixing), there are three peaks, due to the MSW resonance, and the parametric ridge. The positions of all these and other structures are determined by the generalized phase and amplitude conditions.
In the case of normal mass hierarchy, the resonance peaks induced by the 1-3 mixing are in the neutrino channels. For inverted mass hierarchy they are in the antineutrino channels. This is the foundation for determining the neutrino mass hierarchy. The resonance structures due to the 1-2 mixing are always in the neutrino channels, since the sign of the small mass square difference has been fixed. 10. The CP-properties of the oscillograms (their dependence on CP-phase) are determined by the CP-domains: areas in which the CP-violation effect has the same sign. The borders of these domains are approximately determined by the grids of the magic lines (solar and atmospheric magic lines) and the lines where the oscillation phase condition is fulfilled. 11. Measurements of matter effects in neutrino oscillations provides a good opportunity to determine the neutrino mass hierarchy. The 1-2 ordering has been determined due to the matter effect of solar neutrinos. The 1-3 ordering can be identified by studying the matter effects in accelerator and atmospheric neutrino experiments.
There is a good chance that future studies of the atmospheric neutrinos with multi-megaton underwater (ice) detectors will be able to establish the mass hierarchy. With a threshold of a few GeV, these detectors will be sensitive to the resonance region (∼ 6 − 10) GeV, where the difference of probabilities for the normal and inverted mass hierarchies is maximal.
The challenges here are the accuracy of reconstruction of the neutrino energies and directions. Integration over the energy and angle, as well summation of neutrino and antineutrino signals, diminish the sensitivity to the hierarchy. Another problem is the degeneracy of the hierarchy effects with the effects of other neutrino parameters, in particular with ∆m 2 32 and θ 32 . 12. In accelerator experiments, many of the problems mentioned above are absent. However, existing and proposed accelerator experiments will cover only periferal regions of oscillograms where enhancement of oscillations is very weak and oscillatory structures are rather poor. As a consequence the problem of degeneracy here is even more severe.