A Event numbers

Exploring the Earth matter effect with atmospheric neutrinos in ice



Exploring the Earth matter effect with

atmospheric neutrinos in ice


Sanjib Kumar Agarwalla1, Tracey Li2, Olga Mena3 and Sergio Palomares-Ruiz4

567 Instituto de Física Corpuscular, CSIC-Universitat de València,

Apartado de Correos 22085, E-46071 Valencia, Spain

8 Centro de Física Teórica de Partículas (CFTP), Instituto Superior Técnico,

Universidade Técnica de Lisboa, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal


We study the possibility to perform neutrino oscillation tomography and to determine the neutrino mass hierarchy in kilometer-scale ice Čerenkov detectors by means of the -driven matter effects which occur during the propagation of atmospheric neutrinos deep through the Earth. We consider the ongoing IceCube/DeepCore neutrino observatory and future planned extensions, such as the PINGU detector, which has a lower energy threshold. Our simulations include the impact of marginalization over the neutrino oscillation parameters and a fully correlated systematic uncertainty on the total number of events. For the current best-fit value of the mixing angle , the DeepCore detector, due to its relatively high-energy threshold, could only be sensitive to fluctuations on the normalization of the Earth’s density of at  CL after 10 years in the case of a true normal hierarchy. For the two PINGU configurations we consider, overall density fluctuations of () could be measured at the  CL after 10 years, also in the case of a normal mass hierarchy. We also compare the prospects to determine the neutrino mass hierarchy in these three configurations and find that this could be achieved at the  CL, for both hierarchies, after 5 years in DeepCore and about 1 year in PINGU. This clearly shows the importance of lowering the energy threshold below 10 GeV so that detectors are fully sensitive to the resonant matter effects.

1 Introduction

Cosmic ray interactions in the atmosphere are a natural source of neutrinos with baselines that span three orders of magnitude and energies in the range from MeV to well above TeV. This implies that different scales and effects could be accessible by studying atmospheric neutrinos. In fact, the most important result of the latest years in neutrino physics was the discovery in 1998 of neutrino oscillations at the Super-Kamiokande detector by using atmospheric neutrino data in the GeV energy range [1].

Among the potential effects that could be explored are the resonant effects in the propagation of GeV atmospheric neutrinos through the Earth that mainly affect the subleading and (or and ) transitions. When atmospheric neutrinos deeply cross the Earth’s mantle, the Mikheyev-Smirnov-Wolfenstein (MSW) resonance [2, 3] could be in action [4, 5, 6]. On the other hand, when they also traverse the core, resonant effects, different from the MSW resonance, could show up [7, 8, 9, 10, 11, 12, 13, 14, 15]. These effects are well known and their physics potential in iron-magnetized calorimeter, water-Čerenkov and liquid argon detectors has been thoroughly studied in the literature [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37].

It is very well known that the propagation of atmospheric neutrinos with energies in the range of a few GeV deeply through the Earth is sensitive to the currently unknown neutrino mass hierarchy [4, 16, 18, 21] (or neutrino mass ordering), i.e., whether the third neutrino mass eigenstate is the lightest or the heaviest mass eigenstate. In the latter case, the hierarchy is called a normal hierarchy (NH) and in the former case, inverted hierarchy (IH). However, the ability to measure the neutrino mass hierarchy exploiting the resonant oscillation phenomena of atmospheric neutrinos inside the Earth crucially relies on the value of the mixing angle . Recent measurements of this angle from the Daya Bay [38], RENO [39] and Double Chooz [40] reactor experiments, in addition to the long-baseline T2K experiment [41], seem to indicate that  [42, 43, 44]. This large value of opens up the possibility of using the atmospheric neutrino fluxes not only to determine the neutrino mass hierarchy, but also to study features of the Earth’s matter density profile, by exploiting the -driven matter effects.

There are three different ways that have been proposed to infer some information about the internal structure of the Earth by exploiting the different effects of neutrino propagation in matter: neutrino absorption tomography, neutrino oscillation tomography and neutrino diffraction. The idea of using the absorption of very high energy neutrinos to explore the Earth’s interior dates back to 1974 [45] and is analogous to X–ray tomography, but instead exploits the attenuation of the neutrino flux for energies  TeV [46]. There are numerous studies which have considered neutrinos of different origins, such as man-made neutrinos [45, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56], extraterrestrial neutrinos9 [51, 57, 58, 59, 60] and atmospheric neutrinos [61, 62, 63, 64]. On the other hand, neutrino oscillation tomography relies on the matter effects in neutrino oscillations and it has been considered by studying man-made beams [65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75], solar [76, 77] and supernova neutrinos [78, 77]. The third possibility is based on studying the diffraction pattern of crystalline matter in the interior of the Earth caused by coherent neutrino scattering, but it is technologically unfeasible [79].

In the case of atmospheric neutrinos, only neutrino absorption tomography has been considered within the context of kilometer-scale detectors [61, 62, 63, 64]. However, resonant effects are very sensitive to the matter density that neutrinos traverse, so one could also think of doing neutrino oscillation tomography with atmospheric neutrinos with energies in the GeV range10. This is one of the main goals of this work and represents a completely different technique to determine the Earth’s density distribution to those used in geophysics, mainly based on the analysis of seismic waves. Although, in principle, geophysics can obtain more precise results, its inferences are based on numerous assumptions because the velocities of seismic waves not only depend on the density, but also on the composition, temperature, pressure and elastic properties of the medium. Thus, the results from atmospheric neutrino tomography would represent an independent and complementary assessment of the Earth’s internal structure.

