Earth tomography with atmospheric neutrino oscillations

The study of the flux of atmospheric neutrino crossing the Earth can provide useful information not only on the matter density of the different layers that make up the planet but also on their chemical composition. The key phenomenon that makes this possible is flavor oscillations and their dependence on the electron density along the neutrino baseline. To extract the relevant information, we simulate the energy and azimuth angle distribution of events produced in a generic neutrino telescope by atmospheric neutrinos passing through the deepest parts of the Earth. Changes in the densities of the outer core and the mantle are implemented by varying the location of the boundary between these layers so that the restrictions on the mass of and the moment of inertia of the Earth are both satisfied. This allows us to examine the effect of simultaneous changes in composition and density of the outer core, unlikely other works on the subject, where only one of these quantities was varied.


Introduction
To understand the internal dynamics and evolution of the Earth we need to know its internal structure and composition. These are essential ingredients for a full comprehension of basic geological phenomena, such as volcanology, earthquakes, plate tectonics and mountain formation [1,2]. Current knowledge has led to symmetric spherical models consisting of several concentric shells that differ in composition and mechanical behavior [3], of which the main ones are: crust, upper mantle, lower mantle, outer core, and inner core. Characterizing their properties has not been an easy a e-mail: dolivo@nucleares.unam.mx b e-mail: fis_pp@ciencias.unam.mx (corresponding author) c e-mail: ismaelromero_@hotmail.com d e-mail: sampayo@mdp.edu.ar e e-mail: zapatagabriel@hotmail.com task. Except in the crust, technological difficulties prevent inspecting the deeper layers by drilling holes to take samples. Instead, several type of observations on the surface need to be combined to probe the Earth's inner parts. Most of our knowledge comes from geophysics, in particular examining how the direction and/or speed of seismic waves generated by earthquakes are affected when they travel through different layers [4][5][6]. Valuable information has also been obtained from measurements of the magnetic and gravitational fields, observations of the moment of inertia and the precession motion of the planet, laboratory experiments at high pressures and temperatures, and physical, chemical and mineralogical analysis of meteorites and xenoliths.
While the distribution of the matter density can be inferred from seismological observations, the compositional structure of the Earth is more difficult to determine. This is particularly true for the lower mantle and the core, whose compositions remain quite uncertain despite significant progress in recent years. According to the most popular model, the core mainly consists of an iron-nickel alloy, with Ni/Fe ∼ 0.06, and is divided into inner and outer regions distinguished by a great density difference at a depth of approximately 5100 km [7]. The inner core is solid, while the outer core is liquid as indicated by the lower density and lack of s-wave propagation in it [3]. The density deficit of 4.5% in the outer core is too large to be due only to the solid-liquid phase transition and requires elements of lower atomic weight, such as Si, O, S, C, and/or H, to be present in about 5-10 wt% (weight percent) [8][9][10][11]. The nature and content of the light elements has important implications for convection in the outer core, which in turn is closely related to the geodynamo that generates the planetary magnetic field. They are also linked to the scenario of Earth's differentiation, the rate of core cooling, and the way the core and the mantle interact. It has been challenging (and not free of controversy) to ascertain the identity and concentrations of the light elements. In order to discrim-inate between the variety of plausible models, the most direct method should be to compare laboratory measurements with the seismically observed density and acoustic velocity of the core. These experiments are difficult to conduct for liquid iron alloys samples at the extreme temperature and pressure conditions of the core [12]; as consequence, the available accurate experimental data are still insufficient to make a direct comparisons with the seismological data. On the other hand, in a recent work, the chemical composition of liquid iron alloys at outer-core pressure and temperature conditions has been constrained by comparing first principle calculations of the density and bulk sound velocity with the PREM profiles [13]. The primary light component was found to be hydrogen when the inner-core boundary temperature T I C B is not high (4800-5400 K) or oxygen if T I C B = 6000 K. Hence, as the authors conclude, it is still difficult to specify the chemical composition of the outer core based solely on comparing theoretical calculations and seismological observations of the density and speed of sound.
From the above, it is clearly important to develop alternative methods that can provide complementary and independent information about the deep interior of the Earth. In this regard, the use of neutrinos seems to be a promising and viable option. In general, there are two possible ways to perform geophysical investigations with neutrinos: geoneutrinos and neutrino tomography. Geoneutrinos are electron antineutrinos from the β-decay of long-lived natural isotopes within the Earth. Measurements of their flux provide valuable evidence of the amount and distribution of the radio active elements internally heating the planet [14,15]. Two experiments, Borexino [16] and KamLAND [17], have currently sensitivity to geoneutrinos. In respect to neutrino tomography, the underlying idea is that the propagation of neutrinos in a medium is affected by their interactions with matter [18]. In neutrino absorption tomography the density profile can be reconstructed from the attenuation of the flux of very high energy ( 10 TeV) neutrinos, which are absorbed when passing through the Earth. Measuring at large neutrino telescopes the isotropic flux of ultra-high energy cosmic neutrinos seems to be a promising tool to make an independent determination of the Earth's density profile [19][20][21][22][23][24][25]. On the other hand, neutrino oscillation tomography takes advantage of the matter effect on neutrino oscillations [26][27][28] to probe the Earth's interior with lower energy (MeV to GeV) neutrinos. In matter the flavor transition probabilities depend on the number density of electrons n e along the neutrino baseline, which is proportional to the product of the matter density ρ times the average ratio of the atomic number Z to the mass number A [18,[29][30][31][32][33][34][35].
Since the Z /A ratio is a function of the chemical and isotopic composition of the medium, from the determination of the densities of matter and electrons, it might be possible in principle to restrict the different composition models of the Earth. With this in mind, in this work we reexamine the feasibility of studying the internal structure of the planet by means of oscillation tomography with atmospheric neutrinos, namely, those produced in the atmosphere by the interactions of cosmic rays with air nuclei. We pay particular attention to how changes in the content of light elements in the outer core could influence the number density of electrons. At the same time, we allow some variation in the location of the core-mantle boundary (CMB) as an effective way to account for the existence of a transition zone including the D" layer [36,37]. This is the lowermost ∼ 250 km thick portion of the mantle directly above the CMB and gaining a better picture of it can help to better understand the dynamics of the entire mantle, as well as how the core and mantle interact. Furthermore, requiring variations in the CMB radius to be consistent with the well-measured Earth's mass and moment of inertia allows us to change the densities of the outer core and mantle without violating such important constraints.
The paper is organized as follows. In Sect. 2 we present a simplified model of the Earth's structure with four density layers. In Sect. 3 we briefly review the formalism of neutrino oscillations in matter and present the formulas for the transition probabilities. In Sect. 4 we determine the number of μ-like neutrino events in a generic detector and the effects that changes in the composition and radius of the outer core have on it. Sect. 5 contains our results and final comments. There, we present a Monte Carlo simulation of the number of μ-like events and the application of this observable to test different composition models of the outer core. We included an appendix with details of the derivation of the transition amplitudes between flavor neutrinos in a medium with a symmetric density profile.

Model of the Earth's structure
We want to explore the possibility that neutrino oscillations in matter can provide information on the CMB and its near zone, contributing to a better understanding of some of the mechanism responsible for the phenomena described above. For this purpose, we consider the simplified spherical model for the Earth's interior shown in Fig. 1, which consists of four layers of constant densities delimited by concentric spheres of different radii: inner core, outer core, mantle, and crust. 1 The radii of the inner core R ic , the mantle R m , and the Earth R ⊕ are assumed to have the following fixed values: R ic = 1221.5 km, R m = 5600 km, and R ⊕ = 6371 km (Fig. 2). The primary information on the Earth's density as a function of the position comes from the total mass of the Earth M ⊕ = 5.9724 × 10 27 g and its mean moment of inertia about the polar axis I ⊕ = 8.025 × 10 44 gcm 2 [38]. In terms of the radial density distribution they are given by [39]: These important conditions alone requires that there be a concentration of mass towards the centre of the planet. Unlike previous works on the subject, we incorporate both constraints into our model. This is done by introducing a scheme for density variations, in which they are caused by changes in the radius of the outer core R oc around the value R oc = 3480 km. According to the above considerations we require that the sum of the masses and the sum of the inertial moments of the four spherical shells have to be equal to M ⊕ and I ⊕ , respectively: withR oc = R oc , where the numerical variable takes values close to one. For the densities of the inner core ρ ic and the crust ρ c we take the average values of the Preliminary Reference Earth Model (PREM) [3]: ρ ic = 12.916 g/cm 3 and ρ c = 3.671 g/cm 3 . On the other hand, the (average) densities of the mantle ρ m and the outer core ρ oc are variables that depends on . To determine them, let us rewrite Eq. (3) as with Y ic,oc,m = R ic,oc,m /R ⊕ and where ρ M and ρ I denote the average mass density of the Earth calculated by means of the total mass and the total moment of inertia, respectively: By solving the linear system given in Eq. (4) we express the densities of the outer core and the mantle as functions of with The coefficients of the polynomials F ( ), G ( ) and D( ) are evaluated in terms of the radii and the average densities given above. The resulting formulas for ρ oc ( ) and ρ m ( ) has been plotted in Fig. 3 in the interval of relevant to us. In particular, for = 1 (i.e.,R oc = R oc ) we get ρ oc = 11.526 g/cm 3 and ρ m = 4.742 g/cm 3 . Later, we will allow the radius and the composition of the outer core to vary simultaneously to examine their join effects on the flavor transformations and, in particular, possible cancellations among them that could leave the flow of neutrinos in the detector unaltered.

Atmospheric neutrino oscillations
The discovery of flavor oscillations with atmospheric and solar neutrinos, and their subsequent verification using reactor and accelerator neutrinos, has firmly established that neutrinos are massive mixed particles [40]. Except for some unconfirmed anomalies (see Refs. [41,42] for recent reviews) the current data set can be interpreted in terms of the minimal extension of the Standard Model, where the known flavor states |ν α (α = e, μ, τ ) are linear combinations of the states The coefficients U αi are elements of the unitary mixing matrix U that appears in the leptonic charged current, which can be parametrized as 2 where O i j are orthogonal matrices representing rotations by angles θ i j ∈ [0, π/2] in the respective planes and = diag(1, 1, e iδ ) is a diagonal matrix with δ ∈ [0, 2π ]. Values 2 The expression in Eq. (9) is for Dirac neutrinos. If they are Majorana particles U has to be multiplied by the right by another diagonal matrix The two additional physical phases do not play a role in neutrino oscillations and can be omitted in the analysis of the phenomenon [43,44]. of the phase δ = 0, π imply C P violation in the leptonic sector. Besides these quantities, the oscillations between three neutrinos are parametrized by two mass-squared differences: The sign of Δm 2 31 is unknown and will be considered in our study. The two possibilities correspond to different array of neutrino masses: the normal ordering (NO) and the inverted ordering (IO), depending on whether m 3 is greater or less than m 1,2 .
Neutrinos are produced and detected in flavor eigenstates. Due to slight mass differences, the phases of the mass eigenstate components of the original flavor state change at different rates and, as a consequence, the flavor content of the neutrino beam oscillates along the trajectory. The relevant quantity in connection with neutrino oscillations is the probability P ν α →ν β (L) that a neutrino born as a ν α can be found in a flavor ν β at a distance L from the source: where the probability amplitude U βα (L) = ν β |Û (L)|ν α is an element of the 3 × 3 unitary matrix U (L) representing the evolution operator in the flavor basis. When neutrinos propagate in a medium, their dispersion relations are modified due to the coherent forward scattering with electrons and nucleons. This effect can be incorporated by means of a refraction index different for ν e and ν μ,τ and, after subtracting an unobservable overall phase, results in a additional contribution to the Hamiltonian H (L) that governs the evolution of flavor amplitudes. Besides the vacuum part, it now contains a potential term V (L) = ∓v(L) diag(1, 0, 0), with v(L) = √ 2G F n e (L). Here, G F is the Fermi constant and the signs refer to neutrinos (minus) and antineutrinos (plus). The electron number density n e (L) is related to the mater density ρ(L) according to where m u = 931.494 MeV is the atomic mass unit and Z /A = λ r λ (Z /A) λ . The summation runs over all the elements present in the medium and (Z /A) λ denotes the ratio between the atomic number Z λ and the atomic mass A λ of the element that contributes a fraction r λ to the total mass. Under proper conditions, the pattern of neutrino oscillations in matter can be significantly modified compared to the oscillations in vacuum. New resonance enhancement effects appear, which are sensitive to the density and composition of the medium. Matter effects in long baseline oscillations inside the Earth are strongly dependent on sin 2 θ 13 , which drives the transitions ν e ↔ ν μ . Since the enhancement features depends on the sign of Δm 2 23 , the large measured value of θ 13 has opened the possibility to distinguish among the two mass ordering by using atmospheric neutrinos going deeply through the Earth. Atmospheric neutrinos are generated as decay products in hadronic showers resulting from collisions of cosmic rays with nuclei in the upper atmosphere. Due to the isotropic nature of the primary cosmic ray flux, these neutrinos are produced around the world, providing a continuous source of neutrinos spanning a very wide range of energies and travelled distances before detection. On this work, we concentrate on those with long baseline and energies in the range of 2-10 GeV, which also offer the possibility to determine the unknown neutrino mass ordering and also to conduct an oscillation tomographic study of the Earth's interior. According to Eq. (10), to calculate P ν α →ν β (L) we have to determine U (L) subject to the condition U (0) = I at L = 0. For three flavors the problem cannot be solved analytically in arbitrary density profiles, but a solution can be given for a constant density in terms of the energy eigenvalues in matter. Therefore, for a model like ours in which the planet consist of concentric spherical layers with constant densities, the evolution operator in the Earth U ⊕ (L) can be expressed as the ordered product of the respective operators in the consecutive layers traversed by neutrinos. For the calculations we take advantage of the fact that the density profile is symmetrical with respect to the midpoint of the entire neutrino trajectory. The details are presented in the Appendix. We need the probabilities for the transitions ν μ (ν μ ) → ν μ (ν μ ) and ν e (ν e ) → ν μ (ν μ ), which are those that are involved in the computation of the μ-like events produced in a detector by 'upward' atmospheric neutrinos. From the results derived in the Appendix, we have where L ⊕ = 2R ⊕ cos η denotes the length of the neutrino baseline within the Earth and The quantities in the right hand side of these equations with a tilde over them are elements of the symmetric matrix in Eq.
(A.5). The expressions for antineutrinos are similar but with the replacements δ → −δ and v → −v. The electron number density depends on both the matter density and Z /A, which is a function of the chemical and isotopic composition of the medium. Then, for neutrinos going through the deepest part of the Earth, P ν μ (ν e )→ν μ (L ⊕ ) become sensitive to changes in the composition of the outer core. In this work, additional effects are introduced by allowing variations in the position of the boundary between the outer core and the mantle, which lead to modifications in the matter densities of these layers.
The Earth's outer core is composed mostly of liquid iron, or iron alloyed with nickel, along with a significant fraction of lighter elements (up to 15%) necessary to meet the density required by seismic wave velocities. Except for hydrogen ((Z /A)) H = 1), for the components of this alloy Z /A ∼ = 0.5 because they have almost equal numbers of protons and neutrons. To appreciate the effect that changes in composition have on the flavor oscillations probabilities we will take a pure iron core as the standard model with which to compare models that contain a small percentage of hydrogen. In Figs. 5 and 4 we show the ν μ → ν μ and ν e → ν μ transition probabilities as a function of the neutrino energy, for different values of and weight percent of hydrogen.The values of the oscillations parameters are those given in Table 1. All figures were made both for normal and inverse ordering and θ 23 in the first octant. The values of the radii and the average densities of the layers are those given in Sect. 2, which were adapted from the PREM model. In the next section, we use these probabilities in the calculation of the μ-like events produced by atmospheric neutrinos after crossing the internal regions of the Earth.

Neutrino events and test of Earth's composition
To examine the feasibility of a tomographic study of the deeper layers of the Earth we look at the mark that variations of the density and composition leave on the flux of atmospheric neutrinos arriving to the detector after experience the effect of flavor oscillations in the terrestrial matter. We concentrate on the outer core, which is believe to play a fundamental role in the creation and behavior of the geomagnetic field, as well as in other important geophysical phenomena, such as plumes and hotspots. For our analysis, we consider a generic detector with a mass of water or ice containing a number of nucleons n N as target. A detector of this type, like for example IceCube, detect efficiently the Cherenkov radiation emitted along the trajectory of the μ ± produced by the interactions of ν μ andν μ with the nucleons in the instrumented volume or near it.
We implement a Monte Carlo simulation of neutrino propagation in order to have a data set that simulates real observa- tions, with the inevitable statistical fluctuations. To study the possibility of distinguishing between the different compositions and densities of the outer core, we applied the statistical methods for hypothesis testing considering a large number of pseudo experimental outcomes. To identify the most sensitive regions, we first do a scan dividing the cone under the detector into a series of angular and energy bins and calculate the number of μ-type events N μ within given angular and energy intervals. This can be expressed as where T is the detection time. The flux of atmospheric neutrinos and antineutrinos have been taken from Ref. [48], while the charged-current cross sections for the ν μ (ν μ )nucleon scattering, in the range of neutrino energies considered (1 − 10 GeV), are give approximately by [49] σ cc In Eq. (14) it is understood that oscillations probabilities are evaluated at L ⊕ . The dependence of N μ on the density and composition of the medium is incorporated through these probabilities, which are calculated from the expressions given in Sect. 3 (see also the Appendix).
To find the angular and energy intervals where N μ is more sensitive to changes in the upper core composition we intro- duce the quantity which gives the percentage difference between the number of events for the standard composition N 0 μ , without hydrogen in the outer core, and the number of events for a different composition N H μ . The values for the radii and the average densities of the layers are those given in Sect. 2, which were adapted from the PREM model.
In Fig. 6a and b we represent the level surfaces for ϒ in the (E, η) plane with a hydrogen proportion of 1wt%, both for the normal and inverted hierarchies. From the pictures, it is apparent that in all the cases the most sensitive region corresponds to energies around the 5 GeV and angles compatible with the shadow of the outer core. Deviations in the number of events are considerably more pronounced for NO, mainly because the effects of matter for IO occur for antineutrinos, whose charged cross section is approximately twice as small as for neutrinos. We also examine the change in the radii of the outer core and find that the sensitivity zone do not change appreciably.
Our observable is the number of μ-like events and in what follows we use it to test different hypothesis about the composition and density of the outer core. To this end, we consider events with 4 GeV < E < 6 GeV and 10 • < η < 30 • , and divide each of these intervals into 200 bins. In order to perform the Monte Carlo simulation we consider that each pseudo experiment consists of tossing, in each square bin of the grid, a number of Poisson distributed events with the mean value as given by Eq. (14). Thus, each pseudo experiment consist in 200 × 200 numbers corresponding to events, one for each bin. This sample of events are then distributed in angle and energy. We call it the true events and suppose that they are distributed according to the probability distribution function (pdf) f where the integral is over the energy and angle intervals of the bin (i, j).
To obtain a realistic distribution of events we must allow for a limited resolution of the detector, both in energy and angle. The net effect is the redistribution of the true events ("smearing") in the grid bins, which is implemented by folding the true distribution with a resolution function S (E o , η o |E, η). We assume a Gaussian smearing and write with the detector characterized by the angular and energy resolutions Δη(E) = α η / √ E/GeV and ΔE(E) = α E E, respectively. For the values of the parameters α η and α E we consider different situations according to the discussion in Ref. [32]. To keep our analysis as simple as possible we assume a detection efficiency of 100%.
When convoluted with the true events the kernel in Eq. (18) gives us what we call the observed events. That is, S(E o , η o |E, η) represents the conditional pdf for the measured values (E o , η o ) given that the true values were (E, η) and, since the event is observed somewhere, it is normalized such that In terms of the resolution function, the number of observed events O i ex p m,n in the bin (m, n) for the i ex p -th experiment is given by As an illustration, in Fig. 7  for an alternative Earth with a given fraction of hydrogen in the outer core and a ratio =R oc /R oc for each of the n ex p pseudo experiment. Thus, for each bin of the observed events we have a sample of the events and determinate the corresponding mean valuē where α label the two-dimensional bins (m, n). For Poisson distributed events and for the likelihood function L we construct the negative log-likelihood ratio function as The same number of bines has been used to make the figures, but in the calculations a much small number of bins is used for the observed events According to Wilks' theorem [50] the χ 2 λ distribution can be approximated by the χ 2 distribution, and from it, the goodness of fit can be established. The statistical significance of the χ 2 test is, as usual, given by the p-value: where n dof is the number of degrees of freedom and f χ (w, n dof ) is the chi-square distribution. In our case, n dof = 9 × 5 − 2 = 43. In this way, we can study the discrepancy level between the standard Earth and different hypothesis about the composition and density of the outer core. This is done in the next section where we present our results and final comments.

Results and final comments
As discussed in Sect. 2, we allow some variation in the location of the CMB that, to meet the constraints on the total mass and moment of inertia of the Earth, results in the dependence of the outer core and mantle densities on (Fig. 3).
As the quantities to be fitted, in what follows we take the fraction of hydrogen in the outer core and the relative change of its density Δρ oc /ρ oc , where Δρ oc = ρ oc ( ) − ρ oc and ρ oc = 11.526 g/cm 3 is the density for = 1. In order to evaluate the effect that changes on both of these quantities have on our observable, we construct the regions in the (Δ ρ oc /ρ oc , H) plane where the statistical significance, given by the p-value, for the discrepancy between the standard and the modified Earth is less than 1σ . In Fig. 8 we show these regions for the oscillation parameters given in Table 1, for 10 years operation of a 10 and 100 Mton detector, with resolution parameters α η = α E = 0.1. As can be seen, the regions for NO are more restricted than for IO. In both cases, a negative correlation is observed between Δρ oc /ρ oc and H, which indicates that an increase in the percentage of hydrogen in the outer core can be compensated by a decrease in its density.
We also examine the effect that variations of the detector resolutions have on the bounds at 1 σ confidence level. Fig. 9 left (right) panel shows the bound for Δρ oc /ρ oc as a function of the resolution parameters, when no hydrogen is present in the outer-core (H = 0), for normal (inverted) hierarchy. Likewise, the panels in Fig.10 show the corresponding bounds on the hydrogen fraction in the outer-core, when the matter density is that given by the PREM (Δρ oc /ρ oc = 0).
In summary, we have studied the possibility of conducting an oscillation tomography of the Earth based on the matter effects on the flavor oscillations of atmospheric neutrinos propagating through the depths of the planet. Using the μlike events in a generic large Cherenkov detector as physical observables and making a Monte Carlo simulation of the energy and azimuthal angle distribution of these event, we tested possible variants with respect to a reference geophysical model with the average densities as given by PREM and no hydrogen in the outer core. Unlike previous studies, the procedure we followed in this work allowed us to simultaneously vary the composition and density of the outer core. When one of these quantities was fixed in the value of the reference geophysical model, our results, shown in Figs. 9 and 10, are compatible with those obtained by others authors in an uncorrelated way [33,35]. The remote located interface between the rocky mantle and iron core is physically the most significant in the Earth's interior. The phenomena that happens in the region around this interface play an important role in several process that impact he whole Earth and its complete understanding requires a combined effort of different geophysical disciplines. This task could greatly benefit from the use of phenomena of the type described here, normally studied in particle physics, and the arrival of new and improved neutrino telescopes, such as KM3NeT, ORCA, PINGO, and Hyper-Kamiokande [51][52][53][54]. Even though our results were derived considering a simplified PREM model with four layer of constant densities, and without taking into account the systematics related to the experimental device, they exemplify the potentialities of neutrino oscillations as a novel techniques to explore the interior of the Earth.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data sharing not applicable to this article as no datasets were generated or analysed during the current study.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Flavor neutrino evolution in the Earth
In this appendix we describe in detail the calculation of the flavor transition amplitudes for neutrinos produced in the atmosphere and detected after they go through the Earth. Let denote by L D the distance between the detector and the source, then for the U (L D ) = U ⊕ (L ⊕ ) U A (L A ), where U ⊕ and U A are the matrix representations of the evolution operators in the Earth and the atmosphere, respectively. We have denoted by L ⊕ and L A the lengths of the neutrino path en each of these media, with L D = L ⊕ + L A . The flavor evolution along the atmosphere is well described assuming that neutrinos propagate in vacuum. Thus, in the flavor basis we have where U is the mixing matrix in Eq. (9) and Here, E i = |p| 2 + m 2 i are the energies of the neutrino mass eigenstates described as plane waves with the same momentum p. Discarding a global phase, for relativistic neutrinos we have , and E ∼ = |p|. Typically, L A ∼ 10 km and, for neutrinos with energies E 1 GeV, both ϕ 21 and ϕ 31 are small quantities. Therefore, D is almost equal to the identity matrix and U A (L A ) UU † = I . Consequently, the probabilities P ν α →ν β in the detection volume can be calculated by means of the matrix elements of the evolution operator within the Earth. In what follows we will focus on deriving them. The 3 × 3 matrix U ⊕ (L) is the solution of the differential equation where L is the distance traveled by neutrinos along the Earth and is In one-dimensional models of the type we use here, which represent the average properties of the Earth, the density profile seen by neutrinos passing through the terrestrial matter is symmetrical with respect to the midpoint of their entire path. Hence, the effective potential energy satisfies the equality To compute the matrix in Eq. (A.5) we keep in mind that, as described in Sect. (2), we model the Earth as an onion like sphere made of four layers with different constant densities. Consequently, the full evolution operator can be expressed in terms of the product of the evolution operators in each of the successive layers through which the neutrinos pass on their way to the detector. Thus, for neutrinos traversing the inner core we havẽ It should be noted that, in general, the matrices of the different layers do not commute between them and the factors must go in the order prescribed. Here, L is equal to half of the whole distance that neutrinos travel in layer , so L =L. These distances can be easily determined as functions of the nadir angle η by applying elemental trigonometry (see Fig. 1): with 0 ≤ η 11 • . If neutrinos goes trough the outer core but not the inner core, then = oc, c, m and 11 • η 33 • . In this case, L ic = 0 and L oc = R ⊕ Y 2 oc − sin 2 η, while the formulas for L m and L c remain unchanged.
Since the layers are assumed to have a constant electron number density n e , for each factor in Eq. (A.6) we havẽ U (L ) = exp(−iH L ), where inH the effective potential takes the fixed value v = √ 2 G F n e . The exponential of a constant matrix can be evaluated by means of the Putzer's algorithm [55]. For our 3 × 3 matrices, discarding in material phases, the results are [56]  with ω 21 = E 2 −E 1 and similarly for ω 31 and ω 32 . The quantities E k (k = 1, 2, 3) are the energy eigenvalues in matter and are given by the real roots of the characteristic (cubic) equation They are the eigenvalues ofH (which coincide with those of H ), do not depend on θ 23 or δ, are positive and unequal, and can be expressed as In this way, we have completed the determination of all the ingredients necessary to calculate each factor in Eq. (A.6) using formula (A.8). OnceŨ ⊕ (L) has been obtained by multiplying the required factors, it is immediate to findŨ ⊕ (L ⊕ ) by means of Eq. (A.5). The resulting matrix, when substituted into Eq. (A.4), allows us to express the elements of U ⊕ (L ⊕ ) in terms of those of the matrix in (A.5).