Analytical description of CP violation in oscillations of atmospheric neutrinos traversing the Earth

Flavour oscillations of sub-GeV atmospheric neutrinos, traversing different distances inside the Earth, are a promising source of information on the leptonic CP phase $\delta$. In that energy range, the oscillations are very fast, far beyond the resolution of modern neutrino detectors. However, the necessary averaging over the experimentally typical energy and azimuthal angle bins does not wash out the CP violation effects. In this paper we derive very accurate analytic compact expressions for the averaged oscillations probabilities. Assuming spherically symmetric Earth, the averaged oscillation probabilities are described in terms of two analytically calculable effective parameters. Based on those expressions, we estimate maximal magnitude of CP-violation effects in such measurements and propose optimal observables best suited to determine the value of the CP phase in the PMNS mixing matrix.


Introduction
Determination of the leptonic CP phase by measuring neutrino oscillations is a challenging issue [1][2][3][4][5][6][7][8][9]. It is well known that sensitivity of oscillations to the CP phase δ generically decreases with the increasing neutrino energy. Matter effects may be helpful in measuring δ but they also fade away when the neutrino energy increases [10]. Thus the oscillations of low energy atmospheric neutrinos, hitting a detector at different angles, after traversing different distance inside the Earth, look as a particularly promising source of information on the leptonic δ [9]. However, in that case the limitations come from the difficulties with precise determination of the neutrino energy and the angle it hits a detector. Therefore, analytical understanding of the oscillation probabilities for low energy (say, below O(1)GeV ) atmospheric neutrinos, as a function of their energy and the number of layers they traverse in the Earth, would be very useful for optimising measurements of the leptonic CP phase in realistic experimental setups. This is the purpose of the present paper.
Oscillations of sub-GeV neutrinos differ substantially from those of higher energy neutrinos [3,[11][12][13][14][15][16][17]. First of all, they are very fast in energy, far beyond the energy resolution of modern neutrino detectors, [18,19], because they are affected by both solar and atmospheric mass splittings. Thus, the relevant "observables" carrying the physical information are the oscillation probabilities averaged over typical experimental energy and angle bins. In addition, the patterns of matter effects also change with energy [20][21][22][23]. A useful insight can be obtained from the description of oscillation probabilities in matter in the conventional parametric form as in the vacuum but with effective mixing angles and mass eigenvalues [24,25]. The main effect resides in the energy dependence of effective mixing angles θ m 12 and θ m 13 . At sub-GeV energies and the matter densities typical for the Earth structure, θ m 13 is close to, and θ m 12 is significantly different from their vacuum values, whereas the opposite is true at higher energies (in that parametrization θ m 23 remains to be the vacuum angle).
In this paper we derive analytical parametrization of the averaged oscillation probabilities for sub-GeV neutrinos, after traversing arbitrary number of Earth layers, each with a constant matter density. Assuming spherically symmetric Earth, the averaged oscillation probabilities are described in terms of two effective parameters. Based on those expressions, we estimate maximal magnitude of CP-violation effects in such measurements. We also propose optimal observables best suited to determine the value of the CP-phase in the PMNS mixing matrix.
Our article is organised as follows. In Section 2 we derive the exact formulae for the transition matrix for neutrinos traversing the Earth, divided into layers of constant matter density. In Section 3 we propose the approximations which can be done for the considered neutrino energy range and assuming symmetric Earth layout and we derive simple and accurate analytical formulae for the averaged oscillation probabilities. Section 4 is devoted to the discussion of the optimal experimental setup and choice of observables best suited to measure the leptonic CP-phase. In Section 5 we discuss the issue of averaging oscillation probabilities over the experimental bins in energy and azimuthal angle. We conclude in Section 6. In Appendix A for completeness we collect the formulae for the neutrino track lengths in the Earth layers. Finally in Appendix B we discuss the numerical quality of the approximations done when deriving the analytical formulae.
2 Oscillation probabilities for neutrinos traversing the Earth .
We consider neutrino oscillations when traversing the Earth. Our main focus is on sub-GeV atmospheric neutrinos but the framework we develop in this Section is a general one. In the next Section we shall discuss the approximations appropriate for low energy neutrinos.
In order to estimate possible effects of the CP phase in the PMNS mixing matrix on the transition probabilities, we calculate them analytically assuming the Earth structure based on the PREM model [26]. In such an approximation the Earth is divided into a finite number of layers, each having a constant density of matter. Although our analysis can be applied to any number of layers, for numerical estimates we use 5-layer pattern of the structure of our planet -starting from the center, one has inner core, outer core, lower mantle, upper mantle and crust. Our schematic setup is illustrated in Fig. 1, and the layer radii and densities are collected in Table 1. Depending on the azimuthal angle θ, neutrinos can traverse 1,3,5,7 or 9 Earth layers. Other numerical inputs used throughout the paper are collected in Table 2 Figure 1: Schematic picture of the Earth structure (not to scale) and definition of the azimuthal angle θ. The detector, marked by black blob, is located at depth h below the Earth surface, the average atmosphere width is denoted by a.
The neutrino oscillation probabilities are determined by the S-matrix elements (α, β = e, µ, τ ): where U denotes the neutrino mixing matrix in the vacuum, with the parametrization: The V (x) is the neutrino weak interaction potential energy in matter. As shown in ref. [24], for a constant V , to a good approximation the Hamiltonian H can be diagonalized by the matrix U m where only the θ 12 and θ 13 angles have values modified compared to vacuum: and approximate explicit formulae for O m 13 , O m 12 are given in ref. [24]. It is convenient to work in a new basis, rotated by the matrix so that the rotated Hamiltonian has the form The transition matrix S can written as the time-ordered product of the transition matrices in the Earth layers, where within the i-th layer of constant density the matrix S i is simply given by Using the rotated basis defined above, one can easily show that (up to an unimportant overall phase denoted as e iξ ) the matrix S can be expressed as where we have defined Formula (2.12) is general and does not involve any approximation yet (other than the "layered Earth" model). In the next Section we introduce analytical approximations appropriate for the oscillation probabilities of the sub-GeV atmospheric neutrinos. 3 Analytical approximations for sub-GeV atmospheric neutrinos

Averaging of probabilities over energy bins
The transition probabilities for sub-GeV atmospheric neutrinos oscillate quickly with neutrino energy and azimuthal angle. This can be traced back to the fact that small variations of both quantities can significantly change the ratio ∆m 2 a L(θ)/E. Realistically, the oscillation probabilities have to be averaged over bins in energy and angle corresponding to the relevant experimental resolutions. As long as the period of the neutrino oscillation frequency is far smaller than the experimental resolution significant simplifications can be performed in calculating analytically the averaged oscillation probabilities.
First, we observe that in the product (2.12) the following structure repeats itself: with the most inner multiplication matrix depending on the differences of the θ m 13 mixing angle between the neighbouring layers: Contrary to high energy neutrinos, like for instance in the Dune experiment [24], for the neutrino energies below E < O(1) GeV and typical values of the Earth density, θ m i13 angle in matter vary only very slightly, as illustrated in Fig. 2 between the layers are typically of the order of 0.01 radian, even less for the lower neutrino energies. Therefore, to a good approximation products of O mT i13 O m (i+1)13 can be replaced by the unit matrices. In contrast, the dependence of the θ m i12 on the matter density is stronger for this energy range.
Then, neglecting the overall phase, the time-ordered product on the RHS of eq. (2.12) takes the form First layer on the neutrino track is the atmosphere, so that O m 13−f irst ≈ O 13 in vacuum. The last layer is the Earth crust around the detector, so that O m 13−last = O m 13−crust . The inner product in eq. (3.3) contain only O m i12 mixing matrices, thus the result has the structure: Approximating again the product O mT 13−crust O 13 by the unit matrix we arrive at the following expression for the transition matrix S (up to an unimportant overall phase factor): The matrix A = A(E, θ) can be calculated by the numerical diagonalization of the Hamiltonian given in eq. (2.2) and the matrix B is constant, given by the vacuum mixing angles. Oscillation probabilities are given by In eq. (3.6) only the quantity Π i (E i ) 33 (being a pure phase) depends on the larger neutrino mass splitting and varies quickly with energy and azimuthal angle. When averaged over bins in energy ∆E and angle ∆θ larger than the period of the oscillation frequency, that term vanishes and we get This is the first important result of the paper -the averaging over energy and azimuthal angle can be now done using some standard 2-dimensional numerical integration techniques, expected to be quickly converging and accurate as the numerically most difficult and CPU-time consuming averaging over fast oscillations of probabilities has been done analytically while obtaining the formulae 3.7. However, as we show below, one can also derive for the matrix A, and thus for the integrand in eq. (3.7, an excellent analytical approximation in terms of of only two effective parameters. Furthermore, for neutrino energies larger than 300 − 400 MeV, they are very accurately calculable analytically.

Analytical results for the matrix A
We begin with the discussion of the properties of the matrix X. The full 2 × 2 matrix X defined in eq. (3.3) is a time-ordered product of matrices of the form for each Earth layer). With the phase conventions chosen in eq. (2.13), each of matrices X i is unitary, symmetric and have determinant equal to 1. Any such matrix has only 2 free real parameters and can be expressed as: direct calculations lead to the formulae so that comparing with eq. (3.8) one has Let us note that excluding azimuthal angles close to π/2 or bigger (when the length of the neutrino track in the atmosphere and the asymmetric position of the detector under the Earth surface cannot be neglected) our setup is symmetric with respect to the Earth center. Thus, the full matrix X is to a good approximation given by a symmetric product of X i and has the same symmetry properties as each of them separately: The quality of this approximation turns out to be very good, as discussed in Appendix B.2.
In the next step, one can obtain analytical formulae for the angles φ X , α X . We observe that (using the formulae from ref. [24]) in the limit of large E V product the quantity sin 2θ m 12 can be expanded as sin 2θ m 12 = cos θ 13 sin 2θ 12 2 cos 2 2θ 13 Therefore, for increasing energy sin 2θ m 12 is suppressed approximately by 1/(E V ) factor (as illustrated in Fig. 3). For sufficiently high values of E one can expand the product (3.12) using the expression for X i matrices given by eq. (3.10) and keeping at most the terms linear in i ≡ sin 2θ mi 12 . Direct multiplication leads then to the remarkably compact formulae Such an approximation works well even for relatively low energies and large values of sin 2θ mi 12 . This can be attributed to the fact that the neglected higher order terms in eq. (3.17) are suppressed by additional 2 i factors. Eq. (3.17) has some remarkable properties. Firstly, it shows that the overall neutrino oscillation phase is just a direct sum of phases in all layers. In addition, using the approximate formulae from ref. [24] (working very well for the considered energy and neutrino potential ranges), one sees that with growing value of neutrino energy the difference between the eigenvalues of the Hamiltonian defined in eq. (2.2) becomes constant: Already for E > 300 − 400 MeV corrections to the constant term are small and the phases ν i in eq. (3.17) (thus also the overall phase φ X ) depend to a good approximation only on the azimuthal angle and the Earth layers density: where the explicit formulae for the oscillation lengths ∆x i (θ) are given in Appendix A. Secondly, since i ≡ sin 2θ mi 12 ∼ 1/E and the phases ν i become energy independent, for E > 300 − 400 MeV to a good approximation one has sin α X = f (θ)/E, where f (θ) is some function of the azimuthal angle only.
The dependence of φ X and sin α X on energy and azimuthal angle and the comparison of numerical fitting (see Appendix B.2) and analytical approximate formula (3.17) for both parameters is illustrated in Fig. 4 (where the normal neutrino mass ordering is assumed). As can be seen, for E > 300 − 400 MeV numerical and analytical results agree very well. The dependence of the φ X and sin α X on the angle θ becomes universal with energy, up to an overall 1/E scaling of the sin α X amplitude. The dependence of φ X and sin α X on the azimuthal angle for the inverse mass ordering is almost identical, with small differences of the order of few % appearing only for small values of θ.

Energy and the angular dependence of the oscillation probabilities
The arguments presented in Section 3.2 show that, for the neutrino energy larger than 300 − 400 MeV, the angle φ X and the product E × sin α X should be to a good approximation functions of the azimuthal angle only. These properties are well confirmed in Fig. 5. As already mentioned in the previous Section, the dependence of φ X and E × sin α X on the azimuthal angle is almost identical for the normal and inverted neutrino mass ordering, thus the problem of the CP-phase determination is not affected by the assumption of the mass hierarchy. Using the parametrization of eq. (3.15), for energies larger than 300−400 MeV one then obtains very compact expressions for the averaged oscillation probabilities. Assuming the central values for the measured vacuum mixing angles (and the normal mass ordering), for the most important ν µ → ν e transition probability we get the simple but very accurate approximation: I µe ≈ 0.024 + 0.450 sin 2 α X − 0.0724 sin 2α X sin(δ + φ X ) (3.21) As an immediate consequence of eq. (3.21) we observe that the variation ofP µe with the CP phase δ cannot be larger than cos 2 θ 13 sin θ 13 sin 2θ 23 ≈ 0.15. Furthermore, the effects of CP violation, which as can be seen from eq. (3.15) are proportional to sin 2α X , decrease approximately like 1/E. For any given energy E and the azimuthal angle θ the parameters φ X and sin α X are calculable using either the numerical diagonalization of neutrino Hamiltonian and fitting procedure described in Appendix B.2 or, for sufficiently large E, the approximate formulae of eq. (3.17) and they depend only on the assumed Earth density profile. Therefore, as follows from eq. (3.15), one can subtract the theoretically known CP-independent terms from the experimentally measured ν µ → ν e transition probability and obtain directly the constraints on the combination δ + φ X (E, θ). Performing a fit to many bins in energy and azimuthal angle one can determine the value of phase δ itself. 4 Optimal observables for the CP-phase detection 4

.1 Optimal azimuthal angles
The experimental chances of measuring the CP-phase in ν µ → ν e transitions are best when the coefficient of sin(δ + φ X ) is maximal. This happens when the sin 2α X reaches maximal or minimal value. In Fig. 6 we plot the dependence of sin 2α X (E, θ) as a function of the azimuthal angle for few chosen values of neutrino energy. As can be seen, independently of the neutrino energy, the extreme values of sin 2α X (E, θ) are reached for three values of the azimuthal angle θ 1 = 0.12π, θ 2 = 0.18π, θ 3 = 0.39π and these three values give the best chance for a successful measurement. As the extreme values of sin 2α X (E, θ) and sin 2 α X (E, θ) (the latter deciding about the magnitude of the CP-independent term in eq. (3.15)) occur for approximately similar values of the azimuthal angle, maximal CP effects are correlated with the maximal values of the totalP µe probability. This is illustrated in Fig. 7, where we plot the dependence of the quantity I µe on the azimuthal angle for E = 400 MeV and several values of the phase δ assuming the normal neutrino mass ordering (corresponding plot for the inverse neutrino mass ordering is almost identical.) For the angles θ 1 and θ 3 , corresponding to maxima of sin 2α X , one has φ X (θ 1 ) ≈ φ X (θ 3 ) and sin α X (θ 1 ) ≈ sin α X (θ 3 ). Therefore, measurements done for θ 1 and θ 3 provide information on almost the same combination of δ + φ X . For the angle θ 2 we get sin α X (θ 2 ) ≈ − sin α X (θ 1 ) but different phase φ X , thus combining measurements for all three azimuthal angles gives a chance for determining the phase δ itself. Discussed effects are illustrated in Fig. 8, where we plot the dependence of I µe as a function o the CP-phase for E = 400 MeV and optimal angles θ 1 , θ 2 , θ 3 . As expected, in this case the variation of I µe due to the phase dependence reaches maximal allowed value of 0.15. The maxima and minima of the line corresponding to angle θ 2 are shifted compared to other two lines due to different value of φ X (θ 2 ).  Figure 8: I µe as a function of the CP-phase for E = 400 MeV and optimal angles θ 1 , θ 2 , θ 3 (blue, yellow and green line, respectively). Normal mass ordering is assumed.

Observables optimised for the CP-phase detection
As discussed in the previous Section, eq. (3.15) can be used to determine the phase δ by subtracting from the experimentally measured transition probability the theoretically calculated CP-independent terms. The analytical understanding of φ X and sin α X behaviour allows us to design the alternative observable best suited to measure the CP violating phase, based on subtracting experimentally measured quantities and a well know combination of neutrino vacuum mixing angles. Inspection of the expression for ν µ → ν e probability shows that I µe in eq. (3.15) consist of a constant term depending only on the vacuum oscillation angles, term proportional to sin 2 α X but not depending on the phase δ (which for E > 300 − 400 MeV scales to a very good accuracy like 1/E 2 ) and a term proportional to sin(δ + φ X ) which scales approximately like 1/E. Therefore, for E 1 , E 2 > 300 − 400 MeV and for any azimuthal angle θ the quantity is to a good approximation proportional solely to the sine of the CP-violating phase shifted by φ X . To maximise ∆I µe (E 1 , E 2 , θ), one can choose θ equal or close to the values maximising | sin 2α X |, as described in the previous Section, and large splitting between E 1 and E 2 , like e.g. E 1 = 400 MeV, E 2 = 1000 MeV.
To illustrate the dependence of ∆I µe (E 1 , E 2 , θ) on the phase δ, in Fig. 9 we plot it for chosen values of energy and azimuthal angles. As expected, the dependence on δ resembles pure sine function symmetric with respect to the horizontal axis.

Measurements with the finite energy and angular resolution
Using the approximate explicit energy scaling properties holding well in the energy range above 300 − 400 MeV: we can estimate the effect of averaging over hypothetical experimental bins. The integral over energy in eq. (3.13) can be calculated analytically: For the energy resolutions achievable at Dune [18] or HyperK [19] experiments the term is small and can be neglected. Integration over the azimuthal angle leads then to the following formula: P µe = 2 cos 2 θ 13 sin 2 θ 13 sin 2 θ 23 + cos 2 θ 13 (cos 2 θ 23 − sin 2 θ 13 sin 2 θ 23 ) η(E, θ) − 1 2 cos 2 θ 13 sin θ 13 sin 2θ 23 (ξ 1 (E, θ) sin δ + ξ 2 (E, θ) cos δ) where the function η(E, θ), ξ 1 (E, θ) and ξ 2 (E, θ) are given by the integrals Since in the presented formalism α X and φ X are known and regular functions of the neutrino energy and the azimuthal angle, for any value of E and θ the coefficients η, ξ 1 , ξ 2 can be easily calculated by simple 1-dimensional numerical integration. Therefore measurements done for different angular momentum bins can provide information on different (but known) combinations of sin δ and cos δ, ultimately giving a good chance to measure the CP-phase itself. Obviously, if necessary one can evaluate them also assuming more complicated Earth structure models including more internal layers, like the full PREM model [26]. The same procedure can be applied to averaging of the quantity ∆I µe , defined in , the same cancellations between terms as in eq. (4.1) occur for barred probabilities and we can define observable ∆P µe averaged over the energy and angular bin as where the functions ρ 1 (E 1 , E 2 , θ), ρ 2 (E 1 , E 2 , θ) are defined as

Summary
We have investigated flavour oscillations of neutrinos created in the atmosphere by cosmic ray interactions with the air and traversing the Earth. We have focused on sub-GeV neutrinos (E < O(1) GeV) where CP violation effects are large but the oscillation probabilities vary very fast with neutrino energy and its azimuthal angle, far beyond the typical experimental resolution. Therefore, the "observables", carrying the physical information, are the averaged probabilities, where the oscillation pattern is washed out. Using the Earth model with layers of constant matter density, we have derived very simple analytic formulae for those averaged probabilities. There are three main formulae summarising our results. Equation (3.7) is the most general expression suitable for fast numerical calculations of the oscillations probabilities averaged over any experimental bins larges than the oscillation periods. Equations (3.13-3.15) give very accurate approximation to the averaged probabilities, where all matter effects are encoded in two effective parameters. And finally, eqs. (3.17,3.19) provides for the neutrino energies larger than 300 − 400 MeV approximate simple analytical expression for these effective parameters.
The obtained analytical parametrization is very accurate when compared with the exact numerical calculations. It opens up the possibility of better understanding the dependence of the averaged flavour oscillations of sub-GeV atmospheric neutrinos as a function of their energy and the azimuthal angle with which they hit the detector. In turn, our results can be useful in optimising the experimental measurements of the leptonic CP phase in oscillations of sub-GeV atmospheric neutrinos. We have made several suggestions in that direction, such as the best choice of the azimuthal angles or taking combinations of the data that are directly measuring the CP phase.

A Oscillation lengths
.
For completeness we include expressions for the length of the neutrino tracks in Earth layers and in the atmosphere. We consider the latter because despite the fact that neutrinos passing through the atmosphere only do not have time to oscillate, they can be important for azimuthal angle θ ≈ π/2 since they come from full 360 degree plane, while those passing through Earth core come only from the small cone. Thus atmospheric-only neutrinos may produce serious background.
Calculating track lengths is a straightforward exercise in trigonometry. We assume setup defined in Fig. 1, with detector at distance h below Earth surface (it is 1600m for Dune and 650m for HyperK) and atmosphere width denoted by a. Obviously h, a r i , R thus, in all expressions below we neglect quadratic terms h 2 , a 2 . Let's consider 3 cases: 1) Neutrino track length in the atmosphere.
Let's define θ i as angles for which neutrino track is tangent to i-th layer: Then for θ ≥ θ 2 neutrino has in 1st layer single undivided track with the length For θ ≤ θ 2 track has 2 parts, next to detector and on the opposite side of Earth: 3) Neutrino track length in inner layers.

B Quality of analytical approximations B.1 Averaged oscillation probability
In order to test the quality of approximation of eq. (3.7), we employ the following procedure.
1. We numerically diagonalize Hamiltonian of eq. (2.2) in each Earth layer and calculate the full transition matrix without any approximations, multiplying layer transition matrices as in eq. (2.10). Resulting transition probability, P (E, θ) of eq. (3.6), is exact but exhibits fast variations with neutrino energy and with the azimuthal angle.
2. We average P (E, θ) over energy with the use of numerical integration, using the formulaP Averaging is done approximately over 4 periods ∆E of "fast" oscillations in energy, which (in vacuum) are given by where L(θ) is the total neutrino track length in Earth for a given azimuthal angle. Actual value of ∆E in matter differ from the vacuum, but tests show that the result of numerical averaging is stable against variations of ∆E as long as it has the correct order of magnitude and we integrate over several (here 4) periods ∆E.
3. For each Earth layer we diagonalize numerically Hamiltonian H and calculate relevant transition matrix S i in rotated basis of eq. (2.8). We assume the upper 2 × 2 sub-block of S i to be matrix X i , as defined in eq. (3.4). Further, we evaluate full matrix X as a time-ordered product of X i (see eq. (3.12)). Finally, knowing matrix X and hence also the matrices A, B defined in eq. (3.5), we calculate the quantity I (see eq. (3.14)). Finally, for the analytically averaged oscillation probability we use the approximationP (E, θ) ≈ I(E, θ), as discussed in Sec. 4.
The comparison of P ,P andP is illustrated in Fig. 10. In general, analytical average of eq. (3.7) works very well, some differences betweenP andP can be attributed more to the inaccuracies in numerical integration rather then in the approximations used when deriving the formula (3.7).
B.2 Numerical fits for α X (θ, E) and φ X (θ, E) angles Matrix X obtained numerically as a 2 × 2 sub-block of full 3 × 3 transition matrix (as described in point 3 of the previous Section) is only approximately unitary and symmetric and has determinant slightly different from unity. We obtain best values of angles φ X , α X minimising the difference between the symmetric form on the RHS of eq. (3.8) and the X Such a procedure reproduces very well X matrix derived by numerical diagonalization. Fig. 11 shows the relative error of a fit as a function of E and θ. The error is defined as with X(α X , φ X ) defined in eq. (3.12). Figure 11: Relative error ε X of α X , φ X fit plotted as a function of energy varied from 200 to 1000 MeV and azimuthal angle varied from 0 to π/2.
As one can see, only for low energies E < 400 and azimuthal angle close to π/2, where the asymmetry of underground detector position and neutrino track in atmosphere becomes relevant, the error can reach 3-4%. For smaller θ angles it is always small, confirming the assumed analytical symmetry properties of X matrix and justifying the approximations done in derivation of eq. (3.12).