Neutrino telescopes such as AMANDA [80], IceCube [81] and Antares [82] have already collected a large number of atmospheric neutrino events [83, 84, 85, 86, 87, 88, 89], even in spite of their high energy threshold ( 100 GeV) and the steeply falling atmospheric neutrino spectrum (). Whereas this is a very high threshold for studying neutrino oscillations with atmospheric neutrinos, a new generation of neutrino telescopes with lower energy thresholds ( 10 GeV), such as the DeepCore extension of the Icecube detector, is currently taking data successfully [90, 91, 92] and further natural extensions of this are being planned. The Precision IceCube Next Generation Upgrade (PINGU) has been proposed in order to reduce the detection threshold down to a few GeV [93]. Although reaching those energy thresholds in these multi-Mton scale neutrino telescopes is very challenging, if successful, the atmospheric neutrino events detected at them would offer a great opportunity for detailed oscillation studies, including the determination of the neutrino mass hierarchy [94, 95], tau neutrino appearance searches [96], and precise measurements of the atmospheric neutrino oscillation parameters [97]. Given the large value of measured by the different reactor [38, 39, 40] and the T2K [41] experiments, the prospects look promising. It therefore seems timely to assess their ability to infer some information about the Earth’s matter density profile and how this could also affect the determination of the neutrino mass hierarchy from the atmospheric neutrino events. In this work we study the future prospects for DeepCore and PINGU by considering the experimental capabilities of the two detectors and present the expected sensitivities to both of these observables.

The structure of the paper is as follows. We start in Sec. 2 by revisiting the atmospheric neutrino fluxes, the Earth’s internal structure and the matter effects which affect the propagation of GeV neutrinos inside the Earth. Sec. 3 contains a description of the detectors considered here, as well as of the analysis methods used to calculate the sensitivities to the Earth’s matter density and neutrino mass hierarchy, which are presented in Sec. 4. Finally, in Sec. 5 we draw our conclusions.

2 Atmospheric neutrino fluxes and oscillations in the Earth

2.1 Atmospheric neutrino flux

When cosmic rays traverse the Earth’s atmosphere, hadronic showers are produced by their interactions with the particles that constitute it. Depending on the energy of the primary cosmic ray, different secondaries can be produced. Below  100 GeV, the neutrino flux is dominated by the pion decay chain, whereas above these energies, kaon decays dominate neutrino production. The decay chain is the following:


Thus, when the conditions for the decay of all these particles are satisfied, we expect flavor ratios for the atmospheric neutrino fluxes of and . This flux is usually referred to as the conventional atmospheric neutrino flux and the energies of these atmospheric neutrinos range from a few MeV up to hundreds of TeV.

The gross features of the spectrum are easy to understand (see e.g., Refs. [98, 99]). For energies , with being the zenith angle, all pions decay and the neutrino spectrum has approximately the same power as the primary cosmic-ray spectrum, which is . Above these energies, the expected power-law neutrino spectrum (from pion decays) is asymptotically one power steeper. This has to do with the extra factor of in the ratio of the pion decay length to interaction length. The fact that the flux is higher in the horizontal direction (it is proportional to ) is explained by the longer decay path near the horizon: a pion traveling through the atmosphere horizontally () has a higher probability of decaying in the atmosphere than a pion traversing the atmosphere vertically () and therefore the flux is higher for the horizontal component. A similar behavior occurs for kaons, but shifted to energies about one order of magnitude higher. On the other hand, the muon spectrum is very similar to the parent meson spectrum, but for energies above  1 GeV, many muons reach the ground before decaying (at  100 GeV, only  15% of atmospheric muon neutrinos come from muon decay), and hence the neutrino flavor ratio grows with increasing energy. Similarly to the angular dependence of neutrinos from pion and kaon decays, this results in the steepening of the electron neutrino spectrum occurring at lower energies for vertical directions and at higher energies for propagation close to the horizon. Existing analytical approximations to the atmospheric neutrino fluxes versus energy and zenith angle at high energies also allow a good qualitative understanding of all these features [98, 99].

The huge range of energies and baselines provided by the atmospheric neutrinos opens up the possibility of exploring many different and exciting physics topics, from the measurement of neutrino mixing parameters [100, 101] to more exotic ones [102, 103, 104, 105, 106, 107, 108, 109, 110, 84, 111, 112, 113, 114]. At very high energies, above 10 TeV, neutrino interaction cross sections become high enough for neutrinos going through the Earth to start becoming attenuated. This effect is sensitive to the neutrino interaction cross sections and to the density profile of the Earth and thus, a measurement of the atmospheric neutrino-induced event rate at these energies could provide information about the Earth’s internal structure [61, 62, 63, 64]. The intermediate-energy region, between  50 GeV and  1 TeV, provides information about the atmospheric neutrino flux normalization [84]. Finally, in the low-energy region which we exploit here (below about 50 GeV), neutrino oscillation effects could be significant [4, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 21], in particular due to the large value of the mixing angle recently measured [38, 39, 40, 41]. Resonant matter effects in the propagation of neutrinos inside the Earth could be very important for  GeV and could strongly enhance or decrease the oscillation probabilities, for neutrinos in the case of normal neutrino mass hierarchy and for antineutrinos if the neutrino mass hierarchy is inverted.

The atmospheric neutrino fluxes we use in this work are those from Ref. [115] (see also Refs. [116, 117, 118] for other atmospheric neutrino flux calculations11). The absolute electron and muon atmospheric (anti)neutrino fluxes are uncertain at the level of 10% - 20% in the energy region of interest here, which is related to our ignorance regarding hadron production models. The uncertainties quoted above are largely canceled for the muon neutrino-antineutrino flavor ratio as well as for the up-down ratio, the former expected to be of order above 1 GeV [119]. Accelerator data from the HARP, MIPP and NA61 experiments on particle multi-production are also expected to improve the predictions of the absolute atmospheric neutrino fluxes [120]. However, as we will see when discussing the effects of correlated systematic uncertainties in the analysis, the uncertainties related to the normalization of the atmospheric fluxes are alleviated when higher energy bins, as well as different angular bins, in which oscillation effects are negligible, are also used.

2.2 The Earth’s internal structure

Most of the knowledge we have about the internal structure of the Earth and the physical properties of its different layers comes from geophysics and, in particular, from the data we obtain from seismic waves. Other information, from geomagnetic and geodynamical data, solid state theory and high temperature-pressure experimental results is also used. It turns out that a reliable estimate of the density of the Earth is essential to solve a number of important problems in geophysics, such as the dynamics of the core and mantle, the gravity field, the mechanism of the geomagnetic dynamo, the bulk composition of the Earth, etc. (see, e.g., Ref. [121] and references therein).

