Locating the Core-Mantle Boundary using Oscillations of Atmospheric Neutrinos

Atmospheric neutrinos provide a unique avenue to explore the internal structure of Earth based on weak interactions, which is complementary to seismic studies and gravitational measurements. In this work, we demonstrate that the atmospheric neutrino oscillations in the presence of Earth matter can serve as an important tool to locate the core-mantle boundary (CMB). An atmospheric neutrino detector like the proposed 50 kt magnetized ICAL at INO can observe the core-passing neutrinos efficiently. These neutrinos would have experienced the MSW resonance and the parametric or neutrino oscillation length resonance. The net effect of these resonances on neutrino flavor conversions depends upon the location of CMB and the density jump at that radius. We quantify the capability of ICAL to measure the location of CMB in the context of multiple three-layered models of Earth. For the model where the density and the radius of core are kept flexible while the mass and radius of Earth as well as the densities of outer and inner mantle are fixed, ICAL can determine the location of CMB with a 1$\sigma$ precision of about 250 km with an exposure of 1000 kt$\cdot$yr. With the 81-layered PREM profile, this $1\sigma$ precision would be about 350 km. The charge identification capability of ICAL plays an important role in achieving this precision.


Introduction and Motivation
The exact nature of deep interior of Earth has been a long-standing puzzle. The internal regions of Earth are inaccessible due to extreme temperatures and pressures. Direct measurements are not feasible after a depth of a few km. Most of the information available about the internal structure of Earth has been obtained using indirect probes such as the gravitational measurements [1][2][3][4][5][6] and the seismic studies [7][8][9][10][11][12][13][14][15]. For example, gravitational measurements provide information on the mass [1][2][3][4][5], and the moment of inertia of Earth [5,6]. These tell us that the average density of Earth deep below must be much larger than the typical rock density at the surface.
Seismic waves are generated when the tectonics plates in the outermost layers of Earth slide against each other. The origin of earthquake, called epicenter, is typically located up to a depth of about 700 km [9,12]. These seismic waves travel in all directions through the bulk of Earth and can reach the surface all over the globe. They are measured by seismometers which provide the information about the time, location, and intensity of the earthquake. While passing though Earth, these waves can get reflected or refracted on encountering matter with different densities and composition [9]. The features observed in the seismic waves are used to infer the internal structure of Earth.
From seismic studies [7][8][9][10][11][12][13][14][15], we know that the Earth consists of a layered structure in the form of concentric shells. As we go from outside to deep inside the Earth, the layers are arranged in the order of increasing densities as crust, outer mantle, inner mantle, outer core and inner core. The density changes sharply at the core-mantle boundary (CMB) which has been infered to be 3483 ± 5 km using seismological measurements [7,16]. Even though seismic studies and gravitational measurements have provided us a huge amount of data and revealed crucial information about the internal structure of Earth, there are still many open issues [13]. For example, the radius, mass and chemical composition of the core are not very well known. The density jump between outer core and inner core is not known precisely. There are many unexplained structures and heterogeneities observed in the lowermost mantle, beneath Africa and Pacific, that show lower-than-average seismic wave speeds which are known as large low-shear-velocity provinces (LLSVPs) and ultra low velocity zones (ULVZs) [14,17,18]. Active research is being pursued to determine the composition of the bulk silicate Earth [13,19]. If we talk about the presence of light elements inside Earth, we don't know how much H 2 O is present in the mantle and how much H is present in the core [15,20,21]. The precise measurement of radioactive power produced in the mantle and core is important to understand the thermal dynamics inside Earth [22][23][24][25].
Complementary information using additional probes such as neutrinos can improve our understanding about the structure of Earth. Neutrinos interact only via weak interactions, which enables most of the neutrinos to pass through Earth without getting absorbed. The cross section for interaction of neutrino with nucleons increases with energy. At energies higher then a few TeV, the interaction cross-section is large enough to have sizable absorption of neutrinos inside the Earth [26,27]. The idea of using the attenuation in the flux of neutrinos to probe the internal structure of Earth was proposed in Ref. [28,29], and detailed studies using neutrinos from various sources, such as man-made neutrinos [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42], extraterrestrial neutrinos [35,[43][44][45][46], and atmospheric neutrinos [47][48][49][50] have been performed. The analysis in Ref. [47] made a forecast that using the absorption of high-energy atmospheric neutrinos (in the range of 10 to 100 TeV), IceCube can reject the homogeneity of the Earth at 3.4σ confidence level in 10 years with conservative assumptions on the theoretical and systematic uncertainties. The analysis of the one-year data of multi-TeV atmospheric muon neutrinos at IceCube, carried out in Ref. [51], estimated the densities of various Earth layers using the absorption of high-energy neutrinos for the first time. They determined the mass of the Earth to be M ν E = 6.0 +1.6 −1.3 ×10 24 kg, which is in good agreement with the gravitational measurement. Another way of exploring the interior of Earth can be the use of diffraction patterns produced by the coherent scattering of neutrinos with the Earth's matter, but this is not feasible with the present technology [52].
Further information about the chemical composition of the Earth can also be obtained using the geoneutrinos which are produced during the decay of radioactive elements such as Uranium, Thorium and Potassium [22-25, 99, 100]. Geoneutrinos may shed light on the radiogenic contribution to the heat budget of Earth. Since neutrinos use weak interactions to probe the internal structure of Earth, which is complementary to seismic studies based on electromagnetic interactions and gravitational measurements based on gravitational interactions, this would pave the way for "multi-messenger tomography of Earth".
The proposed 50 kt ICAL detector at the upcoming India-based Neutrino Observatory (INO) [101] would be able to detect atmospheric neutrinos and antineutrinos in the multi-GeV range of energies and over a wide range of baselines. Thanks to the presence of magnetic field of about 1.5 T, ICAL would be able to distinguish between neutrino (by observing µ − ) and antineutrino (by observing µ + ) events separately. Further, due to its good directional resolution, ICAL would be able to observe and identify the neutrinos passing through core and mantle separately. The ICAL detector would be sensitive to Mikheyev-Smirnov-Wolfenstein (MSW) resonance [102][103][104], which occurs around the energies of 6 to 10 GeV for mantle-passing neutrinos. It can also observe neutrino oscillation length resonance (NOLR) [105][106][107][108][109] or parametric resonance (PR) [110,111], which occurs around the energies of 3 to 6 GeV for some of the core-passing neutrinos. The NOLR/ parametric resonance occurs due to a sudden jump in the density when we go from mantle to core. The pattern of NOLR/ parametric resonance depends upon the amount of density jump and the location of CMB. In the present work, we explore the impact of modifying the CMB radius on neutrino oscillations and hence, on the reconstructed muon events at the ICAL detector. We demonstrate that the location of core-mantle boundary can be measured using the matter effects in atmospheric neutrino oscillations with the help of the ICAL detector.
In section 2, we discuss the internal structure of Earth and the propagation of seismic waves through different regions of Earth. In section 3.1, we describe our methodology for exploring various modified-CMB scenarios with a simple three-layered density profile of Earth. The impact of CMB modification on neutrino oscillation probabilities is discussed in section 3.2 and 3.3. The neutrino event simulation at ICAL is described in section 4. The method of statistical analysis to quantify the impact of CMB modification is explained in section 5. In section 6, we evaluate the sensitivity of ICAL for determining the location of CMB. Finally, we present the summary and concluding remarks of this study in section 7. In appendix A, we demonstrate the possible origin of asymmetric and non-monotonic behavior of sensitivities for modified-CMB scenarios with respect to the standard CMB in threelayered profile at the probability level. In appendix B, we present the sensitivities for determining the location of CMB using the 81-layered PREM profile.

A Brief Review of the Internal Structure of Earth
The surface of Earth consists of soil, sand, rocks, mountains, rivers, oceans, ice, and lava etc, which differ from each other geologically as well as in terms of chemical composition. The layers beneath these structures are mostly solid and consist of various types of rocks. The direct observations of layers inside Earth have been possible only up to a few kms because the deepest hole in the world till today, with current technology, is only about 12 km which was drilled on the rocks of Kola Peninsula in the Murmansk region of Russia [112,113]. Below the depth of 6 km, a temperature gradient of 20 • C per km was observed instead of the expected gradient of 16 • C per km. Eventually, the drilling had to stop due to extreme temperature of about 180 • C. This exploration clearly shows that if we want to know more about the structure of Earth at the deeper locations, then we need to use indirect methods, for example, gravitational measurements [1][2][3][4][5][6], studies of earthquakes or seismic waves [7][8][9][10][11][12][13][14][15], neutrino tomography [114], etc.
Seismic waves are generated when the tectonic plates inside Earth slide and the energy is released in the form of vibrations. They get modified while propagating through Earth. The velocity and timing measurements of seismic waves are used to unravel the internal structure of Earth [7,8,10,11]. Seismic studies have revealed that the Earth consists of the layers in the form of concentric shells which can be broadly divided into core and mantle. The radius of core is about half the radius of Earth (R E ) which is about 6371 km. The contribution of core and mantle to the mass of Earth are about 32% and 68%, respectively.
The information about Earth is also obtained using gravitational measurements. The mass of Earth (M E ) is gravitationally measured to be about (5.9722 ± 0.0006)×10 24 kg [1][2][3][4] whereas its moment of inertia is about (8.01736 ± 0.00097) × 10 37 kg m 2 [5,6]. Since the measured moment of inertia is less than the expected one (2/5M E R 2 E ∼ 9.7 × 10 37 kg m 2 ) for a uniform density distribution inside Earth, the density should be more as we go deeper inside the Earth. In other words, the matter should be concentrated more towards the center of Earth. Using M E and R E , one can obtain the average density of Earth to be Earth as a function of radial distance from its center. 3 , which is larger than the density of rocks (∼2.8 g/cm 3 ) present on the surface. This observation also points towards the presence of regions of higher densities deep inside the Earth. Now, let us try to understand, how the interiors of Earth are probed using seismology. Seismic waves can be classified into two categories -shear (S) waves and pressure (P) waves. The S-waves result into the vibration of rocks perpendicular to the direction of propagation whereas the P-waves force the particles to vibrate along the direction of propagation. The S-waves and P-waves modify differently while passing through the layers of Earth. The P-wave can travel through both solid as well as liquid layers but the velocity of the P-wave decreases while passing through liquid layers. As far as the S-waves are concerned, they cannot travel through the liquid layers.

g/cm
The seismic waves originating from the center of the earthquake travel through the Earth and land at various positions on the surface. During their propagation through Earth, S and P-waves may get reflected and refracted at multiple density discontinuities inside the Earth, splitting into different segments as shown in the left panel of Fig. 1. The nomenclature of these segments is as follows: • S and P indicate the S-wave and P-wave traveling through the mantle.
• K and I stand for the P-waves passing through the outer and inner core, respectively.
• J indicates the segment of the S-wave traveling through the inner core.
• s and p indicate those S and P-waves, respectively, which initially travel upward from the center of the earthquake, and then get reflected downward from the Earth's outer surface.
• c represents an upward reflection at the core-mantle boundary.
The seismic studies indicate that the mantle consists of hot rocks of silicate [8]. The rocks in mantle are not molten. However, they are plastic in nature which allows them to change their shape over long timescales. Beneath the mantle, the Earth consists of a high density core which is composed of metals like iron and nickel. The inability of S-waves to pass through the outer core, and the decrease in the velocity of P-waves therein, indicate that the outer core is expected be liquid [8]. In 1936, I. Lehmann discovered the inner core by its higher P-wave velocity [115]. Inner core is expected to be solid [116]. The only possible reason why iron and other heavy metals may be solid at high temperatures inside the inner core would be because their melting points increase significantly at the tremendous pressure present there [116].
The data on the velocities of seismic waves are used to develop the Preliminary Reference Earth Model (PREM) [117] profile where the density of any layer is given as a one-dimensional function 1 of the radial distance of the layer from the center of Earth. The PREM profile is mainly based on two empirical equations which relate the velocities of S and P-waves with the densities of the layers inside Earth. The first one is called the Birch's law [123] and is valid for the outer mantle whereas the second one is known as Adams-Williamson equation [124] that is suitable for the inner mantle and the core. The parameters of these empirical relations depend upon temperature, pressure, composition, and elastic properties of the layers of Earth, which give rise to uncertainties in density distribution. The density of the mantle has a uncertainty of about 5% whereas for the core, it is significantly larger [125][126][127].
The sharp change in the densities of layers results in the partial reflection and refraction of seismic waves. Such a sudden rise in density is also observed at the CMB. In the present work, we explore whether the matter effects in atmospheric neutrino oscillations can be used to probe the location of CMB. We consider a simple three-layered model of Earth which consists of core, inner mantle, and outer mantle as shown in the right panel of Fig. 1. The core is present up to a radial distance of 3480 km from the center of Earth, the inner mantle spans from 3480 km to 5701 km, and the outer mantle is from 5701 km to the radius of Earth. The densities of core, inner mantle and outer mantle are taken as 11.37, 5, and 3.3 g/cm 3 . Note that we have not considered the crust which has much smaller thickness as compared to the other layers.

Three-layered Models with Modified CMB Locations
In the present work, we use atmospheric neutrinos to probe the location of CMB where we explore how the matter effects in neutrino oscillations change in different modified-CMB scenarios. The standard CMB radius is taken to be R CMB (standard)= 3480 km. For illustration, we show the modification of the CMB radius by ∆R CMB = −500 km (smaller core or SC) and by ∆R CMB = +500 km (larger core or LC) in Fig 2. In the standard, SC, and LC scenarios, the core spans the zenith angle range cos θ ν ≤ −0.84, cos θ ν ≤ −0.88, and cos θ ν ≤ −0.78, respectively. We further consider three different ways of modifying the density profile of the Earth. We term these as three "cases": • Case-I: The densities of all layers of Earth are kept fixed. Note that the mass of Earth is not constrained in this case.
• Case-II: The densities of inner and outer mantle are kept fixed, while the density of core is modified to keep the mass of Earth invariant.
• Case-III: The density of outer mantle is kept fixed, and the densities of core and inner mantle are modified, while keeping their individual masses invariant. It automatically  ensures that the mass of Earth is unchanged.
Note that the Case-II is the most realistic one, since the density of the mantle is quite well measured by the combination of gravitational and seismological studies. However, we also analyze two dummy cases, Case-I and Case-III. In Case-I, we remove the constraint from the mass of Earth while imposing the constraint on the density of the core, compared to Case-II. In Case-III, we add the constraint on the mass of the mantle but remove the constraint on its density, compared to Case-II. This would allow us to obtain insights on the role of these constraints. In this study, we have not taken into account the constraint from the moment of inertia of Earth. Table 1 presents the modified boundaries between various layers and their densities in the SC and LC scenarios. The densities of core and mantle may change depending upon the cases as explained above. We discuss the effect of R CMB modification on neutrino oscillations in the next section.

Effect of Modified CMB Radius on Oscillation Probabilities
Atmospheric neutrinos and antineutrinos can be observed at ICAL separately in the mulit-GeV range of energies over a wide range of baselines from about 15 km to 12757 km. The upward-going neutrinos pass through the Earth and undergo charged-current interactions with the ambient electrons. This coherent forward scattering results in a matter potential V CC given by where G F is Fermi coupling constant and N e is the ambient electron number density. Further, Y e = N e /(N p +N n ) is the relative number density of electron inside the matter having density ρ where N p and N n denote the number densities of protons and neutrons. In this analysis, we assume that Earth is neutral and isoscalar, which implies N n ≈ N p = N e and Y e = 0.5. The positive and negative signs correspond to neutrinos and antineutrinos, respectively. Due to these opposite signs, the matter effects modify the oscillation patterns for neutrinos and antineutrinos differently. The charge identification capability of ICAL plays a crucial role in observing these different matter effects in neutrinos and antineutrinos separately, which enhances the sensitivity of ICAL towards locating the CMB as we demonstrate in our results. In the present work, we use the benchmark values of oscillation parameters given in Table 2. The normal mass ordering (NO) represents m 1 < m 2 < m 3 whereas for inverted mass ordering (IO), m 3 < m 1 < m 2 . We implement NO vs. IO by taking opposite signs of the effective atmospheric mass-squared difference 2 ∆m 2 eff . The matter effects result in a resonant enhancement in the effective value of the smallest lepton mixing angle θ 13 , which can become as large as 45 • . This adiabatic resonance due to θ 13 -driven matter effects is known as the Mikheyev-Smirnov-Wolfenstein (MSW) resonance [102][103][104]. For NO (IO), the MSW resonance occurs in neutrinos (antineutrinos). For a given matter density ρ, the energy at which the MSW resonance occurs can be given as Now, if we consider the average density of mantle to be around 4.5 g/cm 3 , then the MSW resonance occurs at approximately 7 GeV: Therefore, the mantle-passing neutrinos feel the MSW resonance in the energy range of about 5 to 10 GeV. On the other hand, if we consider the density of core to be about 11.3 g/cm 3 , the MSW resonance in the core occurs at approximately 2.8 GeV: The effective atmospheric mass-squared difference can be given in terms of ∆m 2 31 and ∆m 2 21 as follows [128,129] ∆m 2 eff = ∆m 2 31 − ∆m 2 21 (cos 2 θ12 − cos δCP sin θ13 sin 2θ12 tan θ23).  Table 2, where we assume NO and sin 2 θ23 = 0.5.
As we go along the trajectory of a core-passing neutrino, first the density increases from mantle to core and then it decrease from core to mantle. The quasi-periodic nature of densities along the neutrino path, combined with a sharp contrast in the densities of core and mantle at CMB, may result in specific phase relationships between neutrino oscillation amplitudes in the core and the mantle for particular path lengths. This phenomenon is known as the neutrino oscillation length resonance (NOLR) [105][106][107][108][109] or parametric resonance (PR) [110,111], which may enhance neutrino flavor conversions.
The muon neutrino events at ICAL are contributed by both ν µ → ν µ survival as well as ν e → ν µ appearance channels. However, since the contribution of the survival channel is more than 98%, here we present the effect of R CMB modification only on the oscillation probabilities for ν µ → ν µ survival channel. Figure 3 shows the three-flavor ν µ survival probability P (ν µ → ν µ ) as a function of energy for NO. The top, middle, and bottom panels are for Case-I, Case-II, and Case-III, respectively. We consider the baseline (antineutrino). We use the three-flavor neutrino oscillation parameters given in Table 2, where we assume NO and sin 2 θ23 = 0.5.
of 11000 km (11500 km) in the left (right) panels. For the baseline of 11500 km, neutrinos pass through the core for all the three scenarios -standard core (R CMB = 3480 km), SC (∆R CMB = − 500 km), and LC (∆R CMB = + 500 km). However, for the baseline of 11000 km, neutrinos pass through the core in the scenarios where we have standard core and the larger core (LC). But in the scenario where we consider smaller core (SC), neutrinos do not pass through the core for the 11000 km baseline. This feature gives rise to additional modifications in neutrino oscillation probabilities for baselines around 11000 km and thus, contributes to the sensitivity for locating CMB. The θ 13 -driven matter resonances inside Earth at atmospheric frequency result in significant modification of ν µ survival probability, P (ν µ → ν µ ). The curves corresponding to SC and LC scenarios have different patterns with respect to each other for a given case and baseline. The oscillation patterns are also quite different in the three cases. The deviations in oscillation patterns mostly occur in the energy range of 3 to 6 GeV, which corresponds to the NOLR/parametric resonance. These observations imply that the neutrinos are highly sensitive to the modification of CMB as well as the way of modifying the density profile. Further, a comparison of the left and right panels shows that the modifications in neutrino oscillation patterns due to R CMB modification also depend upon the baselines through which neutrinos have passed.

Effect of Modified CMB Radius on Oscillograms
Since the modifications in neutrino oscillation probabilities depend upon the baselines, let us discuss the effects of modification of R CMB on neutrino oscillation probabilities in the two-dimensional plane of energy and direction. First, we consider the standard scenario of R CMB = 3480 km in Fig. 4, where we present the three-flavor ν µ survival probability oscillograms in the plane of (E ν , cos θ ν ) for NO. The standard R CMB is denoted by a vertical dotted-white line. The left and right panels show the survival probabilities P (ν µ → ν µ ) and P (ν µ →ν µ ), respectively. In each panel, we consider the neutrino energy range of 1 to 25 GeV and neutrino zenith angle (cos θ ν ) range of -1 to 0. Here, cos θ ν = -1 (cos θ ν = 1) represents the upward-going (downward-going) neutrinos. The dark-blue diagonal band which starts from (E ν = 1 GeV, cos θ ν = 0) and ends at (E ν = 25 GeV, cos θ ν = -1), corresponds to the first oscillation minimum, which is also known as "oscillation valley" [130,131]. In the left panel, the red patch around -0.8 < cos θ ν < -0.5 and 6 GeV < E ν < 10 GeV corresponds to the MSW resonance whereas the yellow patches around cos θ ν < -0.8 and 3 GeV < E ν < 6 GeV are found to be due to the NOLR/ parametric resonance. With NO, both these resonances are experienced only by neutrinos. As far as antineutrinos are concerned, we do not see the MSW and the NOLR/ parametric resonance in the right panel for normal ordering. On the other hand, for inverted ordering, this trend is opposite, i.e., antineutrinos experience a significant amount of matter effects whereas neutrinos do not. In this section, we will only focus on neutrino survival channel and NO. Now, we describe how do the neutrino oscillograms modify during the modification of CMB radius. The top, middle, and bottom rows in Fig. 5 present the ν µ survival probability oscillograms for the Case-I, Case-II, and Case-III, respectively. The left (right) panels correspond to the SC (LC) scenario. The observations are as follows: • Comparing the oscillograms in Fig. 5 with those in Fig. 4, we observe that, for both Case-I and Case-II, significant effect on the probabilities occurs in the NOLR/ parametric resonance region.
• In the standard core scenario (left panel of Fig. 4), there are three yellow patches in the NOLR/ parametric resonance region. In both Case-I and Case-II, the SC scenario shows that one of the yellow patch is missing and the other two are shrunk. In the LC scenario, these patches are stretched towards lower baselines.
• However, for Case-III, we see that the NOLR/ parametric resonance as well as the MSW resonance region both are affected in the SC and LC scenarios, the effect being more prominent for LC.
For a clear picture of the energies and baselines where the probability is affected significantly, we present Fig. 6 which shows the difference between ν µ survival probability for the standard R CMB and SC/LC scenario. The left and right panels show the values of ∆P SC and ∆P LC , respectively, which are defined as It may be observed that the probability differences are nonzero only at the core regions for Case-I (top panels) and Case-II (middle panels). ∆P SC/LC vanishes for the mantle region  Table 2, where we assume NO and sin 2 θ23 = 0.5.
because in both these cases, the density of the mantle remains the same (see Table 1). In contrast, for the Case-III (bottom panels), the differences are significant in the core as well as mantle.  Table 2, where we assume NO and sin 2 θ23 = 0.5.
From Figs. 4, 5, and 6, we can infer that the effect of R CMB modification manifests itself mainly at the higher baselines and lower energies. The dependence of this effect on the zenith angle (cos θ ν ) implies that we need a detector like ICAL with high directional resolution to probe the position of the CMB.

Event Generation at ICAL
The 50 kton magnetized iron calorimeter (ICAL) detector at the proposed India-based Neutrino Observatory (INO) [101] is going to detect the atmospheric neutrinos and antineutrinos separately in multi-GeV energy range and over a wide range of baselines. ICAL, with a total size of 48 m × 16 m × 14.5 m, will have about 151 alternative layers of iron having a thickness of 5.6 cm, with glass Resistive Plate Chambers (RPCs) [132][133][134] sandwiched between the iron layers. The iron plates act as passive detector elements for the neutrino interactions, whereas RPCs act as active detector elements. In other words, the iron plates provide the target mass for neutrino interactions, and RPCs detect the secondary particles like muons and hadrons that are produced during the charged-current (CC) interactions of neutrinos with the iron nuclei. The charged particles deposit their energies in the form of hits in the RPCs during their propagation inside the detector. The X and Y coordinates of the hits are given by the pickup strips using the produced electronic signals. At the same time, the RPC layer number provides the Z coordinate of the hit.
A muon in the multi-GeV energy range is a minimum ionizing particle. Hence, it can pass through many layers, leaving a hit in each layer. These hits produced by muons form a track-like event. The magnetic field of about 1.5 T [135] enables ICAL to distinguish between the atmospheric neutrinos and antineutrinos by identifying the opposite curvature of µ − and µ + tracks. Further the nanosecond-level time resolution of RPCs [136][137][138] helps ICAL to separately identify the upward-going and downward-going muon events. In the multi-GeV energy range, the resonance and deep inelastic scatterings (DIS) give rise to the production of hadrons which can deposit energy in the form of multiple hits in the same layer of RPC, resulting into shower-like events. During neutrino interactions, a large fraction of neutrino energy is carried away by the hadrons, which is quantified as In this work, we simulate the unoscillated neutrino events using the NUANCE [139] Monte Carlo (MC) neutrino event generator with the ICAL detector geometry as a target. The atmospheric neutrino flux at the proposed INO site [140,141] at Theni district of Tamil Nadu, India, is used as an input to NUANCE. The solar modulation effect on atmospheric neutrino flux is incorporated by considering the flux with high solar activity (solar maximum) for half exposure and low solar activity (solar minimum) for another half exposure. A mountain coverage of about 1 km (3800 m water equivalent) at the INO site acts as a filter to reduce the downward-going cosmic muon background by a factor of ∼ 10 6 [142]. Further, the ICAL analysis considers the events having vertices far from the edges and completely inside the detector to exclude the events entering from outside [101]. Therefore, the background due to the downward-going cosmic muons is expected to be negligible. In our analysis, the statistical uncertainties are minimized by generating MC unoscillated neutrino events for a large exposure of 1000 years for the ICAL detector. The three-flavor neutrino oscillations in the presence of Earth's matter effects are incorporated using the reweighting algorithm [143][144][145].
In the present analysis, we incorporate the detector response for muons and hadrons as described in refs. [146,147]. These detector responses are obtained by ICAL collaboration SC  LC  SC  LC  Case-I  8842 8858 4032 4032  Case-II 8844 8862 4030 4034  Case-III 8850 8876 4032 4032   Table 3: The total number of reconstructed µ − and µ + events expected at the 50 kt ICAL detector in 20 years for the cases I, II, and III. The SC (LC) scenario corresponds to ∆RCMB = −500 km (+500 km). The number of reconstructed events for the standard CMB is 8850 (4032) for µ − (µ + ). We use the three-flavor neutrino oscillation parameters given in Table 2, where we assume NO and sin 2 θ23 = 0.5.
after performing a detailed GEANT4 [148] simulation of the ICAL detector [146,147]. The procedure for obtaining the detector response for muons has been discussed in detail in ref. [146]. The muon detector response is simulated by passing a large number of muons through the ICAL detector. The muon forms a track-like event while passing through various RPC layers. These tracks are fitted using a Kalman filter technique to obtain the energy, direction, and charge of reconstructed muons [149]. The ICAL reconstruction algorithm requires at least 8 to 10 muon hits to reconstruct a track. Since muon deposits an energy of about 100 MeV in each layer of iron, the energy threshold of ICAL is about 1 GeV. Reference [146] provides the reconstruction efficiency, energy resolution, angular resolution, and charge identification (CID) efficiency of ICAL for muons (figures 13, 11, 6 and 14 therein, respectively).
Along with the reconstruction of muons, ICAL can also retrieve information about the hadron energy using the total number of hits in the shower-like events. The details about the hadron energy resolution of the ICAL detector is given in ref. [147]. Since the hadron energy resolution is much poorer than the muon energy resolution, we do not add them to obtain neutrino energy. Instead, to exploit the measured four-momentum of muon and hadron energy on an event-by-event basis, we employ a binning scheme having reconstructed muon energy (E rec µ ), muon direction (cos θ rec µ ), and hadron energy (E ′ rec had ) as three independent observables. The detector properties are folded in following the procedure mentioned in [143][144][145]. For the analysis, the reconstructed events are scaled from 1000-yr MC to 20-yr MC. For 1 Mt·yr exposure of ICAL, we would get about 8850 reconstructed µ − and 4032 reconstructed µ + events using the three-flavor neutrino oscillation for propagation through the Earth's matter considering the three-layered profile of Earth with the standard core if the mass ordering is normal. In Table 3, we present the number of reconstructed µ − and µ + events for NO for all the three cases for SC and LC scenarios for the exposure of 1 Mt·yr at the ICAL detector. It is important to note that the total event rate in Table 3 is almost identical for all the three cases. However, after binning these events in the above three observables, the three cases would look quite different. This is possible because ICAL would have good energy and directional resolutions for reconstructed muons. Now we discuss the effect of R CMB modification on the distributions of reconstructed muon events at the ICAL detector. In the left (right) panels of Fig. 7, we present the distributions of differences between events of the standard core and that of smaller (larger) core in the plane of (E rec µ − , cos θ rec µ − ) for 20 years (1 Mt·yr) of exposure considering NO. For demonstrating the difference of reconstructed event distributions with 1 Mt·yr exposure, we use a binning scheme 3 where we have 10 bins in E rec µ and 20 bins in cos θ rec µ . Since for NO, significant amount of matter effects are present for neutrinos but not antineutrinos, we choose to present the distributions of event difference in Fig. 7 for µ − only. The event difference is quite small for µ + , and hence, not shown here. For Case-I and Case-II, event differences mainly occur in the core region, although some differences can also be observed in the mantle region due to the angular smearing originating from the difference between the directions of incoming neutrinos and reconstructed muons. On the other hand, for Case-III, we see significant event differences which span both the core as well as in the mantle regions. In the next section, we describe the numerical procedure used to estimate the median sensitivity of ICAL to probe the CMB.

Binning Scheme for Analysis
In this section, we present the binning scheme used for numerical analysis. From Figs. 6 and 7, it is clear that R CMB modification manifests itself mainly through the regions of higher baselines (−1 < cos θ < −0.7) and lower energies (1 < E ν < 6 GeV). Therefore, we use finer bins for this region. The optimized binning scheme is shown in Table 4. We have total 16 bins for E rec µ spanning in the range of 1 to 25 GeV, 39 bins for cos θ rec µ from -1 to 1, and 4 bins for E ′rec had in the range of 0 to 25 GeV. For our analysis, only upward-going neutrinos are relevant because they have experienced the matter effects, but we have also included downward-going neutrinos (0 < cos θ < 1) because they help in reducing the impact of normalization errors in atmospheric neutrino events. This also includes those upwardgoing (near horizon) neutrino events that result in the downward-going reconstructed muon events due to angular smearing during the interaction of neutrinos as well as reconstruction. We have considered the same binning scheme for µ − and µ + .

Numerical Analysis
To estimate the sensitivity of ICAL, a χ 2 analysis is performed which is expected to give median sensitivity in the frequentist approach [150]. For this analysis, we define the following Poissonian χ 2 − [151] for reconstructed µ − events as considered in ref. [143]: 3 For E rec µ , we take 5 bins of 1 GeV in (1 − 5) GeV, 1 bin of 2 GeV in (5 − 7) GeV, 1 bin of 3 GeV in (7 − 10) GeV, and 3 bins of 5 GeV in (10 − 25) GeV. On the other hand, for cos θ rec µ , we choose uniform bins with a width of 0.1 in the range of −1 to 1. Note that for analysis, we will use a finer binning scheme as described in Sec. 5.1. scenario with ∆RCMB = −500 km (+500 km). The dotted-red and dotted-cyan curves represent the CMB radius for the smaller core and larger core, respectively. Top, middle, and bottom rows correspond to the Case-I, Case-II, and Case-III, respectively. We use the three-flavor neutrino oscillation parameters given in Table 2, where we assume NO and sin 2 θ23 = 0.5.

Range
Bin width Number of bins  with N theory Here N theory ijk and N data ijk correspond to the expected and observed number of reconstructed µ − events in a given (E rec µ , cos θ rec µ , E ′rec had ) bin, respectively. The quantity N 0 ijk stands for the number of expected events without considering systematic uncertainties. From the binning scheme mentioned in Table 4, N E rec µ = 16, N cos θ rec µ = 39, and N E ′rec had = 4. For our analysis, we use the well-known method of pulls [152][153][154] to incorporate the following five systematic uncertainties [144,145]: (i) 20% uncertainty on flux normalization, (ii) 10% uncertainty on cross section, (iii) 5% energy dependent tilt error in flux, (iv) 5% zenith angle dependent tilt error in flux, and (v) 5% overall systematics. The pull variables of systematic uncertainties are denoted by ξ l in Eqs. 5.1 and 5.2.
The χ 2 + for reconstructed µ + events is also obtained by following the same procedure. The separate contributions of both χ 2 − and χ 2 + are added to get the resultant sensitivity of the ICAL detector, which is define as χ 2 : To simulate the MC data for our analysis, we use the benchmark values of oscillation parameters given in Table 2 as true parameters. In the fit, we minimize χ 2 over the pull variables ξ l and the relevant oscillation parameters. We vary the atmospheric mixing angle sin 2 θ 23 in the range (0.36 − 0.66) and atmospheric mass-squared difference |∆m 2 eff | in the range (2.1 − 2.6) × 10 −3 eV 2 . We also minimize over both the choices of neutrino mass orderings, NO and IO. During the fit, we do not vary solar oscillation parameters sin 2 2θ 12 and ∆m 2 21 ; they are kept fixed at their true values given in Table 2. For the reactor mixing angle, which is already very well measured [53][54][55][56], we take a fixed value of sin 2 2θ 13 =  CMB sensitivities for the modification of RCMB by ± 500 km for all three cases with an exposure of 500 kt·yr (10 years) and 1 Mt·yr (20 years). For the numbers without parentheses, the ∆χ 2 CMB minimization has been performed by varying test values of sin 2 θ23 and ∆m 2 eff over their current uncertainties, and taking both mass orderings NO and IO. The numbers given in parentheses correspond to a fixed-parameter scenario where we do not marginalize over oscillation parameters in the fit. The results in the second and fourth (third and fifth) columns are with (without) charge identification. For simulating MC data, we use three-flavor neutrino oscillation parameters given in Table 2, where we assume NO and sin 2 θ23 = 0.5. 0.0875 both in MC data and theory. We kept δ CP = 0 throughout our analysis in both data and theory.

Results
In this section, we present the statistical significance of the ICAL detector to probe the location of the core-mantle boundary. For numerical analysis, we simulate the MC data assuming the three-layered profile as a true profile of the Earth with standard core. We quantify the statistical significance of the analysis for measuring the position of CMB in the following way: where χ 2 (modified R CMB ) and χ 2 (standard R CMB ) is calculated by performing a fit to the MC data with the modified and standard R CMB , respectively. Here, we calculate the Asimov sensitivity representing the median ∆χ 2 CMB in the frequentist approach where the statistical fluctuations are suppressed such that χ 2 (standard R CMB ) ≈ 0.
In Table 5, we present the ∆χ 2 CMB sensitivity with an exposure of 500 kt·yr and 1 Mt·yr for all three cases in the SC (R CMB = 2980 km) and LC (R CMB = 3980 km) scenarios. It is apparent from Table 5 that the ICAL detector is sensitive to modification of R CMB by ± 500 km with a statistical significance of more than 1σ for all the three cases with an exposure of 1 Mt·yr and the CID capability of ICAL detector plays a crucial role in achieving these sensitivities. Comparing the numbers which are given with and without parentheses in Table 5, we learn that the impact of marginalization over oscillation parameters in the fit is negligible in our study.   Table 2, where we assume NO and sin 2 θ23 = 0.5.
In Fig. 8, we present the sensitivities of the ICAL detector while modifying R CMB in smaller steps of 50 -100 km up to ∆R CMB = ± 1000 km with respect to the standard R CMB of 3480 km. The figure shows that the ICAL detector would be able to measure the location of CMB at 1σ with a precision of about ± 380 km, ± 250 km, and ± 120 km for Case-I, Case-II, and Case-III, respectively. We can observe that the sensitivities are in general increasing as we go from Case-I to Case-II and Case-III. Note that Case-II is the most realistic one as mentioned earlier. Needless to mention that the CID capability of ICAL plays a crucial role to achieve this precision in locating R CMB . For an example, in the absence of CID, the 1σ precision for Case-II would be around ± 330 km.
The most striking feature is the asymmetry and non-monotonic behavior observed in ∆χ 2 CMB about the standard R CMB value in Case-II and Case-III. This turns out to be the effect of NOLR/ parametric resonance as explained in appendix A. We also evaluate the sensitivity for the 81-layered PREM profile for Case-II in appendix B. We find that the 1σ precision on R CMB is about ± 350 km.

Summary and Concluding Remarks
Understanding the detailed internal structure of Earth is an active field of research. This quest has been pursued traditionally using gravitational measurements and seismic studies. In recent times, neutrinos have emerged as important messengers to achieve this goal of multi-messenger tomography of Earth. For example, geoneutrinos can shed light on the composition and energy budget of Earth. The oscillations of atmospheric neutrinos at GeV energies and their absorption at TeV-PeV energies deep inside the Earth can provide complementary information on density profile and composition of Earth. We start by discussing how the information about the interior of Earth is obtained indirectly using gravitational measurements and seismic studies, and focus on how neutrinos can provide complementary information through their weak interactions.
While passing through different regions inside Earth, the multi-GeV atmospheric neutrinos experience the Earth's matter effects due to interactions with ambient electrons. These matter effects depend upon the density of electrons and alter the neutrino oscillation probabilities. In particular, neutrinos with energies 6-10 GeV experience MSW resonance while passing through mantle. The core-passing neutrinos with energies 3-6 GeV may also experience the neutrino parametric/ oscillation length resonance (NOLR) which depends upon the density jump at the core-mantle boundary (CMB) and its location. These effects make the neutrino signal at atmospheric neutrino detectors sensitive to the density profile of Earth. In the present work, we explore the effect of changing the radius of the CMB with respect to its standard value, R CMB = 3480 km, on the neutrino oscillation probabilities, and calculate the sensitivity of the ICAL detector at INO for measuring the CMB location.
We perform our analysis in the approximation of three-layered profile of Earth: core, inner mantle and outer mantle. While exploring the effect of changing CMB location, we consider the modification of the density profile of Earth in three ways. In Case-I, we modify R CMB while keeping the density in each layer constant; the total mass of Earth is not constrained. In Case-II, we modify R CMB and allow the density of core to modify such that the Earth's mass remains the same. In Case-III, we modify R CMB and the densities of core and inner mantle such that the masses of core and inner mantle, and hence that of Earth, are not affected. For each case, we consider two scenarios: SC (smaller core with ∆R CMB = −500 km) and LC (larger core with ∆R CMB = +500 km). For these scenarios, we observe that the effect of modification in R CMB manifests itself mainly for core-passing neutrinos in the NOLR/parametric resonance region. Moreover, for Case-III, these effects can also be seen in the MSW resonance region. The good direction and energy resolution of ICAL enables it to preserve these features at the level of reconstructed muons.
We estimate the sensitivity of ICAL for measuring the position of CMB in the three cases described above. To determine the precision on the location of CMB, we simulate the prospective data with the standard R CMB and calculate ∆χ 2 when the data are fitted with modified R CMB values. We observe that the sensitivity is significantly enhanced by the charge identification capability of ICAL, but not affected much by the uncertainties in oscillation parameters. The ICAL detector would be able to measure the location of CMB at 1σ with a precision of about ± 380 km, ± 250 km, and ± 120 km for Case-I, Case-II, and Case-III, respectively. For Case-II, we also perform the sensitivity study for locating CMB using the 81-layered PREM profile, and find the 1σ precision of about ± 350 km. Note that Case-II is the most realistic one, since it is based on the assumption that the density of mantle is known from seismic studies and the total mass of Earth is well measured from gravitational probes.
Also, we find that the sensitivity is not symmetric about the standard R CMB when we go to smaller or larger core radii. At large values of |∆R CMB | in Case-II and Case-III, the value of ∆χ 2 is not even monotonic. This a real physical effect, the origin of which may be traced to interesting interference effects in the NOLR/ parametric resonance region. As a result of this, larger core radii would be more constrained as compared to smaller core radii using atmospheric neutrino oscillation data.
In this study, we have found that the NOLR/ parametric resonance effects play a crucial role in measuring the radius of the CMB. Note that the NOLR/ parametric resonance comes into picture because of the contrast in the density of core and mantle, and specific relation between oscillation phases gained by neutrinos during their travel in mantle and core. Such quantum mechanical effects are possible with neutrino oscillations where details of phase change along the neutrino path are important, but not with neutrino absorption where only the integrated matter profile encountered along the path is relevant. Therefore, neutrino oscillations have a great potential for probing the internal structure of Earth, owing to the current precision of neutrino oscillation parameters.
The large amount of data from the next-generation atmospheric neutrino experiments like ORCA, IceCube/DeepCore/Upgrade, Hyper-K, DUNE, and P-ONE will significantly enhance the prospects of neutrino oscillation tomography of Earth. Over time, this data will augment the gravitational measurements and seismic studies to give us a clearer picture of the internal structure of Earth. Since neutrino oscillations are sensitive to the electron number densities, as opposed to the gravitational and seismic measurements that are sensitive to the baryonic number densities and material properties of the Earth's interior, neutrino data will provide independent and complementary information to these traditional probes.

A Asymmetric Effects with Smaller and Larger Core
While showing our results in Sec. 6, we have observed in Fig. 8 that ∆χ 2 CMB is asymmetric about the standard R CMB value in Case-II and Case-III. It also shows non-monotonic behavior at smaller test values of R CMB . In this appendix, we demonstrate the possible origin of these effects at the probability level in the context of Case-II using Figs. 9 and 10, which show the possible impact of smaller and larger core on three-flavor ν µ → ν µ oscillograms.
In Fig. 9, we present ν µ → ν µ oscillograms for smaller core radii with ∆R CMB = (−200, −400, −600, −800, −1000) km. Note that as R CMB decreases (core becomes smaller), the oscillograms show changes in the energy range of ∼ 3 to 15 GeV for core-passing neutrinos. We observe that the oscillation patterns for core-passing neutrinos get compressed towards cos θ ν = −1. Moreover, the oscillations in the NOLR/ parametric resonance region becomes more rapid as the core becomes smaller. This oscillatory behavior explains the non-monotonic nature of ∆χ 2 CMB for smaller core. In Fig. 10, we present ν µ → ν µ oscillograms for larger core radii with ∆R CMB = (+200, +400, +600, +800, +1000) km. Here also, as as R CMB increases (core becomes larger), we see changes in the oscillograms in the energy range of ∼ 3 to 15 GeV for core-passing neutrinos. However, here, the oscillation patterns for core-passing neutrinos get stretched towards smaller | cos θ ν | values. For ∆R CMB ≳ 400 km, the oscillatory behavior in the NOLR/ parametric resonance region becomes smooth, which accounts for the monotonic behavior of ∆χ 2 CMB at larger core radii. As a result of the above-mentioned differences between oscillation patterns for smaller and larger core radii, the values of ∆χ 2 CMB are not symmetric about ∆R CMB = 0 km, i.e., about the standard R CMB . The same observations described above for Case-II also hold for Case-III.

B Sensitivity Study using the 81-layered PREM Profile
Throughout the paper, we have considered a simple three-layered density profile of Earth to estimate the sensitivity of an atmospheric neutrino experiment like ICAL for locating R CMB . In this appendix, we explore the effect of going to a more detailed density profile on this sensitivity. We take the Earth density to be according to the 81-layered PREM  Table 2, where we assume NO and sin 2 θ23 = 0.5. Figure 10: Possible impact of larger core radii on three-flavor νµ → νµ oscillograms for the threelayered profile of Earth. The oscillograms are shown for RCMB = 3480 km + ∆RCMB where ∆RCMB = (0, +200, +400, +600, +800, +1000) km, with the density profile modified according to Case-II. We use three-flavor neutrino oscillation parameters given in Table 2, where we assume NO and sin 2 θ23 = 0.5.   Table 2, where we assume NO and sin 2 θ23 = 0.5. profile, where the density is interpolated/extrapolated using polynomial functions [117]. For the case of modified R CMB , the densities of all the layers in the core are scaled by the same fraction, such that the mass of the Earth remains invariant, just like in the Case-II which is discussed in section 3.1. The density distribution as a function of radial distance for the 81-layered PREM profile is shown in the left panel of Fig. 11, for ∆R CMB = ± 500 km.
The right panel of Fig. 11 shows the median sensitivities of the ICAL detector in terms of ∆χ 2 CMB as functions of the location of R CMB for Case-II. The value of R CMB is modified in small steps of 50 -100 km up to ∆R CMB = ± 1000 km with respect to the standard R CMB of 3480 km. Since the effects of the uncertainties of oscillation parameters are observed to be negligible in the fit (as can be seen from Table 5), we keep all the oscillation parameters fixed in the fit while evaluating ∆χ 2 CMB . The figure shows that, with the 81-layered PREM profile, ICAL would be able to measure the location of R CMB at 1σ confidence level with a precision of about ± 350 km. With the simple three-layered profile, this precision was about ± 250 km. One of the major reasons for this is the decrease in the density jump at R CMB in the 81-layered PREM profile.