The Earth is conventionally divided into three main (approximately) spherical concentric shells: the crust, the mantle and the core, each of them further divided into subshells, with different properties (see, e.g., Refs. [122, 123]). Whereas the crust only represents about 0.4% of the Earth’s mass, the mantle and the core constitute about 68% and 32%, respectively. On the other hand, the core has a radius almost half that of the Earth and it is about twice as dense as the mantle.

The divisions of the Earth’s interior are based on the reflection and refraction of compressional and shear waves, i.e., P- and S-waves, respectively. The crust is the brittle outer shell, with a large proportion of incompatible elements. It is composed of three different crustal types: continental (thick and composed primarily of granite), transitional (defined by the islands, island arcs and continental margins) and oceanic (thinner and composed primarily of basalt). The crust is separated from the mantle by a boundary discovered by Mohorovicic in 1909. The mantle is composed mainly of silicates and oxides, being a poor conductor of electricity and heat and with very high temperature-dependent viscosity. Its temperature increases with depth creating a geothermal gradient which is responsible for the different rock behaviors that determine the mantle subdivisions. Whereas in the upper mantle rocks are cool and brittle, in the lower mantle they are hot and soft, but not molten. In 1891 E. Wiechert inferred the existence of a completely different region beneath the mantle from the mean density and moment of inertia of the Earth, but it was not until 1914 that B. Gutenberg made the first accurate determination of the location of the core-mantle boundary [124] from the R. D. Oldham measurement in 1906 of a sharp decrease in the velocities of P- and S-waves deep in the Earth. Soon after this, the first models of the radial structure of elastic velocities started to be developed. This boundary presents a density contrast of a factor of  2 and a viscosity contrast of  20-24 orders of magnitude, which is what ultimately causes the strongly differing time and length scales for convection motions in these two layers. The Earth’s core, the Earth’s source of internal heat, is thought to be composed mainly of an alloy of iron and nickel, with a small admixture of lighter elements, which is a good conductor of electricity and heat. Its radius is approximately 2900 km and it is divided into the outer (liquid) and the inner (solid) cores.

Although the Earth is not exactly spherically symmetric and is irregular, the use of mean models, taking the Earth’s physical parameters as symmetric, is very useful as a reference framework for different studies [125]. During many years, the most widely used radially symmetric model has been the Preliminary Reference Earth Model (PREM) [126], which included anisotropy in the upper mantle, and it was a good fit to free oscillation center frequency measurements, surface-wave dispersion observations, travel-time data for body-waves and some astronomical data such as the Earth’s radius, mass and moment of inertia. However, there was no discussion about the sources of error or covariances. A more recent 1-D model, AK135, which combines travel times and free oscillations, is also available [127, 128]. Although there are observed deviations at the few per cent level that cannot be explained by 1-D models, to first order, the structure of the Earth is determined by its radial properties and these reference models are used as a starting point for more refined studies. In Fig. 1 we show the matter density profile of the Earth (and fluctuations of ) as a function of the distance to the center according to the PREM [126]. Up to a radius of  km, the density increases gradually from the crust to the mantle. Then, in the core–mantle boundary, there is a very sharp transition, followed by a gradually increasing density until the center of the Earth.

Figure 1: The Earth’s density profile according to the PREM [126] (red) and with a overall density fluctuation (green dash-dotted and blue dashed lines), as a function of the distance from the center of the Earth.

For the last half a century the problem of determining the density distribution of the Earth has mainly consisted of perturbation inversion using different seismic data. This is a very demanding problem, which in many cases is non-linear. Moreover, most studies of the Earth’s radial structure make inferences about different properties given a specific density model or based on empirical relations between seismic waves velocities and density. Nevertheless, there are significant trade-offs with other crucial parameters because wave velocities also depend on composition, temperature, pressure and elastic properties, which necessarily imply some modeling and uncertainties in the density estimate. Recently, 3-D models, using 1-D models as a reference, have been developed, although some studies concluded that few robust density features could be constrained with existing data [129, 130, 131, 132] (see, however, Refs. [133, 134, 135, 136, 137, 138]). All in all, density values, averaged over  100 km are thought to be known at the level of a few per cent at all depths, whereas density gradients are less known. On the other hand, linear integral constraints on the spherically symmetric density distribution are known at the level of from estimates of the Earth’s mass and moment of inertia [139, 140, 141] (see also Ref. [138] for a recent analysis of the spectra of the lowest frequency seismic mode with a precision of ). This implies that global variations of the density are constrained to be within those uncertainties. Obviously then, overall density variations at the level of a few per cent (as those shown in Fig. 1) are already excluded. Nevertheless, we will not impose these linear integral constraints and evaluate, using the PREM as our reference model, how sensitive atmospheric neutrino oscillation tomography is to changes in the density.

2.3 Oscillation probabilities

For neutrino energies in the range of a few GeV, the transition probabilities () and () of atmospheric neutrinos in their propagation through the Earth are relevant if genuine 3-flavor neutrino mixing takes place, i.e., if the mixing angle is different from zero [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Moreover, in this energy range and for these baselines ( km), CP-violation effects are very small and can be safely neglected. Likewise, effects due to the 1-2 sector are also subdominant and, as a first approximation, can also be neglected. In this context, the calculation of the transition probabilities effectively reduces to a 2-neutrino problem, with and playing the role of the relevant 2-neutrino oscillation parameters. Within these approximations, the 3-neutrino oscillation probabilities of interest for atmospheric having energy and crossing the Earth along a trajectory characterized by a zenith angle12 , have the following form [7, 11, 12, 13] (see also Refs. [8, 10, 14, 15]):


where is the 2-neutrino probability describing transitions, where , and and are the phase and the 2-neutrino transition probability amplitude. For antineutrinos the oscillation probabilities are analogous to those for neutrinos: they can be obtained formally from Eqs. (2) - (6) by changing the sign of the matter potential (or equivalently, by ). It is interesting to note that, within the approximation , the probabilities for neutrinos and NH (IH) are the same as those for antineutrinos and IH (NH).

Therefore, the magnitude of the matter effects depends on the 2-neutrino oscillation probability . In case of oscillations in vacuum, , so this probability is small. However, matter effects can strongly enhance and thereby greatly modify the 3-neutrino probabilities. On the other hand, if , then and , and hence matter effects are absent. If this were the situation, these probabilities would get reduced to the case of 2-neutrino oscillations, so for NH and IH they would be equal and identical to the case of vacuum oscillations.

Figure 2: Oscillation probabilities. Upper panels: as a function of , for  GeV (left panel) and  GeV (right panel), for NH (red regions) and IH (green regions limited by blue lines). Lower panels: Same but for . The widths of the bands correspond to varying the matter density by . Similar results are obtained for the case of antineutrinos by exchanging the curves for NH (IH) by those for IH (NH). We have used the current best-fit values for the oscillation parameters [42].

In Fig. 2 we show the (upper panels) and (lower panels) transition probabilities as a function of for two different energies ( GeV in left panels and  GeV in right panels) for NH (red regions) and IH (green regions limited by blue lines). The bands correspond to how the probabilities change if the density (according to the PREM) varies up to . We can see that matter effects, that are very sensitive to the value of the density, tend to greatly enhance the transitions and reduce the survival probability, and also shift the positions of the maxima and minima with respect to the case of negligible matter effects. Moreover, a change in the matter density profile would shift the location of both the resonance energy and the baseline at which the resonant behavior is maximized, modifying the number of events in the different angular and energy bins. For propagation in the mantle only, the resonance is, to a good approximation, the MSW resonance. However, in the case of propagation through the core (), non-trivial resonant effects show up, although for the energies considered in Fig. 2 they are not maximal. Notice that matter effects are larger for  GeV and more baselines are affected than for  GeV, the reason being that resonances are closer to  GeV (see below).

As we discussed above, the Earth is composed, in addition to the crust, of two major density structures that constitute almost all the mass of the planet: the mantle and the core. Although in our calculations we compute the full 3-neutrino evolution and use the detailed PREM for the density distribution13, as a first approximation and in order to gain insight about the physics of neutrino propagation in the Earth, one can consider a simplified model for the density profile, divided into two shells with different densities, and , with the core-mantle boundary at a radius of  2900 km. All the interesting features of the atmospheric neutrino oscillations in the Earth can be understood quite accurately within this framework [7, 8, 9, 10, 11, 12, 13, 14, 15]. There are analytical solutions for the transition probabilities for neutrinos crossing the Earth [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18], but they reduce to the case of neutrino propagation in a medium of constant density for trajectories such that , i.e., for neutrinos which propagate only in the mantle and could experience MSW resonant effects [4]. In this case, Eqs. (2) - (5) read14




and the matter parameter is , with the Fermi constant and the electron number density. The expressions for antineutrinos can be obtained by replacing . Thus, the resonant behavior, when maximal mixing occurs, happens for the case of NH (IH) in the neutrino (antineutrino) channel at the resonant energy,


On the other hand, the first oscillation maximum of the oscillation term happens when the oscillation phase is . Hence, the baseline at which both the conditions for the resonance and the first oscillation maximum are satisfied, is [4]

Figure 3: Line average of the Earth’s density and resonance energy. Upper panel: Line average density of the Earth for a given trajectory with total baseline . Lower panel: Resonance energy for propagation of neutrinos in the Earth assuming the density is constant and equal to the line average density. We have assumed ,  eV and the density distribution according to the PREM (red lines) and with a overall density fluctuation (green dash-dotted and blue dashed lines).

For ,  km and  GeV for densities approximately corresponding to the average value of the mantle. However, strong matter effects are also present around the resonant energy even if the baseline does not correspond exactly to the first oscillation maximum, i.e., even for baselines  [4]. For illustration purposes, we show in Fig. 3 the line average density for a given trajectory in the Earth (upper panel) and the corresponding resonance energy in the case of propagation in a medium of constant density with that value (lower panel).

Using Eqs. (2) - (6) it is straightforward to see that the fluxes of atmospheric of energy which reach the detector after crossing the Earth along a given trajectory, , are given by the following expressions in the case of the 3-neutrino oscillations under discussion [9, 10, 11, 18]:


where is the flux at the production point in the atmosphere and


For neutrino trajectories with , i.e.,  km, the predicted ratio of to atmospheric neutrino fluxes is for sub-GeV (anti)neutrinos while for multi-GeV (anti)neutrinos [115, 116, 117, 118]. Thus, even in the case of resonant matter effects taking place, the effects of -driven transitions are suppressed for the sub-GeV sample due to the term , especially if (the current confidence level (CL) allowed range is  [42]). On the other hand, for multi-GeV neutrinos, , so actually this factor can amplify the matter effects. Nevertheless, neutrino telescopes are Čerenkov detectors, so they have no charge-identification capabilities and therefore cannot distinguish neutrinos from antineutrinos. As the resonances take place for NH for neutrinos and IH for antineutrinos, the -driven matter effects would consequently be smeared and the sensitivity to them reduced as compared to the case when measurements of the neutrino-induced rates and antineutrino-induced rates can be performed separately15. However, notice that the fact that neutrino and antineutrino cross sections are different could allow us to distinguish NH from IH. Unfortunately, neutrino telescopes have very poor angular resolution for cascades, so the electron neutrino-induced event rates, where these resonant effects are expected to be larger (see Fig. 2), would give little information. Therefore, muon tracks are the main signal in this study and hence, the muon neutrino-induced event rates.

Following the discussion above, from Eqs. (2) - (5), Eqs. (17) - (18) and the analogous equations for antineutrinos, we see that -driven matter effects increase the (and ) transitions with , leading to an increase of the electron neutrino-induced event rate and a decrease of the muon neutrino-induced event rate at the detector. The rates are larger for multi-GeV events than for sub-GeV events.

3 Analysis

3.1 Detectors set-up

The IceCube/DeepCore detector [90] is a densely instrumented region located at the bottom center of the IceCube detector at a depth of between 2100 m and 2450 m, in such a way that it avoids a horizontal layer with a high content of dust and therefore poor optical properties. The detector consists of six additional strings instrumented with phototubes of higher quantum efficiency with respect to IceCube. Thus, in general, DeepCore has multiple advantages when compared to Icecube: at the depth at which DeepCore is located, ice is on average twice as clear as the ice above, allowing for a lower neutrino energy threshold. Additionally, the larger amount of photosensors, separated by 7 m instead of 17 m for IceCube, and the higher quantum efficiency, lead to a significant gain in sensitivity of up to a factor of 6, especially for low energy neutrinos. The rest of the IceCube detector, along with a horizontal region with additional instrumentation at a depth of 1800 m, could be used as an active veto for downgoing atmospheric muons, allowing the study also of downgoing atmospheric neutrinos.

At the energies of interest here, the angular reconstruction capabilities of DeepCore for showers is not optimal. Therefore, for the present study, we only consider upgoing muon-like events with the effective volume for the 86-string configuration (IC86) at trigger (SMT3) and online filter level in the 10-15 GeV energy region as shown in Tab. 1 [144]. We follow a very conservative approach by only considering this single energy bin ( GeV) and discarding lower energy bins which in principle have a lower effective volume and worse detection capabilities. The IceCube Collaboration aims to achieve a signal efficiency of for contained and partially contained events [90].

There is also a recent proposal to further upgrade the IceCube detector, with the planned PINGU [93]. It would consist of 20 additional strings within the DeepCore volume to lower the neutrino energy threshold down to  GeV energies and the detector would also benefit from the DeepCore strings.

In this work, we consider two possible PINGU scenarios: a conservative scenario which we refer to as PINGU-0 with 5 GeV energy threshold and with two 5 GeV energy bins ( GeV and  GeV), and a less conservative scenario which we call PINGU-I, with four 2.5 GeV energy bins ( GeV,  GeV,  GeV and  GeV). The effective masses for these configurations (Tab. 1) are computed assuming a trigger setting of 3 digital optical modules hit in 2.5  that are in local coincidence, a containment criterium which is implemented by a cut on the -position that matches the DeepCore fiducial volume ( m) and a tight radius cut from string 36 ( m) which is the center of DeepCore/PINGU [144].

Detector Energy bins (GeV) (Mton)
8.9 - 10.875 3.91
DeepCore 10.875 - 12.85 4.46
12.85 - 14.825 5.69
14.825 - 16.8 6.53
5 - 6 4.21
6 - 7 4.87
7 - 8 4.91
8 - 9 5.79
PINGU 9 - 10 7.54
10 - 11 6.73
11 - 12 6.77
12 - 13 7.58
13 - 14 7.86
14 - 15 8.77
Table 1: Detectors set-up. Effective mass for DeepCore and PINGU for the relevant energy bins which are considered in the current study (see text for details) [144]. In all cases, we add a post-trigger detection efficiency of (not included in these numbers). For DeepCore we consider a single energy bin,  GeV, and for PINGU we consider events in the interval  GeV but we use two different sets of energy bins: for PINGU-0,  GeV and  GeV; and for PINGU-I,  GeV,  GeV,  GeV and  GeV. The size of the angular bins is always in and we consider events in the interval .

Throughout this work we assume that the determination of the neutrino energy and direction could be achieved with some degree of accuracy. This implies that in addition to the characterization of the muon tracks induced by the atmospheric muon neutrino interactions in the detector, the accompanying hadronic cascades must also be detected. This is a challenging task, so in this respect, the results presented here should be taken as optimistic. Nevertheless, we have assumed very wide energy bins, so that a neutrino energy resolution at the level of is not expected to significantly affect our results. On the other hand, if the angular resolution is not better than the typical root mean square value of the scattering angle between the incoming neutrino and the produced muon, , some smearing is expected and hence, the sensitivity would be slightly worse. Let us also note that the effective masses used in this work represent preliminary MonteCarlo results based on the outcome of ongoing simulations by the IceCube Collaboration and these numbers might be further refined in the near future.

In summary, for the DeepCore detector simulations we assume a single energy bin of  GeV. For the PINGU-0 scenario, we assume two energy bins:  GeV and  GeV. For PINGU-I, we assume four energy bins:  GeV,  GeV,  GeV and  GeV. The effective masses we consider are shown in Tab. 1. In all configurations we assume that the post-trigger efficiency of the detector, for all energy bins, is and consider an angular binning in of width 0.1 for .

3.2 Number of muon-like events

Čerenkov detectors do not have charge-identification capabilities, so neutrinos cannot be distinguished from antineutrinos. Hence, these two types of events must be added together in the analysis. The number of contained atmospheric muon-like events (neutrino- plus antineutrino-induced samples) at the detector in a given angular () and energy bin (), is given by


where is the Avogadro number, is the time of data taking, is the post-trigger efficiency, is the density of ice, is the effective volume (the effective mass, , is given in Tab. 1 for the different configurations we consider) and and are the deep inelastic scattering cross sections for neutrinos and antineutrinos off an isoscalar target16,


The number of atmospheric neutrino- and antineutrino-induced contained events that are expected in the PINGU-0 set-up after an exposure of 10 years is shown in Fig. 4 as a function of . The left panel corresponds to the 5-10 GeV bin and the right panel to the 10-15 GeV energy bin. The number of events for both neutrino mass hierarchies, NH and IH, are shown for neutrinos and antineutrinos. We also depict two cases where we increase the Earth’s density by 10% with respect to the PREM: neutrino-induced events for NH and antineutrino-induced events for IH, where resonant matter effects are at play. We show the neutrino and antineutrino rates independently for the sake of illustrating the matter effect taking place for one case or the other, depending on the mass hierarchy, although in the analysis they are added.

These plots show the dependence of the numbers of muon-like events on the mass hierarchy (compare the red solid and the blue dotted lines for neutrinos with the magenta dotted and the black dotted lines for antineutrinos) and on changes of the normalization of the matter profile (compare the red solid with the green dotted lines for neutrinos and NH and the magenta dotted with the brown dotted lines for antineutrinos and IH). We only show the results with a change in the density with respect to the PREM for the neutrino-NH events and the antineutrino-IH events because matter effects take place inside the Earth in the case of NH for antineutrinos and of IH for antineutrinos. The effect is larger for the neutrino channel mainly due to the larger cross section and due to the fact that the flux of neutrinos is slightly higher than that of antineutrinos, and , for the energies of interest here [115, 116, 117, 118].

Figure 4: Number of muon-like contained events in PINGU-0 after 10 years for NH and IH in the case of a PREM matter profile, in the 5-10 GeV bin (left) and 10-15 GeV (right) as a function of (also shown is the corresponding baseline on the upper axis). We also show the event numbers in the case of a PREM matter profile with an overall fluctuation of for neutrino with NH and antineutrino with IH. Note that the vertical axis is different in each panel. Although shown independently, in the analyses we add together the neutrino and antineutrino events, due to the lack of charge-identification of the detectors.

A general feature in both panels is the absence of significant matter effects for ( km), as expected (see Sec. 2.3). For neutrinos crossing the Earth deeply, resonant -driven matter effects are very important and there are clear differences in the number of events between neutrino-NH and neutrino-IH (antineutrino-NH and antineutrino-IH), and between neutrino-NH (antineutrino-IH) for a density profile given by the PREM and for the case of increasing the density by 10%. In particular, notice that there are important differences in the total number of muon-like events in the angular interval ( km) for the two cases of the depicted Earth’s density, being always more important for NH as it is in this case when the neutrino channel is mostly affected by matter effects. The larger the density the lower the resonant energy and the shorter the baseline of the first oscillation maximum. Hence, the crossing of the number of events of these two cases (e.g., the red solid and green dotted lines) for is due to the interplay between these two effects. On the other hand, the matter effect on the propagation of core-crossing neutrinos is smaller than for just-mantle-crossing neutrinos. This is due to the resonance energy for the former being outside the energy range we consider here, and hence the effect is due to non-resonant matter effects that modify the propagation in a less significant manner.

The matter effect is larger in the lower part of the energy range we consider, as is clear from both panels and was discussed above. In general, the resonance energy for just-mantle-crossing neutrinos lies within the 5-10 GeV energy bin, whereas in the 10-15 GeV energy bin the effects are non-resonant for all trajectories. The minimum at for 10-15 GeV (right panel) is explained by the minimum in the probability (upper-right plot in Fig. 2), which moves to smaller values of (longer baselines) for higher energies. Likewise, the minimum at for 5-10 GeV is explained by the first minimum in the probability (upper-left and upper-right plots in Fig. 2). The fact that the second minimum in this probability does not translate into a minimum in the number of muon-like events in the energy range 5-10 GeV is because the first maximum contributes significantly as the full energy range is integrated.

Overall, the best sensitivity to the neutrino mass hierarchy (difference between the red solid and blue dotted lines, and between the magenta dotted and black dotted lines) and to the matter density (difference between the red solid and green dotted lines, and between magenta dotted and brown dotted lines) is achieved in the 5-10 GeV energy range and for neutrino trajectories such that , although some sensitivity remains for other angular and energy bins. In this case, the resonant matter effects take place, whereas for higher energies or for neutrinos traversing the core, non-resonant effects occur. As expected, it is clear from both figures that the sensitivity to the mass hierarchy is better than that to changes in the Earth’s density.

3.3 Statistical analysis

True Values Marginalization Range External error
[0.019,  0.030]
[0.38,   0.66]
[-0.1,   0.1]
[-1,   1]
Table 2: Benchmark parameters used in the analyses. The first column shows the central values of the oscillation parameters used in this work as true values. The second column shows the ranges of the parameters over which the marginalizations have been performed in the fit ( CL ranges [42]). Estimated future  CL errors on the oscillation parameters are given in the last column. In the last two rows, is the correction to the normalization of the Earth’s matter density with respect to the PREM profile and represents the nuisance parameter used to incorporate a fully correlated systematic error () in the normalization of the number of events.

We have generated our simulated data for the true (central) values of the oscillation parameters given in the first column of Tab. 2. The numbers of simulated muon-like events in each bin for 10 years of data taking in DeepCore and PINGU-0 are given in Tabs. 4 and 5 in Appendix A. We also show in these tables the results after varying the normalization of the matter density.

The values of the oscillation parameters we use to estimate the sensitivities lie within the allowed ranges at  CL obtained from global fits of the current neutrino data, except for the value of  [42, 43, 44]. Recent global analyses prefer a non-maximal value of at about the  CL which is mainly due to the recent MINOS disappearance data [148], although the new Super-Kamiokande atmospheric neutrino data still favors maximal 2-3 mixing as the best-fit value in a 2-neutrino oscillation analysis [149]. In our analysis, we choose to use as the reference value. For the atmospheric neutrino mass splitting, we take an effective mass-squared difference of as measured by the accelerator experiments in the disappearance channel [148]. The () sign refers to the case where NH (IH) is the true hierarchy. This effective mass-squared difference is related to the and mass-squared differences through the expression [150, 151]


where . Since the degenerate solutions occur at , which is slightly different from , it is important to take into account this feature when evaluating the sensitivity to the neutrino mass hierarchy.

In all fits, we marginalize over , , and within their presently allowed ranges (see the second column of Tab. 2). We impose a prior (Gaussian error at ) of 5% on considering the expected precision by the Daya Bay experiment by 2016 [152]. We also take expected priors of 2% and 4% on and , respectively [153]. We keep , and fixed both in data and in theory and no marginalization has been performed for these oscillation parameters since they would affect our results in a negligible way. Finally, we also add a fully correlated systematic error in the normalization of the number of events, which could be due to errors on the normalization of the atmospheric neutrino flux, the detector effective mass, the cross section or the efficiency. We set this error to .

The treatment of correlated systematic errors is performed by the Lagrange multiplier method or pull approach [154, 155, 156, 157]. Thus, we consider nuisance parameters that describe the systematic error of the normalization of the number of events () and the errors of (), () and (). The variation of these nuisance parameters in the fit is constrained by adding a quadratic penalty to the corresponding function without systematic errors.

Let the number of muon-like events detected in the -th angular and -th energy bin be and the expected number of events in that bin, where


with and a global normalization factor for the Earth’s density so that , with being the PREM density profile. We estimate the difference between data and theory with the following expression [154, 155, 156, 157]:


where min indicates the minimum with respect to those parameters. The function is given by


and is defined as


In order to estimate the sensitivity to the global normalization of the Earth’s matter density, in the simulated experimental data we set , and in the fit we vary in the range . Then, we marginalize over the two possibilities for the neutrino mass hierarchy (), so that we obtain . The sensitivity is calculated as


We present our results for both NH () and IH () as the true mass hierarchy.

In the case of the determination of the neutrino mass hierarchy, we define


for true NH () and true IH (), where


with () for NH (IH) defined as


Then we marginalize over in the range to incorporate the impact of possible fluctuations in the global normalization of the density of all the layers of the PREM profile (but without imposing any prior) to obtain . The sensitivity in this case is just given by


Like for the sensitivity to the Earth’s density, we also study both cases for the neutrino mass hierarchy, when data is generated assuming NH (IH) as the true neutrino mass hierarchy and we compare with the theoretical expectation when considering IH (NH).

Notice that for both analyses, only the few bins where the -driven matter effects are present (see Fig. 4 and Tabs. 4 and 5 in Appendix A) would contribute to the . Thus, correlated systematic errors are expected to have a much smaller effect than uncorrelated systematic errors (see Appendix B), as the values of the parameters related to them get fixed by the bins where the matter effect is negligible. We have explicitly checked this for both analyses.

4 Results

In this section we show the sensitivities to the Earth’s density distribution and to the neutrino mass hierarchy by using atmospheric neutrino data for the DeepCore, PINGU-0 and PINGU-I configurations. We also comment on the impact of systematic errors on the analyses, as well as the effect of the marginalization over the oscillation parameters , and .

4.1 Sensitivity to the Earth’s matter density

In Fig. 5 we show our results for the sensitivity to the Earth’s matter density, for both NH (left panel) and IH (right panel) as the true hierarchy. We show the sensitivities for the three detector configurations detailed in Tab. 1 for an exposure of 10 years. All the results shown in this section include the effects of the marginalization, within the current  CL ranges [42], over , and , on which we impose a prior as discussed in the previous section. We also marginalize over the neutrino mass hierarchy. Finally, we also include a 20% fully correlated systematic uncertainty in the total number of muon-like events. As mentioned above, this systematic error does not have a large impact on the sensitivities due mainly to the presence of the higher (closer to the horizon and thus very little affected by matter effects) angular bins in the analysis, that help to fix the normalization of the number of events in those bins.

Figure 5: Sensitivity to fluctuations on the overall normalization of the Earth’s density with respect to the PREM profile, for the DeepCore (red solid lines), PINGU-0 (blue dashed lines) and PINGU-I (brown dash-dotted lines) configurations for the case of NH as the true hierarchy (left panel) and for IH as the true hierarchy (right panel).

For the DeepCore configuration, fluctuations in the normalization of the Earth’s density of can be detected at the  CL ( CL) in the case where NH (IH) is the true hierarchy. The prospects are significantly improved for the PINGU-0 and PINGU-I configurations, due to the information contained in the lower energy bins (matter effects are maximal at the resonant energies  GeV, which are unaccessible to the considered DeepCore configuration). For NH (IH), PINGU-0 would be able to measure matter fluctuations of () at the  CL. PINGU-I would improve these numbers to () for NH (IH).

The results are always significantly better for NH. This can be understood from the fact that the matter effects are present in the neutrino channel, which has more statistics. For IH the larger number of events in the neutrino channel (not affected by matter effects in this case) substantially swamp the -driven matter effects in the antineutrino channel. As we have mentioned, these have to be added together, which reduces the sensitivity to density fluctuations, specially in the case of IH. On the other hand, it is interesting to note that PINGU-0 and PINGU-I are slightly more sensitive to negative fluctuations of the density than to positive ones. This is because the higher the density, the lower the resonance energy (see lower panel of Fig. 3). Thus, with an energy threshold of 5 GeV, some of the events affected by the matter effects would lie outside the energy range we consider. However, for lower densities, more events around the resonance region are included in the data set, which compensates the loss of sensitivity due the marginalizations. This has a smaller effect in DeepCore since the resonant effects do not take place in the energy bin we consider. We have also explicitly checked that, for the assumed energy thresholds, these detectors are basically sensitive to fluctuations of the normalization of the Earth’s mantle ( km), as the resonance for cross-crossing neutrinos occurs at lower energies.

Finally, let us comment on the impact of the marginalization over all the different parameters, which is small. Indeed, for negative density fluctuations, it is negligible for PINGU-0 and PINGU-I. For positive fluctuations, the marginalization over has a small but non-negligible effect for the three detector configurations. The value of sets the amplitude of the matter effects (see Eqs. (2) - (6) and Eqs. (7) - (11)) and therefore, in principle, one expects its marginalization to have some impact. As a consequence of what was just discussed above, the effect would be larger for the case of positive fluctuations for which the sensitivity is slightly worse. The configuration for which the largest effects of the marginalizations occurs is PINGU-0 (at the level of ). DeepCore does not include the resonance region and hence, the effects of the marginalizations are similar for negative and positive fluctuations. On the other hand, PINGU-I has more energy bins, which helps to distinguish density effects from changes in the mixing parameters.

Notice that the measurement of the matter density profile is much more difficult than the extraction of the neutrino mass hierarchy, which we consider next. This can be explained by the fact that while just measuring matter effects alone can serve to determine the neutrino mass spectrum, for the Earth’s matter density one also has to be able to measure perturbations to the matter effects. These could be small and difficult to disentangle unless the energy and angular information of the signal is reasonably good.

4.2 Sensitivity to the neutrino mass hierarchy

Detector configuration Exposure [years] CL [] NH(true) CL [] IH(true)
1 2.90 2.63
DeepCore 5 6.17 5.36
10 8.38 6.92
1 4.98 4.91
PINGU-0 5 10.2 8.49
10 14.0 11.0
1 6.50 6.22
PINGU-I 5 13.3 10.8
10 18.4 14.5
Table 3: Sensitivity to the neutrino mass hierarchy for the three detector configurations explored here for different exposure times. We show the results for NH (third column) and IH (fourth column) as the true hierarchy.

The sensitivity to the neutrino mass hierarchy by exploiting the number of muon-like events in the three possible neutrino configurations considered in this work is shown in Tab. 3, which contains the results of the analyses for different exposure times (1, 5 and 10 years). We show the results for the two possible choices of the true neutrino mass hierarchy, NH and IH. As in the case of the study of the matter density fluctuations, for the analysis of the determination of the neutrino mass hierarchy, we also marginalize over the current allowed range of the oscillation parameters , and , imposing the same priors expected from future measurements, and add a 20% fully correlated systematic error related to the overall normalization of the number of events. In addition, here we also marginalize over the Earth’s matter density distribution, by allowing its normalization to vary freely within .

From the results in Tab. 3 we can see that the prospects of measuring the neutrino mass hierarchy with these three detector configurations are very promising. The ongoing DeepCore detector with the signal efficiency and the energy reconstruction capabilities assumed here could provide a () measurement of the neutrino mass hierarchy after slightly more than 1 year (less than 5 years) if nature has chosen NH and after less than 2 years (5 years) for IH. As we have already mentioned, the results are worse for the IH scenario because in this case the resonant behavior occurs in the antineutrino channel, which is statistically suppressed due to the smaller antineutrino cross sections and initial fluxes. On the other hand, with detector configurations such as those considered here for PINGU-0 and PINGU-I, the mass hierarchy could be measured with an astonishing precision, and after 1 year for both NH and IH, respectively.

We have checked that our results agree with previous findings in the literature [94, 95]. However, note that those analyses differ on a number of important points from the one presented here. In particular, in the recent Ref. [95], a larger effective volume for PINGU was used and no post-trigger efficiency was considered. Overall this amounts to about a factor of four more statistics in their case. In that work, a much finer binning, both in and , was also considered. This is mainly the reason why they get a sensitivity to the mass hierarchy for 5 years of  CL for NH for the case of no systematic errors and no smearing, whereas for PINGU-0 (PINGU-I), we get  CL ( CL). On the other hand, whereas we consider large bins for the neutrino energy and zenith angle to approximately take into account the uncertainties in their reconstruction which smear the signal, in Ref. [95] these uncertainties are studied more explicitly. Even with the much larger bins we consider, we would expect our results to slightly worsen when adding this information in a more accurate way. Another important difference is the way in which systematic errors are added. Whereas they add uncorrelated errors on the normalization of the number of events, we add fully correlated uncertainties. As briefly discussed in Appendix B and checked numerically, this has important consequences when only a few out of all the considered bins contribute to the , as is the case here. In the case of uncorrelated systematic errors in the normalization, the sensitivity gets reduced as the systematic uncertainty increases, whereas for correlated systematic errors the sensitivity is not much affected with respect to the case of no systematic errors, thanks to the use of bins where this error gets fixed. Hence, their results get worse the larger the systematic error is, whereas the impact of the systematic error on our results is very small. Finally, we also use different central values for the oscillation parameters, although this has a small effect.

Finally, note that the neutrino mass hierarchy sensitivities computed here are almost independent of the value of the CP violating phase, . Therefore the mass hierarchy determination for these three detector configurations would be a clean measurement, free of degeneracies, and could serve as an input for other experiments which focus on CP violation searches, such as the T2K [41] or NO[158] long-baseline neutrino experiments (see Ref. [159] for a recent analysis combining the capabilities of both experiments).

5 Conclusions

Atmospheric neutrinos produced by the interactions of cosmic rays in the Earth’s atmosphere cover a huge range of energies and baselines, offering the opportunity to address many open physics questions. In this work we have exploited the low energy region,  GeV, where matter effects could affect the propagation of these atmospheric neutrinos while they pass deep through the Earth. If is sufficiently large, the transition probabilities () and () of atmospheric neutrinos become relevant and resonant matter effects would affect the neutrino (antineutrino) oscillation channels if nature has chosen NH (IH). Hence, the recent measurement of from several reactor neutrino oscillation experiments [38, 39, 40], as well as from the accelerator based T2K experiment [41], brings the opportunity to study these -driven matter effects. This could allow us not only to determine the neutrino mass hierarchy, but also to infer some general features of the Earth’s density profile.

Indeed, the idea of exploiting matter effects in atmospheric neutrino oscillations to distinguish the type of neutrino mass hierarchy has been extensively explored in the literature [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 94, 95]. On the other hand, three different ways to study the density profile of the Earth by using neutrinos have been considered: neutrino absorption tomography [45, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64], neutrino oscillation tomography [65, 66, 67, 68, 69, 70, 71, 72, 74, 75, 76, 77, 78] and neutrino diffraction [79]. These techniques have been discussed for accelerator, extra-terrestrial, atmospheric, solar and supernova neutrinos. With regards to exploiting the atmospheric neutrino flux, only neutrino absorption tomography had been studied in some detail. Here we have considered neutrino oscillation tomography with atmospheric neutrinos.

Obtaining a reliable estimate of the density of the Earth is essential to solve a number of important problems in geophysics. In this work we have presented a completely different technique to determine the Earth’s density profile to those used in geophysics, which are mainly based on the analysis of the velocities of seismic waves and on empirical relations between these velocities and density. Although, in principle, geophysics can obtain more precise results, its inferences are based on numerous assumptions as there are significant trade-offs with other parameters. Thus, the results from atmospheric neutrino tomography would represent an independent and complementary assessment of the Earth’s internal structure.

The goal of this work is the study of the Earth matter effects taking place for atmospheric neutrinos in the GeV range. We have explored the possibility to do neutrino oscillation tomography to infer information on the Earth’s matter density and have studied future sensitivities to determine the neutrino mass hierarchy. In order to do so, we have considered kilometer-scale ice Čerenkov detectors such as the ongoing DeepCore extension of IceCube [90, 91, 92], and PINGU [93], which is a further proposed extension to IceCube. Although these multi-Mton detectors have no charge-identification capabilities that would allow us to distinguish neutrinos from antineutrinos and hence increase the sensitivity to matter effects, the different cross sections for neutrinos and antineutrinos (and also the slightly different fluxes), the large statistics which can be accumulated by these detectors, and the large value measured for the mixing angle , make the study of these matter effects feasible and make possible the determination of the neutrino mass hierarchy in relatively short time scales.

In this work, three different experimental scenarios have been examined (see Tab. 1). For the simulations of the DeepCore set-up we have assumed a single energy bin of  GeV and an effective mass of 5 Mton. For the proposed PINGU experiment, we have considered two possible configurations. For the PINGU-0 set-up, we have assumed two energy bins:  GeV and  GeV and an effective mass of  Mton. For the PINGU-I set-up we have assumed four energy bins: